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:


*

*How does this eigen decomposition helps inverting $(X'\Sigma^{-1}X)$?

*Is there any other viable alternatives (that is robust and stable)? (e.g. qr.solve or chol2inv?)

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