LOBPCG

Lokalnie optymalny blok wstępnie uwarunkowany gradient sprzężony ( LOBPCG ) to bezmacierzowa metoda znajdowania największych (lub najmniejszych) wartości własnych i odpowiadających im wektorów własnych symetrycznego uogólnionego problemu wartości własnych

${\ Displaystyle Ax = \ lambda Bx,}$

$dla$ danej pary $hermitowskich$ lub rzeczywistych symetrycznych macierzy macierz jest również dodatnio .

Tło

Kantorowicz w $=Ax-\lambda (x)x}$ roku zaproponował obliczenie najmniejszej wartości własnej $macierzy$ symetrycznej przez strome zejście przy użyciu kierunku $)$${\$ skalowanego gradientu ilorazu Rayleigha ${\ Displaystyle \ lambda (x) = (x, Ax) / (x, x)}$ w iloczynie skalarnym ${\ Displaystyle (x, y) ) = x'y}$ , z wielkością kroku obliczoną przez zminimalizowanie ilorazu Rayleigha w $.$$i$ wektorów , tj w lokalnie optymalny sposób. Samokish zaproponował zastosowanie a uwarunkowanie wstępne $resztkowego$$T}$ w $asymptotycznego$ wygenerowania kierunku uwarunkowanego $T {$ , gdy zbliża się do wektora własnego , współczynnik zbieżności miedza. D'yakonov zasugerował spektralnie równoważne warunkowanie wstępne i pochodne nieasymptotyczne granice tempa konwergencji. Blok lokalnie optymalny wieloetapowy najbardziej stromy spadek dla problemów z wartością własną został opisany w. Lokalna minimalizacja ilorazu Rayleigha na podprzestrzeni rozpiętej przez bieżące przybliżenie, bieżące przybliżenie resztkowe i poprzednie przybliżenie, a także jego wersję blokową, pojawiła się w. Wstępnie uwarunkowane wersja została przeanalizowana w i.

Główne cechy

• Bezmacierzowy , tj. nie wymaga jawnego przechowywania macierzy współczynników, ale może uzyskać dostęp do macierzy poprzez ocenę iloczynów macierzowo-wektorowych.
• faktoryzacji , tj. nie wymaga dekompozycji macierzy nawet dla uogólnionego problemu wartości własnych .
• Koszty na iterację i wykorzystanie pamięci są konkurencyjne w stosunku do metody Lanczosa , obliczającej pojedynczą skrajną parę własną macierzy symetrycznej.
• Zbieżność liniowa jest teoretycznie gwarantowana i praktycznie obserwowana.
• Przyspieszona zbieżność dzięki warunkowaniu bezpośredniemu , w przeciwieństwie do metody Lanczosa , obejmującemu warunkowanie zmienne i niesymetryczne oraz stałe i dodatnio określone .
• Umożliwia trywialne włączenie wydajnej dekompozycji domen i technik wielosiatkowych poprzez wstępne uwarunkowanie.
• Ciepły start i oblicza przybliżenie wektora własnego w każdej iteracji.
• Bardziej stabilny numerycznie w porównaniu z metodą Lanczosa i może działać w arytmetyce komputerowej o niskiej precyzji.
• Łatwy do wdrożenia, z wieloma wersjami już się pojawił.
• Blokowanie pozwala na wykorzystanie wysoce wydajnych operacji macierz-macierz, np. BLAS 3.
• Rozmiar bloku można dostosować, aby zrównoważyć szybkość zbieżności z kosztami komputerowymi ortogonalizacji i metodą Rayleigha-Ritza w każdej iteracji.

Algorytm

Wersja jednowektorowa

Wstępy: Zejście gradientowe dla problemów z wartościami własnymi

Metoda wykonuje iteracyjną maksymalizację (lub minimalizację) uogólnionego ilorazu Rayleigha

${\ Displaystyle \ rho (x): = \ rho (A, B; x): = {\ Frac { x^{T}Topór}{x^{T}Bx}},}$

co skutkuje znalezieniem największych (lub najmniejszych) par własnych ${\ Displaystyle Ax = \ lambda Bx.}$

Kierunek najbardziej stromego wzniesienia, czyli gradient , uogólnionego ilorazu Rayleigha jest dodatnio proporcjonalny do wektora

${\ Displaystyle r: = Ax- \ rho (x) Bx,}$

resztą wektora własnego . Jeśli warunek wstępny $stosowany$ do reszty i daje wektor

${\ Displaystyle w: = Tr,}$

zwaną wstępnie uwarunkowaną resztą. Bez uwarunkowania wstępnego ustawiamy ${\ displaystyle T: = I}$ i tak ${\ displaystyle w: = r}$ . Metoda iteracyjna

${\ Displaystyle x ^ {i + 1}: = x ^ {i} + \ alfa ^ {i }T(Ax^{i}-\rho (x^{i})Bx^{i}),}$

lub, w skrócie,

${\ Displaystyle x ^ {i + 1}: = x ^ {i} + \ alfa ^ {i} w ^ {i}, \,}$

${\ Displaystyle w ^ {i}: = Tr ^ {i}, \,}$

${\ Displaystyle r ^ {i}: =Ax^{i}-\rho (x^{i})Bx^{i},}$

jest znany jako wstępnie uwarunkowane $najbardziej$ strome wzniesienie ( zejście), gdzie skalar nazywany rozmiarem kroku. Optymalną wielkość kroku można określić maksymalizując iloraz Rayleigha, tj.

${\ Displaystyle x ^ {i + 1}: = \ arg \ max _ {y \ w rozpiętości \ { x^{i},w^{i}\}}\rho (y)}$

(lub ${\ displaystyle \ arg \ min}$ w przypadku minimalizacji), w takim przypadku metoda nazywana jest lokalnie optymalną.

Nawrót trzykrotny

Aby radykalnie przyspieszyć zbieżność lokalnie optymalnego, wstępnie uwarunkowanego, najbardziej stromego wzniesienia (lub zejścia), do dwuczłonowej relacji powtarzalności można dodać jeden dodatkowy wektor , aby uczynić ją trójczłonową:

${\ Displaystyle x ^ {i + 1}: = \ arg \ max _ {y \in span\{x^{i},w^{i},x^{i-1}\}}\rho (y)}$

(użyj ${\ displaystyle \ arg \ min}$ w przypadku minimalizacji). Maksymalizacja/minimalizacja ilorazu Rayleigha w trójwymiarowej podprzestrzeni może być przeprowadzona numerycznie metodą Rayleigha-Ritza . Dodanie większej liczby wektorów, patrz np. ekstrapolacja Richardsona , nie powoduje znacznego przyspieszenia, ale zwiększa koszty obliczeń, więc ogólnie nie jest zalecane.

Numeryczne ulepszenia stabilności

$się$ miarę zbiegania się iteracji wektory $–$ stają liniowo zależne powoduje utratę precyzji i powoduje, że Ritz metoda numerycznie niestabilna w obecności błędów zaokrągleń. Utraty precyzji można uniknąć, zastępując wektor wektorem $x ^$ ja ${i-1} }$$displaystyle span \ { x^{i},w^{i},x^{i-1}\}}$ mogą być dalej od , na podstawie trójwymiarowej podprzestrzeni $s$ , zachowując niezmienioną podprzestrzeń i unikając ortogonalizacji lub innych dodatkowych operacji. Ponadto ortogonalizacja podstawy trójwymiarowej podprzestrzeni może być potrzebna w przypadku źle uwarunkowanych problemy z wartościami własnymi w celu poprawy stabilności i osiągalnej dokładności.

Analogi podprzestrzeni Kryłowa

Jest to jednowektorowa wersja metody LOBPCG — jedna z możliwych uogólnień wstępnie uwarunkowanych liniowych solwerów gradientu sprzężonego na przypadek symetrycznych problemów z wartościami własnymi . $trywialnym$ przypadku $ja {\ displaystyle$ b $I}$ wynikowe przybliżenie z będzie inne niż uzyskane Lanczosa , chociaż oba przybliżenia będą należeć do tej samej podprzestrzeni Kryłowa .

Praktyczne scenariusze użycia

Ekstremalna prostota i wysoka wydajność jednowektorowej wersji LOBPCG sprawiają, że jest ona atrakcyjna dla aplikacji związanych z wartościami własnymi przy poważnych ograniczeniach sprzętowych, począwszy od wykrywania anomalii w czasie rzeczywistym w oparciu o klastry widmowe, poprzez partycjonowanie wykresów na wbudowanym ASIC lub FPGA , po modelowanie fizycznych zjawisk złożoności obliczeniowej na eksaskalowych superkomputerach TOP500 .

Wersja blokowa

Streszczenie

Kolejne pary własne można obliczać jeden po drugim za pomocą pojedynczego wektora LOBPCG uzupełnionego ortogonalną deflacją lub jednocześnie jako blok. W pierwszym podejściu niedokładności już obliczonych przybliżonych wektorów własnych dodatkowo wpływają na dokładność później obliczonych wektorów własnych, zwiększając w ten sposób błąd przy każdym nowym obliczeniu. Iteracja kilku przybliżonych wektorów własnych razem w bloku w lokalnie optymalny sposób w blokowej wersji LOBPCG. umożliwia szybkie, dokładne i niezawodne obliczanie wektorów własnych, w tym odpowiadających prawie wielokrotnym wartościom własnym, w których pojedynczy wektor LOBPCG cierpi z powodu powolnej zbieżności. Rozmiar bloku można dostosować, aby zrównoważyć stabilność numeryczną a szybkość zbieżności w stosunku do kosztów komputerowych ortogonalizacji i metody Rayleigha-Ritza w każdej iteracji.

Rdzeń projektu

Podejście blokowe $_$$_$ _ $}$ macierze $^ {i$ np . każda kolumna $X$$styl wyświetlania X^{i}}$ przybliża jeden z wektorów własnych. Wszystkie kolumny są iterowane jednocześnie $jest$ a następna macierz przybliżonych wektorów własnych określana metodą Rayleigha-Ritza obejmującej wszystkie kolumny ${\ Displaystyle X ^ {i}, \, W ^ {i},} i$$\ Displaystyle P ^ {i}}$ { . Każda kolumna ${\ Displaystyle W ^ {i}}$ jest obliczana po prostu jako wstępnie uwarunkowana reszta dla każdej kolumny ${\ Displaystyle X ^ {i}.}$$ja$$} , \, P ^ {i}]}$ ] { i ${\ Displaystyle [X ^ {i}, \, X ^ {i-1}]}$ są takie same.

Stabilność numeryczna a wydajność

Wynik $i$ Rayleigha-Ritza jest określony przez podprzestrzeń obejmującą wszystkie kolumny macierzy $\ Displaystyle X ^ {i}, \, W ^ {i},$ } , gdzie baza podprzestrzeni może być teoretycznie dowolna. Jednak w niedokładnej arytmetyce komputerowej metoda Rayleigha-Ritza staje się numerycznie niestabilny, jeśli niektóre wektory bazowe są w przybliżeniu liniowo zależne. Niestabilności liczbowe zwykle występują, np. jeśli niektóre wektory własne w bloku iteracyjnym osiągają już osiągalną dokładność dla danej precyzji komputera i są szczególnie widoczne w przypadku niskiej precyzji, np. pojedynczej precyzji .

Sztuka wielu różnych implementacji LOBPCG polega na zapewnieniu numerycznej stabilności metody Rayleigha-Ritza przy minimalnych kosztach obliczeniowych poprzez wybór dobrej bazy podprzestrzeni. Prawdopodobnie najbardziej stabilne podejście polegające na uczynieniu wektorów bazowych ortogonalnymi, np. za pomocą procesu Grama-Schmidta , jest również najbardziej kosztowne obliczeniowo. Na przykład implementacje $macierzach$ ale wydajną dekompozycję Cholesky'ego normalnej macierzy , która jest wykonywana tylko na poszczególnych i ${\ displaystyle P ^ {i}}$ , a nie na całej podprzestrzeni. Stale rosnąca ilość pamięci komputera pozwala obecnie na typowe rozmiary bloków w $,$ którym procent czasu obliczeniowego spędzanego na ortogonalizacjach i metoda zaczyna dominować.

Blokowanie wcześniej zbieżnych wektorów własnych

Metody blokowe dla problemów z wartościami własnymi, które iterują podprzestrzenie, zwykle mają niektóre iteracyjne wektory własne zbieżne szybciej niż inne, co motywuje blokowanie już zbieżnych wektorów własnych, tj. usuwanie ich z pętli iteracyjnej, w celu wyeliminowania niepotrzebnych obliczeń i poprawy stabilności numerycznej. Proste usunięcie wektora własnego może prawdopodobnie spowodować utworzenie jego duplikatu w wciąż powtarzających się wektorach. Fakt, że wektory własne symetrycznych problemów z wartościami własnymi są ortogonalne parami, sugeruje utrzymywanie wszystkich wektorów iteracyjnych ortogonalnych względem zablokowanych wektorów.

Blokowanie można zaimplementować w różny sposób, zachowując dokładność i stabilność numeryczną, jednocześnie minimalizując koszty obliczeniowe. Na przykład implementacje LOBPCG oddzielają twarde blokowanie, tj. deflację przez ograniczenie, w której zablokowane wektory własne służą jako wejście kodu i nie zmieniają się, od miękkiego blokowania, w którym zablokowane wektory nie uczestniczą w typowo najdroższym kroku iteracyjnym obliczania reszt, jednak w pełni uczestniczą w metodzie Rayleigha-Ritza, a zatem mogą być zmieniane metodą Rayleigha-Ritza.

Modyfikacje, LOBPCG II

LOBPCG obejmuje wszystkie kolumny macierzy w metodę Rayleigha-Ritza, $}$ wyniku czego ${\ Displaystyle X ^ {i}, \, W ^ {i$ } do ${\ displaystyle 3k}$ -by- $problem$$obliczenia$ wartością własną potrzebny do rozwiązania i do przy każdej iteracji, ${\ displaystyle k}$ oznacza rozmiar bloku - liczbę kolumn. W przypadku $dużych$ rozmiarów bloków to dominować w kosztach obliczeniowych i we / wy oraz ograniczać równoległość, w której wiele urządzeń obliczeniowych działa jednocześnie.

Oryginalny artykuł LOBPCG opisuje modyfikację, zwaną LOBPCG II, mającą na celu rozwiązanie takiego problemu, uruchamiającą jednowektorową wersję metody LOBPCG dla każdej pożądanej pary własnej z rozwiązywaniem procedury Rayleigha-Ritza 3 na 3 k {\ displaystyle $}$ problemy z przewidywaną wartością własną. Globalna procedura Rayleigha-Ritza dla wszystkich $własnych$ dotyczy każdej iteracji, ale tylko na kolumnach macierzy , zmniejszając $.$ liczbę niezbędnych iloczynów skalarnych do ${\ Displaystyle k ^ {2} + 3 k}$ od ${\ Displaystyle 9k ^ {2}}$ i wielkości globalnego przewidywanego problemu z wartością własną do ${\ displaystyle k}$ -by- ${\ displaystyle k }$ od ${\ displaystyle 3k}$ -by- ${\ displaystyle 3k}$ przy każdej iteracji. Odniesienie idzie dalej, stosując algorytm LOBPCG do każdego przybliżonego wektora własnego oddzielnie, tj. uruchamiając niezablokowaną wersję metody LOBPCG dla każdej pożądanej pary własnej dla ustalonej liczby iteracji. Procedury Rayleigha-Ritza w tych przebiegach wymagają jedynie rozwiązania zestawu problemów z przewidywaną wartością własną 3 × 3. Globalna procedura Rayleigha-Ritza dla wszystkich pożądanych par własnych jest stosowana tylko okresowo na końcu ustalonej liczby niezablokowanych iteracji LOBPCG.

Takie modyfikacje mogą być mniej niezawodne w porównaniu z oryginalnym LOBPCG. Indywidualnie działające gałęzie LOBPCG pojedynczego wektora mogą nie podążać ciągłymi iteracyjnymi ścieżkami, zamiast tego odwracać i tworzyć zduplikowane przybliżenia do tego samego wektora własnego. Pojedynczy wektor LOBPCG może być nieodpowiedni dla klastrowanych wartości własnych, ale oddzielne przebiegi LOBPCG małych bloków wymagają automatycznego określenia ich rozmiarów bloków podczas procesu iteracji, ponieważ liczba klastrów wartości własnych i ich rozmiary mogą być nieznane a priori.

Teoria i praktyka konwergencji

Konstrukcja LOBPCG gwarantuje minimalizację ilorazu Rayleigha nie wolniej niż blok o najbardziej stromym spadku gradientu , który ma wszechstronną teorię zbieżności. Każdy wektor własny jest stacjonarnym punktem ilorazu Rayleigha , w którym gradient znika. Zatem opadanie gradientu może zwolnić w pobliżu dowolnego wektora własnego , jednak zagwarantowane jest albo zbieżność do wektora własnego z liniową szybkością zbieżności, albo, jeśli ten wektor własny jest punktem siodłowym , iteracja iloraz Rayleigha spadnie poniżej odpowiedniej wartości własnej i zacznie zbiegać się liniowo do następnej wartości własnej poniżej. Najgorsza wartość współczynnika konwergencji liniowej została określona i zależy od względnej luki między wartością własną a resztą widma macierzy oraz jakości kondycjonera wstępnego , jeśli jest obecny.

W przypadku ogólnej macierzy najwyraźniej nie ma sposobu na przewidzenie wektorów własnych, a tym samym wygenerowanie początkowych przybliżeń, które zawsze działają dobrze. Rozwiązanie iteracyjne LOBPCG może być wrażliwe na początkowe przybliżenia wektorów własnych, np. zbieżność trwa dłużej, zwalniając jako przechodzące pośrednie pary własne. Co więcej, teoretycznie nie można zagwarantować zbieżności koniecznie do najmniejszej pary własnej, chociaż prawdopodobieństwo pomyłki wynosi zero. Dobrej jakości losowa funkcja Gaussa ze średnią zerową jest zwykle wartością domyślną w LOBPCG do generowania początkowych przybliżeń. Aby ustalić początkowe przybliżenia, można wybrać stałe ziarno dla generatora liczb losowych .

W przeciwieństwie do metody Lanczosa , LOBPCG rzadko wykazuje w praktyce asymptotyczną zbieżność superliniową .

Analiza częściowych głównych składowych (PCA) i dekompozycja wartości osobliwych (SVD)

LOBPCG można w prosty sposób zaadaptować do obliczania kilku największych wartości osobliwych i odpowiadających im wektorów osobliwych (częściowe SVD), np. do iteracyjnego obliczania PCA , dla macierzy danych D ze średnią zerową, bez jawnego obliczania macierzy kowariancji D ^T D , tj. w moda bez matrixa . Głównym obliczeniem jest obliczenie funkcji iloczynu D ^T (DX) macierzy kowariancji D ^T D i wektora blokowego X , który iteracyjnie przybliża pożądane wektory osobliwe. PCA wymaga największych wartości własnych macierzy kowariancji, podczas gdy LOBPCG jest zwykle implementowany do obliczania najmniejszych. Prostym obejściem jest zanegowanie funkcji, zastąpienie -D ^T (DX) zamiast DT (DX) , a ^tym samym odwrócenie kolejności wartości własnych, ponieważ LOBPCG nie dba o to, czy macierz problemu wartości własnej jest dodatnio określona, ​​czy nie .

LOBPCG dla PCA i SVD jest zaimplementowany w SciPy od wersji 1.4.0

Implementacje oprogramowania ogólnego

Wynalazca LOBPCG, Andrew Knyazev , opublikował referencyjną implementację o nazwie Block Locally Optimal Pre Condition Eigenvalue Xolvers (BLOPEX) z interfejsami do PETSc , hypre i Parallel Hierarchical Adaptive MultiLevel method (PHAML). Inne implementacje są dostępne np. w GNU Octave , MATLAB (w tym dla tablic rozproszonych lub kafelkowych), Java , Anasazi ( Trilinos ), SLEPc , SciPy , Julia , MAGMA, Pytorch , Rust , OpenMP i OpenACC , CuPy ( biblioteka tablic kompatybilna z NumPy , akcelerowana przez CUDA ), Google JAX i NVIDIA AMGX. LOBPCG jest zaimplementowany, ale nie jest uwzględniony w TensorFlow .

Aplikacje

Eksploracja danych

Pakiety oprogramowania scikit-learn i Megaman wykorzystują LOBPCG do skalowania klastrów widmowych i wielorakiego uczenia się za pomocą Laplace'owskich map własnych do dużych zbiorów danych. NVIDIA zaimplementowała LOBPCG w swojej bibliotece nvGRAPH wprowadzonej w CUDA 8. Sphynx, hybrydowy, równoległy partycjoner grafów z obsługą rozproszonej i współdzielonej pamięci - pierwsze narzędzie do partycjonowania grafów, które działa na procesorach graficznych w ustawieniach pamięci rozproszonej - wykorzystuje grupowanie spektralne do partycjonowania grafów , obliczanie wektorów własnych na Laplace'owska macierz wykresu przy użyciu LOBPCG z pakietu Anasazi .

Nauki o materiałach

LOBPCG jest zaimplementowany w ABINIT (w tym w wersji CUDA ) i Octopus . Został użyty do wielomiliardowych macierzy przez Nagrody Gordona Bella na superkomputerze Earth Simulator w Japonii. Model Hubbarda dla silnie skorelowanych układów elektronów w celu zrozumienia mechanizmu stojącego za nadprzewodnictwem wykorzystuje LOBPCG do obliczenia stanu podstawowego hamiltonianu na komputerze K i systemach z wieloma procesorami graficznymi.

Istnieją wersje MATLAB i Julia LOBPCG dla równań Kohna-Shama i teorii funkcjonału gęstości (DFT) wykorzystujące podstawę fali zwykłej. Ostatnie wdrożenia obejmują TTPY, Platypus-QM, MFDn, ACE-Molecule, LACONIC.

Mechanika i płyny

LOBPCG firmy BLOPEX jest używany do konfiguracji uwarunkowania wstępnego w bibliotece solwera Multilevel Balancing Domain Decomposition by Constraints (BDDC) BDDCML, która jest włączona do OpenFTL (Open Finite Element Library) i Flow123d symulatora przepływu wody podziemnej, substancji rozpuszczonej i transportu ciepła w spękanych porowatych mediach . LOBPCG został zaimplementowany w LS-DYNA .

Równania Maxwella

LOBPCG jest jednym z podstawowych narzędzi do rozwiązywania wartości własnych w PYFEMax i wysokowydajnym oprogramowaniu do multifizyki elementów skończonych Netgen/NGSolve. LOBPCG firmy hypre jest włączony do lekkiej, skalowalnej biblioteki C++ typu open source dla metod elementów skończonych MFEM , która jest używana w wielu projektach, w tym BLAST , XBraid, VisIt , xSDK, instytucie FASTMath w SciDAC oraz współprojektowaniu Centrum efektywnych dyskrecji eksaskalowych (CEED) w obliczeniach eksaskalowych Projekt.

Odszumianie

filtr dolnoprzepustowy oparty na LOBPCG może być użyty do odszumiania ; patrz np. przyspieszenie odszumiania całkowitej zmienności .

Segmentacja obrazu

Segmentacja obrazu za pomocą grupowania widmowego polega na osadzeniu niskowymiarowym przy użyciu macierzy powinowactwa między pikselami, po czym następuje grupowanie składowych wektorów własnych w przestrzeni niskowymiarowej, np. przy użyciu wykresu Laplaciana dla filtra dwustronnego . Segmentacja obrazu poprzez partycjonowanie wykresów widmowych przez LOBPCG z kondycjonowaniem wstępnym wielu siatek został po raz pierwszy zaproponowany i faktycznie przetestowany w i . To drugie podejście zostało później zaimplementowane w Pythonie scikit-learn , który wykorzystuje LOBPCG z SciPy z algebraicznymi warunkami wstępnymi multigrid do rozwiązania problemu wartości własnej dla grafu Laplaciana.

Linki zewnętrzne

• LOBPCG w MATLABIE
• LOBPCG w oktawie
• LOBPCG w SciPy
• LOBPCG w Javie w Google Code
• LOBPCG w bloku lokalnie optymalnych wstępnie uwarunkowanych wartości własnych Xolvers (BLOPEX) w GitHub i zarchiwizowanych w Google Code