Jak operator odwrotnego ukośnika MATLAB rozwiązuje dla macierzy kwadratowych?

Porównywałem kilka moich kodów z „zapasowymi” kodami MATLAB. Jestem zaskoczony wynikami.

Uruchomiłem przykładowy kod (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;

Wyniki:

    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.

W przypadku gęstej matrycy:

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;

Wyniki:

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.

Jak do cholery jest tak niesamowicie?

W Matlabie komenda „\” wywołuje algorytm, który zależy od struktury macierzy A i obejmuje sprawdzenie (mały narzut) właściwości A.

1. Jeśli A jest rzadki i zwinięty, zastosuj solver zwężający się.
2. Jeśli A jest górną lub dolną trójkątną matrycą, zastosuj algorytm podstawiania wstecznego.
3. Jeśli A jest symetryczny i ma rzeczywiste dodatnie elementy ukośne, spróbuj dokonać faktoryzacji Choleskiego. Jeśli A jest rzadkie, najpierw zastosuj zmianę kolejności, aby zminimalizować wypełnienie.
4. Jeśli żadne z powyższych kryteriów nie jest spełnione, wykonaj ogólną trójkątną faktoryzację, stosując eliminację Gaussa z częściowym przestawieniem.
5. Jeśli A jest rzadkie, zastosuj bibliotekę UMFPACK.
6. Jeśli A nie jest kwadratowe, zastosuj algorytmy oparte na faktoryzacji QR dla nieokreślonych systemów.

Aby zmniejszyć koszty ogólne, można użyć polecenia linsolve w Matlabie i wybrać odpowiedni solver spośród tych opcji.

Jeśli chcesz zobaczyć, co a\brobi dla twojej konkretnej matrycy, możesz ustawić spparms('spumoni',1)i dokładnie określić, jakim algorytmem byłeś pod wrażeniem. Na przykład:

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

wyjdzie

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.

więc widzę, że w tym przypadku „\” użyło „CHOLMOD”.

W przypadku rzadkich macierzy Matlab używa UMFPACK do operacji „ \”, która w twoim przykładzie zasadniczo przegląda wartości a, odwraca je i mnoży je przez wartości b. Jednak w tym przykładzie powinieneś użyć b./diag(a), który jest znacznie szybszy.

W przypadku gęstych systemów operator odwrotnego ukośnika jest nieco bardziej skomplikowany. Krótki opis tego, co się dzieje, gdy jest podany tutaj . Zgodnie z tym opisem, w twoim przykładzie Matlab rozwiązałby a\bza pomocą podstawienia wstecznego. W przypadku ogólnych macierzy kwadratowych stosuje się rozkład LU.

Zasadniczo, jeśli masz rzadką matrycę o rozsądnej złożoności (tj. Nie musi to być szablon 5-punktowy, ale w rzeczywistości może być dyskretyzacją równań Stokesa, dla których liczba niezerowych w rzędzie wynosi znacznie większy niż 5), wówczas rzadki bezpośredni solver, taki jak UMFPACK, zazwyczaj pokonuje iteracyjny solver Kryłowa, jeśli problem nie jest większy niż około 100 000 niewiadomych.

Innymi słowy, dla większości rzadkich macierzy wynikających z dyskretyzacji 2D, najszybszym rozwiązaniem jest bezpośredni solver. Tylko w przypadku problemów 3D, w których szybko dostajesz ponad 100 000 niewiadomych, konieczne jest użycie iteracyjnego solvera.