Interpolative Matrixzerlegung (scipy.linalg.interpolative)#

Hinzugefügt in Version 0.13.

Geändert in Version 1.15.0: Die zugrunde liegenden Algorithmen wurden aus dem ursprünglichen Fortran77-Code nach Python portiert. Weitere Details finden Sie in den unten stehenden Referenzen.

Eine interpolative Zerlegung (ID) einer Matrix \(A \in \mathbb{C}^{m \times n}\) vom Rang \(k \leq \min \{ m, n \}\) ist eine Faktorisierung

\[A \Pi = \begin{bmatrix} A \Pi_{1} & A \Pi_{2} \end{bmatrix} = A \Pi_{1} \begin{bmatrix} I & T \end{bmatrix},\]

wobei \(\Pi = [\Pi_{1}, \Pi_{2}]\) eine Permutationsmatrix mit \(\Pi_{1} \in \{ 0, 1 \}^{n \times k}\) ist, d.h. \(A \Pi_{2} = A \Pi_{1} T\). Dies kann äquivalent geschrieben werden als \(A = BP\), wobei \(B = A \Pi_{1}\) und \(P = [I, T] \Pi^{\mathsf{T}}\) die Skelett- und Interpolationsmatrizen sind.

Wenn \(A\) keinen exakten Rang \(k\) hat, dann existiert eine Annäherung in Form einer ID, so dass \(A = BP + E\), wobei \(\| E \| \sim \sigma_{k + 1}\) von der Ordnung des \((k + 1)\)-ten größten Singulärwerts von \(A\) ist. Beachten Sie, dass \(\sigma_{k + 1}\) der bestmögliche Fehler für eine Rang-\(k\)-Annäherung ist und tatsächlich durch die Singulärwertzerlegung (SVD) \(A \approx U S V^{*}\) erreicht wird, wobei \(U \in \mathbb{C}^{m \times k}\) und \(V \in \mathbb{C}^{n \times k}\) orthonormale Spalten haben und \(S = \mathop{\mathrm{diag}} (\sigma_{i}) \in \mathbb{C}^{k \times k}\) diagonal mit nicht-negativen Einträgen ist. Die Hauptvorteile der Verwendung einer ID gegenüber einer SVD sind:

  • die Konstruktion ist kostengünstiger;

  • die Struktur von \(A\) bleibt erhalten;

  • und die Berechnung ist effizienter aufgrund der Einheitsmatrix von \(P\).

Routinen#

Hauptfunktionalität

interp_decomp(A, eps_or_k[, rand, rng])

Berechnet die ID einer Matrix.

reconstruct_matrix_from_id(B, idx, proj)

Rekonstruiert eine Matrix aus ihrer ID.

reconstruct_interp_matrix(idx, proj)

Rekonstruiert die Interpolationsmatrix aus der ID.

reconstruct_skel_matrix(A, k, idx)

Rekonstruiert die Skelettmatrix aus der ID.

id_to_svd(B, idx, proj)

Konvertiert ID in SVD.

svd(A, eps_or_k[, rand, rng])

Berechnet die SVD einer Matrix über eine ID.

estimate_spectral_norm(A[, its, rng])

Schätzt die Spektralnorm einer Matrix mit der randomisierten Potenzmethode.

estimate_spectral_norm_diff(A, B[, its, rng])

Schätzt die Spektralnorm der Differenz zweier Matrizen mit der randomisierten Potenzmethode.

estimate_rank(A, eps[, rng])

Schätzt den Rang einer Matrix mit einer angegebenen relativen Genauigkeit unter Verwendung von randomisierten Methoden.

Folgende Hilfsfunktionen sind veraltet und werden in SciPy 1.17.0 entfernt

seed([seed])

Diese Funktion diente historisch dazu, den Seed der Randomisierungsalgorithmen zu setzen, die in den Fortran77-geschriebenen Funktionen von scipy.linalg.interpolative verwendet werden.

rand(*shape)

Diese Funktion diente historisch dazu, gleichverteilte Zufallszahlen für die Randomisierungsalgorithmen zu generieren, die in den Fortran77-geschriebenen Funktionen von scipy.linalg.interpolative verwendet werden.

Referenzen#

Dieses Modul verwendet die Algorithmen aus dem ID-Softwarepaket [R5a82238cdab4-1] von Martinsson, Rokhlin, Shkolnisky und Tygert, einer Fortran-Bibliothek zur Berechnung von IDs mit verschiedenen Algorithmen, einschließlich des Rang-aufdeckenden QR-Ansatzes von [R5a82238cdab4-2] und der neueren randomisierten Methoden, die in [R5a82238cdab4-3], [R5a82238cdab4-4] und [R5a82238cdab4-5] beschrieben werden.

Wir empfehlen dem Benutzer, auch die Dokumentation des ID-Pakets zu konsultieren.

[R5a82238cdab4-1]

P.G. Martinsson, V. Rokhlin, Y. Shkolnisky, M. Tygert. “ID: a software package for low-rank approximation of matrices via interpolative decompositions, version 0.2.” http://tygert.com/id_doc.4.pdf.

[R5a82238cdab4-2]

H. Cheng, Z. Gimbutas, P.G. Martinsson, V. Rokhlin. “On the compression of low rank matrices.” SIAM J. Sci. Comput. 26 (4): 1389–1404, 2005. DOI:10.1137/030602678.

[R5a82238cdab4-3]

E. Liberty, F. Woolfe, P.G. Martinsson, V. Rokhlin, M. Tygert. “Randomized algorithms for the low-rank approximation of matrices.” Proc. Natl. Acad. Sci. U.S.A. 104 (51): 20167–20172, 2007. DOI:10.1073/pnas.0709640104.

[R5a82238cdab4-4]

P.G. Martinsson, V. Rokhlin, M. Tygert. “A randomized algorithm for the decomposition of matrices.” Appl. Comput. Harmon. Anal. 30 (1): 47–68, 2011. DOI:10.1016/j.acha.2010.02.003.

[R5a82238cdab4-5]

F. Woolfe, E. Liberty, V. Rokhlin, M. Tygert. “A fast randomized algorithm for the approximation of matrices.” Appl. Comput. Harmon. Anal. 25 (3): 335–366, 2008. DOI:10.1016/j.acha.2007.12.002.

Tutorial#

Initialisierung#

Der erste Schritt ist der Import von scipy.linalg.interpolative durch Ausführen des Befehls

>>> import scipy.linalg.interpolative as sli

Nun erstellen wir eine Matrix. Dazu betrachten wir eine Hilbert-Matrix, die bekanntermaßen Rang-niedrig ist.

>>> from scipy.linalg import hilbert
>>> n = 1000
>>> A = hilbert(n)

Wir können dies auch explizit tun über

>>> import numpy as np
>>> n = 1000
>>> A = np.empty((n, n), order='F')
>>> for j in range(n):
...     for i in range(n):
...         A[i,j] = 1. / (i + j + 1)

Beachten Sie die Verwendung des Flags order='F' in numpy.empty. Dies instanziiert die Matrix in Fortran-kontiguer Reihenfolge und ist wichtig, um Datenkopien beim Übergeben an das Backend zu vermeiden.

Dann definieren wir Multiplikationsroutinen für die Matrix, indem wir sie als scipy.sparse.linalg.LinearOperator betrachten.

>>> from scipy.sparse.linalg import aslinearoperator
>>> L = aslinearoperator(A)

Dadurch werden automatisch Methoden eingerichtet, die die Wirkung der Matrix und ihrer Adjungierten auf einen Vektor beschreiben.

Berechnung einer ID#

Wir haben mehrere Algorithmusoptionen zur Berechnung einer ID. Diese lassen sich weitgehend nach zwei Dichotomien einteilen:

  1. wie die Matrix dargestellt wird, d.h. über ihre Einträge oder über ihre Wirkung auf einen Vektor;

  2. und ob sie auf eine feste relative Genauigkeit oder einen festen Rang angenähert werden soll.

Wir gehen unten schrittweise durch jede Option.

In allen Fällen wird die ID durch drei Parameter dargestellt:

  1. einen Rang k;

  2. ein Index-Array idx; und

  3. Interpolationskoeffizienten proj.

Die ID ist durch die Beziehung np.dot(A[:,idx[:k]], proj) == A[:,idx[k:]] spezifiziert.

Aus Matrixeinträgen#

Wir betrachten zunächst eine Matrix, die in Bezug auf ihre Einträge gegeben ist.

Um eine ID mit fester Genauigkeit zu berechnen, geben Sie ein:

>>> eps = 1e-3
>>> k, idx, proj = sli.interp_decomp(A, eps)

wobei eps < 1 die gewünschte Genauigkeit ist.

Um eine ID mit festem Rang zu berechnen, verwenden Sie:

>>> idx, proj = sli.interp_decomp(A, k)

wobei k >= 1 der gewünschte Rang ist.

Beide Algorithmen verwenden Zufallsstichproben und sind in der Regel schneller als die entsprechenden älteren, deterministischen Algorithmen, die über die Befehle aufgerufen werden können:

>>> k, idx, proj = sli.interp_decomp(A, eps, rand=False)

und

>>> idx, proj = sli.interp_decomp(A, k, rand=False)

beziehungsweise.

Aus Matrixaktion#

Betrachten wir nun eine Matrix, die in Bezug auf ihre Wirkung auf einen Vektor als scipy.sparse.linalg.LinearOperator gegeben ist.

Um eine ID mit fester Genauigkeit zu berechnen, geben Sie ein:

>>> k, idx, proj = sli.interp_decomp(L, eps)

Um eine ID mit festem Rang zu berechnen, verwenden Sie:

>>> idx, proj = sli.interp_decomp(L, k)

Diese Algorithmen sind randomisiert.

Rekonstruktion einer ID#

Die obigen ID-Routinen geben nicht explizit die Skelett- und Interpolationsmatrizen aus, sondern geben die relevanten Informationen in einer kompakteren (und manchmal nützlicheren) Form zurück. Um diese Matrizen zu erstellen, schreiben Sie:

>>> B = sli.reconstruct_skel_matrix(A, k, idx)

für die Skelettmatrix und

>>> P = sli.reconstruct_interp_matrix(idx, proj)

für die Interpolationsmatrix. Die ID-Annäherung kann dann berechnet werden als:

>>> C = np.dot(B, P)

Dies kann auch direkt über folgende Funktion konstruiert werden:

>>> C = sli.reconstruct_matrix_from_id(B, idx, proj)

ohne dass P zuerst berechnet werden muss.

Alternativ kann dies auch explizit mit folgender Funktion erfolgen:

>>> B = A[:,idx[:k]]
>>> P = np.hstack([np.eye(k), proj])[:,np.argsort(idx)]
>>> C = np.dot(B, P)

Berechnung einer SVD#

Eine ID kann über den Befehl in eine SVD umgewandelt werden:

>>> U, S, V = sli.id_to_svd(B, idx, proj)

Die SVD-Annäherung ist dann:

>>> approx = U @ np.diag(S) @ V.conj().T

Die SVD kann auch „frisch“ berechnet werden, indem sowohl die ID- als auch die Konvertierungsschritte zu einem Befehl kombiniert werden. Nach den verschiedenen oben genannten ID-Algorithmen gibt es entsprechend verschiedene SVD-Algorithmen, die verwendet werden können.

Aus Matrixeinträgen#

Wir betrachten zuerst SVD-Algorithmen für eine Matrix, die in Bezug auf ihre Einträge gegeben ist.

Um eine SVD mit fester Genauigkeit zu berechnen, geben Sie ein:

>>> U, S, V = sli.svd(A, eps)

Um eine SVD mit festem Rang zu berechnen, verwenden Sie:

>>> U, S, V = sli.svd(A, k)

Beide Algorithmen verwenden Zufallsstichproben; für die deterministischen Versionen geben Sie wie oben das Schlüsselwort rand=False ein.

Aus Matrixaktion#

Betrachten wir nun eine Matrix, die in Bezug auf ihre Wirkung auf einen Vektor gegeben ist.

Um eine SVD mit fester Genauigkeit zu berechnen, geben Sie ein:

>>> U, S, V = sli.svd(L, eps)

Um eine SVD mit festem Rang zu berechnen, verwenden Sie:

>>> U, S, V = sli.svd(L, k)

Hilfsroutinen#

Es sind auch mehrere Hilfsroutinen verfügbar.

Um die Spektralnorm einer Matrix zu schätzen, verwenden Sie:

>>> snorm = sli.estimate_spectral_norm(A)

Dieser Algorithmus basiert auf der randomisierten Potenzmethode und erfordert daher nur Matrix-Vektor-Produkte. Die Anzahl der Iterationen kann mit dem Schlüsselwort its (Standard: its=20) eingestellt werden. Die Matrix wird als scipy.sparse.linalg.LinearOperator interpretiert, aber es ist auch gültig, sie als numpy.ndarray zu übergeben, in diesem Fall wird sie trivial mit scipy.sparse.linalg.aslinearoperator konvertiert.

Dieselbe Algorithmus kann auch die Spektralnorm der Differenz zweier Matrizen A1 und A2 wie folgt schätzen:

>>> A1, A2 = A**2, A
>>> diff = sli.estimate_spectral_norm_diff(A1, A2)

Dies ist oft nützlich, um die Genauigkeit einer Matrixannäherung zu überprüfen.

Einige Routinen in scipy.linalg.interpolative erfordern auch die Schätzung des Rangs einer Matrix. Dies kann entweder erfolgen mit:

>>> k = sli.estimate_rank(A, eps)

oder

>>> k = sli.estimate_rank(L, eps)

je nach Darstellung. Der Parameter eps steuert die Definition des numerischen Rangs.

Schließlich kann die Zufallszahlengenerierung, die für alle randomisierten Routinen erforderlich ist, gesteuert werden, indem NumPy-Pseudo-Zufallsgeneratoren mit einem festen Seed bereitgestellt werden. Siehe numpy.random.Generator und numpy.random.default_rng für weitere Details.

Anmerkungen#

Die oben genannten Funktionen erkennen automatisch die geeignete Schnittstelle und arbeiten sowohl mit reellen als auch mit komplexen Datentypen, wobei die Eingabeargumente an die entsprechende Backend-Routine weitergegeben werden.