Python SciPy linalg.svds用法及代碼示例

本文簡要介紹 python 語言中 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', random_state=None, options=None)#

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

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

在下麵的說明中，讓 M, N = A.shape 。

參數 ：：

A： ndarray、稀疏矩陣或LinearOperator

用於分解浮點數字數據類型的矩陣。

k： 整數，默認：6

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

ncv： 整數，可選

當 solver='arpack' 時，這是生成的 Lanczos 向量的數量。有關詳細信息，請參閱‘arpack’。當 solver='lobpcg' 或 solver='propack' 時，忽略此參數。

tol： 浮點數，可選

奇異值的公差。零(默認)表示機器精度。

which： {‘LM’, ‘SM’}

要查找哪 k 個奇異值：最大幅度 (‘LM’) 或最小幅度 (‘SM’) 奇異值。

v0： ndarray，可選

迭代的起始向量；有關詳細信息，請參閱 method-specific 文檔 (‘arpack’, ‘lobpcg’) 或 ‘propack’。

maxiter： 整數，可選

最大迭代次數；有關詳細信息，請參閱 method-specific 文檔 (‘arpack’, ‘lobpcg’) 或 ‘propack’。

return_singular_vectors： {對，錯，“u”, “vh”}

奇異值總是被計算並返回；該參數控製奇異向量的計算和返回。

• True：返回奇異向量。

• False：不返回奇異向量。

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

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

如果 solver='propack' ，無論矩陣形狀如何，都會尊重該選項。

solver： {‘arpack’, ‘propack’, ‘lobpcg’}，可選

使用的求解器。‘arpack’,‘lobpcg’， 和‘propack’都支持。默認：‘arpack’.

random_state： {無，整數， numpy.random.Generator ，

numpy.random.RandomState}, optional

用於生成重采樣的偽隨機數生成器狀態。

如果random_state是None(或者np.random)， 這numpy.random.RandomState使用單例。如果random_state是一個 int，一個新的RandomState使用實例，播種random_state.如果random_state已經是一個Generator或者RandomState實例然後使用該實例。

options： 字典，可選

solver-specific 選項的字典。當前不支持solver-specific選項；此參數保留供將來使用。

返回 ：：

u： ndarray，形狀=(M，k)

具有左奇異向量作為列的酉矩陣。

s： ndarray，形狀=(k，)

奇異值。

vh： ndarray，形狀=(k，N)

具有右奇異向量作為行的酉矩陣。

注意：

這是一種簡單的實現，使用 ARPACK 或 LOBPCG 作為矩陣 A.conj().T @ A 或 A @ A.conj().T 上的特征求解器，具體取決於哪個尺寸較小，然後使用 Rayleigh-Ritz 方法作為後處理；請參閱使用法線矩陣，在 Rayleigh-Ritz 方法中，(2022 年，11 月 19 日)，維基百科，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, random_state=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

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

>>> 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.RandomState(0)
>>> X_dense = rng.random(size=(100, 100))
>>> X_dense[:, 2 * np.arange(50)] = 0
>>> X = sparse.csr_matrix(X_dense)
>>> _, singular_values, _ = svds(X, k=5, random_state=rng)
>>> print(singular_values)
[ 4.3293...  4.4491...  4.5420...  4.5987... 35.2410...]

無需顯式構造輸入矩陣的轉置即可調用該函數。

>>> rng = np.random.default_rng()
>>> G = sparse.rand(8, 9, density=0.5, random_state=rng)
>>> Glo = aslinearoperator(G)
>>> _, singular_values_svds, _ = svds(Glo, k=5, random_state=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’是bi-diagonal，也可以直接創建

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

可以看作是關聯矩陣；請參閱具有 5 個頂點和 4 個邊的線性圖的關聯矩陣(2022 年 11 月 19 日)，維基百科，https://w.wiki/5YXU。因此 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]]

是所謂的基於邊的拉普拉斯算子；請參閱拉普拉斯矩陣中的對稱拉普拉斯通過關聯矩陣，(2022 年，11 月 19 日)，維基百科，https://w.wiki/5YXW 。

“LinearOperator”設置需要選項 ‘rmatvec’ 和 ‘rmatmat’ 乘以矩陣轉置 M.T ，但我們希望成為 matrix-free 以節省內存，因此了解 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

現在我們設置matrix-free“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')

奇異值平方和奇異向量是明確已知的；請參閱純狄利克雷邊界條件，在二階導數的特征值和特征向量中，(2022 年，11 月 19 日)，維基百科，https://w.wiki/5YX6，因為 ‘diff’ 對應於一階導數，其較小的大小 n-1 x n-1 法線矩陣M @ M.T 表示具有狄利克雷邊界條件的離散二階導數。我們使用這些分析表達式進行驗證。

>>> 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
>>> print(np.allclose(np.abs(u), np.abs(ue), atol=1e-6))
True