Wie man Matlab's Mldivide (auch bekannt als Backslash-Operator "\") implementiert

Question

Ich versuche gerade eine kleine Matrix-orientierte Math-Bibliothek zu entwickeln (ich verwende Eigen 3 für Matrix Datenstrukturen und Operationen) und ich wollte einige implementieren Praktische Matlab-Funktionen, wie der weit verbreitete Backslash-Operator (entspricht mldivide), um die Lösung linearer Systeme zu berechnen (ausgedrückt in Matrixform).

Gibt es eine gute detaillierte Erklärung, wie dies erreicht werden könnte? (Ich habe implementiert bereits die pseudoinverse pinv Funktion mit einer klassischen SVD Zersetzung, aber ich habe irgendwo gelesen, dass A\b ist nicht immer pinv(A)*b, zumindest Matalb nicht einfach tun)

Dank

Answer 1

Für x = A\b umfasst der Operator backslash eine Anzahl von algorithms, um verschiedene Arten von Eingangsmatrizen zu behandeln. So wird die Matrix A diagnostiziert und ein Ausführungspfad wird entsprechend seiner Eigenschaften ausgewählt.

Die folgenden page beschreibt in Pseudo-Code, wenn A eine vollständige Matrix ist:

if size(A,1) == size(A,2)   % A is square 
    if isequal(A,tril(A))   % A is lower triangular 
     x = A \ b;    % This is a simple forward substitution on b 
    elseif isequal(A,triu(A))  % A is upper triangular 
     x = A \ b;    % This is a simple backward substitution on b 
    else 
     if isequal(A,A')   % A is symmetric 
      [R,p] = chol(A); 
      if (p == 0)   % A is symmetric positive definite 
       x = R \ (R' \ b); % a forward and a backward substitution 
       return 
      end 
     end 
     [L,U,P] = lu(A);   % general, square A 
     x = U \ (L \ (P*b));  % a forward and a backward substitution 
    end 
else        % A is rectangular 
    [Q,R] = qr(A); 
    x = R \ (Q' * b); 
end 

Bei nicht quadratischen Matrizen, QR decomposition verwendet wird. Für quadratische Dreiecksmatrizen führt es eine einfache forward/backward substitution aus. Für quadratische symmetrische positiv definite Matrizen wird Cholesky decomposition verwendet. Ansonsten wird LU decomposition für allgemeine quadratische Matrizen verwendet.

Update: MathWorks die algorithm section in der doc Seite von mldivide mit einigen netten Flussdiagramme aktualisiert hat. Siehe here und here (vollständige und spärliche Fälle).

Alle diese Algorithmen haben entsprechende Verfahren in LAPACK, und es ist wahrscheinlich das, was MATLAB tut (beachten Sie, dass neuere Versionen von MATLAB Schiff mit der optimierten Intel MKL Implementierung) in der Tat. Der Grund dafür, verschiedene Methoden zu haben, besteht darin, dass man versucht, den Algorithmus mit dem spezifischsten zu verwenden, um das System der Gleichungen zu lösen, das alle Eigenschaften der Koeffizientenmatrix ausnutzt (entweder weil es schneller oder numerisch stabiler wäre). Man könnte also durchaus einen allgemeinen Löser verwenden, aber es ist nicht der effizienteste.

In der Tat, wenn Sie wissen, wie A im Voraus ist, können Sie den zusätzlichen Testprozess überspringen, indem Sie linsolve aufrufen und die Optionen direkt angeben.

wenn A rechteckig oder Singular, auch PINV verwenden könnte eine minimale Norm Least-Squares-Lösung (implementiert SVD decomposition) zu finden:

x = pinv(A)*b 

---

Alle oben genannten dichten Matrizen gilt, spärliche Matrizen sind eine ganz andere Geschichte. Normalerweise werden in solchen Fällen iterative solvers verwendet. Ich glaube, MATLAB verwendet UMFPACK und andere verwandte Bibliotheken aus dem SuiteSpase-Paket für direkte Löser.

Wenn mit spärlichen Matrizen arbeiten, können Sie auf Diagnoseinformationen wenden und die Tests durchgeführt und Algorithmen sehen spparms gewählt werden:

spparms('spumoni',2) 
x = A\b; 

---

Was mehr ist, auch der Backslash Operator auf gpuArray ‚s arbeitet, In diesem Fall verlässt man sich auf cuBLAS und MAGMA, um auf der GPU auszuführen.

Es ist auch für distributed arrays implementiert, die in einer verteilten Computerumgebung funktioniert (Arbeit unter einem Cluster von Computern aufgeteilt, wo jeder Worker nur einen Teil des Arrays hat, möglicherweise wo die gesamte Matrix nicht im Speicher auf einmal gespeichert werden kann). Die zugrunde liegende Implementierung verwendet ScaLAPACK.

Das ist eine ziemlich große Herausforderung ist, wenn Sie all das selbst :)