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，其中 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’}，可选

使用的求解器。支持‘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-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 的特征值求解器的朴素实现，具体取决于哪个尺寸较小，然后是瑞利-里兹方法作为后处理；参见 Wikipedia 中使用法矩阵，瑞利-里兹方法，（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, 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

内存效率最高的方案是既不显式构造原始矩阵也不显式构造其转置。我们的示例计算“LinearOperator”的最小奇异值和向量，该运算符由 numpy 函数“np.diff”按列方式构造，以便与作用于列的“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 日），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]]

是所谓的基于边的拉普拉斯算子；参见拉普拉斯矩阵中通过关联矩阵的对称拉普拉斯算子，（2022 年 11 月 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')

奇异值的平方和奇异向量是显式已知的；参见 Wikipedia 中二阶导数的特征值和特征向量中的纯狄利克雷边界条件，（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