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)[源代码]#
稀疏矩阵的部分奇异值分解。
计算稀疏矩阵 A 的最大或最小的 k 个奇异值和对应的奇异向量。返回奇异值的顺序不保证。
在下面的描述中，令
M, N = A.shape。- 参数:
- Andarray、稀疏矩阵或 LinearOperator
要分解的浮点数值类型矩阵。
- kint，默认值: 6
要计算的奇异值和奇异向量的数量。必须满足
1 <= k <= kmax，其中kmax=min(M, N)对于solver='propack'，否则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’}, 可选
- 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 是一个 int，则使用一个新的
RandomState实例，并使用 random_state 进行种子设置。如果 random_state 已经是
Generator或RandomState实例，则使用该实例。
在版本 1.15.0 中更改: 作为从使用
numpy.random.RandomState过渡到numpy.random.Generator的 SPEC-007 转换的一部分，此关键字从 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（取决于哪个尺寸较小）的特征求解器的简单实现，然后使用瑞利-里兹方法进行后处理；请参阅维基百科中关于使用正规矩阵，瑞利-里兹方法（2022 年 11 月 19 日）， https://w.wiki/4zms。或者，可以调用 PROPACK 求解器。
输入矩阵 A 的数值类型 (numeric dtype) 选择可能受到限制。只有
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
与这些向量所跨越的二维不变子空间之间的角度相反，右奇异向量的角度很小
>>> 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]]
可以视为关联矩阵；请参阅关联矩阵，（2022 年 11 月 19 日），维基百科，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]]
是所谓的基于边的拉普拉斯算子；请参阅拉普拉斯矩阵中通过关联矩阵的对称拉普拉斯算子，（2022 年 11 月 19 日），维基百科，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')
奇异值的平方和奇异向量是明确已知的；请参阅二阶导数的特征值和特征向量中的纯狄利克雷边界条件，（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 >>> np.allclose(np.abs(u), np.abs(ue), atol=1e-6) True