Verstehen, wie Numpy SVD macht

13

Ich habe verschiedene Methoden verwendet, um sowohl den Rang einer Matrix als auch die Lösung eines Matrixgleichungssystems zu berechnen. Ich bin auf die Funktion linalg.svd gestoßen. Vergleicht man dies mit meiner eigenen Anstrengung, das System mit der Gaußschen Eliminierung zu lösen, scheint es sowohl schneller als auch präziser zu sein. Ich versuche zu verstehen, wie das möglich ist.

Soweit ich weiß, verwendet die Funktion linalg.svd einen QR-Algorithmus, um die Eigenwerte meiner Matrix zu berechnen. Ich weiß, wie das mathematisch funktioniert, aber ich weiß nicht, wie Numpy es so schnell und ohne großen Genauigkeitsverlust schafft.

Also meine Frage: Wie funktioniert die Funktion numpy.svd und wie gelingt es, dies schnell und genau zu tun (im Vergleich zur Gauß-Eliminierung)?

RobVerheyen
quelle
2
numpy verwendet die Lapack-Routine dgesddfür echte SVDs. Ihre eigentliche Frage lautet also wahrscheinlich: "Wie funktioniert Lapack Dgesdd?"
Talonmies
Wenn Sie wirklich neugierig sind, würde ich vorschlagen, die LAPACK-Quelle zu untersuchen.
Vielen Dank für Ihre Kommentare, und ich entschuldige mich offtopic.
RobVerheyen
Dieser Beitrag ist ein Cross-Post von Stack Overflow . Cross-Posting wird auf Stack Exchange-Sites im Allgemeinen nicht empfohlen. Das Standardprotokoll zum erneuten Bereitstellen einer Frage auf einer anderen Site besteht darin, den ursprünglichen Beitrag zu schließen, zu löschen oder zu migrieren, bevor versucht wird, auf einer anderen Site erneut bereitzustellen. (Wenn Sie die Frage migrieren, wird sie automatisch erneut veröffentlicht.)
Geoff Oxberry
Es tut mir leid, das Protokoll war mir nicht bekannt. Ich hoffe ich bekomme noch eine Antwort.
RobVerheyen

Antworten:

15

Ihre Frage enthält eine Reihe von Problemen.

Verwenden Sie nicht die Gaußsche Eliminierung (LU-Faktorisierung), um den numerischen Rang einer Matrix zu berechnen. Die LU-Faktorisierung ist für diesen Zweck in der Gleitkommaarithmetik unzuverlässig. Verwenden Sie stattdessen eine rangaufschlussreiche QR-Zerlegung (z. B. xGEQPXoder xGEPQYin LAPACK, wobei x C, D, S oder Z ist, obwohl diese Routinen schwer zu finden sind; siehe die Antwort von JedBrown zu einer verwandten Frage ), oder verwenden Sie eine SVD (Singulärwertzerlegung wie xGESDDoder xGESVD, wobei x wieder C, D, S oder Z ist). Die SVD ist ein genauerer, zuverlässigerer Algorithmus zur Bestimmung des numerischen Ranges, erfordert jedoch mehr Gleitkommaoperationen.

Für die Lösung eines linearen Systems ist die LU-Faktorisierung (mit partiellem Schwenken, die Standardimplementierung in LAPACK) in der Praxis äußerst zuverlässig. Es gibt einige pathologische Fälle, in denen die LU-Faktorisierung mit partiellem Schwenken instabil ist (siehe Vorlesung 22 in Numerical Linear Algebra)von Trefethen und Bau für Details). Die QR-Faktorisierung ist ein stabilerer numerischer Algorithmus zum Lösen linearer Systeme, weshalb Sie wahrscheinlich so genaue Ergebnisse erhalten. Es erfordert jedoch mehr Gleitkommaoperationen als eine LU-Faktorisierung um den Faktor 2 für quadratische Matrizen (ich glaube, JackPoulson kann mich diesbezüglich korrigieren). Für rechteckige Systeme ist die QR-Faktorisierung die bessere Wahl, da sie Lösungen mit den kleinsten Quadraten für überbestimmte lineare Systeme liefert. SVD kann auch zum Lösen linearer Systeme verwendet werden, ist jedoch teurer als die QR-Faktorisierung.

janneb ist richtig, dass numpy.linalg.svd ein Wrapper xGESDDin LAPACK ist. Singulärwertzerlegungen laufen in zwei Stufen ab. Zunächst wird die zu zersetzende Matrix in eine bidiagonale Form reduziert. Der Algorithmus, der in LAPACK zum Reduzieren auf Bidiagonalform verwendet wird, ist wahrscheinlich der Lawson-Hanson-Chan-Algorithmus und verwendet an einem Punkt die QR-Faktorisierung. Die Vorlesung 31 in Numerical Linear Algebra von Trefethen und Bau gibt einen Überblick über diesen Prozess. Dann xGESDDwird ein Teile-und-Herrsche - Algorithmus die Singulärwerte und linke und rechte Singulärvektoren aus der Bidiagonalmatrix zu berechnen. Um Hintergrundinformationen zu diesem Schritt zu erhalten, müssen Sie Matrix Computations von Golub und Van Loan oder Applied Numerical Linear Algebra von Jim Demmel konsultieren .

Schließlich sollten Sie singuläre Werte nicht mit Eigenwerten verwechseln . Diese beiden Mengen sind nicht gleich. Die SVD berechnet die Singularwerte einer Matrix. Cleve Molers Numerical Computing mit MATLAB gibt einen schönen Überblick über die Unterschiede zwischen Singularwerten und Eigenwerten . Im Allgemeinen gibt es keine offensichtliche Beziehung zwischen den Singularwerten einer gegebenen Matrix und ihren Eigenwerten, außer im Fall von normalen Matrizen, bei denen die Singularwerte der Absolutwert der Eigenwerte sind.

Geoff Oxberry
quelle
Ich denke, "nicht verwandt" ist ziemlich stark für die Beziehung zwischen Eigenwerten und Singularwerten. Die Beziehung ist ziemlich dunkel, es sei denn, Sie kennen die vollständige Jordan-Zerlegung Ihrer Matrix, aber Sie können eine verwenden, um Schätzungen der anderen zu erhalten, wenn Sie Informationen über die Jordan-Zerlegung haben (oder bereit sind, Annahmen zu treffen).
Dan
Was würden Sie stattdessen vorschlagen?
Geoff Oxberry
Zunächst einmal vielen Dank für die ausführliche Antwort. Ich fand heraus, dass ich die LU-Zerlegung nicht verwenden kann, um den Matrixrang auf die harte Tour zu bestimmen. Ihre Antwort scheint zu implizieren, dass die QR-Faktorisierung tatsächlich eine schnellere Methode zur Lösung meines Problems wäre, richtig? Gibt es einen deutlichen Vorteil bei der Verwendung von SVD? Mir war klar, dass Singularwerte keine Eigenwerte sind. Ich bezog mich auf die Tatsache, dass Singularwerte als Eigenwerte der Matrix multipliziert mit ihrer Transponierung von links berechnet werden können. Es tut mir leid, dass das nicht klar war.
RobVerheyen
Ich könnte hinzufügen, dass die Matrix, die ich löse, tatsächlich singulär ist. Tatsächlich ist der Matrixrang nur etwa halb so groß wie die Matrix. Vielleicht ist dadurch eine Methode vorzuziehen?
RobVerheyen
1
@RobVerheyen: QR wird langsamer sein als LU, aber wesentlich genauer. SVD ist sogar langsamer als QR, aber SVD wird als die zuverlässigste Methode zur Bestimmung des numerischen Rangs angesehen (MATLAB verwendet beispielsweise die SVD in seiner rankFunktion). Bei beiden Ansätzen ist auch ein gewisses Maß an Diskretion erforderlich. Beim SVD-Ansatz ist der numerische Rang die Anzahl der Singularwerte über einem festgelegten (normalerweise sehr kleinen) Grenzwert. (Der QR-Ansatz ist ähnlich, ersetzt jedoch singuläre Werte durch diagonale Einträge der R-Matrix.)
Geoff Oxberry
8

Aufgrund des Wortlauts Ihrer Frage gehe ich davon aus, dass Ihre Matrix quadratisch ist. Die SVD-Routinen von LAPACK wie zgesvd verlaufen für quadratische Matrizen im Wesentlichen in drei Schritten:

  1. UEINVEINEINB: =UEINHEINVEINUEINVEINBÖ(n3)
  2. {UB,VB,Σ}B=UBΣVBHÖ(n2)Ö(n3)
  3. UEINBVEINH=EINEIN=(UEINUB)Σ(VEINVB)HUEINVEINUBVBÖ(n3)
Jack Poulson
quelle
7

numpy.linalg.svd ist ein Wrapper um {Z, D} GESDD von LAPACK. LAPACK wiederum wurde von einigen der weltweit führenden Experten für numerische lineare Algebra mit größter Sorgfalt geschrieben. Es wäre in der Tat sehr überraschend, wenn es jemandem, der mit dem Gebiet nicht vertraut ist, gelingen würde, LAPACK zu schlagen (entweder in Bezug auf Geschwindigkeit oder Genauigkeit).

Warum QR besser ist als Gaußsche Eliminierung, ist wahrscheinlich besser geeignet für /scicomp//

janneb
quelle
Vielen Dank für die Antwort und den Hinweis. Ich werde es dort versuchen.
RobVerheyen