Explique cómo ayuda el `eigen` a invertir una matriz

Question

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

En una configuración de modelo lineal mixto: $$ Y=X\beta+Zb+e $$ donde $\beta$ y $b$ son los efectos fijos y aleatorios respectivamente. También, $\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, por lo que geoR emplear el truco

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

(puede verse en geoR:::.negloglik.GRF y geoR:::.solve.geoR ) lo que equivale a descomponer $$ (X'\Sigma^{-1}X)=\Lambda D \Lambda^{-1}\\ $$ donde $\Lambda'=\Lambda^{-1}$ y por lo tanto $$ (X'\Sigma^{-1}X)^{-1}=(D^{-1/2}\Lambda^{-1})'(D^{-1/2}\Lambda^{-1}) $$

Dos preguntas:

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

Answer 1

1) La eigendecomposición no ayuda mucho. Ciertamente es más estable numéricamente que una factorización Cholesky, lo cual es útil si tu matriz está mal condicionada/casi singular/tiene un número de condición alto. Así que puedes usar la eigendecomposición y te dará A solución a su problema. Pero hay pocas garantías de que sea la DERECHO solución. Honestamente, una vez que se invierte explícitamente $\Sigma$ El daño ya está hecho. Formando $X^T \Sigma^{-1} X$ sólo empeora las cosas. La eigendecomposición te ayudará a ganar la batalla, pero la guerra está perdida con toda seguridad.

2) Sin conocer los detalles de tu problema, esto es lo que yo haría. En primer lugar, realizar una factorización de Cholesky en $\Sigma$ para que $\Sigma = L L^T$ . A continuación, realice una factorización QR en $L^{-1} X$ para que $L^{-1} X = QR$ . Por favor, asegúrese de calcular $L^{-1} X$ utilizando la sustitución hacia adelante - NO invertir explícitamente $L$ . Así que entonces se obtiene: $$ \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} $$ A partir de aquí, puedes resolver el lado derecho que quieras. Pero, de nuevo, no inviertas explícitamente $R$ (o $R^T R$ ). Utilice las sustituciones hacia adelante y hacia atrás según sea necesario.

Por cierto, tengo curiosidad por la parte derecha de tu ecuación. Has escrito que es $X^T \Sigma Y$ . ¿Estás seguro de que no es $X^T \Sigma^{-1} Y$ ? Porque si fuera así, podrías utilizar un truco similar en el lado derecho: $$ \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} $$ Y luego puedes dar el golpe de gracia cuando vayas a resolver por $\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} $$ Por supuesto, nunca se invertiría explícitamente $R$ para el paso final, ¿verdad? Eso es sólo una sustitución hacia atrás. :-)