Jak wdrożyć Matusbowi mldivide (znany również jako operator odwrotnego ukośnika "\")

Question

Obecnie próbuję stworzyć małą bibliotekę matematyczną (używam Eigen 3 dla struktur danych macierzy i operacji) i chciałem zaimplementować niektóre poręczne funkcje Matlaba, takie jak szeroko stosowany operator ukośnika odwrotnego (który jest odpowiednikiem mldivide) w celu obliczenia rozwiązania układów liniowych (wyrażonych w postaci macierzowej).

Czy są jakieś dobre szczegółowe wyjaśnienia, w jaki sposób można to osiągnąć? (Ja już wdrożył uogólniona macierz odwrotna pinv funkcji z klasycznym rozkładem SVD, ale czytałem gdzieś, że nie zawsze jest A\bpinv(A)*b przynajmniej Matalb nie robi po prostu to zrobić)

Dzięki

Answer 1

Dla x = A\b operator backslash obejmuje liczbę algorithms do obsługi różnych rodzajów matryc wejściowych. Tak więc diagnozuje się matrycę A i wybiera się ścieżkę wykonania zgodnie z jej charakterystyką.

Poniżej page opisano w pseudokodzie gdy A jest pełna macierz:

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 

dla macierzy nie kwadratowych QR decomposition używany. W przypadku kwadratowych macierzy trójkątnych wykonuje prostą forward/backward substitution. W przypadku kwadratowych, symetrycznych macierzy dodatnich o określonym numerze, stosuje się Cholesky decomposition. W przeciwnym razie LU decomposition jest używany do ogólnych kwadratowych matryc.

Aktualizacja: MathWorks zaktualizowała algorithm section na stronie doc mldivide nice niektórych schematów. Zobacz here i here (pełne i rzadkie przypadki).

Wszystkie te algorytmy mają odpowiednie metody w LAPACK, aw rzeczywistości jest to prawdopodobnie co robi MATLAB (zauważ, że najnowsze wersje MATLAB statku ze zoptymalizowaną Intel MKL realizacji).

Powodem różnych metod jest to, że stara się użyć najbardziej konkretnego algorytmu do rozwiązania układu równań, który wykorzystuje wszystkie cechy macierzy współczynników (albo dlatego, że byłby szybszy, albo bardziej stabilny liczbowo). Więc na pewno skorzystasz z ogólnego rozwiązania, ale nie będzie to najbardziej efektywne.

W rzeczywistości, jeśli wiesz wcześniej, jak wygląda A, możesz pominąć dodatkowy proces testowania, dzwoniąc pod numer linsolve i określając opcje bezpośrednio.

jeśli A jest prostokątny lub pojedynczej, można też użyć PINV znaleźć minimalną normę rozwiązanie najmniejszych kwadratów (realizowane z wykorzystaniem SVD decomposition):

x = pinv(A)*b 

---

Wszystkie powyższe stosuje się do gęstych matryc, rzadkie macierze to zupełnie inna historia. Zwykle w takich przypadkach stosuje się iterative solvers. Wierzę, że MATLAB używa UMFPACK i innych powiązanych bibliotek z pakietu SuiteSpase dla bezpośrednich solverów.

Podczas pracy z nielicznych macierzy, można włączyć informacje diagnostyczne i zobaczyć Przeprowadzone testy i algorytmy wybranych użyciu spparms:

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

---

Co więcej, operator odwrotny ukośnik działa również na gpuArray „s, w takim przypadku opiera się na cuBLAS i MAGMA, aby wykonać na GPU.

Jest również zaimplementowany na distributed arrays, który działa w środowisku przetwarzania rozproszonego (praca podzielona między klastry komputerów, gdzie każdy pracownik ma tylko część tablicy, prawdopodobnie tam, gdzie cała macierz nie może być przechowywana w pamięci jednocześnie). Podstawowa implementacja używa ScaLAPACK.

To bardzo wysokie zamówienie, jeśli chcesz to wszystko wdrożyć samemu :)