scipy.sparse.linalg.

lobpcg#

scipy.sparse.linalg.lobpcg(A, X, B=None, M=None, Y=None, tol=None, maxiter=None, largest=True, verbosityLevel=0, retLambdaHistory=False, retResidualNormsHistory=False, restartControl=20)[Quelle]#

Locally Optimal Block Preconditioned Conjugate Gradient Method (LOBPCG).

LOBPCG ist ein vorkonditionierter Eigenwertlöser für große reelle symmetrische und komplexe Hermitesche definite verallgemeinerte Eigenwertprobleme.

Parameter:
A{sparse matrix, ndarray, LinearOperator, callable object}

Der Hermitesche lineare Operator des Problems, üblicherweise gegeben durch eine spärliche Matrix. Oft als "Steifigkeitsmatrix" bezeichnet.

Xndarray, float32 oder float64

Anfängliche Annäherung an die k Eigenvektoren (nicht spärlich). Wenn A die Form shape=(n,n) hat, muss X die Form shape=(n,k) haben.

B{sparse matrix, ndarray, LinearOperator, callable object}

Optional. Standardmäßig ist B = None, was äquivalent zur Identität ist. Der Operator der rechten Seite in einem verallgemeinerten Eigenwertproblem, falls vorhanden. Oft als "Massenmatrix" bezeichnet. Muss Hermitesch positiv definit sein.

M{sparse matrix, ndarray, LinearOperator, callable object}

Optional. Standardmäßig ist M = None, was äquivalent zur Identität ist. Präkonditionierer zur Beschleunigung der Konvergenz.

Yndarray, float32 oder float64, Standard: None

Ein n-by-sizeY ndarray von Nebenbedingungen mit sizeY < n. Die Iterationen werden im B-orthogonalen Komplement des Spaltenraums von Y durchgeführt. Y muss vollen Rang haben, falls vorhanden.

tolSkalar, optional

Standardmäßig tol=n*sqrt(eps). Solver-Toleranz für das Abbruchkriterium.

maxiterint, Standard: 20

Maximale Anzahl von Iterationen.

largestbool, Standard: True

Wenn True, werden die größten Eigenwerte berechnet, andernfalls die kleinsten.

verbosityLevelint, optional

Standardmäßig verbosityLevel=0, keine Ausgabe. Steuert die Standard-/Bildschirmausgabe des Solvers.

retLambdaHistorybool, Standard: False

Ob die iterative Eigenwertverlauf zurückgegeben werden soll.

retResidualNormsHistorybool, Standard: False

Ob der iterative Verlauf der Residuen-Normen zurückgegeben werden soll.

restartControlint, optional.

Iterationen werden neu gestartet, wenn die Residuen 2**restartControl Mal größer sind als die kleinsten aufgezeichneten in retResidualNormsHistory. Der Standardwert ist restartControl=20, was die Neustarts zur Rückwärtskompatibilität selten macht.

Rückgabe:
lambdandarray der Form (k, ).

Array mit k ungefähren Eigenwerten.

vndarray der gleichen Form wie X.shape.

Ein Array mit k ungefähren Eigenvektoren.

lambdaHistoryndarray, optional.

Der Eigenwertverlauf, wenn retLambdaHistory True ist.

ResidualNormsHistoryndarray, optional.

Der Verlauf der Residuen-Normen, wenn retResidualNormsHistory True ist.

Hinweise

Die iterative Schleife läuft höchstens maxit=maxiter (20, wenn maxit=None) Iterationen und endet früher, wenn die Toleranz erreicht ist. LOBPCG bricht die Rückwärtskompatibilität mit der vorherigen Version, da LOBPCG nun den Block der iterativen Vektoren mit der besten Genauigkeit anstelle des zuletzt iterierten zurückgibt, als Heilmittel für mögliche Divergenz.

Wenn X.dtype == np.float32 und benutzerdefinierte Operationen/Multiplikationen mit A, B und M den np.float32 Datentyp beibehalten, erfolgen alle Berechnungen und die Ausgabe in np.float32.

Die Größe des Iterationsverlaufs entspricht der Anzahl der besten (begrenzt durch maxit) Iterationen plus 3: Anfang, Ende und Nachbearbeitung.

Wenn sowohl retLambdaHistory als auch retResidualNormsHistory True sind, hat das Rückgabebündel das folgende Format (lambda, V, lambda history, residual norms history).

Im Folgenden bezeichnet n die Matrixgröße und k die Anzahl der benötigten Eigenwerte (kleinsten oder größten).

Der LOBPCG-Code löst intern Eigenwertprobleme der Größe 3k in jeder Iteration, indem er den dichten Eigenwertlöser eigh aufruft. Wenn k also nicht klein genug im Vergleich zu n ist, macht es keinen Sinn, den LOBPCG-Code aufzurufen. Darüber hinaus, wenn man den LOBPCG-Algorithmus für 5k > n aufruft, würde er intern wahrscheinlich fehlschlagen, sodass der Code stattdessen die Standardfunktion eigh aufruft. Es ist nicht so, dass n groß sein muss, damit LOBPCG funktioniert, sondern vielmehr, dass das Verhältnis n / k groß sein sollte. Wenn Sie LOBPCG mit k=1 und n=10 aufrufen, funktioniert es, obwohl n klein ist. Die Methode ist für extrem große n / k gedacht.

Die Konvergenzgeschwindigkeit hängt im Wesentlichen von drei Faktoren ab

  1. Qualität der anfänglichen Annäherungen X an die gesuchten Eigenvektoren. Zufällig um den Ursprung verteilte Vektoren funktionieren gut, wenn keine bessere Wahl bekannt ist.

  2. Relative Trennung der gewünschten Eigenwerte von den restlichen Eigenwerten. Man kann k variieren, um die Trennung zu verbessern.

  3. Geeignete Vorkonditionierung zur Verringerung der spektralen Ausdehnung. Zum Beispiel ist ein Testproblem für die Vibration einer Stange (unter Verzeichnis tests) für große n schlecht konditioniert, sodass die Konvergenz langsam sein wird, es sei denn, es wird eine effiziente Vorkonditionierung verwendet. Für dieses spezielle Problem wäre eine gute einfache Präkonditionierungsfunktion eine lineare Lösung für A, die einfach zu codieren ist, da A tridiagonal ist.

Referenzen

[1]

A. V. Knyazev (2001), Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method. SIAM Journal on Scientific Computing 23, Nr. 2, S. 517-541. DOI:10.1137/S1064827500366124

[2]

A. V. Knyazev, I. Lashuk, M. E. Argentati, und E. Ovchinnikov (2007), Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) in hypre und PETSc. arXiv:0705.2626

[3]

A. V. Knyazevs C- und MATLAB-Implementierungen: lobpcg/blopex

Beispiele

Unser erstes Beispiel ist minimalistisch - finden Sie den größten Eigenwert einer Diagonalmatrix, indem Sie das nicht-verallgemeinerte Eigenwertproblem A x = lambda x ohne Nebenbedingungen oder Vorkonditionierung lösen.

>>> import numpy as np
>>> from scipy.sparse import diags_array
>>> from scipy.sparse.linalg import LinearOperator, aslinearoperator
>>> from scipy.sparse.linalg import lobpcg

Die Größe der quadratischen Matrix ist

>>> n = 100

und ihre Diagonaleinträge sind 1, ..., 100, definiert durch

>>> vals = np.arange(1, n + 1).astype(np.int16)

Der erste obligatorische Eingabeparameter in diesem Test ist die spärliche Diagonalmatrix A des Eigenwertproblems A x = lambda x, das gelöst werden soll.

>>> A = diags_array(vals, offsets=0, shape=(n, n))
>>> A = A.astype(np.int16)
>>> A.toarray()
array([[  1,   0,   0, ...,   0,   0,   0],
       [  0,   2,   0, ...,   0,   0,   0],
       [  0,   0,   3, ...,   0,   0,   0],
       ...,
       [  0,   0,   0, ...,  98,   0,   0],
       [  0,   0,   0, ...,   0,  99,   0],
       [  0,   0,   0, ...,   0,   0, 100]], shape=(100, 100), dtype=int16)

Der zweite obligatorische Eingabeparameter X ist ein 2D-Array, dessen Zeilendimension die Anzahl der angeforderten Eigenwerte bestimmt. X ist eine anfängliche Schätzung für die Ziel-Eigenvektoren. X muss linear unabhängige Spalten haben. Wenn keine anfänglichen Annäherungen verfügbar sind, funktionieren zufällig orientierte Vektoren normalerweise am besten, z. B. mit Komponenten, die normal um Null verteilt sind oder gleichmäßig auf dem Intervall [-1 1] verteilt sind. Das Festlegen der anfänglichen Annäherungen auf den Datentyp np.float32 erzwingt, dass alle iterativen Werte den Datentyp np.float32 haben, was die Ausführung beschleunigt und dennoch genaue Eigenwertberechnungen ermöglicht.

>>> k = 1
>>> rng = np.random.default_rng()
>>> X = rng.normal(size=(n, k))
>>> X = X.astype(np.float32)

>>> eigenvalues, _ = lobpcg(A, X, maxiter=60)
>>> eigenvalues
array([100.], dtype=float32)

lobpcg benötigt nur Zugriff auf das Matrixprodukt mit A und nicht auf die Matrix selbst. Da die Matrix A in diesem Beispiel diagonal ist, kann man eine Funktion des Matrixprodukts A @ X unter Verwendung der Diagonalwerte vals schreiben, z. B. durch elementweise Multiplikation mit Broadcasting in der Lambda-Funktion

>>> A_lambda = lambda X: vals[:, np.newaxis] * X

oder der regulären Funktion

>>> def A_matmat(X):
...     return vals[:, np.newaxis] * X

und den Handle zu einer dieser aufrufbaren Funktionen als Eingabe verwenden

>>> eigenvalues, _ = lobpcg(A_lambda, X, maxiter=60)
>>> eigenvalues
array([100.], dtype=float32)
>>> eigenvalues, _ = lobpcg(A_matmat, X, maxiter=60)
>>> eigenvalues
array([100.], dtype=float32)

Der traditionelle aufrufbare LinearOperator ist nicht mehr notwendig, wird aber immer noch als Eingabe für lobpcg unterstützt. Die explizite Angabe von matmat=A_matmat verbessert die Leistung.

>>> A_lo = LinearOperator((n, n), matvec=A_matmat, matmat=A_matmat, dtype=np.int16)
>>> eigenvalues, _ = lobpcg(A_lo, X, maxiter=80)
>>> eigenvalues
array([100.], dtype=float32)

Die am wenigsten effiziente aufrufbare Option ist aslinearoperator

>>> eigenvalues, _ = lobpcg(aslinearoperator(A), X, maxiter=80)
>>> eigenvalues
array([100.], dtype=float32)

Wir wechseln nun zur Berechnung der drei kleinsten Eigenwerte und geben

>>> k = 3
>>> X = np.random.default_rng().normal(size=(n, k))

und den Parameter largest=False an

>>> eigenvalues, _ = lobpcg(A, X, largest=False, maxiter=90)
>>> print(eigenvalues)  
[1. 2. 3.]

Das nächste Beispiel illustriert die Berechnung von 3 kleinsten Eigenwerten derselben Matrix A, die durch den Funktionshandle A_matmat gegeben ist, jedoch mit Nebenbedingungen und Vorkonditionierung.

Nebenbedingungen - ein optionaler Eingabeparameter ist ein 2D-Array, das Spaltenvektoren enthält, zu denen die Eigenvektoren orthogonal sein müssen

>>> Y = np.eye(n, 3)

Der Präkonditionierer wirkt in diesem Beispiel als Invers von A, jedoch in reduzierter Präzision np.float32, obwohl das anfängliche X und damit alle Iterationen und die Ausgabe in vollem np.float64 erfolgen.

>>> inv_vals = 1./vals
>>> inv_vals = inv_vals.astype(np.float32)
>>> M = lambda X: inv_vals[:, np.newaxis] * X

Lassen Sie uns nun das Eigenwertproblem für die Matrix A zuerst ohne Vorkonditionierung lösen und 80 Iterationen anfordern

>>> eigenvalues, _ = lobpcg(A_matmat, X, Y=Y, largest=False, maxiter=80)
>>> eigenvalues
array([4., 5., 6.])
>>> eigenvalues.dtype
dtype('float64')

Mit Vorkonditionierung benötigen wir nur 20 Iterationen von demselben X

>>> eigenvalues, _ = lobpcg(A_matmat, X, Y=Y, M=M, largest=False, maxiter=20)
>>> eigenvalues
array([4., 5., 6.])

Beachten Sie, dass die in Y übergebenen Vektoren die Eigenvektoren der 3 kleinsten Eigenwerte sind. Die oben zurückgegebenen Ergebnisse sind orthogonal zu diesen.

Die primäre Matrix A kann indefinit sein, z. B. nachdem vals um 50 von 1, ..., 100 auf -49, ..., 50 verschoben wurde, können wir immer noch die 3 kleinsten oder größten Eigenwerte berechnen.

>>> vals = vals - 50
>>> X = rng.normal(size=(n, k))
>>> eigenvalues, _ = lobpcg(A_matmat, X, largest=False, maxiter=99)
>>> eigenvalues
array([-49., -48., -47.])
>>> eigenvalues, _ = lobpcg(A_matmat, X, largest=True, maxiter=99)
>>> eigenvalues
array([50., 49., 48.])
Auf dieser Seite