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!


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.


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 .


4

An×n Axi=bii=1…mm

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

A−1LUm>125

Kilka notatek

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

m≪n

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.


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.


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.


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
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.