Wann übertreffen orthogonale Transformationen die Gaußsche Elimination?

22

Wie wir wissen, sind orthogonale Transformationsmethoden (Givens-Rotationen und Housholder-Reflexionen) für lineare Gleichungssysteme teurer als die Gauß-Elimination, haben jedoch theoretisch bessere Stabilitätseigenschaften in dem Sinne, dass sie die Bedingungszahl des Systems nicht ändern. Obwohl ich nur ein akademisches Beispiel für eine Matrix kenne, die durch Gaußsche Eliminierung mit partiellem Schwenken zerstört wird. Und es gibt die gängige Meinung, dass es sehr unwahrscheinlich ist, dass ein solches Verhalten in der Praxis erreicht wird (siehe diese Vorlesungsunterlagen [pdf] ).

Wo sollen wir also nach der Antwort zum Thema suchen? Parallele Implementierungen? Aktualisierung?..

Faleichik
quelle

Antworten:

24

Richtigkeit

Trefethen und Schreiber haben eine exzellente Arbeit geschrieben, Average-Case Stability of Gaussian Elimination , in der die Genauigkeitsseite Ihrer Frage erörtert wird. Hier sind einige seiner Schlussfolgerungen:

  1. „Für QR - Faktorisierung mit oder ohne Säule schwenkbar, das durchschnittliche maximale Element der Restmatrix ist O(n1/2) , während für Gaußsche Eliminations ist O(n) . Dieser Vergleich zeigt , dass die Gaußsche Eliminations leicht instabil ist, aber die Instabilität kann nur bei sehr großen Matrixproblemen festgestellt werden, die mit geringer Genauigkeit gelöst werden. Bei den meisten praktischen Problemen ist die Gaußsche Elimination im Durchschnitt sehr stabil. "(Schwerpunkt Mine)

  2. "Nach den ersten Schritten der Gaußschen Eliminierung sind die verbleibenden Matrixelemente ungefähr normal verteilt, unabhängig davon, ob sie auf diese Weise begonnen haben."

Es gibt viel mehr in dem Artikel, das ich hier nicht erfassen kann, einschließlich der Diskussion der Worst-Case-Matrix, die Sie erwähnt haben. Ich empfehle Ihnen daher nachdrücklich, sie zu lesen.

Performance

Für quadratische reelle Matrizen erfordert LU mit teilweisem Schwenken ungefähr Flops, während Householder-basierte QR grob erfordert 4 / 3 n 3 - Flop. Für einigermaßen große Quadratmatrizen ist die QR-Faktorisierung daher nur etwa doppelt so teuer wie die LU-Faktorisierung.2/3n34/3n3

Für - Matrizen, wobei m n , LU mit teilweisem Verschwenkung erfordert m n 2 - n 3 / 3 - Flop, im Vergleich zu der QR 2 m n 2 - 2 n 3 / 3 (das immer noch doppelt so groß ist , dass die LU - Faktorisierung). Allerdings ist es überraschend häufig für Anwendungen sehr hohe dünne Matrizen (zur Herstellung von m » n ) und Demmel et al. Ich habe ein nettes Referat, Kommunikationsvermeidende parallele und sequentielle QR-Faktorisierungm×nmnmn2-n3/32mn2-2n3/3mn, in dem (in Abschnitt 4) ein cleverer Algorithmus erörtert wird, der nur das Senden von Nachrichten erfordert , wenn p Prozessoren verwendet werden, im Vergleich zu den n Protokoll- p- Nachrichten herkömmlicher Ansätze. Der Aufwand besteht darin, dass O ( n 3 log p ) zusätzliche Flops durchgeführt werden, aber für sehr kleine n wird dies oft den Latenzkosten für das Senden von mehr Nachrichten vorgezogen (zumindest wenn nur eine einzige QR-Faktorisierung durchgeführt werden muss).logppnlogpO(n3logp)n

Jack Poulson
quelle
10

Ich bin überrascht, dass niemand Probleme mit linearen kleinsten Quadraten erwähnt hat , die beim wissenschaftlichen Rechnen häufig auftreten. Wenn Sie die Gaußsche Elimination verwenden möchten, müssen Sie die normalen Gleichungen bilden und lösen, die wie folgt aussehen:

ATAx=ATb,

woher eine Matrix von Datenpunkten ist, die Beobachtungen von unabhängigen Variablen entsprechen, x ein Vektor von zu findenden Parametern ist und b ein Vektor von Datenpunkten ist, die Beobachtungen einer abhängigen Variablen entsprechen.Axb

Wie Jack Poulson häufig betont, ist die Zustandsnummer von das Quadrat der Bedingungszahl von A , so dass die normalen Gleichungen katastrophal schlecht konditioniert werden können. In solchen Fällen sind QR- und SVD-basierte Ansätze zwar langsamer, liefern jedoch viel genauere Ergebnisse.ATAA

Geoff Oxberry
quelle
2
Upvoted, aber QR soll eigentlich auf dem Niveau sein mit LU , wenn Sie den unnötigen betrachten Operationen erforderlich bilden A H A (QR erfordert nur 2 / 3 n 3 mehr Flops als LU). Der SVD-Ansatz sollte jedoch langsamer sein (man kann sich seine Kosten als ungefähr 6 n 3 vorstellen ). n3AHA2/3n36n3
Jack Poulson
1
Neben der durch die Verwendung von orthogonalen Transformationen garantierten Stabilität besteht der große Vorteil der SVD darin, dass die Zerlegung eine eigene Zustandsüberprüfung bietet, da das Verhältnis des größten zum kleinsten Singularwert genau die (2-Norm-) Zustandszahl ist. Für die anderen Zerlegungen ist die Verwendung eines Zustandsschätzers (z. B. Hager-Higham), obwohl nicht so teuer wie die eigentliche Zerlegung, etwas "angeheftet".
JM
1
@JackPoulson Hast du aus Neugier eine Referenz für deinen Flop Count für SVD? Nach dem, was ich aus einem kurzen Blick auf Golub & Van Loan (S. 254, 3. Auflage) ersehen kann, scheint die Konstante für die Verwendung der SVD bei der Lösung von Least-Squares-Problemen höher zu sein, aber ich könnte mich täuschen. Danke im Voraus.
OscarB
1
@OscarB: Es war eine sehr grobe Zahl von oben auf meinem Kopf, die niedriger ist als die Bildung der vollständigen SVD (weil wir Rücktransformationskosten vermeiden können). Arbeit für die Reduktion zu bidiagonal Form benötigt wird (beispielsweise A = F B G H ), eine gewisse Menge an Arbeit, sagen wir C , wird für die bidiagonal SVD (benötigt B = U Σ V H ), und dann x : = ( G ( V ( i n v ( Σ ) ( U H (8/3n3A=FBGHCB=UΣVH , was O ( n 2 ) Arbeiterfordern sollte. Es ist also alles eine Frage derGrößevon C ... wenn MRRR jemals hier funktioniert, ist es O ( n 2 ) , aber bis dahin ist es kubisch und problemabhängig. x:=(G(V(inv(Σ)(UH(FHb)))))O(n2)CO(n2)
Jack Poulson
1
@JM Beachten Sie jedoch, dass die Bedingungsnummer des Problems der kleinsten Quadrate nicht die "klassische" Bedingungsnummer eine Matrix; es ist eine kompliziertere Menge. σ1σn
Federico Poloni
3

Wie messen Sie die Leistung? Geschwindigkeit? Richtigkeit? Stabilität? Ein schneller Test in Matlab ergibt Folgendes:

>> N = 100;
>> A = randn(N); b = randn(N,1);
>> tic, for k=1:10000, [L,U,p] = lu(A,'vector'); x = U\(L\b(p)); end; norm(A*x-b), toc
ans =
   1.4303e-13
Elapsed time is 2.232487 seconds.
>> tic, for k=1:10000, [Q,R] = qr(A); x = R\(Q'*b); end; norm(A*x-b), toc             
ans =
   5.0311e-14
Elapsed time is 7.563242 seconds.

Das Lösen eines einzelnen Systems mit einer LU-Zerlegung ist also etwa dreimal so schnell wie das Lösen mit einer QR-Zerlegung, und kostet eine halbe Dezimalstelle Genauigkeit (in diesem Beispiel!).

Pedro
quelle
Alle von Ihnen vorgeschlagenen Vorteile sind willkommen.
Faleichik
3

Der Artikel, den Sie zitieren, verteidigt die Gaußsche Eliminierung, indem er sagt, dass, obwohl es numerisch instabil ist, es bei Zufallsmatrizen gut abschneidet, und da die meisten Matrizen, die man sich vorstellen kann, wie Zufallsmatrizen sind, sollten wir in Ordnung sein. Dieselbe Aussage kann für viele numerisch instabile Methoden getroffen werden.

Betrachten Sie den Raum aller Matrizen. Diese Methoden funktionieren fast überall. Das sind 99,999 ...% aller Matrizen, die man erstellen kann, haben keine Probleme mit instabilen Methoden. Es gibt nur einen sehr kleinen Bruchteil von Matrizen, für die GE und andere Schwierigkeiten haben werden.

Die Probleme, die den Forschern am Herzen liegen, liegen in der Regel in dieser kleinen Fraktion.

Wir konstruieren Matrizen nicht zufällig. Wir konstruieren Matrizen mit ganz besonderen Eigenschaften, die ganz besonderen, nicht zufälligen Systemen entsprechen. Diese Matrizen sind oft schlecht konditioniert.

Geometrisch können Sie den linearen Raum aller Matrizen berücksichtigen. Es gibt einen Subraum mit einem Volumen / Maß von Null von singulären Matrizen, die diesen Raum durchschneiden. Viele Probleme, die wir konstruieren, konzentrieren sich auf diesen Unterraum. Sie werden nicht zufällig verteilt.

Betrachten Sie als Beispiel die Wärmegleichung oder -dispersion. Diese Systeme neigen dazu, Informationen aus dem System zu entfernen (alle Anfangszustände tendieren zu einem einzigen Endzustand), und als Ergebnis sind Matrizen, die diese Gleichungen beschreiben, enorm singulär. Dieser Prozess ist in einer zufälligen Situation sehr unwahrscheinlich und in physischen Systemen allgegenwärtig.

MRocklin
quelle
2
Wenn das lineare System anfangs schlecht konditioniert ist, spielt es keine Rolle, welche Methode Sie verwenden: Sowohl die LU- als auch die QR-Zerlegung führen zu ungenauen Ergebnissen. QR kann nur gewinnen, wenn der Prozess der Gaußschen Eliminierung eine gute Matrix "verdirbt". Das Hauptproblem ist, dass praktische Fälle eines solchen Verhaltens nicht bekannt sind.
Faleichik
Für die meisten wissenschaftlichen Anwendungen erhalten wir im Allgemeinen spärliche, symmetrische, positiv definite und / oder diagonal dominante Matrizen. Mit wenigen Ausnahmen gibt es eine Struktur in der Matrix, die es uns ermöglicht, bestimmte Techniken gegenüber der traditionellen Gaußschen Eliminierung auszunutzen.
Paul
@Paul: Auf der anderen Seite wird bei der multifrontalen Methode für spärliche unsymmetrische Matrizen die meiste Zeit für die dichte Gaußsche Eliminierung aufgewendet.
Jack Poulson
6
@Paul Es ist einfach nicht wahr, dass "die meisten Anwendungen SPD / diagonal dominante Matrizen erzeugen". Ja, es gibt normalerweise eine Art ausnutzbare Struktur, aber unsymmetrische und unbestimmte Probleme sind äußerst häufig.
Jed Brown
4
"In 50 Jahren Computerarbeit sind unter natürlichen Umständen keine Matrixprobleme aufgetreten, die eine explosive Instabilität hervorrufen." - LN Trefethen und D. Bau In ihrem Buch geben sie eine interessante probabilistische Analyse.
JM