Mit dem Inversen einer positiv bestimmten symmetrischen (Kovarianz-) Matrix umgehen?

In der Statistik und ihren verschiedenen Anwendungen berechnen wir häufig die Kovarianzmatrix , die (in den betrachteten Fällen) positiv bestimmt und für verschiedene Verwendungen symmetrisch ist. Manchmal benötigen wir die Inverse dieser Matrix für verschiedene Berechnungen (quadratische Formen mit dieser Inverse als (einzige) Mittelmatrix zum Beispiel). Angesichts der Eigenschaften dieser Matrix und des Verwendungszwecks frage ich mich:

Was ist in Bezug auf die numerische Stabilität der beste Weg, um diese Inverse zu berechnen oder zu verwenden (z. B. für quadratische Formen oder Matrix-Vektor-Multiplikation im Allgemeinen)? Eine Faktorisierung, die nützlich sein kann?

Die Cholesky-Faktorisierung $C=R^TR$ führt zu einer Cholesky-ähnlichen Faktorisierung des Inversen $C^{-1}=SS^T$ mit der oberen Dreiecksmatrix $S=R^{-1}$ .

In der Praxis ist es am besten, den inversen Faktor beizubehalten. Wenn $R$ dünn ist, ist es normalerweise noch besser, $S$ implizit zu halten , da Matrixvektorprodukte $y=C^{-1}x$ durch Lösen der beiden Dreieckssysteme $R^Tz=x$ und berechnet werden können $Ry=z$.

Eine Cholesky-Faktorisierung ist für die beste Stabilität und Geschwindigkeit am sinnvollsten, wenn Sie mit einer Kovarianzmatrix arbeiten, da die Kovarianzmatrix eine positive semi-definite symmetrische Matrix ist. Cholesky ist hier selbstverständlich. ABER...

Wenn Sie beabsichtigen, eine Cholesky-Faktorisierung zu berechnen, bevor Sie jemals die Kovarianzmatrix berechnen, tun Sie sich selbst einen Gefallen. Machen Sie das Problem maximal stabil, indem Sie eine QR-Faktorisierung Ihrer Matrix berechnen. (Ein QR ist auch schnell.) Das heißt, wenn Sie die Kovarianzmatrix als berechnen würden

wo die Spaltenmittel entfernt hatte, dann sehen , dass , wenn Sie bilden C , es Quadrate die Konditionszahl. Es ist also besser, die QR-Faktoren von A zu bilden, als explizit eine Cholesky-Faktorisierung von A T A zu berechnen .$A$$C$$A$$A^{T}A$

Da Q orthogonal ist,

So erhalten wir den Cholesky - Faktor direkt von der QR - Faktorisierung, in Form von . Wenn eine Q- freie QR-Faktorisierung verfügbar ist, ist dies sogar noch besser, da Sie Q nicht benötigen . Ein QR ohne Q ist schnell zu berechnen, da Q nie generiert wird. Es wird nur eine Folge von Verwandlungen der Hausbesitzer. (Eine Säule drehte sich, $R^{T}$$Q$$Q$$Q$$Q$$Q$ -less QR, wäre logischerweise noch stabiler, was zusätzliche Arbeit bei der Auswahl der Pivots kostet.)

Die große Tugend der Verwendung des QR hier ist, dass er bei unangenehmen Problemen hoch numerisch stabil ist. Dies liegt wiederum daran, dass wir die Kovarianzmatrix nie direkt bilden mussten, um den Cholesky-Faktor zu berechnen. Sobald Sie das Produkt $A^{T}A$ , quadrieren Sie die Bedingungsnummer der Matrix. Tatsächlich verlieren Sie Informationen in den Teilen dieser Matrix, in denen Sie ursprünglich nur sehr wenige Informationen hatten.

Schließlich müssen Sie, wie eine andere Antwort hervorhebt, nicht einmal die Inverse berechnen und speichern, sondern sie implizit in Form von Backsolves auf Dreieckssystemen verwenden.

Ich habe dies kürzlich zum ersten Mal mit Vorschlägen von mathSE gemacht.

SVD wurde von den meisten empfohlen, aber ich habe mich für die Einfachheit von Cholesky entschieden:

Wenn die Matrix , zerlege ich M mit Cholesky in eine dreieckige Matrix L , so dass M = L L ⊤ ist . Ich verwende dann die Rückwärtssubstitution oder die Vorwärtssubstitution (abhängig davon, ob ich L als oberes oder unteres Dreieck wähle), um L zu invertieren , so dass ich L - 1 habe . Daraus kann ich schnell M - 1 = ( L L ⊤ ) - 1 = L - ⊤ L - 1 berechnen .$M = A A^\top$$M$$L$$M = L L^\top$$L$$L^{-1}$$M^{-1} = \left(L L^\top\right)^{-1} = L^{-\top}L^{-1}$

Beginnen mit:

, wobei$M = A A^\top$$M$ bekannt und implizit symmetrisch und auch positiv-definit ist.

Cholesky-Faktorisierung:

, wobei$M \rightarrow L L^\top$$L$ quadratisch und nicht singulär ist

Zurück Substitution:

, wahrscheinlich der schnellste Weg, um L zu invertieren$L \rightarrow L^{-1}$$L$ (zitiere mich aber nicht dazu)

Multiplikation:

$M^{-1} = \left(L L^\top\right)^{-1} = L^{-\top} L^{-1}$

Notation verwendet: Untere Indizes sind Zeilen, obere Indizes sind Spalten und ist die Transponierte von L - $L^{-\top}$$L^{-1}$

Mein Cholesky-Algorithmus (wahrscheinlich aus Numerical Recipes oder Wikipedia)

$L_i^j = \frac{M_i^j - M_i \cdot M_j}{M_i^i - M_i \cdot M_i}$

Dies kann fast vor Ort erfolgen (Sie benötigen nur einen temporären Speicher für die diagonalen Elemente, einen Akkumulator und einige ganzzahlige Iteratoren).

Mein Back-Substitution-Algorithmus (siehe Numerische Rezepte, überprüfen Sie deren Version, da ich möglicherweise einen Fehler beim LaTeX-Markup gemacht habe)

$\left(L^{-1}\right)_i^j = \left\{\begin{array}{11} 1 / {L_i^i} & \mbox{if } i = j\\ \left(-L_i \cdot \left(L^{-T}\right)_j\right) / L_i^i & \mbox{otherwise} \end{array}\right.$

$L^{-T}$

Wenn Sie wissen, dass die Matrix eine Inverse hat (dh, wenn sie tatsächlich positiv definit ist) und wenn sie nicht zu groß ist, bietet die Cholesky-Zerlegung eine geeignete Möglichkeit, die Inverse einer Matrix zu charakterisieren.