Zrozumienie, jak Numpy robi SVD

Stosowałem różne metody do obliczania zarówno rangi macierzy, jak i rozwiązania układu równań macierzowych. Natknąłem się na funkcję linalg.svd. Porównując to do własnych wysiłków związanych z rozwiązaniem systemu z eliminacją Gaussa, wydaje się on zarówno szybszy, jak i bardziej precyzyjny. Próbuję zrozumieć, jak to jest możliwe.

O ile mi wiadomo, funkcja linalg.svd wykorzystuje algorytm QR do obliczania wartości własnych mojej macierzy. Wiem, jak to działa matematycznie, ale nie wiem, jak Numpy udało się to zrobić tak szybko i bez utraty precyzji.

Więc moje pytanie: jak działa funkcja numpy.svd, a dokładniej, w jaki sposób potrafi to zrobić szybko i dokładnie (w porównaniu z eliminacją gaussa)?

W twoim pytaniu jest wiele problemów.

Nie używaj eliminacji Gaussa (faktoryzacja LU), aby obliczyć liczbową pozycję macierzy. Faktoryzacja LU jest w tym celu zawodna w arytmetyki zmiennoprzecinkowej. Zamiast tego zastosuj rozkład QR ujawniający rangę (taki jak xGEQPXlub xGEPQYw LAPACK, gdzie x to C, D, S lub Z, chociaż te procedury są trudne do wyśledzenia; zobacz odpowiedź JedBrown na powiązane pytanie ) lub użyj SVD (rozkład wartości w liczbie pojedynczej, taki jak xGESDDlub xGESVD, gdzie x oznacza ponownie C, D, S lub Z). SVD jest bardziej dokładnym, niezawodnym algorytmem do określania rangi numerycznej, ale wymaga większej liczby operacji zmiennoprzecinkowych.

Jednak w przypadku rozwiązania układu liniowego faktoryzacja LU (z częściowym obrotem, co jest standardową implementacją w LAPACK) jest wyjątkowo niezawodna w praktyce. Istnieją pewne przypadki patologiczne, w których rozkład na czynniki pierwsze z częściowym przestawieniem jest niestabilny (patrz Wykład 22 w Numerycznej algebrze liniowej)Trefethen i Bau w celu uzyskania szczegółowych informacji). Faktoryzacja QR jest bardziej stabilnym algorytmem numerycznym do rozwiązywania układów liniowych i prawdopodobnie dlatego daje tak precyzyjne wyniki. Wymaga to jednak więcej operacji zmiennoprzecinkowych niż faktoryzacja LU o współczynnik 2 dla macierzy kwadratowych (wierzę; JackPoulson może mnie w tym poprawić). W przypadku systemów prostokątnych lepszym wyborem jest faktoryzacja QR, ponieważ pozwoli ona uzyskać rozwiązania najmniejszych kwadratów w przypadku przesadnie określonych systemów liniowych. SVD może być również stosowany do rozwiązywania układów liniowych, ale będzie droższy niż faktoryzacja QR.

janneb ma rację, że numpy.linalg.svd jest opakowaniem xGESDDw LAPACK. Rozkłady poszczególnych wartości przebiegają w dwóch etapach. Po pierwsze, matryca, która ma zostać rozłożona, jest redukowana do postaci dwu- kątnej. Algorytmem zastosowanym w LAPACK jest redukcja do postaci dwukierunkowej, prawdopodobnie jest to algorytm Lawsona-Hansona-Chana, który w pewnym momencie wykorzystuje faktoryzację QR. Wykład 31 w Numerycznej algebrze liniowej autorstwa Trefethena i Baua zawiera przegląd tego procesu. Następnie xGESDDużywa algorytmu dziel i zwyciężaj, aby obliczyć wartości osobliwe oraz lewe i prawe wektory osobliwe na podstawie macierzy dwukierunkowej. Aby uzyskać informacje na temat tego kroku, musisz skonsultować obliczenia macierzowe Goluba i Van Loana lub stosowaną numeryczną algebrę liniową autorstwa Jima Demmela.

Wreszcie, nie należy mylić wartości pojedynczych z wartościami własnymi . Te dwa zestawy wielkości nie są takie same. SVD oblicza pojedyncze wartości macierzy. Obliczenia numeryczne Cleve Moler za pomocą MATLAB dają ładny przegląd różnic między wartościami osobliwymi a wartościami własnymi . Zasadniczo nie ma oczywistej zależności między wartościami osobliwymi danej macierzy a jej wartościami własnymi, z wyjątkiem macierzy normalnych , gdzie wartości osobliwe są wartością bezwzględną wartości własnych.

Ze względu na brzmienie twojego pytania zakładam, że macierz jest kwadratowa. Procedury SVD LAPACK, takie jak zgesvd , zasadniczo przebiegają w trzech etapach dla macierzy kwadratowych:

1. obliczanie macierzy jednostkowych i , jako przekształceń Householdera, tak że ogólna macierz jest zredukowana do rzeczywistej, górnej macierzy dwukierunkowej . Po wyjściu z tego podprogramu wektory Householdera (znormalizowane tak, że ich pierwszy wpis jest domyślnie jeden) dla i są odpowiednio przechowywane w częściach poniżej i na prawo od głównej i super-przekątnej. Ten krok wymaga pracy .$U_A$$V_A$$A$$B := U_A^H A V_A$$U_A$$V_A$$B$$O(n^3)$
2. Odmiana algorytmu do obliczania rozkładu wartości własnej rzeczywistej symetrycznej macierzy tridiagonalnej jest wykorzystywana do obliczania dwuwymiarowego SVD. Chociaż najbardziej znany algorytm dla prawdziwych symetryczny tridiagonal EVP'S (MWSZ) nie jest jeszcze, według mojej wiedzy, stabilne dla bidiagonal SVD, znajduje się interesująca dyskusja tutaj . LAPACK obecnie używa dziel i rządź podejścia do bidiagonal SVD. W bidiagonal wydajności SVD tak, że . Ten krok będzie wymagał pracy gdy MRRR jest stabilny dla SVD, ale obecnie wymaga tyle samo pracy .$\{U_B,V_B,\Sigma\}$$B = U_B \Sigma V_B^H$$O(n^2)$$O(n^3)$
3. Ten krok jest trywialny. Ponieważ , , aby SVD można było jawnie utworzyć po zastosowaniu transformacji Householdera, które domyślnie definiują i dla i . Ten krok wymagał pracy .A = ( U A U B ) Σ ( V A V B ) H U A V A U B V B O ( n 3$U_A B V_A^H = A$$A = (U_A U_B) \Sigma (V_A V_B)^H$$U_A$$V_A$$U_B$$V_B$$O(n^3)$

numpy.linalg.svd to opakowanie wokół {Z, D} GESDD z LAPACK. Z kolei LAPACK jest bardzo starannie napisany przez czołowych światowych ekspertów w dziedzinie numerycznej algebry liniowej. Rzeczywiście, byłoby bardzo zaskakujące, gdyby komuś niezbyt dobrze zaznajomionemu z polem udało się pokonać LAPACK (pod względem szybkości lub dokładności).

Jeśli chodzi o to, dlaczego QR jest lepszy niż eliminacja Gaussa, jest to prawdopodobnie bardziej odpowiednie dla /scicomp//