# Explain how eigen helps inverting a matrix

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?)

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