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)[source]#

稀疏矩陣的部分奇異值分解。

計算稀疏矩陣 A 的最大或最小 k 個奇異值和對應的奇異向量。不保證返回奇異值的順序。

在以下描述中，令 M, N = A.shape。

參數:

Andarray、稀疏矩陣或 LinearOperator

要分解的矩陣，具有浮點數值資料類型。

kint，預設值：6

要計算的奇異值和奇異向量的數量。必須滿足 1 <= k <= kmax，其中當 solver='propack' 時 kmax=min(M, N)，否則 kmax=min(M, N) - 1。

ncvint，選填

當 solver='arpack' 時，這是產生的 Lanczos 向量的數量。詳情請參閱 ‘arpack’。當 solver='lobpcg' 或 solver='propack' 時，此參數會被忽略。

tolfloat，選填

奇異值的容忍度。零（預設值）表示機器精度。

which{‘LM’, ‘SM’}

尋找哪 k 個奇異值：最大量值 ('LM') 或最小量值 ('SM') 奇異值。

v0ndarray，選填

迭代的起始向量；詳情請參閱特定方法的文檔 (‘arpack’、‘lobpcg’) 或 ‘propack’。

maxiterint，選填

最大迭代次數；詳情請參閱特定方法的文檔 (‘arpack’、‘lobpcg’) 或 ‘propack’。

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

奇異值始終會被計算和返回；此參數控制奇異向量的計算和返回。

True：返回奇異向量。
False：不返回奇異向量。
"u"：如果 M <= N，則僅計算左奇異向量，並為右奇異向量返回 None。否則，計算所有奇異向量。
"vh"：如果 M > N，則僅計算右奇異向量，並為左奇異向量返回 None。否則，計算所有奇異向量。

如果 solver='propack'，則無論矩陣形狀如何，此選項都有效。

solver{‘arpack’, ‘propack’, ‘lobpcg’}，選填

使用的求解器。支援 ‘arpack’、‘lobpcg’ 和 ‘propack’。預設值：‘arpack’。

rng{None, int, numpy.random.Generator}，選填

如果 rng 通過關鍵字傳遞，則 numpy.random.Generator 以外的類型會傳遞給 numpy.random.default_rng 以實例化 Generator。如果 rng 已經是 Generator 實例，則使用提供的實例。指定 rng 以獲得可重複的函數行為。

如果此參數通過位置傳遞，或 random_state 通過關鍵字傳遞，則參數 random_state 的舊版行為適用

如果 random_state 為 None（或 numpy.random），則使用 numpy.random.RandomState 單例。
如果 random_state 是一個整數，則會使用一個新的 RandomState 實例，並以 random_state 作為種子。
如果 random_state 已經是 Generator 或 RandomState 實例，則使用該實例。

在 1.15.0 版本中變更：作為從使用 numpy.random.RandomState 過渡到 numpy.random.Generator 的 SPEC-0007 的一部分，此關鍵字已從 random_state 更改為 rng。在過渡期間，這兩個關鍵字將繼續有效，但一次只能指定一個。在過渡期之後，使用 random_state 關鍵字的函數調用將發出警告。上面概述了 random_state 和 rng 的行為，但在新程式碼中應僅使用 rng 關鍵字。

optionsdict，選填

求解器特定選項的字典。目前不支援任何求解器特定選項；此參數保留供未來使用。

返回值:

undarray，形狀=(M, k): 具有左奇異向量作為列的么正矩陣。
sndarray，形狀=(k,): 奇異值。
vhndarray，形狀=(k, N): 具有右奇異向量作為行的么正矩陣。

註解

這是一個樸素的實現，它使用 ARPACK 或 LOBPCG 作為矩陣 A.conj().T @ A 或 A @ A.conj().T 的特徵值求解器（取決於哪個尺寸較小），然後使用 Rayleigh-Ritz 方法作為後處理；請參閱 Using the normal matrix，在 Rayleigh-Ritz method，(2022, Nov. 19)，Wikipedia，https://w.wiki/4zms。

或者，可以調用 PROPACK 求解器。

輸入矩陣 A 數值資料類型的選擇可能受到限制。只有 solver="lobpcg" 支援所有浮點資料類型：實數：‘np.float32’、‘np.float64’、‘np.longdouble’ 和複數：‘np.complex64’、‘np.complex128’、‘np.clongdouble’。solver="arpack" 僅支援 ‘np.float32’、‘np.float64’ 和 ‘np.complex128’。

範例

從奇異值和奇異向量建構矩陣 A。

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

從奇異值和奇異向量建構密集矩陣 A。

>>> 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

僅使用四個奇異值/向量，SVD 即可近似原始矩陣。

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

使用所有五個非零奇異值/向量，我們可以更準確地重現原始矩陣。

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

奇異值與預期的奇異值相符。

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

由於在此範例中奇異值彼此不接近，因此每個奇異向量都如預期般匹配，僅在符號上存在差異。

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

奇異向量也是正交的。

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

如果有（幾乎）多個奇異值，則對應的個別奇異向量可能不穩定，但是包含所有此類奇異向量的整個不變子空間會被準確計算，這可以通過子空間之間的角度（通過 'subspace_angles'）來衡量。

>>> 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

個別精確和計算的奇異向量之間的角度可能不是那麼小。要檢查，請使用

>>> (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

與這些向量跨越的 2 維不變子空間之間的角度相反，對於右奇異向量而言，這些角度很小

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

以及對於左奇異向量。

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

下一個範例遵循 ‘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...]

可以調用該函數，而無需顯式建構輸入矩陣的轉置。

>>> 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

最節省記憶體的情況是原始矩陣及其轉置都沒有顯式建構。我們的範例計算從 numpy 函數 ‘np.diff’ 建構的 ‘LinearOperator’ 的最小奇異值和向量，該函數按列使用，以與 ‘LinearOperator’ 在列上的操作保持一致。

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

讓我們從 ‘diff0’ 建立矩陣，僅用於驗證。

>>> 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]]

矩陣 ‘M_from_diff0’ 是雙對角矩陣，也可以通過以下方式直接建立：

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

它的轉置

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

可以看作是關聯矩陣；請參閱 Incidence matrix，(2022, Nov. 19)，Wikipedia，https://w.wiki/5YXU，具有 5 個頂點和 4 條邊的線性圖。因此，5x5 正規矩陣 M.T @ M 是

>>> 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]]

圖拉普拉斯算子，而實際在 ‘svds’ 中使用的較小尺寸 4x4 正規矩陣 M @ M.T 是

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

所謂的邊緣拉普拉斯算子；請參閱 Symmetric Laplacian via the incidence matrix，在 Laplacian matrix，(2022, Nov. 19)，Wikipedia，https://w.wiki/5YXW。

‘LinearOperator’ 設置需要通過矩陣轉置 M.T 進行乘法的選項 ‘rmatvec’ 和 ‘rmatmat’，但我們希望無矩陣化以節省記憶體，因此在知道 M.T 的外觀後，我們手動建構以下函數以在 rmatmat=diff0t 中使用。

>>> 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

我們檢查我們的矩陣轉置函數 ‘diff0t’ 是否有效。

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

現在我們設置我們的無矩陣 ‘LinearOperator’，稱為 ‘diff0_func_aslo’，並為了驗證，設置基於矩陣的 ‘diff0_matrix_aslo’。

>>> 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)

並在 ‘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

驗證了 ‘LinearOperator’ 設置後，我們運行求解器。

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

奇異值的平方和奇異向量是顯式已知的；請參閱 Pure Dirichlet boundary conditions，在 Eigenvalues and eigenvectors of the second derivative，(2022, Nov. 19)，Wikipedia，https://w.wiki/5YX6，因為 ‘diff’ 對應於一階導數，並且其較小尺寸 n-1 x n-1 正規矩陣 M @ M.T 表示具有 Dirichlet 邊界條件的離散二階導數。我們使用這些解析表達式進行驗證。

>>> 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