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.
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.
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.
