Explicar cómo `eigen` ayuda a invertir una matriz

Mi pregunta se refiere a una técnica de cálculo explotada en geoR:::.negloglik.GRFo geoR:::solve.geoR.

En una configuración de modelo mixto lineal: donde y son los efectos fijos y aleatorios, respectivamente. Además,

$$Y=X\beta+Zb+e$$ $\beta$$b$$\Sigma=\text{cov}(Y)$

Al estimar los efectos, es necesario calcular

$$(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1} Y$$ que normalmente se puede hacer usando algo como solve(XtS_invX,XtS_invY), pero a veces $(X'\Sigma^{-1}X)$ es casi no invertible, así que geoRemplea el truco

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

(puede verse en geoR:::.negloglik.GRFy geoR:::.solve.geoR) lo que equivale a descomponerse where y, por lo tanto,

$$(X'\Sigma^{-1}X)=\Lambda D \Lambda^{-1}\\ $$ ( X ′ Σ - 1 X ) - 1 = ( D $\Lambda'=\Lambda^{-1}$ $$(X'\Sigma^{-1}X)^{-1}=(D^{-1/2}\Lambda^{-1})'(D^{-1/2}\Lambda^{-1})$$

Dos preguntas:

1. ¿Cómo esta descomposición propia ayuda a invertir ?$(X'\Sigma^{-1}X)$
2. ¿Hay alguna otra alternativa viable (que sea robusta y estable)? (p. ej. qr.solveo chol2inv?)

1) La descomposición propia realmente no ayuda mucho. Ciertamente es más estable numéricamente que una factorización de Cholesky, lo que es útil si su matriz está mal acondicionada / casi singular / tiene un número de condición alto. Por lo tanto, puede usar la descomposición propia y le dará una solución a su problema. Pero hay pocas garantías de que será la solución CORRECTA . Honestamente, una vez que inviertes explícitamente , el daño ya está hecho. Formar simplemente empeora las cosas. La descomposición propia lo ayudará a ganar la batalla, pero la guerra ciertamente se pierde.X T Σ - 1$\Sigma$$X^T \Sigma^{-1} X$

2) Sin saber los detalles de su problema, esto es lo que haría. En primer lugar, realizar una factorización de Cholesky de de manera que . Luego realice una factorización QR en para que . Por favor asegúrese de cómputo mediante la sustitución hacia delante - NO explícitamente invertido . Entonces obtienes: Desde aquí, puede resolver cualquier lado derecho que desee. Pero otra vez,Σ = L L T L - 1 X L - 1 X = Q R L - 1 X L X T Σ - 1 X = X T ( L L T ) - 1$\Sigma$$\Sigma = L L^T$$L^{-1} X$$L^{-1} X = QR$$L^{-1} X$$L$ RRT

$$\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}$$ $R$(o ). Use sustituciones hacia adelante y hacia atrás según sea necesario.$R^T R$

Por cierto, tengo curiosidad sobre el lado derecho de su ecuación. Usted escribió que es . ¿Estás seguro de que no es ? Porque si lo fuera, podría usar un truco similar en el lado derecho: Y luego puede entregar el golpe de gracia cuando vaya a resolver : X T Σ - 1 Y X T Σ - 1 Y = X T ( L L T ) - 1$X^T \Sigma Y$$X^T \Sigma^{-1} Y$

$$\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}$$ $\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}$$ $R$para el paso final, ¿verdad? Eso es solo una sustitución hacia atrás. :-)