13
$\begingroup$

My question relates to a computation technique exploited in geoR:::.negloglik.GRF or geoR:::solve.geoR.

In a linear mixed model setup: $$ Y=X\beta+Zb+e $$ where $\beta$ and $b$ are the fixed and random effects respectively. Also, $\Sigma=\text{cov}(Y)$

When estimating the effects, there is a need to compute $$ (X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1} Y $$ which can normally be done using something like solve(XtS_invX,XtS_invY), but sometimes $(X'\Sigma^{-1}X)$ is almost non-invertible, so geoR employ the trick

t.ei=eigen(XtS_invX)
crossprod(t(t.ei$vec)/sqrt(t.ei$val))%*%XtS_invY

(can be seen in geoR:::.negloglik.GRF and geoR:::.solve.geoR) which amounts to decomposing $$ (X'\Sigma^{-1}X)=\Lambda D \Lambda^{-1}\\ $$ where $\Lambda'=\Lambda^{-1}$ and therefore $$ (X'\Sigma^{-1}X)^{-1}=(D^{-1/2}\Lambda^{-1})'(D^{-1/2}\Lambda^{-1}) $$

Two questions:

  1. How does this eigen decomposition helps inverting $(X'\Sigma^{-1}X)$?
  2. Is there any other viable alternatives (that is robust and stable)? (e.g. qr.solve or chol2inv?)
$\endgroup$
15
+50
$\begingroup$

1) The eigendecomposition doesn't really help that much. It is certainly more numerically stable than a Cholesky factorization, which is helpful if your matrix is ill-conditioned/nearly singular/has a high condition number. So you can use the eigendecomposition and it will give you A solution to your problem. But there's little guarantee that it will be the RIGHT solution. Honestly, once you explicitly invert $\Sigma$, the damage is already done. Forming $X^T \Sigma^{-1} X$ just makes matters worse. The eigendecomposition will help you win the battle, but the war is most certainly lost.

2) Without knowing the specifics of your problem, this is what I would do. First, perform a Cholesky factorization on $\Sigma$ so that $\Sigma = L L^T$. Then perform a QR factorization on $L^{-1} X$ so that $L^{-1} X = QR$. Please be sure to compute $L^{-1} X$ using forward substitution - DO NOT explicitly invert $L$. So then you get: $$ \begin{array}{} X^T \Sigma^{-1} X & = & X^T (L L^T)^{-1} X \\ & = & X^T L^{-T} L^{-1} X \\ & = & (L^{-1} X)^T (L^{-1} X) \\ & = & (Q R)^T Q R \\ & = & R^T Q^T Q T \\ & = & R^T R \end{array} $$ From here, you can solve any right hand side you want. But again, please do not explicitly invert $R$ (or $R^T R$). Use forward and backward substitutions as necessary.

BTW, I'm curious about the right hand side of your equation. You wrote that it's $X^T \Sigma Y$. Are you sure it's not $X^T \Sigma^{-1} Y$? Because if it were, you could use a similar trick on the right hand side: $$ \begin{array}{} X^T \Sigma^{-1} Y & = & X^T (L L^T)^{-1} Y \\ & = & X^T L^{-T} L^{-1} Y \\ & = & (L^{-1} X)^T L^{-1} Y \\ & = & (Q R)^T L^{-1} Y \\ & = & R^T Q^T L^{-1} Y \end{array} $$ And then you can deliver the coup de grâce when you go to solve for $\beta$: $$ \begin{array}{} X^T \Sigma^{-1} X \beta & = & X^T \Sigma^{-1} Y \\ R^T R \beta & = & R^T Q^T L^{-1} Y \\ R \beta & = & Q^T L^{-1} Y \\ \beta & = & R^{-1} Q^T L^{-1} Y \end{array} $$ Of course, you would never explicitly invert $R$ for the final step, right? That's just a backward substitution. :-)

$\endgroup$
  • $\begingroup$ Thanks. this is a helpful response. Just to be explicit, is your chol/qr alternative going to help win the war? or just winning the game better than what eigen does? $\endgroup$ – qoheleth Oct 30 '15 at 4:45
  • $\begingroup$ That's a tough question to answer. I am confident that combining the Cholesky and QR factorizations will give you a better answer (and a faster answer). Whether it's the right answer really depends on the source of the problem. In this case, there are 2 potential sources. Either the columns of $X$ are nearly linearly dependent or $\Sigma$ is approaching singular. When you form $X^T \Sigma^{-1} X$, these problems amplify each other. The Cholesky+QR approach does not (and can not) mitigate either of these problems, but it does prevent the situation from getting worse. $\endgroup$ – Bill Woessner Oct 30 '15 at 14:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.