Kiedy transformacje ortogonalne przewyższają eliminację Gaussa?

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? ..

Precyzja

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

1. „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)$

2. „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$

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:

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$

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!).

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.