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$– fedjaAug 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$– DaeyoungAug 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. RughAug 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$– DaeyoungAug 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. RughAug 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$– DaeyoungAug 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$– DaeyoungAug 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$– fedjaAug 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$– DaeyoungAug 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$– fedjaAug 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$– DaeyoungAug 26, 2016 at 18:05
-
$\begingroup$ @DaeyoungLim Sorry, I forgot to define them. $\endgroup$ 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$– DaeyoungAug 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$ Aug 26, 2016 at 19:31