Solving a Hadamard + Linear System

Question

I need to solve a matrix equation which contains a Hadamard product and a standard matrix multiplication: $$A\odot X + BX = C$$ where $A, C, X \in \mathbb{R}^{m \times n}$, $B \in \mathbb{R}^{m \times m}$, and $\odot$ denotes the element-wise (Hadamard) matrix product.

Does this sort of system have a name? I've checked the standard resources (Golub and van Loan, Horn and Johnson), but I didn't find any discussion of this type of problem.

If we let $a = \text{vec}(A)$, $x = \text{vec}(X)$, $c = \text{vec}(C)$ and $\tilde{a} = \text{diag}(a)$, then this becomes: $$a \odot x + (I \otimes B)x = \tilde{a}x + (I \otimes B)x = c$$ which has solution $$x = (\tilde{a} + (I \otimes B))^+c \implies X = \text{unvec}\left[\left(\text{diag}(\text{vec}(A)) + I\otimes B\right)^+\text{vec}(C)\right]$$ where $\otimes$ is the Kronecker product and $(\cdot)^+$ is the (Moore-Penrose) pseudo-inverse.

While it's not too expensive to use an iterative solver (since the $\tilde{a} + (I \otimes B)$ operator can be evaluated more-or-less easily), but I'd like to avoid that if possible. Is it possible to derive an analytical solution for this class of problems?

A similar problem is discussed at https://math.stackexchange.com/questions/2583719/how-to-solve-the-linear-equation-a-circ-xb-cx-d, but they stop at the "vec'd-out" solution.

Answer 1

If you multiply a Hadamard product (from the left or right) by a vector from the standard basis, it distributes over the terms in the product like so $$\eqalign{ e_k^T(A\odot B) &= (e_k^TA)\odot(e_k^TB) \cr (A\odot B)e_k &= (Ae_k)\odot(Be_k) \cr }$$ We can use this to solve for the $X$ matrix, one column at a time $$\eqalign{ (A\odot X)e_k + (BX)e_k &= Ce_k \cr (Ae_k)\odot(Xe_k) + B(Xe_k) &=Ce_k \cr \Big({\rm Diag}(a_k)+B\Big)\,x_k &= c_k \cr x_k &= \Big({\rm Diag}(a_k)+B\Big)^+c_k \cr }$$ Here is some Julia code (version 1.1.0) which demonstrates the technique.

using LinearAlgebra

m,n = 30,50;
X,A,C,B = zeros(m,n), randn(m,n), randn(m,n), randn(m,m);

for k = 1:n
  X[:,k] = (Diagonal(A[:,k]) + B) \ C[:,k];
end

norm(A.*X + B*X - C)
1.3032451036976518e-12