Expliquez comment `eigen` aide à inverser une matrice

13

Ma question concerne une technique de calcul exploitée dans geoR:::.negloglik.GRFou geoR:::solve.geoR.

Dans une configuration de modèle mixte linéaire: où et sont respectivement les effets fixes et aléatoires. En outre,

Y=Xβ+Zb+e
b Σ = cov ( Y )βbΣ=cov(Y)

Lors de l'estimation des effets, il est nécessaire de calculer qui peut normalement être fait en utilisant quelque chose comme , mais parfois est presque non inversible, alors employez l'astuce( X Σ - 1 X )

(XΣ1X)1XΣ1Y
solve(XtS_invX,XtS_invY)(XΣ1X)geoR
t.ei=eigen(XtS_invX)
crossprod(t(t.ei$vec)/sqrt(t.ei$val))%*%XtS_invY

(peut être vu dans geoR:::.negloglik.GRFet geoR:::.solve.geoR) ce qui revient à se décomposer où et donc

(XΣ1X)=ΛDΛ1
( X Σ -Λ=Λ1
(XΣ1X)1=(D1/2Λ1)(D1/2Λ1)

Deux questions:

  1. Comment cette décomposition propre aide-t-elle à inverser ?(XΣ1X)
  2. Existe-t-il d'autres alternatives viables (robustes et stables)? (par exemple qr.solveou chol2inv?)
qoheleth
la source

Réponses:

15

1) La composition eigendec n'aide pas vraiment beaucoup. Elle est certainement plus stable numériquement qu'une factorisation de Cholesky, ce qui est utile si votre matrice est mal conditionnée / presque singulière / a un nombre élevé de conditions. Vous pouvez donc utiliser la composition eigendec et cela vous donnera une solution à votre problème. Mais rien ne garantit que ce sera la bonne solution. Honnêtement, une fois que vous inversez explicitement , le mal est déjà fait. La formation de X T Σ - 1 X ne fait qu'empirer les choses. La composition par eigend vous aidera à gagner la bataille, mais la guerre est certainement perdue.ΣXTΣ-1X

2) Sans connaître les spécificités de votre problème, voici ce que je ferais. Tout d' abord, effectuer une factorisation de Cholesky sur pour que Σ = L L T . Puis effectuer une factorisation QR sur L - 1 X de telle sorte que L - 1 X = Q R . Veuillez vous assurer de calculer L - 1 XΣΣ=LLTL-1XL-1X=QRL-1X en utilisant la substitution de l' avant - NE PAS explicitement inverti . Alors vous obtenez: X T Σ - 1 X = X T ( LL De là, vous pouvez résoudre n'importe quel côté droit que vous voulez. Mais encore une fois, veuillez ne pas inverser explicitementR(ouRTR). Utilisez des substitutions avant et arrière si nécessaire.

XTΣ-1X=XT(LLT)-1X=XTL-TL-1X=(L-1X)T(L-1X)=(QR)TQR=RTQTQT=RTR
RRTR

BTW, je suis curieux de savoir le côté droit de votre équation. Vous avez écrit que c'est . Êtes-vous sûr que ce n'est pas X T Σ - 1 Y ? Parce que si c'était le cas, vous pourriez utiliser une astuce similaire sur le côté droit: X T Σ - 1XTΣOuiXTΣ-1Oui Et puis vous pouvez livrer le coup de grâce quand vous allez résoudre pourβ: X T Σ - 1 X β = X T Σ - 1 Y R T R β = R T Q T L - 1 Y R β = Q T L -

XTΣ-1Oui=XT(LLT)-1Oui=XTL-TL-1Oui=(L-1X)TL-1Oui=(QR)TL-1Oui=RTQTL-1Oui
β Bien sûr, vous n'inverseriez jamais explicitementRpour la dernière étape, non? C'est juste une substitution en arrière. :-)
XTΣ-1Xβ=XTΣ-1OuiRTRβ=RTQTL-1OuiRβ=QTL-1Ouiβ=R-1QTL-1Oui
R
Bill Woessner
la source
Merci. c'est une réponse utile. Juste pour être explicite, votre alternative chol / qr va-t-elle aider à gagner la guerre? ou simplement gagner le jeu mieux que ce que fait eigen?
qoheleth
XΣXTΣ-1X