Come implementare la mldivide di Matlab (ovvero l'operatore backslash "\")

Question

Attualmente sto cercando di sviluppare una piccola libreria matematica orientata alla matrice (sto usando Eigen 3 per le strutture e le operazioni di dati matriciali) e volevo implementare alcuni pratiche funzioni di Matlab, come l'operatore backslash ampiamente utilizzato (che è equivalente a mldivide) per calcolare la soluzione di sistemi lineari (espressi in forma di matrice).

Esiste una buona spiegazione dettagliata su come ciò potrebbe essere realizzato? (Ho già implementato la pseudo-inversa pinv funzione con una decomposizione classica SVD, ma ho letto da qualche parte che A\b non è sempre pinv(A)*b, almeno Matalb non si limita a farlo)

Grazie

Answer 1

Per x = A\b, l'operatore backslash comprende un numero di algorithms per gestire diversi tipi di matrici di input. Quindi viene diagnosticata la matrice A e viene selezionato un percorso di esecuzione in base alle sue caratteristiche.

Il seguente page descrive in pseudo-codice quando A è una matrice completa:

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 

Per matrici non quadrati, QR decomposition viene utilizzato. Per matrici triangolari quadrate, esegue un semplice forward/backward substitution. Per le matrici simmetriche con segno positivo definito, viene utilizzato Cholesky decomposition. Altrimenti LU decomposition viene utilizzato per le matrici quadrate generali.

Aggiornamento: MathWorks ha aggiornato il algorithm section nella pagina doc di mldivide con alcuni bei diagrammi di flusso. Vedere here e here (casi pieni e sparsi).

Tutti questi algoritmi sono metodi corrispondenti a LAPACK, e in effetti è probabilmente quello che sta facendo MATLAB (si noti che le versioni recenti di MATLAB nave con l'ottimizzato Intel MKL implementazione).

La ragione per avere metodi diversi è che tenta di utilizzare l'algoritmo più specifico per risolvere il sistema di equazioni che sfrutta tutte le caratteristiche della matrice dei coefficienti (o perché sarebbe più veloce o più numericamente stabile). Quindi potresti certamente usare un risolutore generale, ma non sarà il più efficiente.

Infatti, se si conosce come A è come prima, è possibile saltare il processo di test extra chiamando linsolve e specificando direttamente le opzioni.

se A è rettangolare o singolari, si potrebbe anche usare PINV per trovare una norma minima soluzione di minimi quadrati (implementata utilizzando SVD decomposition):

x = pinv(A)*b 

---

Tutto quanto sopra vale per le matrici dense, matrici sparse sono una storia completamente diversa. Solitamente iterative solvers vengono utilizzati in questi casi. Credo che MATLAB usi UMFPACK e altre librerie correlate dal pacchetto SuiteSpase per i risolutori diretti.

Quando si lavora con matrici sparse, è possibile attivare le informazioni di diagnostica e di vedere i test eseguiti e gli algoritmi scelti utilizzando spparms:

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

---

Cosa c'è di più, l'operatore backslash funziona anche su gpuArray 's, nel qual caso si basa su cuBLAS e MAGMA da eseguire sulla GPU.

E 'anche implementato per distributed arrays che funziona in un ambiente di elaborazione distribuito (lavoro diviso in un cluster di computer in cui ogni lavoratore ha solo una parte dell'array, probabilmente dove l'intera matrice non può essere memorizzata in memoria tutto in una volta). L'implementazione sottostante utilizza ScaLAPACK.

Questo è un ordine abbastanza alto se si desidera implementare tutto questo da soli :)