Suppose we have two matrices $A$ and $B$ (we can assume they're symmetric; if absolutely necessary I think they may be positive definite). Then, is there any technique for solving $$(A\circ B)x=b,$$ where $\circ$ denotes the Hadamard (elementwise) product?

I am seeking a direct solver if possible rather than iterative methods like CG.

• Which field does this Hadamard system come from ? If we can implement the (A o B)x fast and efficiently, I think that the iterative solvers for this type of linear systems is not bad. Sep 24, 2014 at 20:45
• The systems come from a finite elements discretization of a geometry problem. The Hadamard product falls out of a weirdly-shaped least-squares problem. I can implement the product somewhat efficiently but will have to solve many instances of this system with different right-hand sides -- I'd prefer to find some kind of prefactorization. Sep 24, 2014 at 22:55
• OK, I've added a second question showing where this system comes up: scicomp.stackexchange.com/questions/14713/… Sep 26, 2014 at 18:35

Ultimately, it depends on the sparsity of $A$ and $B$ and the symmetry of the resulting hadamard product.

The hadamard product will aggregate the sparsity structure of $A$ and $B$. So if one or the other is sparse, the product is also sparse (and may be more so if $A$ and $B$ have difference sparsity structures). So any sparse direct solver would be appropriate. However, if the matrix is extremely large, memory becomes an issue and sparse iterative methods are just about the only choice.

Even if $A$ and $B$ are dense, the hadamard product is an $O(n^2)$ operation, which is insignificant compared to a direct linear solve of $O(n^3)$ operations. So even if A and B are extremely large dense matrices, simply computing the hadamard product before solving the resulting system is not that much more work compare to solving the system.

In the dense case, the key question becomes whether the resulting matrix is hermitian or not. In the worst case (non-hermitian hadamard product), gaussian elimination is pretty much the only thing you can do. If it is hermitian, a Cholesky (definite ) or Bunch-Kaufman (indefinite) algorithm would be appropriate.

• Sadly the matrices aren't sparse, but their blocks have nice decompositions (e.g. as $BB^\top$ for a "tall" matrix $B$). I was hoping there was some way to think about $(A\circ B)^{-1}$ without just putting it all together, e.g. through the FFT. Sep 24, 2014 at 19:56
• You'll need to tell us more about the structure you wish to preserve then. A $O(n\log n)$ FFT will not suffice in general, for sure: at the very least because your input size is $O(n^2)$. Sep 25, 2014 at 12:05
• Hmm, it looks like a generic trick won't work here! I'll make a second question on stack exchange with more detail. Thanks. Sep 26, 2014 at 18:22

The nonzero elements of the product are exactly those locations that are nonzero in both matrices. In other words, the sparsity pattern of the product is a subset of the sparsity patterns of the two matrices, and a cheap way to form the product is to take the sparser of the two matrices and for each one of its elements replace it by the product of the two matrices. This can be done in at most ${\cal O}(nnz_A + nnz_B)$ operations ($nnz$ = number of nonzero entries) which is at the very worst ${\cal O}(N^2)$ but in practice for sparse matrices is of course much less and typically only ${\cal O}(N)$.

You'd then use your usual solver to solve the linear problem. This will in general cost much more than forming the product in place.