13

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?)
qoheleth
  • 1,210
  • 1
  • 12
  • 25

1 Answers1

15

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. :-)

Bill Woessner
  • 1,168
  • 11
  • 15
  • 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? – qoheleth Oct 30 '15 at 04:45
  • 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. – Bill Woessner Oct 30 '15 at 14:47