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. – Hsien-Ming Ku Sep 24 '14 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. – Justin Solomon Sep 24 '14 at 22:55
• OK, I've added a second question showing where this system comes up: scicomp.stackexchange.com/questions/14713/… – Justin Solomon Sep 26 '14 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. – Justin Solomon Sep 24 '14 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)$. – Federico Poloni Sep 25 '14 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. – Justin Solomon Sep 26 '14 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.