SPARSE Eigenwertprobleme mit ARPACK#
Einleitung#
ARPACK [1] ist ein Fortran-Paket, das Routinen zum schnellen Finden einiger Eigenwerte/Eigenvektoren großer dünnbesetzter Matrizen bereitstellt. Um diese Lösungen zu finden, ist nur eine Links-Multiplikation mit der fraglichen Matrix erforderlich. Diese Operation wird über eine *Reverse-Communication*-Schnittstelle durchgeführt. Das Ergebnis dieser Struktur ist, dass ARPACK Eigenwerte und Eigenvektoren jeder linearen Funktion finden kann, die einen Vektor auf einen Vektor abbildet.
Die gesamte Funktionalität von ARPACK ist in den beiden High-Level-Schnittstellen scipy.sparse.linalg.eigs und scipy.sparse.linalg.eigsh enthalten. eigs bietet Schnittstellen zum Finden der Eigenwerte/Eigenvektoren von reellen oder komplexen nicht-symmetrischen Quadratmatrizen, während eigsh Schnittstellen für reelle symmetrische oder komplexe hermitesche Matrizen bietet.
Grundlegende Funktionalität#
ARPACK kann entweder Standard-Eigenwertprobleme der Form
oder allgemeine Eigenwertprobleme der Form
Die Stärke von ARPACK liegt darin, dass es nur eine ausgewählte Teilmenge von Eigenwert-/Eigenvektorpaaren berechnen kann. Dies wird durch das Schlüsselwort which erreicht. Die folgenden Werte für which sind verfügbar
which = 'LM': Eigenwerte mit größter Magnitude (eigs,eigsh), d.h. größte Eigenwerte in der euklidischen Norm von komplexen Zahlen.which = 'SM': Eigenwerte mit kleinster Magnitude (eigs,eigsh), d.h. kleinste Eigenwerte in der euklidischen Norm von komplexen Zahlen.which = 'LR': Eigenwerte mit größtem Realteil (eigs).which = 'SR': Eigenwerte mit kleinstem Realteil (eigs).which = 'LI': Eigenwerte mit größtem Imaginärteil (eigs).which = 'SI': Eigenwerte mit kleinstem Imaginärteil (eigs).which = 'LA': Eigenwerte mit größtem algebraischem Wert (eigsh), d.h. größte Eigenwerte, einschließlich negativer Vorzeichen.which = 'SA': Eigenwerte mit kleinstem algebraischem Wert (eigsh), d.h. kleinste Eigenwerte, einschließlich negativer Vorzeichen.which = 'BE': Eigenwerte von beiden Enden des Spektrums (eigsh).
Beachten Sie, dass ARPACK im Allgemeinen besser darin ist, Extremaleigenwerte zu finden, d.h. Eigenwerte mit großer Magnitude. Insbesondere die Verwendung von which = 'SM' kann zu langsamer Ausführungszeit und/oder anomalen Ergebnissen führen. Ein besserer Ansatz ist die Verwendung des *Shift-Invert-Modus*.
Shift-Invert-Modus#
Der Shift-Invert-Modus basiert auf der folgenden Beobachtung. Für das verallgemeinerte Eigenwertproblem
kann gezeigt werden, dass
wo
Beispiele#
Stellen Sie sich vor, Sie möchten die kleinsten und größten Eigenwerte und die entsprechenden Eigenvektoren für eine große Matrix finden. ARPACK kann viele Eingabeformen verarbeiten: dichte Matrizen wie numpy.ndarray-Instanzen, dünnbesetzte Matrizen wie scipy.sparse.csr_matrix oder einen allgemeinen linearen Operator, der von scipy.sparse.linalg.LinearOperator abgeleitet ist. Für dieses Beispiel erstellen wir vereinfachend eine symmetrische, positiv definite Matrix.
>>> import numpy as np
>>> from scipy.linalg import eig, eigh
>>> from scipy.sparse.linalg import eigs, eigsh
>>> np.set_printoptions(suppress=True)
>>> rng = np.random.default_rng()
>>>
>>> X = rng.random((100, 100)) - 0.5
>>> X = np.dot(X, X.T) # create a symmetric matrix
Wir haben nun eine symmetrische Matrix X, mit der wir die Routinen testen können. Berechnen Sie zunächst eine Standard-Eigenwertzerlegung mit eigh
>>> evals_all, evecs_all = eigh(X)
Mit wachsender Dimension von X wird diese Routine sehr langsam. Insbesondere wenn nur wenige Eigenvektoren und Eigenwerte benötigt werden, kann ARPACK eine bessere Option sein. Berechnen wir zunächst die größten Eigenwerte (which = 'LM') von X und vergleichen sie mit den bekannten Ergebnissen
>>> evals_large, evecs_large = eigsh(X, 3, which='LM')
>>> print(evals_all[-3:])
[29.22435321 30.05590784 30.58591252]
>>> print(evals_large)
[29.22435321 30.05590784 30.58591252]
>>> print(np.dot(evecs_large.T, evecs_all[:,-3:]))
array([[-1. 0. 0.], # may vary (signs)
[ 0. 1. 0.],
[-0. 0. -1.]])
Die Ergebnisse sind wie erwartet. ARPACK liefert die gewünschten Eigenwerte und sie stimmen mit den zuvor bekannten Ergebnissen überein. Darüber hinaus sind die Eigenvektoren orthogonal, wie wir es erwarten würden. Versuchen wir nun, die Eigenwerte mit der kleinsten Magnitude zu lösen
>>> evals_small, evecs_small = eigsh(X, 3, which='SM')
Traceback (most recent call last): # may vary (convergence)
...
scipy.sparse.linalg._eigen.arpack.arpack.ArpackNoConvergence:
ARPACK error -1: No convergence (1001 iterations, 0/3 eigenvectors converged)
Ups. Wir sehen, dass ARPACK, wie oben erwähnt, nicht ganz so gut darin ist, kleine Eigenwerte zu finden. Es gibt ein paar Möglichkeiten, dieses Problem zu lösen. Wir könnten die Toleranz (tol) erhöhen, um eine schnellere Konvergenz zu erreichen
>>> evals_small, evecs_small = eigsh(X, 3, which='SM', tol=1E-2)
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 0.99999999 0.00000024 -0.00000049], # may vary (signs)
[-0.00000023 0.99999999 0.00000056],
[ 0.00000031 -0.00000037 0.99999852]])
Das funktioniert, aber wir verlieren die Präzision der Ergebnisse. Eine andere Möglichkeit ist, die maximale Anzahl von Iterationen (maxiter) von 1000 auf 5000 zu erhöhen
>>> evals_small, evecs_small = eigsh(X, 3, which='SM', maxiter=5000)
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 1. 0. 0.], # may vary (signs)
[-0. 1. 0.],
[ 0. 0. -1.]])
Wir erhalten die erhofften Ergebnisse, aber die Rechenzeit ist viel länger. Glücklicherweise enthält ARPACK einen Modus, der eine schnelle Bestimmung von nicht-extremen Eigenwerten ermöglicht: den *Shift-Invert-Modus*. Wie oben erwähnt, beinhaltet dieser Modus die Transformation des Eigenwertproblems in ein äquivalentes Problem mit anderen Eigenwerten. In diesem Fall hoffen wir, Eigenwerte nahe Null zu finden, also wählen wir sigma = 0. Die transformierten Eigenwerte erfüllen dann \(\nu = 1/(\lambda - \sigma) = 1/\lambda\), sodass unsere kleinen Eigenwerte \(\lambda\) zu großen Eigenwerten \(\nu\) werden.
>>> evals_small, evecs_small = eigsh(X, 3, sigma=0, which='LM')
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 1. 0. 0.], # may vary (signs)
[ 0. -1. -0.],
[-0. -0. 1.]])
Wir erhalten die erhofften Ergebnisse mit deutlich geringerem Rechenaufwand. Beachten Sie, dass die Transformation von \(\nu \to \lambda\) vollständig im Hintergrund stattfindet. Der Benutzer muss sich nicht um die Details kümmern.
Der Shift-Invert-Modus bietet mehr als nur eine schnelle Möglichkeit, einige kleine Eigenwerte zu erhalten. Angenommen, Sie möchten innere Eigenwerte und Eigenvektoren finden, z.B. diejenigen, die \(\lambda = 1\) am nächsten liegen. Setzen Sie einfach sigma = 1 und ARPACK erledigt den Rest
>>> evals_mid, evecs_mid = eigsh(X, 3, sigma=1, which='LM')
>>> i_sort = np.argsort(abs(1. / (1 - evals_all)))[-3:]
>>> evals_all[i_sort]
array([0.94164107, 1.05464515, 0.99090277])
>>> evals_mid
array([0.94164107, 0.99090277, 1.05464515])
>>> print(np.dot(evecs_mid.T, evecs_all[:,i_sort]))
array([[-0. 1. 0.], # may vary (signs)
[-0. -0. 1.],
[ 1. 0. 0.]]
Die Eigenwerte kommen in einer anderen Reihenfolge, aber sie sind alle da. Beachten Sie, dass der Shift-Invert-Modus die interne Lösung einer Matrixinverse erfordert. Dies wird automatisch von eigsh und eigs übernommen, aber die Operation kann auch vom Benutzer spezifiziert werden. Weitere Einzelheiten finden Sie in der Docstring von scipy.sparse.linalg.eigsh und scipy.sparse.linalg.eigs.
Verwendung von LinearOperator#
Wir betrachten nun den Fall, in dem Sie vermeiden möchten, eine dichte Matrix zu erstellen, und stattdessen scipy.sparse.linalg.LinearOperator verwenden möchten. Unser erster linearer Operator wendet eine elementweise Multiplikation zwischen dem Eingangsvektor und einem vom Benutzer an den Operator selbst übergebenen Vektor \(\mathbf{d}\) an. Dieser Operator imitiert eine Diagonalmatrix mit den Elementen von \(\mathbf{d}\) entlang der Hauptdiagonale und hat den Hauptvorteil, dass die Vorwärts- und Adjungiertoperationen einfache elementweise Multiplikationen und keine Matrix-Vektor-Multiplikationen sind. Für eine Diagonalmatrix erwarten wir, dass die Eigenwerte gleich den Elementen entlang der Hauptdiagonale sind, in diesem Fall \(\mathbf{d}\). Die mit eigsh erhaltenen Eigenwerte und Eigenvektoren werden mit denen verglichen, die durch die Anwendung von eigh auf die dichte Matrix erhalten werden
>>> from scipy.sparse.linalg import LinearOperator
>>> class Diagonal(LinearOperator):
... def __init__(self, diag, dtype='float32'):
... self.diag = diag
... self.shape = (len(self.diag), len(self.diag))
... self.dtype = np.dtype(dtype)
... def _matvec(self, x):
... return self.diag*x
... def _rmatvec(self, x):
... return self.diag*x
>>> N = 100
>>> rng = np.random.default_rng()
>>> d = rng.normal(0, 1, N).astype(np.float64)
>>> D = np.diag(d)
>>> Dop = Diagonal(d, dtype=np.float64)
>>> evals_all, evecs_all = eigh(D)
>>> evals_large, evecs_large = eigsh(Dop, 3, which='LA', maxiter=1e3)
>>> evals_all[-3:]
array([1.53092498, 1.77243671, 2.00582508])
>>> evals_large
array([1.53092498, 1.77243671, 2.00582508])
>>> print(np.dot(evecs_large.T, evecs_all[:,-3:]))
array([[-1. 0. 0.], # may vary (signs)
[-0. -1. 0.],
[ 0. 0. -1.]]
In diesem Fall haben wir einen schnellen und einfachen Diagonal-Operator erstellt. Die externe Bibliothek PyLops bietet ähnliche Fähigkeiten im Diagonal-Operator sowie mehrere andere Operatoren.
Schließlich betrachten wir einen linearen Operator, der die Anwendung eines Differenzenquotienten erster Ordnung imitiert. In diesem Fall ist der Operator äquivalent zu einer reellen nicht-symmetrischen Matrix. Wiederum vergleichen wir die geschätzten Eigenwerte und Eigenvektoren mit denen einer dichten Matrix, die dieselbe erste Ableitung auf ein Eingangssignal anwendet
>>> class FirstDerivative(LinearOperator):
... def __init__(self, N, dtype='float32'):
... self.N = N
... self.shape = (self.N, self.N)
... self.dtype = np.dtype(dtype)
... def _matvec(self, x):
... y = np.zeros(self.N, self.dtype)
... y[1:-1] = (0.5*x[2:]-0.5*x[0:-2])
... return y
... def _rmatvec(self, x):
... y = np.zeros(self.N, self.dtype)
... y[0:-2] = y[0:-2] - (0.5*x[1:-1])
... y[2:] = y[2:] + (0.5*x[1:-1])
... return y
>>> N = 21
>>> D = np.diag(0.5*np.ones(N-1), k=1) - np.diag(0.5*np.ones(N-1), k=-1)
>>> D[0] = D[-1] = 0 # take away edge effects
>>> Dop = FirstDerivative(N, dtype=np.float64)
>>> evals_all, evecs_all = eig(D)
>>> evals_large, evecs_large = eigs(Dop, 4, which='LI')
>>> evals_all_imag = evals_all.imag
>>> isort_imag = np.argsort(np.abs(evals_all_imag))
>>> evals_all_imag = evals_all_imag[isort_imag]
>>> evals_large_imag = evals_large.imag
>>> isort_imag = np.argsort(np.abs(evals_large_imag))
>>> evals_large_imag = evals_large_imag[isort_imag]
>>> evals_all_imag[-4:]
array([-0.95105652, 0.95105652, -0.98768834, 0.98768834])
>>> evals_large_imag
array([0.95105652, -0.95105652, 0.98768834, -0.98768834]) # may vary
Beachten Sie, dass die Eigenwerte dieses Operators alle imaginär sind. Darüber hinaus erzeugt das Schlüsselwort which='LI' von scipy.sparse.linalg.eigs die Eigenwerte mit der größten absoluten imaginären Komponente (sowohl positiv als auch negativ). Wiederum ist eine fortgeschrittenere Implementierung des Differenzenquotienten erster Ordnung in der PyLops-Bibliothek unter dem Namen FirstDerivative-Operator verfügbar.