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.
quelle
Antworten:
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:
eigsh( Aflip )
für große Eigenpaare ist schneller als Shift-Invert für kleine, weilA * x
es schneller ist als dassolve()
, was Shift-Invert tun muss. Matlab / Octave könnten dies möglicherweiseAflip
automatisch 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
v0
ist ein Zufallsvektor. Ich benutzev0 = np.ones(n)
, was für manche schrecklich sein mag,A
aber reproduzierbar ist :)Diese
A
Matrix 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
k
undtol
aus Protokolldateien unter gist.github :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 ähnlicheA
, kleinere oder sparsamere generieren ?Hoffe das hilft.
quelle
tol
ist etwas chaotisch für kleine Eigenwerte - fragen Sie eine neue Frage auf , dass , wenn Sie möchten, lassen Sie es mich wissen.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) ).
quelle
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 entsprichtk=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.
Ich bin mir jedoch nicht sicher, ob dies sinnvoll ist.
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 ersigma=0
bereitgestellt wird. Aber ich bin mir nicht sicher, obsigma=0
das wirklich so nützlich ist?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.
quelle