Wie löst der MATLAB-Backslash-Operator

36

Ich habe einige meiner Codes mit MATLAB-Codes "auf Lager" verglichen. Ich bin überrascht über das Ergebnis.

Ich habe einen Beispielcode ausgeführt (Sparse Matrix)

n = 5000;
a = diag(rand(n,1));
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('Inv(A)*B');
tic;inv(a)*b;toc;

Ergebnisse :

    For a\b
    Elapsed time is 0.052838 seconds.

    For LU
    Elapsed time is 7.441331 seconds.

    For Conj Grad
    Elapsed time is 3.819182 seconds.

    Inv(A)*B
    Elapsed time is 38.511110 seconds.

Für dichte Matrix:

n = 2000;
a = rand(n,n);
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('For INV(A)*B');
tic;inv(a)*b;toc;

Ergebnisse:

For a\b
Elapsed time is 0.575926 seconds.

For LU
Elapsed time is 0.654287 seconds.

For Conj Grad
Elapsed time is 9.875896 seconds.

Inv(A)*B
Elapsed time is 1.648074 seconds.

Wie zum Teufel ist a \ b so toll?

Anfrage
quelle
1
Der eingebaute Backslash von MATLAB, mit anderen Worten der direkte Löser für ein lineares Gleichungssystem, verwendet die Multifrontal-Methode für eine spärliche Matrix, weshalb A \ B so großartig ist.
Shuhao Cao
1
Es wird der Code von Tim Davis verwendet, der unter cise.ufl.edu/research/sparse verfügbar ist . Auch die Attraktivität verschwindet, wenn Sie ein nicht-triviales Problem haben
stali
1
Was ist "LULU"? Warum ist es Ihrer Meinung nach eine gute Implementierung einer LU-Faktorisierung und einer anschließenden direkten Lösung?
Jed Brown
3
@Nunoxic Welche Implementierung? Hast du es selbst geschrieben? Hochleistungsdichte lineare Algebra ist zwar normalerweise algorithmisch gut verstanden, aber auf moderner Hardware nicht einfach effizient zu implementieren. Die besten BLAS / Lapack-Implementierungen sollten für eine Matrix dieser Größe nahe am Peak sein. Außerdem habe ich aus Ihren Kommentaren den Eindruck, dass LU und Gaußsche Elimination unterschiedliche Algorithmen sind.
Jed Brown
1
Es ruft einen Fortran-Code auf, der mit Intel MKL geschrieben wurde.
Anfrage

Antworten:

37

In Matlab ruft der Befehl '\' einen Algorithmus auf, der von der Struktur der Matrix A abhängt und Überprüfungen (geringer Overhead) der Eigenschaften von A enthält.

  1. Wenn A dünn und gebändert ist, verwenden Sie einen gebänderten Solver.
  2. Wenn A eine obere oder untere Dreiecksmatrix ist, verwenden Sie einen Rückwärtssubstitutionsalgorithmus.
  3. Wenn A symmetrisch ist und echte positive diagonale Elemente aufweist, versuchen Sie eine Cholesky-Faktorisierung. Wenn A spärlich ist, führen Sie zuerst eine Neuordnung durch, um das Ausfüllen zu minimieren.
  4. Wenn keines der oben genannten Kriterien erfüllt ist, führen Sie eine allgemeine Dreiecksfaktorisierung unter Verwendung der Gaußschen Eliminierung mit teilweisem Schwenken durch.
  5. Wenn A dünn ist, verwenden Sie die UMFPACK-Bibliothek.
  6. Wenn A nicht quadratisch ist, wenden Sie Algorithmen an, die auf der QR-Faktorisierung für unbestimmte Systeme basieren.

Um den Overhead zu reduzieren, können Sie den Befehl linsolve in Matlab verwenden und selbst einen geeigneten Löser aus diesen Optionen auswählen.

Allan P. Engsig-Karup
quelle
Angenommen, ich habe es mit einer unstrukturierten dichten Matrix von 10000 x 10000 mit allen Elementen ungleich Null (hohe Dichte) zu tun. Was wäre meine beste Wahl? Ich möchte diesen 1-Algorithmus isolieren, der für dichte Matrizen funktioniert. Ist es LU, QR oder Gaußsche Eliminierung?
Anfrage
1
Klingt nach einem Schritt 4, bei dem die Gaußsche Eliminierung aufgerufen wird, was dem allgemeinsten Fall entspricht, bei dem keine Struktur von A zur Leistungssteigerung ausgenutzt werden kann. Dies ist also im Grunde genommen eine LU-Faktorisierung und eine nachfolgende Vorwärts-Faktorisierung, gefolgt von einem Rückwärts-Substitutionsschritt.
Allan P. Engsig-Karup
Vielen Dank! Ich denke, das gibt mir eine Richtung zum Nachdenken. Derzeit ist die Gaußsche Eliminierung das Beste, was wir zur Lösung derartiger unstrukturierter Probleme haben. Ist das richtig?
Anfrage
37

Wenn Sie sehen möchten, was a\bfür Ihre spezielle Matrix funktioniert, können Sie genau festlegen spparms('spumoni',1)und herausfinden, von welchem ​​Algorithmus Sie beeindruckt waren. Beispielsweise:

spparms('spumoni',1);
A = delsq(numgrid('B',256));
b = rand(size(A,2),1);
mldivide(A,b);  % another way to write A\b

wird ausgegeben

sp\: bandwidth = 254+1+254.
sp\: is A diagonal? no.
sp\: is band density (0.01) > bandden (0.50) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? no.
sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? yes.
sp\: is CHOLMOD's symbolic Cholesky factorization (with automatic reordering) successful? yes.
sp\: is CHOLMOD's numeric Cholesky factorization successful? yes.
sp\: is CHOLMOD's triangular solve successful? yes.

Ich kann also sehen, dass "\" in diesem Fall "CHOLMOD" verwendet.

dranxo
quelle
3
+1 für neue MATLAB-Einstellungen, von denen ich noch nie gehört habe. Gut gemacht, Sir.
Geoff Oxberry
2
Hey danke! Es ist in help mldivide.
Dranxo
16

Für dünne Matrizen verwendet Matlab UMFPACK für die "\ Operation ", die in Ihrem Beispiel im Wesentlichen die Werte von durchläuft a, sie invertiert und mit den Werten von multipliziert b. Für dieses Beispiel sollten Sie jedoch verwendenb./diag(a) , was viel schneller ist.

Bei dichten Systemen ist der Backslash-Operator etwas komplizierter. Eine kurze Beschreibung dessen, was wann getan wird, finden Sie hier . Gemäß dieser Beschreibung würde Matlab in Ihrem Beispiel a\bdurch Rückwärtssubstitution lösen . Für allgemeine Quadratmatrizen wird eine LU-Zerlegung verwendet.

Pedro
quelle
Regd. Sparsamkeit, das inv einer Diag-Matrix wäre einfach der Kehrwert der diagonalen Elemente, so dass b./diag(a) funktionieren würde, aber a \ b funktioniert auch hervorragend für allgemeine spärliche Matrizen. Warum ist linsolve oder LULU (Meine optimierte Version von LU) in diesem Fall nicht schneller als a \ b für dichte Matrizen?
Anfrage
@Nunoxic Tut Ihr LULU etwas, um die Diagonale oder Dreieckigkeit der dichten Matrix zu erkennen? Wenn nicht, dauert es mit jeder Matrix genauso lange, unabhängig von ihrem Inhalt oder ihrer Struktur.
Pedro
Etwas. Mit den linsolve-OPT-Flags habe ich jedoch alles definiert, was an der Struktur zu definieren ist. Dennoch ist a \ b schneller.
Anfrage
@Nunoxic, jedes Mal, wenn Sie eine Benutzerfunktion aufrufen, verursachen Sie Overheads. Matlab erledigt alles für den Backslash intern, z. B. die Zerlegung und anschließende Anwendung der rechten Seite, mit sehr geringem Overhead und ist somit schneller. Außerdem sollten Sie in Ihren Tests mehr als nur einen Anruf verwenden, um zuverlässige Timings zu erhalten, z tic; for k=1:100, a\b; end; toc.
Pedro
5

Als Faustregel gilt, wenn Sie eine spärliche Matrix mit angemessener Komplexität haben (dh, es muss sich nicht um die 5-Punkt-Schablone handeln, sondern es kann sich tatsächlich um eine Diskretisierung der Stokes-Gleichungen handeln, für die die Anzahl der Nichtzeros pro Zeile gilt viel größer als 5), dann schlägt ein spärlicher direkter Löser wie UMFPACK normalerweise einen iterativen Krylov-Löser, wenn das Problem nicht größer als etwa 100.000 Unbekannte ist.

Mit anderen Worten, für die meisten spärlichen Matrizen, die aus 2D-Diskretisierungen resultieren, ist ein direkter Löser die schnellste Alternative. Nur bei 3D-Problemen, bei denen Sie schnell über 100.000 Unbekannte erreichen, muss ein iterativer Löser verwendet werden.

Wolfgang Bangerth
quelle
3
Mir ist nicht klar, wie dies die Frage beantwortet, aber ich habe auch Probleme mit der Prämisse. Es ist wahr, dass Direktlöser bei 2D-Problemen mit bescheidener Größe in der Regel gut funktionieren, während Direktlöser in 3D in der Regel weit vor 100.000 Unbekannten leiden (die Scheitelpunkttrennzeichen sind viel größer als in 2D). Darüber hinaus behaupte ich, dass iterative Löser in den meisten Fällen (z. B. elliptische Operatoren) auch für mittelgroße 2D-Probleme schneller als direkte Löser erstellt werden können, dies kann jedoch erhebliche Anstrengungen erfordern (z. B. Verwendung der richtigen Zutaten, um die Leistung von Multigrid zu verbessern). .
Jed Brown
1
Wenn Ihre Probleme relativ klein sind und Sie nicht über das Entwerfen impliziter Löser nachdenken möchten, oder wenn Ihre Probleme sehr schwierig sind (wie Hochfrequenz-Maxwell) und Sie Ihre Karriere nicht dem Entwerfen guter Löser widmen möchten, dann ich stimmen zu, dass spärliche direkte Löser eine gute Wahl sind.
Jed Brown
Eigentlich besteht mein Problem darin, einen guten dichten Löser zu entwerfen. Ich habe keine bestimmte Anwendung als solche (daher keine Struktur). Ich wollte einige meiner Codes optimieren und sie mit den derzeit verwendeten konkurrenzfähig machen. Ich war nur neugierig auf a \ b und wie es immer den Mist aus meinem Code schlägt.
Anfrage
@JedBrown: Ja, vielleicht kann man mit erheblichem Aufwand iterative Löser dazu bringen, einen direkten Löser sogar für kleine 2D-Probleme zu schlagen. Aber warum? Bei Problemen mit <100.000 Unbekannten sind spärliche Direktlöser in 2d schnell genug. Was ich sagen wollte ist: Verbringen Sie bei so kleinen Problemen nicht Ihre Zeit damit, die beste Kombination von Parametern zu basteln und herauszufinden - nehmen Sie einfach die Blackbox.
Wolfgang Bangerth
Es gibt bereits eine Größenordnung Laufzeitunterschied zwischen einem spärlichen Cholesky- und einem (matrixbasierten) geometrischen Multigrid für "einfache" 2D-Probleme mit 100.000 Dofs unter Verwendung einer 5-Punkt-Schablone (~ 0,5 Sekunden im Vergleich zu ~ 0,05 Sekunden). Wenn Ihre Schablone zweite Nachbarn verwendet (z. B. Diskretisierung vierter Ordnung, Newton für eine Auswahl an nichtlinearer Rheologie, Stabilisierung usw.), verdoppelt sich die Größe der Scheitelpunkttrenner in etwa, sodass sich die Kosten für die direkte Lösung um das 8-fache erhöhen (die Kosten sind höher) problemabhängig für MG). Bei vielen Zeitschritten oder Optimierungs- / UQ-Schleifen sind diese Unterschiede signifikant.
Jed Brown