Wielokrotnie rozwiązując

12

Korzystam z MATLAB, aby rozwiązać problem polegający na rozwiązywaniu za każdym razem, gdy b zmienia się z czasem. Obecnie realizuję to za pomocą MATLAB-a :Ax=bbmldivide

x = A\b

Mam elastyczność, aby wykonać tyle wstępnych obliczeń, ile potrzeba, więc zastanawiam się, czy istnieje szybsza i / lub dokładniejsza metoda niż mldivide. Co zazwyczaj się tutaj robi? Dziękuje wszystkim!

Wątpić
źródło
1
Czy masz konkretną wiedzę na temat struktury ? Na przykład, czy jest symetryczny? Pozytywne określone? Tridiagonal? Prostokątny? A
Dominique
Macierz jest gęstą macierzą kwadratową. A
Wątpliwości
3
Jeśli nie masz innej wiedzy na temat , faktoryzacja L U, jak opisano w odpowiedzi poniżej, jest najlepszym wyborem. ALU
Dominique

Odpowiedzi:

14

Najbardziej oczywistą rzeczą, jaką możesz zrobić, jest wstępne obliczenie

[L,U] = lu(A) ~ O (n ^ 3)

Następnie po prostu oblicz

x = U \ (L \ b) ~ O (2 n ^ 2)

To ogromnie obniżyłoby koszt i przyspieszyło. Dokładność byłaby taka sama.

Milind R.
źródło
1
Zauważ, że z dokumentacji L nie musi być niższy trójkątny. Ta odpowiedź byłaby prawdopodobnie szybsza niż bezpośrednie rozwiązanie, jednak uważałbym, aby upewnić się, że polecenie L \ b jest wystarczająco inteligentne, aby wiedzieć, że rozwiązywać L w prawidłowej kolejności (prawdopodobnie jest, ale nie jest pewne w dokumentacji).
Godric Seer,
Tak, masz rację, L jest iloczynem macierzy o trójkątnym kształcie dolnym i macierzy permutacji. Ale niech mnie diabli, jeśli nie rozpozna, że ​​jedyne, co musi zrobić, to zastąpić go wstecz L\b. Ponieważ widziałem tę dokładną linię w kodzie o wysokiej wydajności przez tych, których uważam za ekspertów.
Milind R
8
O(n2)
1
A
3
@BrianBorcher O ile mi wiadomo, najlepszym sposobem na śledzenie permutacji jest [L,U,p] = lu(A,'vector'); x = U\(L\b(p));Patrz przykład 3 w lu dokumentacji .
Stefano M
5

Przeprowadziliśmy rozległe laboratoria komputerowe na naszych naukowych kursach komputerowych na ten temat. W przypadku „małych” obliczeń, które tam przeprowadziliśmy, operator odwrotnego ukośnika Matlaba był zawsze szybszy niż cokolwiek innego, nawet po tym, jak zoptymalizowaliśmy nasz kod tak bardzo, jak to możliwe i ponownie uporządkowaliśmy wszystkie macierze (na przykład przy zamawianiu rzadkich matryc przez Reverse Cuthill McKee) .

Możesz sprawdzić jedną z naszych instrukcji laboratoryjnych . Odpowiedź na twoje pytanie znajduje się (wkrótce) na stronie 4.

Dobra książka na ten temat została napisana na przykład przez Cheneya .

seb
źródło
4

An×n Axi=bii=1mm

V = inv(A);
...
x = V*b;

O(n3)inv(A)O(n2)V*bm

>> n = 5000;
>> A = randn(n,n);
>> x = randn(n,1);
>> b = A*x;
>> rcond(A)
ans =
   1.3837e-06
>> tic, xm = A\b; toc
Elapsed time is 1.907102 seconds.
>> tic, [L,U] = lu(A); toc
Elapsed time is 1.818247 seconds.
>> tic, xl = U\(L\b); toc
Elapsed time is 0.399051 seconds.
>> tic, [L,U,p] = lu(A,'vector'); toc
Elapsed time is 1.581756 seconds.
>> tic, xp = U\(L\b(p)); toc
Elapsed time is 0.060203 seconds.
>> tic, V=inv(A); toc
Elapsed time is 7.614582 seconds.
>> tic, xv = V*b; toc     
Elapsed time is 0.011499 seconds.
>> [norm(xm-x), norm(xp-x), norm(xl-x), norm(xv-x)] ./ norm(x)
ans =
   1.0e-11 *
    0.1912    0.1912    0.1912    0.6183

A1LUm>125

Kilka notatek

Dla analizy stabilności i błędów zobacz komentarze do tej innej odpowiedzi , szczególnie tej autorstwa VictorLiu.

mn

Pomiar czasu przeprowadzono za pomocą Matlaba R2011b na 12-rdzeniowym komputerze ze dość stałą średnią wartością obciążenia UNIX wynoszącą 5; najlepszy tic, tocczas z trzech sond.

Stefano M.
źródło
Rzeczywiście, istnieje wiele więcej równoległości w mnożeniu wektora macierzy niż trójkątny solver, więc powinno to być jeszcze bardziej widoczne, jeśli obliczenia są wykonywane równolegle (wielordzeniowy / GPU / itp.) W jakikolwiek sposób.
Aron Ahmadia
@AronAhmadia Zgadzam się: szacunki progu rentowności oparte tylko na liczbie operacji mają sens tylko w przypadku realizacji szeregowej.
Stefano M
1
Zauważ, że rzeczy będą się znacznie różnić, jeśli macierz A będzie rzadka - odwrotność zwykle będzie dość gęsta, podczas gdy czynniki LU są zwykle dość rzadkie, przechylając rzeczy do tyłu w kierunku LU.
Brian Borchers,
1
A
1
inv(A)Ax=bbBA\B
2

Spójrz na to pytanie , odpowiedzi pokazują, że mldividejest dość sprytny, a także daje sugestie, jak zobaczyć, co Matlab używa do rozwiązania A\b. Może to dać podpowiedź dotyczącą opcji optymalizacji.

Psirus
źródło
0

Użycie odwrotnego ukośnika jest mniej więcej równoważne inv(A)*B, jeśli swobodnie go kodujesz, ten drugi może być bardziej intuicyjny. Są mniej więcej takie same (po prostu różnią się sposobem przeprowadzania obliczeń), chociaż powinieneś sprawdzić dokumentację Matlaba w celu uzyskania wyjaśnień.

Aby odpowiedzieć na twoje pytanie, odwrotny ukośnik jest ogólnie w porządku, ale zależy to od właściwości macierzy masy.

AlanH
źródło
1
Matematyczne inv (A) * b jest takie samo jak \ jednak liczbowo, w rzeczywistości tworzenie odwrotności jest zarówno mniej wydajne, jak i mniej dokładne. Jeśli pracujesz nad nauką algebry liniowej, może to być do przyjęcia, ale argumentowałbym, że potrzebujesz bardzo dobrego powodu, aby utworzyć odwrotność.
Godric Seer,
Ale dlaczego miałbyś kiedykolwiek obliczać, inv(A)skoro to samo jest droższe niż A\b?
Dominique
7
@Godric: Niedawny artykuł omawia „mit”, że inv (A) * b jest mniej dokładny: na ArXiv . Nie mówię, że zwykle jest powód, aby obliczyć rzeczywistą odwrotność, ale po prostu mówię.
Victor Liu
3
@Dominique: Rozwiązania trójkątne są znacznie mniej równoległe niż mnożenie macierzy i wektora, a wyrafinowane, wstępnie przygotowane metody iteracyjne często wykorzystują metody bezpośrednie w poddomenach. Często przydatne jest jednoznaczne uformowanie odwrotności kilku gęstych trójkątnych matryc o niewielkich rozmiarach w celu poprawy równoległości.
Jack Poulson
@VictorLiu: Dziękuję za artykuł. Stoję poprawiony na moim oświadczeniu o dokładności (przynajmniej dla inteligentnych implementacji inv (A)).
Godric Seer