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

14
Spacekiller23 19 grudzień 2019, 22:50
2
1 - Zakładam, że chciałeś umieścić nawiasy 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
19 grudzień 2019, 23:25
1 - Dokładnie tak, edytowałem mój post. 2 - Wypróbowałem wydajność dla niektórych podmacierzy moich danych i wygląda na to, że Octave jest zawsze szybszy, ale dla mniejszych matryc poniżej 5000x5000 jest to tylko współczynnik 2-5 razy, powyżej jest naprawdę brzydki. A ponieważ jego "prawdziwe dane" nie mogę podać skryptu generatora. Czy jest jakiś standardowy sposób na przesłanie przykładu? Ze względu na rzadkość plik npz jest dość mały.
 – 
Spacekiller23
19 grudzień 2019, 23:53
Myślę, że możesz udostępnić link do dowolnego magazynu w chmurze.
 – 
Patol75
20 grudzień 2019, 02:05
Dzięki. W oryginalnym poście umieściłem link do skrzynki i zaktualizowałem kod do działającego przykładu.
 – 
Spacekiller23
20 grudzień 2019, 03:01
1
Juts, aby wzmocnić twój punkt widzenia, przetestowałem na Matlab R2019b i uzyskałem 84 sekundy w porównaniu z 36,5 minutami w Python 3.7, Scipy 1.2.1 (26 razy szybciej).
 – 
Bill
20 grudzień 2019, 04:07

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]
0
Patol75 20 grudzień 2019, 07:30
Dziękuję za wkład i opinie. Próbowałem kilku rzeczy, aby dać przyzwoitą odpowiedź na twoje punkty. 1. Moje zadanie polega na znalezieniu k najmniejszych wartości własnych/wektorów. Dlatego podejście wykorzystujące sigma=0 jest nawet podane w dokumentacji SciPy: docs .scipy.org/doc/scipy/reference/tutorial/arpack.html 2. Wypróbowałem kilka innych opcji, które zredagowałem do oryginalnego pytania. 3. Jak rozumiem filmy dokumentalne Octave i SciPy, eigs(M,6,0) i k=6, simga=0 powinny być takie same.
 – 
Spacekiller23
20 grudzień 2019, 21:06
4. Ponieważ moja macierz jest rzeczywista 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 tutaj na złej drodze? Ogólnie oznacza to więcej pytań, ale nie mam ode mnie prawdziwych odpowiedzi ani rozwiązań.
 – 
Spacekiller23
20 grudzień 2019, 21:06

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

1
Anthony Gatti 17 kwiecień 2020, 06:20
Czy to rzeczywiście zmieniło dla Ciebie wydajność? To byłaby dla mnie niespodzianka. Znałem link, który zamieściłeś, a także domyślnie określiłem, który = „LM” w moim przykładzie. Ponieważ domyślną wartością dla unset, która jest 'LM'. Jednak nadal sprawdzałem, a wydajność jest niezmieniona dla mojego przykładu.
 – 
Spacekiller23
1 maj 2020, 13:33
Rzeczywiście, wydaje się, że ma dla ciebie podobną różnicę jak z Pythona do oktawy. Miałem też dużą macierz, którą próbowałem rozłożyć i skończyłem używając eigsh(macierz, k=7, która='LM', sigma=1e-10). Początkowo błędnie określałem, który = „SM” myślałem, że muszę to zrobić, aby uzyskać najmniejsze wartości własne, a uzyskanie odpowiedzi zajęło mi wieczność. Potem znalazłem ten artykuł i zdałem sobie sprawę, że wystarczy podać go do szybszego „LM” i ustawić k na dowolne, co przyspieszy sprawę. Czy twoja matryca jest rzeczywiście pustelnikiem?
 – 
Anthony Gatti
5 maj 2020, 04:41
Dobra, ale to inna sprawa. W rzeczywistości, jeśli określisz sigma bliską zeru, wtedy wywołasz transformację, w której LM transformowanej macierzy jest SM macierzy oryginalnej. Tak więc, jak powiedziałem, nie ma różnicy w wydajności między domyślną wartością who i which='LM'. ALE określona sigma zmienia znaczenie której. Zatem which='SM' z małą sigma jest po prostu innym zapytaniem, wyjaśniającym różnicę w wydajności, ale z tego wynikają również nieprawidłowe wartości.
 – 
Spacekiller23
11 maj 2020, 23:01

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:

enter image description here


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.

1
denis 19 maj 2020, 15:57
Dziękuję za Twój wkład i przepraszam za (bardzo) spóźnioną odpowiedź. Projekt, do którego tego użyłem, jest już gotowy, 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óra jest również niezależna od wartości tolerancji w zakresie od 1e-5 do 1e-7. Więc jest tu kolejna różnica. Nie sądzisz, że to dziwne, że scipy (łącznie z twoją sugestią) daje różne małe wartości w zależności od liczby zapytanych wartości?
 – 
Spacekiller23
1 maj 2020, 13:25
Przepraszam, niż nie rozumiem. W oryginalnej odpowiedzi podajesz scipy 1.4.1, a liczba bardzo małych wartości własnych zmienia się wraz z k. A wartości własne, które zamieściłem w mojej odpowiedzi powyżej, są w oktawie. Różnią się one oczywiście od tych, które masz w przypadku scipy. Więc problem nadal występuje, czy czegoś tu brakuje?
 – 
Spacekiller23
11 maj 2020, 22:57