Wyjaśnij, w jaki sposób „własny” pomaga odwrócić macierz

Moje pytanie dotyczy techniki obliczeniowej wykorzystywanej w geoR:::.negloglik.GRFlub geoR:::solve.geoR.

W liniowym modelu mieszanym: gdzie i to odpowiednio efekty stałe i losowe. Ponadtoβ b Σ = cov ( Y

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

Podczas szacowania efektów konieczne jest obliczenie które normalnie można wykonać za pomocą czegoś podobnego , ale czasami jest prawie nieodwracalny, więc zastosuj lewę( X ′ Σ - 1 X

$$(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1} Y$$ solve(XtS_invX,XtS_invY)$(X'\Sigma^{-1}X)$geoR

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

(można zobaczyć w geoR:::.negloglik.GRFi geoR:::.solve.geoR), co oznacza rozkład gdzie a zatem

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

Dwa pytania:

1. W jaki sposób rozkład własny pomaga odwrócić ?$(X'\Sigma^{-1}X)$
2. Czy istnieją inne realne alternatywy (które są solidne i stabilne)? (np. qr.solveczy chol2inv?)

1) Skład eigend nie bardzo pomaga. Jest z pewnością bardziej stabilny numerycznie niż rozkład na czynniki Choleskiego, co jest pomocne, jeśli macierz jest źle uwarunkowana / prawie pojedyncza / ma wysoką liczbę warunków. Więc można użyć eigendecomposition i to daje A rozwiązanie problemu. Ale nie ma gwarancji, że będzie to PRAWIDŁOWE rozwiązanie. Szczerze mówiąc, po jawnym odwróceniu obrażenia już zostały wyrządzone. Formowanie tylko pogarsza sprawę. Skład eigend pomoże ci wygrać bitwę, ale wojna z pewnością jest przegrana.X T Σ - 1$\Sigma$$X^T \Sigma^{-1} X$

2) Nie znając specyfiki twojego problemu, właśnie to bym zrobił. Po pierwsze, należy wykonać na czynniki Cholesky'iego o , tak że . Następnie dokonaj faktoryzacji QR na , aby . Pamiętaj, aby obliczyć stosując substytucję przód - NIE wyraźnie inwertowany . Otrzymujesz: Odtąd możesz rozwiązać dowolną prawą stronę. Ale znowu,$\Sigma$$\Sigma = L L^T$$L^{-1} X$$L^{-1} X = QR$$L^{-1} X$$L$

$$\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$(lub ). W razie potrzeby użyj substytucji do przodu i do tyłu.$R^T R$

BTW, jestem ciekawa prawej strony twojego równania. Napisałeś, że to . Czy na pewno nie jest to ? Ponieważ gdyby tak było, możesz użyć podobnej sztuczki po prawej stronie: A potem możesz podać zamach stanu, gdy pójdziesz rozwiązać : $X^T \Sigma Y$$X^T \Sigma^{-1} Y$ β X T Σ - 1 X β = X T Σ - 1 Y R T R β = R T Q T L - 1 Y R β = Q T L - 1 Y β = R - 1 Q T L - 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$na ostatni krok, prawda? To tylko wsteczna zamiana. :-)