【Python】scipy稀疏矩阵的奇异值分解svds

基本原理

当$A$是方阵时，可以很容易地进行特征分解：$A=W\Sigma W^{-1}$，其中$\Sigma$是$A$的特征值组成的对角矩阵。如果$W$由标准正交基组成，则$W^{-1}=W^T$，特征分解可进一步写成$W^T\Sigma W$。

然而，当$A$不是方阵时，情况大不一样了，但仍然可以将$A$表示成$A=U\Sigma V^T$的形式，其中$\Sigma$也是对角矩阵，对角线上的每个元素被称作奇异值。

奇异值的求解过程和特征值息息相关，因为把$A$变成方阵很简单，只要乘以转置就行。故令$L=A{A}^T$，$R=A^{T}A$，则$L,R$都可以求特征值${\lambda }_i$和特征向量，其中$L$的特征向量为$A$的左奇异向量，$R$的特征向量为右奇异向量。对应的奇异值${\sigma }_i=\sqrt{{\lambda }_i}$。

scipy实现

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 奇异值个数，必须在$[k,k_{\max }]$之间， 当solver='propack'时，$k_{max}=\min (M,N)$，否则$k_{max}=\min (M,N)-1$
• ncv solver='arpack'时，此为Lanczos向量个数，否则此项忽略。
• tol 奇异值容忍度，为0表示达到机器的精度
• which 为'LM'时，选取最大的奇异值；'SM'则选取最小奇异值
• v0 迭代初值
• maxiter 迭代次数
• return_singular_vectors 可选4个值
  • True 返回奇异向量
  • False 不返回奇异向量
  • "u": 如果M <= N，只计算左奇异向量
  • "vh": 如果M > N，只计算右奇异向量；如果 solver='propack'，这个选项将忽略矩阵维度
• solver 可选'arpack', 'propack', 'lobpcg'，但比较吊诡的是，似乎并没有关于这三者区别的文档
• random_state 设置随机数状态
• optionsdict 求解器参数

其返回值有三

• u 即$U$
• s 即奇异值数组，也就是$\Sigma$的对角线
• vh 即$V^T$

测试

下面对奇异值分解做个测试

import numpy as np
from scipy.linalg import svd
from scipy.sparse import csc_array
from scipy.sparse.linalg import svds
np.random.seed(42)  # 设置随机数状态
mat = np.random.rand(500,800)
mat[mat<0.9] = 0
csc = csc_array(mat)
u1, s1, vh1 = svds(csc, k=10)
u2, s2, vh2 = svd(mat)

结果是svds得到的结果和svd的前十个值完全相同，只是排序不一样，但也无关紧要。

下面测试一下二者的时间，由于在Windows下用不了propack，所以svds计算的奇异值数最多只能是$M-1$，也就是499，所以只能测试这个和svd返回500个奇异值的结果相比对，结果如下

>>> from timeit import timeit
>>> timeit(lambda : svds(csc, k=499), number=10)
3.651770199999987
>>> timeit(lambda : svd(mat), number=10)
0.47201400000005833

可见，稀疏矩阵在计算上的确是比不上规整的矩阵。