Kiedy transformacje ortogonalne przewyższają eliminację Gaussa?

22

Jak wiemy, metody transformacji ortogonalnej (obroty Givensa i odbicia Housholder'a) dla układów równań liniowych są droższe niż eliminacja Gaussa, ale teoretycznie mają lepsze właściwości stabilności w tym sensie, że nie zmieniają numeru stanu układu. Chociaż znam tylko jeden akademicki przykład matrycy, która psuje Gaussowska eliminacja z częściowym przestawieniem. Istnieje powszechna opinia, że jest bardzo mało prawdopodobne, aby w praktyce spełnić tego rodzaju zachowania (patrz uwagi do wykładu [pdf] ).

Gdzie więc szukać odpowiedzi na ten temat? Równoległe wdrożenia? Aktualizujesz? ..

linear-algebra reference-request

— faleichik
źródło

24

Precyzja

Trefethen i Schreiber napisali znakomity artykuł, „ Stabilność eliminacji Gaussa w średnich przypadkach” , w którym omawiamy dokładność pytania. Oto kilka wniosków:

„W przypadku faktoryzacji QR z obrotem kolumny lub bez niej, średni maksymalny element macierzy resztkowej wynosi , podczas gdy dla eliminacji Gaussa jest to . Porównanie to pokazuje, że eliminacja Gaussa jest nieznacznie niestabilna , ale niestabilność byłaby wykrywalna tylko w przypadku bardzo dużych problemów z matrycą rozwiązanych z małą precyzją. W przypadku większości problemów praktycznych eliminacja Gaussa jest średnio bardzo stabilna. ”(moje podkreślenie) $O(n^{1/2})$ $O(n)$
„Po kilku pierwszych krokach eliminacji Gaussa pozostałe elementy macierzy są w przybliżeniu normalnie rozmieszczone, niezależnie od tego, czy rozpoczęły się w ten sposób”.

W tym dokumencie jest o wiele więcej, czego nie mogę uchwycić, w tym dyskusja na temat matrycy najgorszego przypadku, o której wspomniałeś, więc zdecydowanie zalecamy jej przeczytanie.

Wydajność

W przypadku rzeczywistych macierzy kwadratowych LU z częściowym przestawieniem wymaga około klap, podczas gdy QR oparty na domownikach wymaga około klap. Tak więc, w przypadku dość dużych macierzy kwadratowych, rozkład na czynniki QR będzie tylko około dwa razy droższy niż rozkład na współczynniki LU. $2/3 n^3$ $4/3 n^3$

Dla macierzy , gdzie , LU z częściowym przestawieniem wymaga 3/3 przerzutów, w porównaniu do QR (co jest nadal dwukrotnie wyższe niż rozkład na czynniki LU) . Jednak zaskakująco często aplikacje wytwarzają bardzo wysokie chude matryce ( ), a Demmel i in. mają ładny artykuł, Unikanie komunikacji równoległej i sekwencyjnej faktoryzacji QR , która omawia sprytny algorytm, który wymaga tylko wysyłania wiadomości, gdy używane są procesory , w porównaniu z wiadomościami tradycyjnych podejść . Koszt jest taki $m \times n$ $m \ge n$ $mn^2 - n^3/3$ $2mn^2 - 2n^3/3$ $m \gg n$ $\log p$ $p$ $n \log p$ $O(n^3 \log p)$ wykonywane są dodatkowe klapy, ale dla bardzo małych jest to często preferowane w stosunku do kosztów opóźnienia wysyłania większej liczby wiadomości (przynajmniej wtedy, gdy trzeba wykonać tylko jedną faktoryzację QR). $n$

— Jack Poulson
źródło

10

Dziwi mnie, że nikt nie wspominał o liniowych problemach najmniejszych kwadratów , które często występują w obliczeniach naukowych. Jeśli chcesz użyć eliminacji Gaussa, musisz utworzyć i rozwiązać równania normalne, które wyglądają następująco:

A^{T} A x = A^{T} b,

$A^{T}Ax = A^{T}b,$

gdzie jest macierzą punktów danych odpowiadających obserwacjom zmiennych niezależnych, jest wektorem parametrów do znalezienia, a jest wektorem punktów danych odpowiadających obserwacjom zmiennej zależnej. $A$ $x$ $b$

Jak często zauważa Jack Poulson, liczba warunków jest kwadratem liczby warunków , więc normalne równania mogą być katastrofalnie źle uwarunkowane. W takich przypadkach, chociaż podejścia oparte na QR i SVD są wolniejsze, dają znacznie dokładniejsze wyniki. $A^{T}A$ $A$

— Geoff Oxberry
źródło

2

Pozytywnie oceniany, ale QR powinien faktycznie być na równi z LU, jeśli weźmiesz pod uwagę niepotrzebne operacje potrzebne do utworzenia (QR wymaga tylko więcej flopów niż LU). Podejście SVD powinno być jednak wolniejsze (jego koszt może być mniej więcej ).

n^{3}

$n^3$

A^{H} A

$A^H A$

2 / 3 n^{3}

$2/3 n^3$

6 n^{3}

$6n^3$

— Jack Poulson

1

Oprócz stabilności gwarantowanej przez zastosowanie przekształceń ortogonalnych, wielką zaletą SVD jest to, że rozkład zapewnia własną kontrolę stanu, ponieważ stosunek największej do najmniejszej liczby pojedynczej jest dokładnie liczbą warunków (2-normalnych). W przypadku innych dekompozycji użycie estymatora warunków (np. Hagera-Highama) jest, choć nie tak drogie, jak sam dekompozycja, nieco „przylepione”.

— JM

1

@JackPoulson Czy z ciekawości masz odniesienie do liczby flopa dla SVD? Z tego, co mogę powiedzieć po krótkim spojrzeniu w Golub i Van Loan (str. 254 3. wydanie), stała wydaje się wyższa w przypadku używania SVD do rozwiązywania problemów z najmniejszymi kwadratami, ale mogę się mylić. Z góry dziękuję.

— OscarB

1

@OscarB: To była bardzo szorstka liczba od czubka mojej głowy, która jest niższa niż formowanie pełnego SVD (ponieważ możemy uniknąć kosztów ponownej transformacji). pracy jest potrzebne do redukcji do postaci dwukierunkowej (powiedzmy

8 / 3 n^{3}

$8/3 n^3$

A = F B G^{H}

$A=FBG^H$

C

$C$

B = U Σ V^{H}

$B=U\Sigma V^H$

x := (G (V (i n v (Σ) (U^{H} (F^{H} b)))))

$x := (G (V (\mathrm{inv}(\Sigma) (U^H (F^H b)))))$

O (n^{2})

$O(n^2)$

C

$C$

O (n^{2})

$O(n^2)$

1

\frac{σ_{1}}{σ_{n}}

$\frac{\sigma_1}{\sigma_n}$

3

Jak mierzysz wydajność? Prędkość? Precyzja? Stabilność? Szybki test w Matlabie daje następujące wyniki:

>> N = 100;
>> A = randn(N); b = randn(N,1);
>> tic, for k=1:10000, [L,U,p] = lu(A,'vector'); x = U\(L\b(p)); end; norm(A*x-b), toc
ans =
   1.4303e-13
Elapsed time is 2.232487 seconds.
>> tic, for k=1:10000, [Q,R] = qr(A); x = R\(Q'*b); end; norm(A*x-b), toc             
ans =
   5.0311e-14
Elapsed time is 7.563242 seconds.

Tak więc rozwiązanie pojedynczego systemu z rozkładem LU jest około trzy razy szybsze niż rozwiązanie z rozkładem QR, kosztem połowy cyfry dziesiętnej dokładności (ten przykład!).

— Pedro
źródło

Wszelkie zasugerowane przez Ciebie zalety są mile widziane.

— faleichik

3

Artykuł, który cytujesz, broni eliminacji Gaussa, mówiąc, że chociaż jest niestabilna numerycznie, zwykle dobrze sobie radzi na macierzach losowych, a ponieważ większość macierzy, o których można myśleć, jest jak macierze losowe, powinniśmy być w porządku. To samo stwierdzenie można powiedzieć o wielu metodach niestabilnych numerycznie.

Rozważ przestrzeń wszystkich macierzy. Te metody działają dobrze prawie wszędzie. To 99,999 ...% wszystkich macierzy, które można utworzyć, nie będzie miało problemów z niestabilnymi metodami. Jest tylko bardzo mały ułamek macierzy, dla których GE i inni będą mieli trudności.

Problemy, na które zwracają uwagę naukowcy, są zwykle w tej niewielkiej części.

Nie konstruujemy macierzy losowo. Konstruujemy macierze o bardzo specjalnych właściwościach, które odpowiadają bardzo specjalnym, nieprzypadkowym układom. Te matryce są często źle uwarunkowane.

Geometrycznie możesz wziąć pod uwagę przestrzeń liniową wszystkich macierzy. Przecięcie tej przestrzeni ma podprzestrzeń zerowej objętości / miary pojedynczych macierzy. Wiele problemów, które tworzymy, skupionych jest wokół tej podprzestrzeni. Nie są dystrybuowane losowo.

Jako przykład rozważ równanie cieplne lub dyspersję. Układy te mają tendencję do usuwania informacji z układu (wszystkie stany początkowe sprowadzają się do jednego stanu końcowego), w wyniku czego macierze opisujące te równania są niezwykle osobliwe. Proces ten jest bardzo mało prawdopodobny w przypadkowej sytuacji, ale wszechobecny w systemach fizycznych.

— MRocklin
źródło

2

Jeśli układ liniowy jest początkowo źle uwarunkowany, bez względu na zastosowaną metodę: zarówno rozkład LU, jak i QR da niedokładne wyniki. QR może wygrać tylko w przypadkach, gdy proces eliminacji Gaussa „psuje” dobrą matrycę. Głównym problemem jest to, że praktyczne przypadki takiego zachowania nie są znane.

— faleichik

W przypadku większości zastosowań naukowych zazwyczaj uzyskujemy matryce, które są rzadkie, symetryczne, określone pozytywnie i / lub dominują po przekątnej. Z nielicznymi wyjątkami, w macierzy znajduje się struktura, która pozwala nam wykorzystywać niektóre techniki zamiast tradycyjnej eliminacji gaussowskiej.

— Paweł

@Paul: Z drugiej strony gęsta eliminacja Gaussa jest miejscem, w którym większość czasu spędza się na wielopłaszczyznowej metodzie rzadkich macierzy niesymetrycznych.

— Jack Poulson

6

@Paul Po prostu nie jest prawdą, że „większość aplikacji wytwarza matryce SPD / diagonalnie dominujące”. Tak, zwykle istnieje jakaś struktura nadająca się do wykorzystania, ale problemy niesymetryczne i nieokreślone są niezwykle powszechne.

— Jed Brown

4

„W ciągu pięćdziesięciu lat komputerów nie pojawiły się żadne problemy z matrycą, które wzbudzałyby wybuchową niestabilność w naturalnych okolicznościach”. - LN Trefethen i D. Bau W swojej książce podają interesującą analizę probabilistyczną.

— JM