Próbuję obliczyć kilka (5-500) wektorów własnych odpowiadających najmniejszym wartościom własnym dużych symetrycznych kwadratowych macierzy rzadkich (do 30000 x 30000), przy czym mniej niż 0,1% wartości jest niezerowe.
Obecnie używam scipy.sparse.linalg.eigsh w trybie shift-invert (sigma = 0,0), co do którego doszedłem przez różne posty na ten temat, jest preferowanym rozwiązaniem. Jednak w większości przypadków rozwiązanie problemu zajmuje 1 godzinę. Z drugiej strony funkcja jest bardzo szybka, jeśli poproszę o największe wartości własne (poniżej sekundy w moim systemie), czego oczekiwano z dokumentacji.
Ponieważ jestem bardziej zaznajomiony z Matlabem w pracy, próbowałem rozwiązać problem w Octave, co dało mi ten sam wynik za pomocą eigs (sigma = 0) w ciągu zaledwie kilku sekund (poniżej 10s). Ponieważ chcę przeprowadzić analizę parametrów algorytmu, w tym obliczenia wektora własnego, tego rodzaju zysk czasu byłby świetny również w Pythonie.
Najpierw zmieniłem parametry (szczególnie tolerancję), ale niewiele to zmieniło w czasie.
Używam Anacondy w systemie Windows, ale próbowałem przełączyć LAPACK / BLAS używany przez scipy (co było ogromnym bólem) z mkl (domyślna Anaconda) na OpenBlas (używany przez Octave zgodnie z dokumentacją), ale nie mogłem zobaczyć zmiany w wydajność.
Nie byłem w stanie ustalić, czy było coś do zmiany w używanym ARPACK (i jak)?
Przesłałem walizkę testową dla poniższego kodu do następującego folderu dropbox: https://www.dropbox.com/sh/l6aa6izufzyzqr3/AABqij95hZOvRpnnjRaETQmka?dl=0
W Pythonie
import numpy as np
from scipy.sparse import csr_matrix, csc_matrix, linalg, load_npz
M = load_npz('M.npz')
evals, evecs = linalg.eigsh(M,k=6,sigma=0.0)
W oktawie:
M=dlmread('M.txt');
M=spconvert(M);
[evecs,evals] = eigs(M,6,0);
Każda pomoc jest doceniana!
Kilka dodatkowych opcji, które wypróbowałem na podstawie komentarzy i sugestii:
Oktawa:
eigs(M,6,0)
i eigs(M,6,'sm')
daj mi ten sam wynik:
[1.8725e-05 1.0189e-05 7.5622e-06 7.5420e-07 -1.2239e-18 -2.5674e-16]
podczas gdy eigs(M,6,'sa',struct('tol',2))
zbiega się do
[1.0423 2.7604 6.1548 11.1310 18.0207 25.3933]
znacznie szybciej, ale tylko wtedy, gdy wartości tolerancji przekraczają 2, w przeciwnym razie w ogóle się nie zbiegają, a wartości są bardzo różne.
Python:
eigsh(M,k=6,which='SA')
i eigsh(M,k=6,which='SM')
oba nie są zbieżne (błąd ARPACK przy braku osiągniętej zbieżności). eigsh(M,k=6,sigma=0.0)
Podaje tylko niektóre wartości własne (po prawie godzinie), które różnią się od oktawy dla najmniejszych (znaleziono nawet 1 dodatkową małą wartość):
[3.82923317e-17 3.32269886e-16 2.78039665e-10 7.54202273e-07 7.56251500e-06 1.01893934e-05]
Jeśli tolerancja jest wystarczająco wysoka, otrzymuję również wyniki eigsh(M,k=6,which='SA',tol='1')
, które zbliżają się do innych uzyskanych wartości
[4.28732218e-14 7.54194948e-07 7.56220703e-06 1.01889544e-05, 1.87247350e-05 2.02652719e-05]
ponownie z inną liczbą małych wartości własnych. Czas obliczeń nadal wynosi prawie 30 minut. Chociaż różne bardzo małe wartości mogą być zrozumiałe, ponieważ mogą reprezentować wielokrotności 0, inna mnogość zaskakuje mnie.
Ponadto wydaje się, że istnieją pewne podstawowe różnice w SciPy i Octave, których nie mogę jeszcze zrozumieć.
źródło
Odpowiedzi:
Przypuszczenie i kilka komentarzy, ponieważ nie mam Matlaba / Octave:
Aby znaleźć małe wartości własne macierzy symetrycznych o wartościach własnych> = 0, podobnie jak Twoja, następujące jest o wiele szybsze niż odwrócenie shift:
eigsh( Aflip )
dla dużych par własnych jest szybszy niż shift-invert dla małych, ponieważA * x
jest szybszy niż to,solve()
co musi zrobić inwersja shift . Możliwe, że Matlab / Octave zrobi toAflip
automatycznie, po szybkim teście pozytywnej oceny z Cholesky.Czy umiesz biegać
eigsh( Aflip )
w Matlabie / Oktawie?Inne czynniki, które mogą mieć wpływ na dokładność / prędkość:
Domyślnym Arpack dla wektora początkowego
v0
jest wektor losowy. Używamv0 = np.ones(n)
, co dla niektórych może być okropne,A
ale jest powtarzalne :)Ta
A
matryca jest prawie dokładnie sigularna,A * ones
~ 0.Wielordzeniowy: scipy-arpack z openblas / Lapack używa ~ 3,9 z 4 rdzeni na moim komputerze iMac; czy Matlab / Octave używają wszystkich rdzeni?
Oto wartości własne scipy-Arpack dla kilku
k
itol
, grep z logów pod gist.github :Czy Matlab / Octave jest mniej więcej taki sam? Jeśli nie, wszystkie zakłady są wyłączone - najpierw sprawdź poprawność, a następnie szybkość.
Dlaczego wartości własne tak chwieją się? Małe <0 dla matrycy, która jak się wydaje nie jest ujemna, jest oznaką błędu zaokrąglenia , ale zwykła sztuczka drobnego przesunięcia
A += n * eps * sparse.eye(n)
nie pomaga.Skąd to się
A
bierze, z jakiego obszaru problemowego? Czy możesz wygenerować podobnyA
, mniejszy lub rzadszy?Mam nadzieję że to pomoże.
źródło
tol
Jest również bałagan dla małych wartości własnych - zadaj nowe pytanie, jeśli chcesz, daj mi znać.Wiem, że to już stare, ale miałem ten sam problem. Czy recenzja tutaj ( https://docs.scipy.org/doc/scipy/reference/tutorial/arpack.html )?
Wygląda na to, że gdy ustawisz sigma na niską liczbę (0), powinieneś ustawić, która = „LM”, nawet jeśli chcesz niskich wartości. Wynika to z faktu, że ustawienie sigma przekształca wartości, które chcesz (w tym przypadku niskie), wydają się wysokie, a więc nadal możesz korzystać z metod „LM”, które są znacznie szybsze, aby uzyskać to, czego chcesz (niskie wartości własne ).
źródło
Chcę najpierw powiedzieć, że nie mam pojęcia, dlaczego wyniki zgłoszone przez Ciebie i @Bill są takie, jakie są. Zastanawiam się tylko, czy
eigs(M,6,0)
w Octave odpowiadak=6 & sigma=0
, czy może jest to coś innego?Dzięki scipy, jeśli nie dostarczę sigmy, w ten sposób mogę uzyskać wynik w przyzwoitym czasie.
Nie jestem jednak całkowicie pewien, czy ma to sens.
Jedynym sposobem, w jaki znalazłem użycie sigmy i uzyskanie wyniku w przyzwoitym czasie, jest podanie M jako LinearOperator. Nie jestem zbyt obeznany z tym, ale z tego, co zrozumiałem, moja implementacja reprezentuje matrycę tożsamości, która powinna być M, jeśli nie zostanie podana w wywołaniu. Powodem tego jest to, że zamiast wykonywania bezpośredniego rozwiązania (dekompozycji LU), scipy użyje iteracyjnego solwera, który jest potencjalnie bardziej odpowiedni. Dla porównania, jeśli podasz
M = np.identity(a.shape[0])
, co powinno być dokładnie takie samo, eigsh zajmie wieczność, aby dać wynik. Pamiętaj, że to podejście nie działa, jeślisigma=0
jest zapewnione. Ale nie jestem pewien, czysigma=0
to naprawdę jest przydatne?Znów nie mam pojęcia, czy jest poprawny, ale zdecydowanie różni się od poprzedniego. Byłoby wspaniale mieć wkład kogoś z scipy.
źródło