Matrix equation with Hadamard product and its own inverse involved I know there is an almost exactly same question here but I have further specifications. So my problem is as follows:
$$
\Omega^{-1}=\dfrac{1}{n}\left(\Omega\odot \mathbf{W}+\mathbf{X}'\mathbf{X}+\lambda\mathbf{I}_{p}\right)
$$
where $\odot$ is the Hadamard product, $\Omega$ is a $p\times p$ positive semi-definite matrix, $\mathbf{X}$ is a $n\times p$ matrix, $\mathbf{I}_{p}$ is a $p\times p$ identity matrix, and $n,\lambda$ are scalars. One more condition is that $\mathbf{W}$ is a symmetric matrix with zero diagonal entries. I tried the following fixed-point iteration but it sometimes fails and sometimes alternates between two states...
$$
\Omega^{\left(k+1\right)}=n\left(\Omega^{\left(k\right)}\odot \mathbf{W}+\mathbf{X}'\mathbf{X}+\lambda\mathbf{I}_{p}\right)^{-1}
$$
Are there any some other ways that I can try? It doesn't necessarily have to be an analytical solution..
 A: Removing all unnecessary parameters, we come to the equation $\Omega^{-1}=2 W\odot \Omega + B$ where $B$ is positive definite. We need to find a solution in the cone $M_+$ of positive definite matrices. The solutions are exactly the stationary points of 
$F(\Omega)=\log\det\Omega- \operatorname{Tr} [(W\odot\Omega)\Omega+B\Omega]$ (I hope that is not where your equation came from in the first place). If the off-diagonal entries of $W$ are positive, you have a functional that is concave and tends to $-\infty$ on the boundary of $M_+$ and at $\infty$, so you just use any decent convex optimization algorithm to find the maximum. If not, the story gets way more complicated, so I'll stop here until you tell me that you need that option. 
A: Another standard thing you can try is Newton's method on the inverse-free form
$$
(\Omega \odot W) \Omega + A\Omega -nI = 0, \quad A=X'X+\lambda I
$$
Its Fréchet derivative (evaluated at a test matrix $H$) is
$$
J_\Omega(H) = (H\odot W)\Omega + (\Omega \odot W)H + AH,
$$
which is linear in $H$ and can be vectorized using Kronecker products and vectorizations. 
The iteration is 
$$J_{\Omega_k}(H_{k+1}) = (\Omega_k \odot W) \Omega_k + A\Omega_k -nI
$$
$$
\Omega_{k+1} = \Omega_k - H_{k+1}.
$$
The cost for iteration is going to be $O(p^6)$, though -- how big are your matrices?
