Znajdowanie najmniejszych wektorów własnych dużej rzadkiej macierzy, ponad 100 razy wolniejszych w SciPy niż w Octave

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ć.

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:

# 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ż shift-invert dla małych, ponieważ A * xjest szybszy niż to, solve()co musi zrobić inwersja shift . Możliwe, że Matlab / Octave zrobi to Aflipautomatycznie, 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 v0jest wektor losowy. Używam v0 = np.ones(n), co dla niektórych może być okropne, Aale jest powtarzalne :)

Ta Amatryca 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?

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 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.

Mam nadzieję że to pomoże.

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 ).

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 odpowiada k=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.

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 ma to 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 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śli sigma=0jest zapewnione. Ale nie jestem pewien, czy sigma=0to naprawdę jest 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)))

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.

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]