Finden kleinster Eigenvektoren mit großer Matrix mit geringer Dichte, die in SciPy über 100-mal langsamer sind als in Octave

12

Ich versuche, wenige (5-500) Eigenvektoren zu berechnen, die den kleinsten Eigenwerten großer symmetrischer quadratischer Sparse-Matrizen (bis zu 30000 x 30000) entsprechen, wobei weniger als 0,1% der Werte ungleich Null sind.

Ich verwende derzeit scipy.sparse.linalg.eigsh im Shift-Invert-Modus (Sigma = 0.0), was ich durch verschiedene Beiträge zum Thema herausgefunden habe. Dies ist die bevorzugte Lösung. In den meisten Fällen dauert es jedoch bis zu 1 Stunde, um das Problem zu lösen. Andererseits ist die Funktion sehr schnell, wenn ich nach den größten Eigenwerten (Subsekunden auf meinem System) frage, die aus der Dokumentation erwartet wurden.

Da ich mit Matlab von der Arbeit aus besser vertraut bin, habe ich versucht, das Problem in Octave zu lösen, was mir das gleiche Ergebnis mit Eigs (Sigma = 0) in nur wenigen Sekunden (unter 10 Sekunden) brachte. Da ich einen Parameter-Sweep des Algorithmus einschließlich der Eigenvektorberechnung durchführen möchte, wäre diese Art von Zeitgewinn auch in Python großartig.

Ich habe zuerst die Parameter (insbesondere die Toleranz) geändert, aber das hat sich auf der Zeitskala nicht viel geändert.

Ich verwende Anaconda unter Windows, habe aber versucht, das von scipy verwendete LAPACK / BLAS (was ein großer Schmerz war) von mkl (Standard-Anaconda) auf OpenBlas (laut Dokumentation von Octave verwendet) umzustellen, konnte jedoch keine Änderung in feststellen Performance.

Ich konnte nicht herausfinden, ob (und wie) etwas an dem verwendeten ARPACK zu ändern war.

Ich habe einen Testfall für den folgenden Code in den folgenden Dropbox-Ordner hochgeladen: https://www.dropbox.com/sh/l6aa6izufzyzqr3/AABqij95hZOvRpnnjRaETQmka?dl=0

In Python

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)

In der Oktave:

M=dlmread('M.txt');
M=spconvert(M);
[evecs,evals] = eigs(M,6,0);

Jede Hilfe wird geschätzt!

Einige zusätzliche Optionen, die ich aufgrund der Kommentare und Vorschläge ausprobiert habe:

Oktave: eigs(M,6,0)und eigs(M,6,'sm')gib mir das gleiche Ergebnis:

[1.8725e-05 1.0189e-05 7.5622e-06 7.5420e-07 -1.2239e-18 -2.5674e-16]

während eigs(M,6,'sa',struct('tol',2))konvergiert zu

[1.0423 2.7604 6.1548 11.1310 18.0207 25.3933] 

viel schneller, aber nur, wenn die Toleranzwerte über 2 liegen, sonst konvergiert es überhaupt nicht und die Werte sind stark unterschiedlich.

Python: eigsh(M,k=6,which='SA')und eigsh(M,k=6,which='SM')beide konvergieren nicht (ARPACK-Fehler, wenn keine Konvergenz erreicht wurde). Nur eigsh(M,k=6,sigma=0.0)gibt einige Eigenwerte (nach fast einer Stunde), die für die Kleinsten Oktave verschieden ist (sogar 1 zusätzlicher kleiner Wert gefunden wird ):

[3.82923317e-17 3.32269886e-16 2.78039665e-10 7.54202273e-07 7.56251500e-06 1.01893934e-05]

Wenn die Toleranz hoch genug ist, erhalte ich auch Ergebnisse eigsh(M,k=6,which='SA',tol='1'), die den anderen erhaltenen Werten nahe kommen

[4.28732218e-14 7.54194948e-07 7.56220703e-06 1.01889544e-05, 1.87247350e-05 2.02652719e-05]

wieder mit einer anderen Anzahl kleiner Eigenwerte. Die Rechenzeit beträgt noch fast 30min. Während die verschiedenen sehr kleinen Werte verständlich sein mögen, da sie ein Vielfaches von 0 darstellen könnten, verwirrt mich die unterschiedliche Vielfalt.

Außerdem scheint es bei SciPy und Octave einige grundlegende Unterschiede zu geben, die ich noch nicht herausfinden kann.

Spacekiller23
quelle
2
1 - Ich nehme an, Sie wollten im Oktavcode Klammern um [evals, evecs] setzen? 2 - können Sie ein kleines Beispiel für M einfügen? oder vielleicht ein Generatorskript für eines, wenn das möglich ist?
Nick J
1 - Ja genau, ich habe meinen Beitrag bearbeitet. 2 - Ich habe die Leistung für einige Untermatrizen meiner Daten ausprobiert und es scheint, dass Octave immer schneller ist, aber für kleinere Matrizen unter 5000x5000 ist es nur ein Faktor von 2-5 Mal, darüber hinaus wird es wirklich hässlich. Und da es "echte Daten" sind, kann ich kein Generatorskript geben. Gibt es eine Standardmethode, um ein Beispiel irgendwie hochzuladen? Aufgrund der Sparsamkeit ist eine npz-Datei relativ klein.
Spacekiller23
Ich denke, Sie können einen Link zu jedem Cloud-Speicher freigeben.
Patol75
Danke. Ich habe einen Dropbox-Link in den ursprünglichen Beitrag eingefügt und den Code auf ein funktionierendes Beispiel aktualisiert.
Spacekiller23
1
Um Ihren Standpunkt zu bekräftigen, habe ich auf Matlab R2019b getestet und 84 Sekunden gegenüber 36,5 Minuten in Python 3.7, Scipy 1.2.1 (26-mal schneller) erhalten.
Bill

Antworten:

1

Eine Vermutung und einige Kommentare, da ich Matlab / Octave nicht habe:

Um kleine Eigenwerte von symmetrischen Matrizen mit Eigenwerten> = 0 wie Ihre zu finden, ist Folgendes schneller als Shift-Invert:

# 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 )für große Eigenpaare ist schneller als Shift-Invert für kleine, weil A * xes schneller ist als das solve(), was Shift-Invert tun muss. Matlab / Octave könnten dies möglicherweise Aflipautomatisch tun , nach einem kurzen Test auf Positiv-Definitiv mit Cholesky.
Kannst du eigsh( Aflip )in Matlab / Octave laufen ?

Andere Faktoren, die die Genauigkeit / Geschwindigkeit beeinflussen können:

Arpacks Standard für den Startvektor v0ist ein Zufallsvektor. Ich benutze v0 = np.ones(n), was für manche schrecklich sein mag, Aaber reproduzierbar ist :)

Diese AMatrix ist fast genau sigular, A * ones~ 0.

Multicore: scipy-arpack mit openblas / Lapack verwendet ~ 3,9 der 4 Kerne auf meinem iMac; Verwenden Matlab / Octave alle Kerne?


Hier sind die scipy-Arpack-Eigenwerte für mehrere kund tolaus Protokolldateien unter 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 

Sind Matlab / Octave ungefähr gleich? Wenn nicht, sind alle Wetten geschlossen - überprüfen Sie zuerst die Richtigkeit und dann die Geschwindigkeit.

Warum wackeln die Eigenwerte so stark? Winzig <0 für eine angeblich nicht negativ-definierte Matrix sind ein Zeichen für einen Rundungsfehler , aber der übliche Trick einer winzigen Verschiebung A += n * eps * sparse.eye(n)hilft nicht weiter.


Woher kommt das A, welcher Problembereich? Können Sie ähnliche A, kleinere oder sparsamere generieren ?

Hoffe das hilft.

denis
quelle
Vielen Dank für Ihre Eingabe und entschuldigen Sie die (sehr) späte Antwort. Das Projekt, für das ich das verwendet habe, ist bereits abgeschlossen, aber ich bin immer noch neugierig, also habe ich es überprüft. Leider fallen die Eigenwerte in Ocatve unterschiedlich aus, für k = 10 finde ich [-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], der auch unabhängig vom Toleranzwert im Bereich von 1e-5 bis 1e-7 ist. Hier gibt es also einen weiteren Unterschied. Finden Sie es nicht seltsam, dass scipy (einschließlich Ihres Vorschlags) je nach Anzahl der abgefragten Werte unterschiedliche kleine Werte liefert?
Spacekiller23
@ Spacekiller23, dies war ein Fehler, der jetzt in scipy 1.4.1 behoben wurde (siehe scipy / issue / 11198 ); Könntest du deine Version überprüfen? Auch tolist etwas chaotisch für kleine Eigenwerte - fragen Sie eine neue Frage auf , dass , wenn Sie möchten, lassen Sie es mich wissen.
Denis
1

Ich weiß, dass das jetzt alt ist, aber ich hatte das gleiche Problem. Haben Sie hier eine Bewertung abgegeben ( https://docs.scipy.org/doc/scipy/reference/tutorial/arpack.html )?

Es scheint, als ob Sie, wenn Sie Sigma auf eine niedrige Zahl (0) setzen, welche = 'LM' setzen sollten, obwohl Sie niedrige Werte wollen. Dies liegt daran, dass durch das Setzen von Sigma die gewünschten Werte (in diesem Fall niedrig) als hoch erscheinen und Sie dennoch die LM-Methoden nutzen können, mit denen Sie viel schneller das erreichen, was Sie möchten (die niedrigen Eigenwerte) ).

Anthony Gatti
quelle
Hat dies tatsächlich die Leistung für Sie verändert? Das wäre eine Überraschung für mich. Ich kannte den Link, den Sie postet, und habe implizit angegeben, welcher = 'LM' in meinem Beispiel. Weil der Standardwert für ein Unset 'LM' ist. Ich habe es trotzdem überprüft und die Leistung ist für mein Beispiel unverändert.
Spacekiller23
In der Tat scheint ein ähnlicher Unterschied wie von Python zu Oktave für Sie zu haben. Ich hatte auch eine große Matrix, die ich zu zerlegen versuchte und die schließlich eigsh verwendete (Matrix, k = 7, die = 'LM', Sigma = 1e-10). Ursprünglich habe ich falsch angegeben, welches = 'SM'-Denken ich tun musste, um die kleinsten Eigenwerte zu erhalten, und es dauerte ewig, bis ich eine Antwort erhielt. Dann fand ich diesen Artikel und erkannte, dass Sie ihn nur für das schnellere 'LM' spezifizieren mussten und k so einstellen mussten, wie Sie es wollten, und es würde die Dinge beschleunigen. Ist Ihre Matrix tatsächlich Einsiedler?
Anthony Gatti vor
0

Ich möchte zunächst sagen, dass ich keine Ahnung habe, warum die Ergebnisse, die Sie und @Bill gemeldet haben, so sind, wie sie sind. Ich frage mich einfach, ob eigs(M,6,0)in Octave entspricht k=6 & sigma=0, oder vielleicht ist es etwas anderes?

Wenn ich mit scipy kein Sigma zur Verfügung stelle, kann ich auf diese Weise in angemessener Zeit ein Ergebnis erzielen.

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)

Ich bin mir jedoch nicht sicher, ob dies sinnvoll ist.

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]

Die einzige Möglichkeit, Sigma zu verwenden und in angemessener Zeit ein Ergebnis zu erzielen, besteht darin, M als LinearOperator bereitzustellen. Ich bin mit dieser Sache nicht allzu vertraut, aber nach meinem Verständnis stellt meine Implementierung eine Identitätsmatrix dar, die M sein sollte, wenn sie nicht im Aufruf angegeben wird. Der Grund dafür ist, dass scipy anstelle einer direkten Lösung (LU-Zerlegung) einen iterativen Löser verwendet, der möglicherweise besser geeignet ist. Zum Vergleich: Wenn Sie M = np.identity(a.shape[0])angeben, was genau gleich sein sollte, dauert es ewig, bis eigsh ein Ergebnis liefert. Beachten Sie, dass dieser Ansatz nicht funktioniert, wenn er sigma=0bereitgestellt wird. Aber ich bin mir nicht sicher, ob sigma=0das wirklich so nützlich ist?

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

Wieder keine Ahnung, ob es richtig ist, aber definitiv anders als zuvor. Es wäre großartig, wenn jemand von scipy etwas dazu sagen würde.

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
quelle
Vielen Dank für Ihre Eingabe und Ihr Feedback. Ich habe einige Dinge ausprobiert, um eine anständige Antwort auf Ihre Punkte zu geben. 1. Meine Aufgabe besteht darin, die k kleinsten Eigenwerte / Vektoren zu finden. Daher wird der Ansatz mit sigma = 0 sogar in den SciPy-Dokumenten angegeben: docs.scipy.org/doc/scipy/reference/tutorial/arpack.html 2. Ich habe einige weitere Optionen ausprobiert, die ich in der ursprünglichen Frage bearbeitet habe. 3. Nach meinem Verständnis der Dokumentationen von Octave und SciPy sollten eigs (M, 6,0) und k = 6, simga = 0 gleich sein.
Spacekiller23
4. Da meine Matrix real und quadratisch ist, dachte ich, dass es keinen Unterschied zwischen SA und SM als Option geben sollte, aber dies ist offensichtlich zumindest bei der Berechnung der Fall. Bin ich hier auf einem falschen Weg? Insgesamt bedeutet das mehr Fragen und aber keine wirklichen Antworten oder Lösungen von mir.
Spacekiller23