Próbuję obliczyć kilka (5-500) wektorów własnych odpowiadających najmniejszym wartościom własnym dużych symetrycznych kwadratowych macierzy rzadkich (do 30000x30000) z mniej niż 0,1% wartości niezerowych.
Obecnie używam scipy.sparse.linalg.eigsh w trybie shift-invert (sigma = 0.0), co jest preferowanym rozwiązaniem w różnych postach na ten temat. Jednak w większości przypadków rozwiązanie problemu zajmuje do 1 godziny. Z drugiej strony funkcja jest bardzo szybka, jeśli poproszę o największe wartości własne (subsekundy w moim systemie), czego oczekiwano od dokumentacji.
Ponieważ Matlab jest mi bardziej zaznajomiony z pracy, próbowałem rozwiązać problem w Octave, co dało mi ten sam wynik przy użyciu eigs (sigma = 0) w zaledwie kilka sekund (poniżej 10s). Ponieważ chcę wykonać przemiatanie parametrów algorytmu, w tym obliczenia wektora własnego, ten rodzaj przyrostu czasu byłby świetny również w Pythonie.
Najpierw zmieniłem parametry (zwłaszcza tolerancję), ale nie zmieniło to wiele w skali czasowej.
Używam Anacondy w systemie Windows, ale próbowałem zmienić LAPACK / BLAS używany przez scipy (co było ogromnym bólem) z mkl (domyślny Anaconda) na OpenBlas (używany przez Octave zgodnie z dokumentacją), ale nie widziałem zmiany w wydajność.
Nie udało mi się ustalić, czy było coś do zmiany w używanym ARPACKU (i jak)?
Przesłałem testcase dla poniższego kodu do następującego folderu dropbox: https://www.dropbox.com/sh/l6aa6izufzyzqr3/AABqij95hET=Qmkannj?
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ą 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 są powyżej 2, w przeciwnym razie wcale się nie zbiegają i wartości są mocno różne.
Pyton: eigsh(M,k=6,which='SA') i eigsh(M,k=6,which='SM') oba nie są zbieżne (błąd ARPACK przy braku zbieżności). Tylko eigsh(M,k=6,sigma=0.0) podaje pewne 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 z eigsh(M,k=6,which='SA',tol='1'), które są zbliżone do innych uzyskanych wartości
[4.28732218e-14 7.54194948e-07 7.56220703e-06 1.01889544e-05, 1.87247350e-05 2.02652719e-05]
Znowu 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, różna krotność wprawia mnie w zakłopotanie.
Ponadto wydaje się, że istnieją pewne fundamentalne różnice w SciPy i Octave, których jeszcze nie mogę rozgryźć.
3 odpowiedzi
Chcę najpierw powiedzieć, że nie mam pojęcia, dlaczego wyniki zgłoszone przez Ciebie i @Billa są takie, jakie są. Po prostu zastanawiam się, czy eigs(M,6,0) w oktawie odpowiada k=6 & sigma=0, czy może jest to coś innego?
Z scipy, jeśli nie podam sigmy, w ten sposób mogę uzyskać wynik w przyzwoitym czasie.
import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigsh
from time import perf_counter
M = np.load('M.npz')
a = csr_matrix((M['data'], M['indices'], M['indptr']), shape=M['shape'])
t = perf_counter()
b, c = eigsh(a, k=50, which='SA', tol=1e1)
print(perf_counter() - t)
print(b)
Nie jestem jednak całkowicie pewien, czy to ma sens.
0.4332823531003669
[4.99011753e-03 3.32467891e-02 8.81752215e-02 1.70463893e-01
2.80811313e-01 4.14752072e-01 5.71103821e-01 7.53593653e-01
9.79938915e-01 1.14003837e+00 1.40442848e+00 1.66899183e+00
1.96461415e+00 2.29252666e+00 2.63050114e+00 2.98443218e+00
3.38439528e+00 3.81181747e+00 4.26309942e+00 4.69832271e+00
5.22864462e+00 5.74498014e+00 6.22743988e+00 6.83904055e+00
7.42379697e+00 7.97206446e+00 8.62281827e+00 9.26615266e+00
9.85483434e+00 1.05915030e+01 1.11986296e+01 1.18934953e+01
1.26811461e+01 1.33727614e+01 1.41794599e+01 1.47585155e+01
1.55702295e+01 1.63066947e+01 1.71564622e+01 1.78260727e+01
1.85693454e+01 1.95125277e+01 2.01847294e+01 2.09302671e+01
2.18860389e+01 2.25424795e+01 2.32907153e+01 2.37425085e+01
2.50784800e+01 2.55119112e+01]
Jedynym sposobem, w jaki znalazłem użycie sigmy i uzyskanie wyniku w przyzwoitym czasie, jest zapewnienie M jako LinearOperator. Nie jestem zbyt zaznajomiony z tą rzeczą, ale z tego, co zrozumiałem, moja implementacja reprezentuje macierz tożsamości, która jest tym, czym powinno być M, jeśli nie zostało określone w wywołaniu. Powodem tego jest to, że zamiast wykonywać bezpośrednie rozwiązanie (dekompozycję LU), scipy użyje iteracyjnego solwera, który jest potencjalnie lepiej dopasowany. Dla porównania, jeśli podasz M = np.identity(a.shape[0]), które powinno być dokładnie takie samo, wynik eigsh zajmie wieki. Zauważ, że ta metoda nie działa, jeśli podano sigma=0. Ale nie jestem pewien, czy sigma=0 jest naprawdę przydatne?
import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigs, eigsh, LinearOperator
from time import perf_counter
def mv(v):
return v
M = np.load('M.npz')
a = csr_matrix((M['data'], M['indices'], M['indptr']), shape=M['shape'])
t = perf_counter()
b, c = eigsh(a, M=LinearOperator(shape=a.shape, matvec=mv, dtype=np.float64),
sigma=5, k=50, which='SA', tol=1e1, mode='cayley')
print(perf_counter() - t)
print(np.sort(-5 * (1 + b) / (1 - b)))
Ponownie, nie mam pojęcia, czy jest poprawne, ale zdecydowanie różni się od poprzedniego. Byłoby wspaniale mieć wkład kogoś ze scipy.
1.4079377939924598
[3.34420263 3.47938816 3.53019328 3.57981026 3.60457277 3.63996294
3.66791416 3.68391585 3.69223712 3.7082205 3.7496456 3.76170023
3.76923989 3.80811939 3.81337342 3.82848729 3.84137264 3.85648208
3.88110869 3.91286153 3.9271108 3.94444577 3.97580798 3.98868207
4.01677424 4.04341426 4.05915855 4.08910692 4.12238969 4.15283192
4.16871081 4.1990492 4.21792125 4.24509036 4.26892806 4.29603036
4.32282475 4.35839271 4.37934257 4.40343219 4.42782208 4.4477206
4.47635849 4.51594603 4.54294049 4.56689989 4.58804775 4.59919363
4.63700551 4.66638214]
Wiem, że to już stare, ale miałem ten sam problem. Czy sprawdziłeś tutaj (https://docs.scipy.org/doc /scipy/reference/tutorial/arpack.html)?
Wygląda na to, że kiedy ustawisz sigma na niską liczbę (0), powinieneś ustawić which='LM', nawet jeśli chcesz mieć niskie wartości. Dzieje się tak dlatego, że ustawienie sigma przekształca żądane wartości (w tym przypadku niskie) tak, aby wydawały się wysokie, dzięki czemu nadal możesz korzystać z metod „LM”, które są znacznie szybsze, aby uzyskać to, czego chcesz (niskie wartości własne ).
Dodano 19 maja: Wewnętrzny solwer Cholesky'ego:
dokument dotyczący scipy eigsh mówi
tryb odwrócenia przesunięcia… wymaga od operatora obliczenia rozwiązania układu liniowego
(A - sigma * I) x = b... Jest to obliczane wewnętrznie przez rzadką dekompozycję LU (splu) dla wyraźnej macierzy, lub przez iteracyjny solwer dla ogólnego operatora liniowego.
ARPACK wywołuje ten "wewnętrzny rozwiązywacz" wiele razy, w zależności od tol itd.; oczywiście wolny wewnętrzny solwer => wolny eigs. Dla A poz., sksparse.cholmod jest o wiele szybciej niż splu.
Eig Matlab używa również cholesky:
Jeśli A jest hermitowskie, a B jest hermitowskie, dodatnio określone, to domyślną wartością algorytmu jest „chol”
Fwiw, np.linalg.eigh znajduje wszystkie wartości własne i wektory własne gęstej macierzy 7 GB A.A w niecałą godzinę na moim starym 4-rdzeniowym imac'u -- wow. Jego widmo wygląda tak:
Luty 2020, TL;DR
Domysł i kilka uwag, skoro nie mam Matlaba / Octave:
Aby znaleźć małe wartości własne macierzy symetrycznych o wartościach własnych >= 0, takich jak twoja, poniższa metoda jest o wiele szybsza niż odwrócenie przesunięcia:
# flip eigenvalues e.g.
# A: 0 0 0 ... 200 463
# Aflip: 0 163 ... 463 463 463
maxeval = eigsh( A, k=1 )[0] # biggest, fast
Aflip = maxeval * sparse.eye(n) - A
bigevals, evecs = eigsh( Aflip, which="LM", sigma=None ... ) # biggest, near 463
evals = maxeval - bigevals # flip back, near 463 -> near 0
# evecs are the same
eigsh( Aflip ) dla dużych par własnych jest szybszy niż odwrócenie przesunięcia dla małych, ponieważ A * x jest szybszy niż solve(), który musi wykonać odwrócenie przesunięcia. Matlab / Octave mógłby zrobić to automatycznie Aflip, po szybkim teście na wynik pozytywny z Choleskim.
Czy możesz uruchomić eigsh( Aflip ) w Matlab / Octave?
Inne czynniki, które mogą wpływać na dokładność / szybkość:
Domyślnie Arpacka dla wektora początkowego v0 jest wektorem losowym. Używam v0 = np.ones(n), co może być straszne dla niektórych A, ale można je odtworzyć :)
Ta macierz A jest prawie dokładnie osobliwa, A * ones ~ 0.
Multicore: scipy-arpack z openblas / Lapack używa ~ 3,9 z 4 rdzeni na moim iMacu; czy Matlab/Octawa używają wszystkich rdzeni?
k
tol
Gist.github
k 10 tol 1e-05: 8 sec eigvals [0 8.5e-05 0.00043 0.0014 0.0026 0.0047 0.0071 0.0097 0.013 0.018]
k 10 tol 1e-06: 44 sec eigvals [0 3.4e-06 2.8e-05 8.1e-05 0.00015 0.00025 0.00044 0.00058 0.00079 0.0011]
k 10 tol 1e-07: 348 sec eigvals [0 3e-10 7.5e-07 7.6e-06 1.2e-05 1.9e-05 2.1e-05 4.2e-05 5.7e-05 6.4e-05]
k 20 tol 1e-05: 18 sec eigvals [0 5.1e-06 4.5e-05 0.00014 0.00023 0.00042 0.00056 0.00079 0.0011 0.0015 0.0017 0.0021 0.0026 0.003 0.0037 0.0042 0.0047 0.0054 0.006
k 20 tol 1e-06: 73 sec eigvals [0 5.5e-07 7.4e-06 2e-05 3.5e-05 5.1e-05 6.8e-05 0.00011 0.00014 0.00016 0.0002 0.00025 0.00027 0.0004 0.00045 0.00051 0.00057 0.00066
k 20 tol 1e-07: 267 sec eigvals [-4.8e-11 0 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 5.8e-05 6.4e-05 6.9e-05 8.3e-05 0.00011 0.00012 0.00013 0.00015
k 50 tol 1e-05: 82 sec eigvals [-4e-13 9.7e-07 1e-05 2.8e-05 5.9e-05 0.00011 0.00015 0.00019 0.00026 0.00039 ... 0.0079 0.0083 0.0087 0.0092 0.0096 0.01 0.011 0.011 0.012
k 50 tol 1e-06: 432 sec eigvals [-1.4e-11 -4e-13 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 ... 0.00081 0.00087 0.00089 0.00096 0.001 0.001 0.0011 0.0011
k 50 tol 1e-07: 3711 sec eigvals [-5.2e-10 -4e-13 7.5e-07 7.6e-06 1e-05 1.9e-05 2e-05 2.2e-05 4.2e-05 5.1e-05 ... 0.00058 0.0006 0.00063 0.00066 0.00069 0.00071 0.00075
versions: numpy 1.18.1 scipy 1.4.1 umfpack 0.3.2 python 3.7.6 mac 10.10.5
Czy Matlab / Octave o tym samym? Jeśli nie, wszystkie zakłady są wykluczone – najpierw sprawdź poprawność, a potem szybkość.
Dlaczego wartości własne tak bardzo się chwieją? Tiny < 0 dla rzekomo nieujemnej określonej macierzy to znak błąd zaokrąglenia, ale zwykła sztuczka małej zmiany, A += n * eps * sparse.eye(n), nie pomaga.
A
A
Mam nadzieję że to pomoże.