# Solving a Hadamard + Linear System

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.

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

• Thanks @greg! It looks like the big advantage here is that I get to solve $n$ $m \times m$ systems $\mathcal{O}(nm^3)$ instead of one $nm \times nm$ system $\mathcal{O}(n^3m^3)$ and avoid instantiating the big Kronecker product. For my specific application, I actually need to solve for different values of $C$ and the same values of $A, B$ many times, so using your approach with $n$ cached factorizations should work well. Commented Feb 11, 2019 at 16:36