scipy.sparse.linalg.

svds#

scipy.sparse.linalg.svds(A, k=6, ncv=None, tol=0, which='LM', v0=None, maxiter=None, return_singular_vectors=True, solver='arpack', rng=None, options=None)[Quelle]#

Partielle Singulärwertzerlegung einer dünnbesetzten Matrix.

Berechnet die größten oder kleinsten k Singulärwerte und die zugehörigen Singulärvektoren einer dünnbesetzten Matrix A. Die Reihenfolge, in der die Singulärwerte zurückgegeben werden, ist nicht garantiert.

In den folgenden Beschreibungen sei M, N = A.shape.

Parameter:

Andarray, dünnbesetzte Matrix oder LinearOperator

Zu zerlegende Matrix mit numerischem Gleitkomma-Datentyp.

kint, Standardwert: 6

Anzahl der zu berechnenden Singulärwerte und Singulärvektoren. Muss 1 <= k <= kmax erfüllen, wobei kmax=min(M, N) für solver='propack' und kmax=min(M, N) - 1 andernfalls gilt.

ncvint, optional

Wenn solver='arpack', ist dies die Anzahl der erzeugten Lanczos-Vektoren. Einzelheiten finden Sie unter ‘arpack’. Wenn solver='lobpcg' oder solver='propack', wird dieser Parameter ignoriert.

tolfloat, optional

Toleranz für Singulärwerte. Null (Standardwert) bedeutet Maschinengenauigkeit.

which{‘LM’, ‘SM’}

Welche k Singulärwerte gefunden werden sollen: entweder die mit der größten Magnitude (‘LM’) oder die mit der kleinsten Magnitude (‘SM’) Singulärwerte.

v0ndarray, optional

Der Startvektor für die Iteration; siehe methodenspezifische Dokumentation (‘arpack’, ‘lobpcg’) oder ‘propack’ für Details.

maxiterint, optional

Maximale Anzahl von Iterationen; siehe methodenspezifische Dokumentation (‘arpack’, ‘lobpcg’) oder ‘propack’ für Details.

return_singular_vectors{True, False, “u”, “vh”}

Singulärwerte werden immer berechnet und zurückgegeben; dieser Parameter steuert die Berechnung und Rückgabe von Singulärvektoren.

True: Singulärvektoren zurückgeben.
False: keine Singulärvektoren zurückgeben.
"u": Wenn M <= N, nur die linken Singulärvektoren berechnen und None für die rechten Singulärvektoren zurückgeben. Andernfalls alle Singulärvektoren berechnen.
"vh": Wenn M > N, nur die rechten Singulärvektoren berechnen und None für die linken Singulärvektoren zurückgeben. Andernfalls alle Singulärvektoren berechnen.

Wenn solver='propack', wird die Option unabhängig von der Matrixform berücksichtigt.

solver{‘arpack’, ‘propack’, ‘lobpcg’}, optional

Der verwendete Solver. ‘arpack’, ‘lobpcg’ und ‘propack’ werden unterstützt. Standard: ‘arpack’.

rng{None, int, numpy.random.Generator}, optional

Wenn rng als Schlüsselwort übergeben wird, werden andere Typen als numpy.random.Generator an numpy.random.default_rng übergeben, um einen Generator zu instanziieren. Wenn rng bereits eine Generator-Instanz ist, dann wird die bereitgestellte Instanz verwendet. Geben Sie rng für reproduzierbares Funktionsverhalten an.

Wenn dieses Argument positional übergeben wird oder random_state als Schlüsselwort übergeben wird, gilt das ältere Verhalten für das Argument random_state.

Wenn random_state None ist (oder numpy.random), wird die Singleton-Instanz numpy.random.RandomState verwendet.
Wenn random_state eine Ganzzahl ist, wird eine neue RandomState-Instanz verwendet, die mit random_state initialisiert wurde.
Wenn random_state bereits eine Generator- oder RandomState-Instanz ist, wird diese Instanz verwendet.

Geändert in Version 1.15.0: Im Rahmen des SPEC-007-Übergangs von der Verwendung von numpy.random.RandomState zu numpy.random.Generator wurde dieser Schlüssel von random_state zu rng geändert. Für eine Übergangszeit werden beide Schlüssel weiterhin funktionieren, obwohl nur einer gleichzeitig angegeben werden darf. Nach der Übergangszeit werden Funktionsaufrufe mit dem Schlüssel random_state Warnungen ausgeben. Das Verhalten von random_state und rng wird oben erläutert, aber nur der Schlüssel rng sollte in neuem Code verwendet werden.

optionsdict, optional

Ein Wörterbuch mit Solver-spezifischen Optionen. Derzeit werden keine Solver-spezifischen Optionen unterstützt; dieser Parameter ist für die zukünftige Verwendung reserviert.

Rückgabe:

undarray, Form=(M, k): Unitäre Matrix mit linken Singulärvektoren als Spalten.
sndarray, Form=(k,): Die Singulärwerte.
vhndarray, Form=(k, N): Unitäre Matrix mit rechten Singulärvektoren als Zeilen.

Hinweise

Dies ist eine naive Implementierung, die ARPACK oder LOBPCG als Eigenlöser auf der Matrix A.conj().T @ A oder A @ A.conj().T verwendet, je nachdem, welche die kleinere Größe hat, gefolgt von der Rayleigh-Ritz-Methode als Nachbearbeitung; siehe Using the normal matrix, in Rayleigh-Ritz method, (2022, Nov. 19), Wikipedia, https://w.wiki/4zms.

Alternativ kann der PROPACK-Solver aufgerufen werden.

Die Auswahl der numerischen Datentypen der Eingabematrix A kann begrenzt sein. Nur solver="lobpcg" unterstützt alle Gleitkomma-Datentypen real: ‘np.float32’, ‘np.float64’, ‘np.longdouble’ und komplex: ‘np.complex64’, ‘np.complex128’, ‘np.clongdouble’. solver="arpack" unterstützt nur ‘np.float32’, ‘np.float64’ und ‘np.complex128’.

Beispiele

Konstruieren Sie eine Matrix A aus Singulärwerten und -vektoren.

>>> import numpy as np
>>> from scipy import sparse, linalg, stats
>>> from scipy.sparse.linalg import svds, aslinearoperator, LinearOperator

Konstruieren Sie eine dichte Matrix A aus Singulärwerten und -vektoren.

>>> rng = np.random.default_rng()
>>> orthogonal = stats.ortho_group.rvs(10, random_state=rng)
>>> s = [1e-3, 1, 2, 3, 4]  # non-zero singular values
>>> u = orthogonal[:, :5]         # left singular vectors
>>> vT = orthogonal[:, 5:].T      # right singular vectors
>>> A = u @ np.diag(s) @ vT

Mit nur vier Singulärwerten/-vektoren approximiert die SVD die ursprüngliche Matrix.

>>> u4, s4, vT4 = svds(A, k=4)
>>> A4 = u4 @ np.diag(s4) @ vT4
>>> np.allclose(A4, A, atol=1e-3)
True

Mit allen fünf nicht-null Singulärwerten/-vektoren können wir die ursprüngliche Matrix genauer wiederherstellen.

>>> u5, s5, vT5 = svds(A, k=5)
>>> A5 = u5 @ np.diag(s5) @ vT5
>>> np.allclose(A5, A)
True

Die Singulärwerte entsprechen den erwarteten Singulärwerten.

>>> np.allclose(s5, s)
True

Da die Singulärwerte in diesem Beispiel nicht nahe beieinander liegen, entspricht jeder Singulärvektor bis auf ein Vorzeichen dem erwarteten.

>>> (np.allclose(np.abs(u5), np.abs(u)) and
...  np.allclose(np.abs(vT5), np.abs(vT)))
True

Die Singulärvektoren sind ebenfalls orthogonal.

>>> (np.allclose(u5.T @ u5, np.eye(5)) and
...  np.allclose(vT5 @ vT5.T, np.eye(5)))
True

Wenn es (nahezu) mehrere Singulärwerte gibt, können die entsprechenden einzelnen Singulärvektoren instabil sein, aber der gesamte invariante Unterraum, der alle solchen Singulärvektoren enthält, wird korrekt berechnet, was sich durch Winkel zwischen Unterräumen mittels 'subspace_angles' messen lässt.

>>> rng = np.random.default_rng()
>>> s = [1, 1 + 1e-6]  # non-zero singular values
>>> u, _ = np.linalg.qr(rng.standard_normal((99, 2)))
>>> v, _ = np.linalg.qr(rng.standard_normal((99, 2)))
>>> vT = v.T
>>> A = u @ np.diag(s) @ vT
>>> A = A.astype(np.float32)
>>> u2, s2, vT2 = svds(A, k=2, rng=rng)
>>> np.allclose(s2, s)
True

Die Winkel zwischen den einzelnen exakten und berechneten Singulärvektoren sind möglicherweise nicht so klein. Zur Überprüfung verwenden Sie

>>> (linalg.subspace_angles(u2[:, :1], u[:, :1]) +
...  linalg.subspace_angles(u2[:, 1:], u[:, 1:]))
array([0.06562513])  # may vary
>>> (linalg.subspace_angles(vT2[:1, :].T, vT[:1, :].T) +
...  linalg.subspace_angles(vT2[1:, :].T, vT[1:, :].T))
array([0.06562507])  # may vary

Im Gegensatz zu den Winkeln zwischen den 2-dimensionalen invarianten Unterräumen, die diese Vektoren aufspannen, welche für rechte Singulärvektoren klein sind

>>> linalg.subspace_angles(u2, u).sum() < 1e-6
True

sowie für linke Singulärvektoren.

>>> linalg.subspace_angles(vT2.T, vT.T).sum() < 1e-6
True

Das nächste Beispiel folgt dem von 'sklearn.decomposition.TruncatedSVD'.

>>> rng = np.random.default_rng()
>>> X_dense = rng.random(size=(100, 100))
>>> X_dense[:, 2 * np.arange(50)] = 0
>>> X = sparse.csr_array(X_dense)
>>> _, singular_values, _ = svds(X, k=5, rng=rng)
>>> print(singular_values)
[ 4.3221...  4.4043...  4.4907...  4.5858... 35.4549...]

Die Funktion kann aufgerufen werden, ohne dass die Transponierte der Eingabematrix jemals explizit konstruiert wird.

>>> rng = np.random.default_rng()
>>> G = sparse.random_array((8, 9), density=0.5, rng=rng)
>>> Glo = aslinearoperator(G)
>>> _, singular_values_svds, _ = svds(Glo, k=5, rng=rng)
>>> _, singular_values_svd, _ = linalg.svd(G.toarray())
>>> np.allclose(singular_values_svds, singular_values_svd[-4::-1])
True

Das speichereffizienteste Szenario ist, wenn weder die ursprüngliche Matrix noch ihre Transponierte explizit konstruiert werden. Unser Beispiel berechnet die kleinsten Singulärwerte und Vektoren eines 'LinearOperator', der aus der NumPy-Funktion 'np.diff' spaltenweise konstruiert wird, um mit 'LinearOperator' zu übereinstimmen, das auf Spalten operiert.

>>> diff0 = lambda a: np.diff(a, axis=0)

Lassen Sie uns die Matrix aus 'diff0' zur Validierung erstellen.

>>> n = 5  # The dimension of the space.
>>> M_from_diff0 = diff0(np.eye(n))
>>> print(M_from_diff0.astype(int))
[[-1  1  0  0  0]
 [ 0 -1  1  0  0]
 [ 0  0 -1  1  0]
 [ 0  0  0 -1  1]]

Die Matrix 'M_from_diff0' ist bi-diagonal und könnte alternativ direkt durch

>>> M = - np.eye(n - 1, n, dtype=int)
>>> np.fill_diagonal(M[:,1:], 1)
>>> np.allclose(M, M_from_diff0)
True

Ihre Transponierte

>>> print(M.T)
[[-1  0  0  0]
 [ 1 -1  0  0]
 [ 0  1 -1  0]
 [ 0  0  1 -1]
 [ 0  0  0  1]]

kann als Inzidenzmatrix betrachtet werden; siehe Incidence matrix, (2022, Nov. 19), Wikipedia, https://w.wiki/5YXU, eines linearen Graphen mit 5 Knoten und 4 Kanten. Die 5x5-Normalmatrix M.T @ M ist somit

>>> print(M.T @ M)
[[ 1 -1  0  0  0]
 [-1  2 -1  0  0]
 [ 0 -1  2 -1  0]
 [ 0  0 -1  2 -1]
 [ 0  0  0 -1  1]]

der Graph-Laplace-Operator, während die tatsächlich in 'svds' verwendete kleinere 4x4-Normalmatrix M @ M.T

>>> print(M @ M.T)
[[ 2 -1  0  0]
 [-1  2 -1  0]
 [ 0 -1  2 -1]
 [ 0  0 -1  2]]

der sogenannte kantenbasierte Laplace-Operator ist; siehe Symmetric Laplacian via the incidence matrix, in Laplacian matrix, (2022, Nov. 19), Wikipedia, https://w.wiki/5YXW.

Die 'LinearOperator'-Einrichtung benötigt die Optionen 'rmatvec' und 'rmatmat' für die Multiplikation mit der Matrixtransponierten M.T, aber wir wollen speichereffizient sein, daher kennen wir das Aussehen von M.T und konstruieren manuell die folgende Funktion, die in rmatmat=diff0t verwendet wird.

>>> def diff0t(a):
...     if a.ndim == 1:
...         a = a[:,np.newaxis]  # Turn 1D into 2D array
...     d = np.zeros((a.shape[0] + 1, a.shape[1]), dtype=a.dtype)
...     d[0, :] = - a[0, :]
...     d[1:-1, :] = a[0:-1, :] - a[1:, :]
...     d[-1, :] = a[-1, :]
...     return d

Wir prüfen, ob unsere Funktion 'diff0t' für die Matrixtransponierte gültig ist.

>>> np.allclose(M.T, diff0t(np.eye(n-1)))
True

Nun richten wir unseren speichereffizienten 'LinearOperator' namens 'diff0_func_aslo' und zur Validierung den matrixbasierten 'diff0_matrix_aslo' ein.

>>> def diff0_func_aslo_def(n):
...     return LinearOperator(matvec=diff0,
...                           matmat=diff0,
...                           rmatvec=diff0t,
...                           rmatmat=diff0t,
...                           shape=(n - 1, n))
>>> diff0_func_aslo = diff0_func_aslo_def(n)
>>> diff0_matrix_aslo = aslinearoperator(M_from_diff0)

Und validieren sowohl die Matrix als auch ihre Transponierte im 'LinearOperator'.

>>> np.allclose(diff0_func_aslo(np.eye(n)),
...             diff0_matrix_aslo(np.eye(n)))
True
>>> np.allclose(diff0_func_aslo.T(np.eye(n-1)),
...             diff0_matrix_aslo.T(np.eye(n-1)))
True

Nachdem die 'LinearOperator'-Einrichtung validiert wurde, führen wir den Solver aus.

>>> n = 100
>>> diff0_func_aslo = diff0_func_aslo_def(n)
>>> u, s, vT = svds(diff0_func_aslo, k=3, which='SM')

Die Singulärwerte im Quadrat und die Singulärvektoren sind explizit bekannt; siehe Pure Dirichlet boundary conditions, in Eigenvalues and eigenvectors of the second derivative, (2022, Nov. 19), Wikipedia, https://w.wiki/5YX6, da 'diff' der ersten Ableitung entspricht und seine kleinere n-1 x n-1 Normalmatrix M @ M.T die diskrete zweite Ableitung mit Dirichlet-Randbedingungen darstellt. Wir verwenden diese analytischen Ausdrücke zur Validierung.

>>> se = 2. * np.sin(np.pi * np.arange(1, 4) / (2. * n))
>>> ue = np.sqrt(2 / n) * np.sin(np.pi * np.outer(np.arange(1, n),
...                              np.arange(1, 4)) / n)
>>> np.allclose(s, se, atol=1e-3)
True
>>> np.allclose(np.abs(u), np.abs(ue), atol=1e-6)
True