I'm trying to solve the equation $$ (\Sigma \circ C)^{-1} \circ C = \left[ (\Sigma \circ C)^{-1} S (\Sigma \circ C)^{-1} \right] \circ C $$ for $\Sigma$ or $\Sigma \circ C$, where $\circ$ denotes the Hadamard or Schur element-wise product. $\Sigma$, $S$ and $C$ are symmetric real matrices. $\Sigma$ and $C$ are positive definite, while $S$ may be positive semidefinite. The elements of $C$ are $\in [0, 1]$, with 1s on the diagonal (which makes $\Sigma \circ C$ positive definite, too), and will in general be sparse. If it helps, we can assume that $C_{ij} = c(D_{ij})$, where $c()$ is a scalar decreasing function and $D$ is a Euclidean distance matrix.
If $C$ had no zero elements and $S$ were positive definite, the solution would be $$ (\Sigma \circ C) = S. $$
Is it possible to find a solution in the general case? I don't believe there is a closed-form solution, but maybe an algorithm to compute it?
The problem arises from the 1-taper approximation to the Wishart likelihood introduced by Kaufman, Schervich, and Nychka (J. Am. Stat. Ass. 2008). Other than in that paper, I don't want to assume a structure for $\Sigma$ beyond being symmetric positive definite.
Let's solve for the variable $$X = (C\circ\Sigma)^{-1}$$ In terms of this variable, the kernel of your equation becomes $$XSX-X = 0$$ We are free to add a matrix $M$ to the RHS of the equation as long as it satisfies the constraint $$C\circ M=0$$
If we denote the complement of $C$ by $B=(1-C)$, then one way to ensure that $M$ satisfies the constraint is to write $$M=B\circ A$$ where $A$ is completely unconstrained.
Now all we need to to is to solve a quadratic matrix equation $$XSX - X = M$$ For that, you can use this method, making the substitutions $$\eqalign{ C &= S \cr {\mathcal H} &= \begin{bmatrix} 0 & S\\ M & I \end{bmatrix} \cr }$$ Note that you are free to adjust $M$ until ${\mathcal H}$ has an acceptable set of eigenvalues/eigenvectors.