Diagonalizacja gęstych, uwarunkowanych matryc

Próbuję diagonalizować niektóre gęste, źle uwarunkowane matryce. W precyzji maszynowej wyniki są niedokładne (zwracając ujemne wartości własne, wektory własne nie mają oczekiwanych symetrii). Przełączyłem się na funkcję Eigensystem [] Mathematiki, aby skorzystać z dowolnej precyzji, ale obliczenia są bardzo wolne. Jestem otwarty na dowolną liczbę rozwiązań. Czy istnieją pakiety / algorytmy, które dobrze nadają się do źle uwarunkowanych problemów? Nie jestem ekspertem w zakresie przygotowywania, więc nie jestem pewien, jak bardzo to może pomóc. W przeciwnym razie wszystko, co mogę wymyślić, to równoległe rozwiązania o dowolnej wartości własnej, ale nie znam niczego poza Mathematica, MATLAB i C ++.

Aby przedstawić tło problemu, matryce są duże, ale nie ogromne (maksymalnie 4096 x 4096 do 32768 x 32768). Są rzeczywiste, symetryczne, a wartości własne są ograniczone od 0 do 1 (wyłączne), przy czym wiele wartości własnych jest bardzo bliskich zeru i żaden nie jest zbliżony do 1. Macierz jest zasadniczo operatorem splotu. Nie potrzebuję przekątnej wszystkich moich matryc, ale im większy mogę przejść, tym lepiej. Mam dostęp do klastrów obliczeniowych z wieloma procesorami i możliwościami przetwarzania rozproszonego.

Dziękuję Ci

linear-algebra parallel-computing eigensystem precision dense-matrix Leigh
źródło

Jakiej procedury używasz do przekątnej swoich prawdziwych macierzy symetrycznych? I w jakim sensie rozkład wartości własnej jest niedokładny?

Jack Poulson,

Oto pomysł związany z odpowiedzią Arnolda: wykonaj rozkład Cholesky'ego swojej macierzy SPD, a następnie znajdź osobliwe wartości właśnie uzyskanego trójkąta Cholesky'ego, być może wykorzystując algorytm typu dqd w celu zachowania dokładności.

@JM: Tworzenie rozkładu Cholesky'ego pojedynczej liczbowo dodatniej określonej macierzy jest numerycznie niestabilne za pomocą zwykłej metody, ponieważ prawdopodobnie napotyka się przestawienia ujemne. (Np. Chol Matlaba (A) zwykle zawodzi.) Trzeba by było ustawić je na zero i unicestwić odpowiednie rzędy czynników. W ten sposób można uzyskać niezawodnie numeryczną spację zerową.

Arnold Neumaier

@Arnold, jeśli pamięć służy są adaptacjami Cholesky'iego używające symetryczny obrotowy dla tych przypadków, w których matryca jest pozytywny pół -definite (lub prawie). Może można by ich użyć ...

@JM: Nie trzeba obracać się, aby rozwiązać przypadek półfinału; przepis, który podałem, wystarczy. Chciałem tylko zaznaczyć, że nie można używać standardowych programów w puszkach, ale trzeba je samodzielnie zmodyfikować.

Arnold Neumaier

Odpowiedzi:

Oblicz SVD zamiast rozkładu widmowego. Wyniki są takie same w przypadku dokładnej arytmetyki, ponieważ macierz jest symetryczna, dodatnia, ale w przypadku arytmetyki o skończonej precyzji otrzymacie małe wartości własne z dużo większą dokładnością.

Edycja: Patrz Demmel i Kahan, Dokładne wartości osobliwe macierzy dwudzielnych, SIAM J. Sci. Stat. Comput. 11 (1990) 873–912.
ftp://netlib2.cs.utk.edu/lapack/lawnspdf/lawn03.pdf

Edytuj2; Zauważ, że żadna metoda nie będzie w stanie rozwiązać wartości własnych mniejszych niż normalna razy dokładność zastosowanej maszyny, ponieważ zmiana pojedynczego wpisu o jedno ulp może już tak bardzo zmienić małą wartość własną. Zatem uzyskanie zerowych wartości własnych w miejsce bardzo małych jest właściwe i żadna metoda (oprócz pracy z większą precyzją) nie rozplątuje odpowiadających sobie wektorów własnych, a jedynie zwraca podstawę dla wspólnej numerycznej przestrzeni zerowej.

Arnold Neumaier
źródło

Nie jestem pewien, czy w to wierzę, ponieważ większość implementacji SVD zaczyna się od jednostkowego zmniejszenia do rzeczywistej postaci dwukierunkowej, a następnie zasadniczo oblicza hermitowską EVD powiązanej macierzy, takiej jak , które można łatwo przekształcić w prawdziwą symetryczną formę tridiagonalną. Względna dokładność jest wysoce zależna od tego, która metoda zostanie zastosowana do rozwiązania skondensowanego EVP / SVD, i nie widzę, gdzie SVD ma przewagę ... Jestem pewien, że jest to omówione w jednym lub więcej artykułów Demmela.

[0, B^{T}; B, 0]

$[0, B^T; B, 0]$

Jack Poulson,

@JackPoulson: Chodzi o to, że forma dwukierunkowa znacznie lepiej określa małe liczby pojedyncze. Przynależna symetryczna forma trójosiowa ma zera na przekątnej, które są zachowane przez dwukrotne zmniejszenie do przekątnej, ale nie przez QR zastosowane do trójosiowej.

Arnold Neumaier,

Odniesienie? Metoda Jacobiego jest znana jako bardzo dokładna (choć powolna).

Jack Poulson,

@JackPoulson: Spróbuj i zobacz. Demmel i Kahan, Dokładne wartości osobliwe macierzy dwudzielnych

Arnold Neumaier

Widzę, do czego zmierzasz: algorytm QR dla dwukierunkowego SVD może być względnie dokładny ze względu na zerową przekątną permutacji matryca, podczas gdy ta technika „zerowego przesunięcia” nie działa w przypadku dowolnych macierzy tridiagonalnych, więc trójpątne EVP oparte na algorytmie QR będą mniej dokładne dla małych wartości własnych. Chodzi o to, że to rozumowanie zakłada eigensolver oparty na użyciu algorytmu QR na skondensowanej formie; Jacobi jest godnym uwagi wyjątkiem, ale być może jest wolniejszy niż SVD oparty na algorytmie QR. MRRR może czasem osiągnąć również wysoką dokładność względną.

[0, B^{T}; B, 0]

$[0,B^T; B,0]$

Jack Poulson,

Dziękuję za tę sugestię. Wypróbowałem polecenie SVD Mathematiki, ale nie dostrzegam zauważalnej poprawy (wciąż brakuje odpowiednich symetrii, „wartości własne” są niepoprawnie zerowe, gdy wcześniej nieprawidłowo wychodziły ujemne). Może potrzebowałbym zaimplementować jeden z opisanych powyżej algorytmów zamiast wbudowanej funkcji? Prawdopodobnie chciałbym uniknąć kłopotów z użyciem takiej metody, chyba że z góry wiedziałem, że przyniesie to znaczną poprawę.

@JackPoulson, przejrzałem artykuł na temat metody Jacobiego, do której się odwoływałeś, i wygląda to obiecująco. Czy ty lub ktokolwiek może polecić dobry sposób na wdrożenie metody Jacobiego w poszukiwaniu eigensystemów? Zgaduję, że gdybym sam to kodował (w MATLAB), byłoby to bardzo wolne.

Leigh
źródło

Nie testowałem tego, ale tutaj jest implementacja MATLAB: groups.google.com/forum/?fromgroups#!msg/sci.math.num-analysis/...

Jack Poulson

Zauważ, że żadna metoda nie będzie w stanie rozwiązać wartości własnych mniejszych niż normalna razy dokładność zastosowanej maszyny, ponieważ zmiana pojedynczego wpisu o jedno ulp może już tak bardzo zmienić małą wartość własną. Zatem uzyskanie zerowych wartości własnych w miejsce bardzo małych jest właściwe i żadna metoda (z wyjątkiem pracy z większą precyzją) nie rozplątuje odpowiadających sobie wektorów własnych, ale po prostu zwraca podstawę dla wspólnej numerycznej przestrzeni zerowej. Do czego potrzebujesz wartości własnych?

Arnold Neumaier,

@ArnoldNeumaier: Przeprowadziłem kilka testów w MATLAB z wartościami własnymi w zakresie [0,1], z jedną wartością własną ręcznie ustawioną na wartości takie jak 6.3e-16 i procedurą SVD Octave (opartą na dgesvd, która wykorzystuje redukcję do dwukierunkowej i następnie QR) pobiera te wartości znacznie dokładniej niż eig Octave'a. Połączony kod Jacobi wydaje się być zbyt wolny, aby używać go nawet w przypadku matryc o niewielkich rozmiarach.

Jack Poulson,

@JackPoulson: Tak. Ale Leigh wydaje się narzekać na wiele bardzo małych wartości własnych, a ich wektorami własnymi rzadko są te zaprojektowane, ale swobodnie się mieszają, bez względu na zastosowaną metodę. I dodatnie bardzo małe dodatnie wartości (mniejsze niż 1e-16) będą oczywiście wynosić zero.

Arnold Neumaier,

@ArnoldNeumaier ma rację, że znajduję wiele bardzo małych wartości własnych, które, jak sądzę, zaostrzają problem. Nie zdawałem sobie sprawy (choć z perspektywy czasu jest oczywiste), że wartości własne mniejsze niż 1e-16 będą zerowe w zmiennoprzecinkowych. Wydaje mi się, że chociaż liczba może być przechowywana, błąd zaokrąglania występuje podczas dodawania jej do większej liczby. Wektory własne mówią mi, czy jakiś problem można rozwiązać. Wektor własny umożliwia rozkład problemu na części możliwe do rozwiązania i nierozwiązywalne. Jeśli jestem zasadniczo ograniczony precyzją, to czy możesz polecić jakieś pakiety dla szybszego rozwiązania?

Leigh,