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

12

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

Spacekiller23
źródło
2
1 - Zakładam, że chciałeś umieścić nawiasy kwadratowe wokół [evals, evecs] w kodzie oktawowym? 2 - czy możesz podać mały przykład dla M? a może skrypt generatora dla jednego, jeśli to możliwe?
Nick J
1 - Tak dokładnie, zredagowałem swój post. 2 - Wypróbowałem wydajność niektórych podmacierzy moich danych i wydaje się, że Octave jest zawsze szybszy, ale w przypadku mniejszych macierzy poniżej 5000 x 5000 to tylko 2-5 razy, ponad to robi się naprawdę brzydki. A ponieważ są to „prawdziwe dane”, nie mogę podać skryptu generatora. Czy istnieje jakiś standardowy sposób przesłania przykładu? Z powodu rzadkości plik npz jest dość mały.
Spacekiller23,
Wydaje mi się, że możesz udostępnić link do dowolnego miejsca w chmurze.
Patol75,
Dzięki. W oryginalnym poście umieściłem link Dropbox i zaktualizowałem kod do działającego przykładu.
Spacekiller23
1
Po prostu, aby wzmocnić swój punkt, przetestowałem na Matlabie R2019b i dostałem 84 sekundy w porównaniu z 36,5 minutami w Python 3.7, Scipy 1.2.1 (26 razy szybciej).
Bill

Odpowiedzi:

1

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?


Oto wartości własne scipy-Arpack dla kilku ki tol, grep z logów pod 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 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ę Abierze, z jakiego obszaru problemowego? Czy możesz wygenerować podobny A, mniejszy lub rzadszy?

Mam nadzieję że to pomoże.

denis
źródło
Dziękujemy za wkład i przepraszamy za (bardzo) spóźnioną odpowiedź. Projekt, z którego korzystałem, jest już ukończony, ale nadal jestem ciekawy, więc sprawdziłem. Niestety wartości własne w Ocatve okazują się inne, dla k = 10 znajduję [-2.5673e-16 -1.2239e-18 7.5420e-07 7.5622e-06 1.0189e-05 1.8725e-05 2.0265e-05 2.1568e- 05 4.2458e-05 5.1030e-05], który jest również niezależny od wartości tolerancji w zakresie 1e-5 do 1e-7. Jest więc inna różnica. Czy nie sądzisz, że to dziwne, że scipy (w tym twoja sugestia) dają różne małe wartości zależne od liczby zapytanych wartości?
Spacekiller23
@ Spacekiller23, to był błąd, teraz naprawiony w scipy 1.4.1 (patrz scipy / Issues / 11198 ); czy mógłbyś sprawdzić swoją wersję? tolJest również bałagan dla małych wartości własnych - zadaj nowe pytanie, jeśli chcesz, daj mi znać.
denis
1

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

Anthony Gatti
źródło
Czy to rzeczywiście zmieniło wydajność? To byłaby dla mnie niespodzianka. Znałem link, który podałeś, i domyślnie określiłem, co = „LM” w moim przykładzie. Ponieważ wartość domyślna dla unset, czyli „LM”. Nadal jednak sprawdzałem, a wydajność nie zmieniła się dla mojego przykładu.
Spacekiller23
Rzeczywiście wydaje się, że ma dla ciebie podobną różnicę jak w Pythonie do oktawy. Miałem też dużą macierz, którą próbowałem rozłożyć, i skończyłem na eigsh (macierz, k = 7, co = „LM”, sigma = 1e-10). Początkowo niepoprawnie określałem, które = „SM”, myśląc, że muszę to zrobić, aby uzyskać najmniejsze wartości własne, a udzielenie odpowiedzi zajęło mi wieczność. Potem znalazłem ten artykuł i zdałem sobie sprawę, że po prostu musisz go sprecyzować na szybsze „LM” i ustawić k na dowolne, co chcesz, a to przyspieszy. Czy twoja matryca jest rzeczywiście pustelnikiem?
Anthony Gatti
0

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]
Patol75
źródło
Dziękujemy za Twój wkład i opinie. Próbowałem kilku rzeczy, aby dać przyzwoitą odpowiedź na twoje pytania. 1. Moje zadanie polega na znalezieniu k najmniejszych wartości / wektorów własnych. Dlatego podejście wykorzystujące sigma = 0 podano nawet w dokumentach SciPy: docs.scipy.org/doc/scipy/reference/tutorial/arpack.html 2. Wypróbowałem jeszcze kilka opcji, które zredagowałem w pierwotnym pytaniu. 3. Jak rozumiem filmy dokumentalne Octave i SciPy, eigs (M, 6,0) i k = 6, simga = 0 powinny być takie same.
Spacekiller23,
4. Ponieważ moja macierz jest prawdziwa i kwadratowa, pomyślałem, że nie powinno być różnicy między SA i SM jako opcją, ale oczywiście jest, przynajmniej w obliczeniach. Czy jestem na złym torze? Ogólnie rzecz biorąc, oznacza to więcej pytań, ale nie ma ode mnie prawdziwych odpowiedzi ani rozwiązań.
Spacekiller23,