SVD do znalezienia największej wartości własnej macierzy 50x50 - czy marnuję znaczną ilość czasu?

Mam program, który oblicza największą wartość własną spośród wielu rzeczywistych symetrycznych macierzy 50 x 50, wykonując dekompozycje o liczbie pojedynczej na wszystkich z nich. SVD stanowi wąskie gardło w programie.

Czy istnieją algorytmy, które są znacznie szybsze w znajdowaniu największej wartości własnej, czy też optymalizacja tej części nie przyniosłaby dużego zwrotu z inwestycji?

linear-algebra eigensystem

— Ania
źródło

Czy możesz podać więcej informacji na temat swoich matryc, np. Jeśli wiadomo coś o ich strukturze, zakresie ich wartości własnych lub ich podobieństwie do siebie?

— Pedro

Jest to macierz kowariancji (

). Testy pokazują, że wszystkie oprócz 5 największych wartości własnych są bliskie zeru i że największa wartość własna jest co najmniej ~ 20% większa niż druga co do wielkości. Ponieważ istnieje wiele wartości własnych bliskich zeru, przypuszczam, że zasięg nie jest ważny? Można go przeskalować do dowolnego zakresu. Skala, której obecnie używam, daje mi zakres 150 ~ 200.

X X^{T}

$XX^T$

— Anna

Ponadto matryca nie jest bardzo ściśle pojedyncza, więc problem SVD jest dobrze uwarunkowany.

— Anna

Ponieważ

jest symetryczny i dodatni (pół) określony, można użyć faktoryzacji Cholesky'ego zamiast SVD. Faktoryzacja Cholesky'ego wymaga obliczenia o wiele mniejszej liczby klap niż SVD, ale bycie dokładną metodą nadal wymaga

X X^{T}

$XX^T$

O (n^{3})

$O(n^3)$

— Ken

@Anna: Czy wypróbowałeś którąś z wielu proponowanych tutaj metod? Byłbym bardzo ciekawy, co dla ciebie działało najlepiej w praktyce ...

— Pedro

Odpowiedzi:

W zależności od precyzji wymaganej dla największej wartości własnej, możesz spróbować użyć iteracji mocy .

W twoim konkretnym przykładzie posunąłbym się do tego, aby nie formułować jawnie, ale obliczać w każdej iteracji. Obliczenie wymagałoby operacji , podczas gdy iloczyn macierz-wektor wymaga tylko . $A=XX^\mathsf{T}$ $x \leftarrow X(X^\mathsf{T}x)$ $A$ $\mathcal O(n^3)$ $\mathcal O(n^2)$

Współczynnik konwergencji zależy od podziału między dwiema największymi wartościami własnymi, więc nie zawsze może to być dobre rozwiązanie,

— Pedro
źródło

Jeśli największa wartość własna jest o 20% większa od następnej, iteracja mocy powinna zbiegać się dość szybko (wszystkie inne wartości własne są tłumione współczynnikiem 5/6 w każdej iteracji, więc otrzymujesz jedną cyfrę za każde 13 iteracji.

— Wolfgang Bangerth

Metody podprzestrzeni Kryłowa są zdecydowanie lepsze niż metody mocy, ponieważ zawierają wektor z iteracji mocy o tej samej liczbie iteracji.

— Jack Poulson

@JackPoulson: Tak, ale każda iteracja jest droższa do obliczenia ... Czy naprawdę warto byłoby tak małym problemem?

— Pedro

@Pedro: oczywiście matvecy wymagają kwadratowej pracy, a eigensolve i późniejszy rozwój Rayleigha są trywialne w porównaniu.

— Jack Poulson,

Koszt kodu? Ponieważ @JackPoulson poruszył ten problem, B. Parlett i wsp. (1982) („O szacowaniu największej wartości własnej z algorytmem Lanczosa”) porównują metodę mocy, metodę mocy + przyspieszenie Aitkena i zastosowanie Lanczosa, kierując się na największą wartość własną rzeczywistej symetryczny (lub hermitowski) poz. def. matryca. Stwierdzili, że metoda Lanczosa jest bardziej wydajna, jeśli potrzebna jest nawet niewielka dokładność (pierwszej wartości własnej w stosunku do drugiej), i lepiej w zapobieganiu nieporozumieniom.

— hardmath

$X(X^Tx)$

— Arnold Neumaier
źródło

B = T = I

$B = T = I$

Nie; Miałem na myśli standardowy algorytm Lanczsosa, ale w pośpiechu napisałem CG. Teraz poprawione.

— Arnold Neumaier

$A = XX^T$ $\mu$ $A - \mu I$ $A$

$||Ax||/||x||$ $\lambda_1$ $[0,\frac{5}{6} \lambda_1]$ $A - \frac{5}{12} \lambda_1 I$ $\frac{7}{12} \lambda_1$ $[\frac{-5}{12} \lambda_1, \frac{5}{12} \lambda_1]$

$A - \mu I$ $\lambda_1$ $\mu$

$A - \mu I$ $(A - \mu I)x = X(X^Tx) - \mu x$ $O(n^2)$

— matematyka
źródło

Wydaje się, że wymaga to dobrego wyobrażenia o wielkości drugiej co do wielkości wartości własnej. Jak byś to przybliżył w takim przypadku?

— Pedro

λ_{1}

$\lambda_1$

| λ_{2} | / | λ_{1} |

$|\lambda_2|/|\lambda_1|$

| λ_{2} | / | λ_{1} |

$|\lambda_2|/|\lambda_1|$

λ_{2}

$\lambda_2$

λ_{1}

$\lambda_1$ w razie potrzeby. Sugerowałem, jakie korzyści zobaczysz w przypadku takim, jak Anna opisuje w swoich komentarzach pod pytaniem.

— hardmath