本文简要介绍 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’},可选
- 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
相关用法
- Python SciPy linalg.svd用法及代码示例
- Python SciPy linalg.svdvals用法及代码示例
- Python SciPy linalg.solve_circulant用法及代码示例
- Python SciPy linalg.spsolve用法及代码示例
- Python SciPy linalg.spsolve_triangular用法及代码示例
- Python SciPy linalg.solve_banded用法及代码示例
- Python SciPy linalg.solve_discrete_lyapunov用法及代码示例
- Python SciPy linalg.splu用法及代码示例
- Python SciPy linalg.spilu用法及代码示例
- Python SciPy linalg.solve用法及代码示例
- Python SciPy linalg.solveh_banded用法及代码示例
- Python SciPy linalg.solve_sylvester用法及代码示例
- Python SciPy linalg.solve_toeplitz用法及代码示例
- Python SciPy linalg.sqrtm用法及代码示例
- Python SciPy linalg.sinm用法及代码示例
- Python SciPy linalg.schur用法及代码示例
- Python SciPy linalg.solve_continuous_lyapunov用法及代码示例
- Python SciPy linalg.solve_continuous_are用法及代码示例
- Python SciPy linalg.solve_discrete_are用法及代码示例
- Python SciPy linalg.solve_triangular用法及代码示例
- Python SciPy linalg.sinhm用法及代码示例
- Python SciPy linalg.signm用法及代码示例
- Python SciPy linalg.subspace_angles用法及代码示例
- Python SciPy linalg.eigvalsh_tridiagonal用法及代码示例
- Python SciPy linalg.cdf2rdf用法及代码示例
注:本文由纯净天空筛选整理自scipy.org大神的英文原创作品 scipy.sparse.linalg.svds。非经特殊声明,原始代码版权归原作者所有,本译文未经允许或授权,请勿转载或复制。