Schnelle Berechnung / Schätzung eines linearen Systems mit niedrigem Rang

10

Lineare Gleichungssysteme sind in der Computerstatistik allgegenwärtig. Ein spezielles System, auf das ich gestoßen bin (z. B. in der Faktoranalyse), ist das System

Ax=b

wobei Hier ist eine Diagonalmatrix mit einer streng positiven Diagonale, ist eine (mit ) symmetrische positive semidefinitive Matrix und ist eine beliebige Matrix. Wir werden gebeten, ein diagonales lineares System (einfach) zu lösen, das durch eine niedrigrangige Matrix gestört wurde. Der naive Weg, um das obige Problem zu lösen, besteht darin, mit der Woodbury-Formel zu invertieren . Dies fühlt sich jedoch nicht richtig an, da Cholesky- und QR-Faktorisierungen die Lösung linearer Systeme (und normaler Gleichungen) normalerweise dramatisch beschleunigen können. Ich bin vor kurzem auf die gekommen

A=D+BΩBT
Dn×nΩm×mmnBn×mADas folgende Papier scheint den Cholesky-Ansatz zu verfolgen und erwähnt die numerische Instabilität von Woodburys Inversion. Das Papier scheint jedoch in Entwurfsform zu sein, und ich konnte keine numerischen Experimente oder unterstützende Forschung finden. Wie ist der Stand der Technik, um das von mir beschriebene Problem zu lösen?
gappy
quelle
1
@gappy, haben Sie in Betracht gezogen, die QR-Zerlegung (oder Cholesky-Zerlegung) für die Matrix (der mittlere Term in der Woodburry-Formel) zu verwenden? Die verbleibenden Operationen sind einfache Matrixmultiplikationen. Die Hauptursache für Instabilität ist dann die Berechnung von . Da vermute ich, dass diese Anwendung von QR oder Cholesky in Kombination mit Woodbury auf allen Matrix schneller als QR ist . Dies ist natürlich kein Stand der Technik, sondern nur allgemeine Beobachtungen. Ω1+BD1BTΩ1m<<nA
mpiktas
Ich vermute , dass Matthias Seeger , was innerhalb befürwortet ist des Standes der Technik, er ist ein sehr helles Macker und diese Art von Fragen auftaucht immer wieder in der Art von Modellen untersucht er. Ich benutze Cholesky-basierte Methoden aus den gleichen Gründen. Ich vermute, dass es in "Matrix Computations" von Golub und Van Loan eine Diskussion gibt, die die Standardreferenz für solche Dinge ist (obwohl ich meine Kopie nicht zur Hand habe). ϵ
Dikran Marsupial
Beachten Sie, dass mit Ihr Problem der Lösung des Systems entspricht wobei . Das vereinfacht das Problem genau dort ein bisschen. Wenn wir nun , wissen wir, dass positiv semidefinit mit höchstens positiven Eigenwerten ist. Seit kann das Finden der größten Eigenwerte und der entsprechenden Eigenvektoren auf verschiedene Arten erfolgen. Die Lösung ist dann wobei die Eigenzusammensetzung von ergibtB¯=D1/2B(I+B¯ΩB¯T)x=b¯b¯=D1/2bΣ=B¯ΩB¯TΣmmnmx=Q(I+Λ)1QTb¯Σ=QΛQTΣ.
Kardinal
Kleine Korrekturen: (1) Das äquivalente System ist und (2) Die endgültige Lösung ist . (Ich hatte in beiden Fällen ein vor abgelegt .) Beachten Sie, dass alle Inversen diagonale Matrizen sind und daher trivial sind. (I+B¯ΩB¯T)D1/2x=b¯x=D1/2Q(I+Λ)1QTD1/2bD1/2x
Kardinal
@mpiktas: Ich denke, Sie meinten da in der von Ihnen geschriebenen Version das Matrixprodukt aufgrund einer Dimensionsinkongruenz nicht genau definiert ist. :)Ω1+BTD1B
Kardinal

Antworten:

2

In "Matrix Computations" von Golub & van Loan wird in Kapitel 12.5.1 ausführlich auf die Aktualisierung von QR- und Cholesky-Faktorisierungen nach Rang-P-Aktualisierungen eingegangen.

Fabian Pedregosa
quelle
Ich weiß, und die relevanten Lapack-Funktionen werden sowohl in dem Artikel, den ich verlinkt habe, als auch in dem Buch erwähnt. Ich frage mich jedoch, was die beste Vorgehensweise für das vorliegende Problem ist, nicht für das generische Aktualisierungsproblem.
Gappy