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..
-
$\begingroup$ Is $\lambda\ge 0$ or not necessarily? $\endgroup$– fedjaCommented Aug 26, 2016 at 14:22
-
$\begingroup$ @fedja oh I forgot to mention it. Yes $\lambda$ should be greater than zero, i.e., $\lambda>0$. $\endgroup$– DaeyoungCommented Aug 26, 2016 at 14:23
-
$\begingroup$ 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)$ ? $\endgroup$– H. H. RughCommented Aug 26, 2016 at 17:44
-
1$\begingroup$ @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. $\endgroup$– DaeyoungCommented Aug 26, 2016 at 17:46
-
$\begingroup$ 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. $\endgroup$– H. H. RughCommented Aug 26, 2016 at 17:47
2 Answers
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.
-
$\begingroup$ That was the objective function! And yes the off-diagonal entries of $W$ are positive. $\endgroup$– DaeyoungCommented Aug 26, 2016 at 23:37
-
1$\begingroup$ 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..? $\endgroup$– DaeyoungCommented Aug 27, 2016 at 4:47
-
$\begingroup$ @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. $\endgroup$– fedjaCommented Aug 28, 2016 at 1:08
-
$\begingroup$ 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... $\endgroup$– DaeyoungCommented Aug 29, 2016 at 12:37
-
1$\begingroup$ 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. $\endgroup$– fedjaCommented Aug 29, 2016 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?
-
$\begingroup$ 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? $\endgroup$– DaeyoungCommented Aug 26, 2016 at 18:05
-
$\begingroup$ @DaeyoungLim Sorry, I forgot to define them. $\endgroup$ Commented Aug 26, 2016 at 18:14
-
$\begingroup$ 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? $\endgroup$– DaeyoungCommented Aug 26, 2016 at 18:22
-
1$\begingroup$ 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). $\endgroup$ Commented Aug 26, 2016 at 19:31