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

• Is $\lambda\ge 0$ or not necessarily? – fedja Aug 26 '16 at 14:22
• @fedja oh I forgot to mention it. Yes $\lambda$ should be greater than zero, i.e., $\lambda>0$. – Daeyoung Lim Aug 26 '16 at 14:23
• possibly a silly idea but have you tried to do the iteration for $U=\Omega^{-1}$. So $U_{k+1} = \frac{1}{n} ( U_k^{-1}\odot W + X'X + \lambda I)$ ? – H. H. Rugh Aug 26 '16 at 17:44
• @H.H.Rugh isn't that essentially the same thing as the one I wrote in my question? But for the moment, i haven't tried that version. – Daeyoung Lim Aug 26 '16 at 17:46
• Not really, I just tested on some random example and it seems to work (on that example). It depends on how 'positive' the other terms are. – H. H. Rugh Aug 26 '16 at 17:47

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.

• That was the objective function! And yes the off-diagonal entries of $W$ are positive. – Daeyoung Lim Aug 26 '16 at 23:37
• I was actually trying to optimize a log-likelihood function with respect to the precision matrix $\Omega$ and the OP was the derivative set to zero. So that is exactly where my equation came from... Could you give me some examples of what the decent convex optimization algorithms would be..? – Daeyoung Lim Aug 27 '16 at 4:47
• @DaeyoungLim Gradient descent would be the first obvious thing to try. If it gives a satisfactory result in a decent time, I wouldn't look any further. – fedja Aug 28 '16 at 1:08
• I've tried the gradient ascent method but it seems to violate the positive semi-definiteness along the way. So I checked and adjusted the step-size accordingly to ensure positive-semidefiniteness but it seems like the gradient is giving me the wrong direction. The objective function sometimes decreases... – Daeyoung Lim Aug 29 '16 at 12:37
• If you leave the domain or the function goes down, it merely means that you've made a jump that is too large, so when you detect it, just go half-step (or if it is still too much quarter step, etc.) when you are at a bad point. The gradient just tells you in which direction to go and reduces each step to a one-dimensional problem. The exact optimal amount by which you should move at each step depends on the problem and is a mixture of art and science in the really interesting cases if you want fast convergence. – fedja Aug 29 '16 at 13:50

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?

• I encountered this equation while I was deriving an algorithm for a statistics estimation problem so I'm afraid there isn't a fixed matrix size... Could you clarify what $A$ and $H$ are? – Daeyoung Lim Aug 26 '16 at 18:05
• @DaeyoungLim Sorry, I forgot to define them. – Federico Poloni Aug 26 '16 at 18:14
• I'm sorry, but I don't understand how that iteration came to be. Why did you subtract the (k+1)th test matrix $H_{k+1}$ from $\Omega_{k}$? And why are the two Jacobian forms different? – Daeyoung Lim Aug 26 '16 at 18:22
• The two equations are essentially the same as in en.wikipedia.org/wiki/…, slightly rearranged and with an intermediate variable. In particular, the second formula for $J$ is not an expression for the Jacobian, it's an equation you have to solve for $H_{k+1}$ (by vectorizing and turning it into a linear system). – Federico Poloni Aug 26 '16 at 19:31