网站建设各模块功能简述怎么自己建立一个网站
文章目录
- 基本原理
- scipy实现
- 测试
基本原理
当AAA是方阵时,可以很容易地进行特征分解:A=WΣW−1A=W\Sigma W^{-1}A=WΣW−1,其中Σ\SigmaΣ是AAA的特征值组成的对角矩阵。如果WWW由标准正交基组成,则W−1=WTW^{-1}=W^TW−1=WT,特征分解可进一步写成WTΣWW^T\Sigma WWTΣW。
然而,当AAA不是方阵时,情况大不一样了,但仍然可以将AAA表示成A=UΣVTA=U\Sigma V^TA=UΣVT的形式,其中Σ\SigmaΣ也是对角矩阵,对角线上的每个元素被称作奇异值。
奇异值的求解过程和特征值息息相关,因为把AAA变成方阵很简单,只要乘以转置就行。故令L=AATL=AA^TL=AAT,R=ATAR=A^TAR=ATA,则L,RL, RL,R都可以求特征值λi\lambda_iλi和特征向量,其中LLL的特征向量为AAA的左奇异向量,RRR的特征向量为右奇异向量。对应的奇异值σi=λi\sigma_i=\sqrt{\lambda_i}σi=λ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,kmax][k, k_{\max}][k,kmax]之间, 当solver='propack'
时,kmax=min(M,N)k_{max}=\min(M,N)kmax=min(M,N),否则kmax=min(M,N)−1k_{max}=\min(M,N)-1kmax=min(M,N)−1ncv
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 即UUU
- s 即奇异值数组,也就是Σ\SigmaΣ的对角线
- vh 即VTV^TVT
测试
下面对奇异值分解做个测试
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−1M-1M−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
可见,稀疏矩阵在计算上的确是比不上规整的矩阵。