[数学基础知识] 线代里的svd, numpy 的svd以及sklearn的TruncatedSVD

SVD的全称Singular Value Decomposition，中文名是矩阵的奇异分解。它是一种常见的做矩阵降维处理的算法，在图像压缩和NLP算法中经常被用到。本文是我在编程过程中，对于数学中的SVD，numpy的svd方法，以及sklearn中的TruncatedSVD方法在实际应用中的一些理解和体会。

线性代数里的SVD

A是mxn的实数矩阵, 则A可以分解成以下的形式
A = U S V T A=USV^T A=USVT ，其中U是mxm的矩阵，S是mxn的对角阵，其主对角线上的每一个值被称为奇异值，V是nxn的矩阵。U和V都是酉矩阵（Unity Matrix），即 U T U = I , V T V = I U^TU=I, V^TV=I UTU=I,VTV=I。U称为左奇异矩阵，V称为右奇异矩阵。
线代里求svd分解的方法主要用到了求特征向量和特征值得方法。
对于矩阵的特征值和特征向量性质，我们可知如果方阵 X m x m X_{mxm} Xmxm如果存在m个线性无关的特征向量，可以做如下分解:
X = Q Λ Q − 1 X=Q \varLambda Q^{-1} X=QΛQ−1其中， Q Q Q是由其特征值组成的非奇异矩阵（即可逆）， Λ \varLambda Λ是由其特征值组成的对角阵，其对角线上的每一个值对应于一个特征值。

A A T = U S V T ( U S V T ) T = U S V T V S T U T = U S S T U T AA^T=USV^T(USV^T)^T=USV^TVS^TU^T=USS^TU^T AAT=USVT(USVT)T=USVTVSTUT=USSTUT
同样的：
A T A = ( U S V T ) T U S V T = V S T U T U S V T = V S T S V T A^TA=(USV^T)^TUSV^T=VS^TU^TUSV^T=VS^TSV^T ATA=(USVT)TUSVT=VSTUTUSVT=VSTSVT
说明：
S S T 和 S T S SS^T和S^TS SST和STS都是对角阵，其对角线上的元素是奇异值的平方，但 S S T 和 S T S SS^T和S^TS SST和STS的维度是不同的， S S T SS^T SST是mxm的，而 S T S S^TS STS是nxn的。这是，我们对于A的奇异值个数应该是取m和n中较小的一个，也就是说多出来的值需要被舍弃，因为我们最终需要的 S S S是一个mxn的矩阵，主对角线上的奇异值的个数只能是min(m,n)。

此时，我们发现可以通过求 A A T AA^T AAT的特征值和特征向量的方法来求U和S。（求特征值和特征向量的方法此处略过）
求得 A A T AA^T AAT的特征值，对其求开根号就可以得到我们需要的奇异值。特征向量组成的mxm的矩阵就是左奇异矩阵。

但是，要注意的是当我们已经求得左奇异矩阵U后，再求右奇异矩阵V时不能再使用 A T A A^TA ATA了，因为对于同一个特征值其特征向量是不唯一的，如果 u i u_i ui是 A T A A^TA ATA对应于 λ i \lambda _i λi特征向量，那么 − u i -u_i −ui也是其特征向量。因为：
A ( u i ) = λ i ( u i ) , A ( − u i ) = λ i ( − u i ) A(u_i) = \lambda_i(u_i), A(-u_i) = \lambda_i(-u_i) A(ui)=λi(ui),A(−ui)=λi(−ui)
因此，如果我们使用 A T A A^TA ATA来求右奇异矩阵V，你很可能会遇到一组或几组特征向量的符号是反的，从而导致你求得的 U S V T USV^T USVT不等于 A A A。
那么为了确保得到正确的右奇异矩阵V，我们在求得左奇异矩阵U后，直接利用 A = U S V T , V T = S − 1 U − 1 A A=USV^T, V^T=S^{-1}U^{-1}A A=USVT,VT=S−1U−1A来得到对应的 V T V^T VT。此处， S − 1 U − 1 S^{-1}U^{-1} S−1U−1是伪逆阵，因为S和U都不是方阵，理论上不存在逆矩阵，但是伪逆阵依旧可以求得。

Numpy里的svd方法

在Numpy的linalg库里提供了现成的svd方法，可以直接分解出U，S，VT。
linalg.svd方法使用起来很简单，只要输入一个矩阵，三个返回值分别是左奇异矩阵，对应的奇异值列表,和右奇异矩阵的转置。而且奇异值是按照大小顺序排列好的。这也是推荐使用的方法，效率比较高。
python代码如下：

import numpy as np
A=np.asarray([[6,1,9],[4,9,9],[4,0,6],[7,4,7]],dtype=np.float32)
U, S, VT=np.linalg.svd(A,full_matrices=1)
print('A:')
print(A)
print('U:')
print(U)
print('S:')
print(S)
print('VT:')
print(VT)
output=========
A:
[[6. 1. 9.][4. 9. 9.][4. 0. 6.][7. 4. 7.]]
U：
[[-0.50626906 -0.5132768  -0.33109968][-0.61136554  0.73364755 -0.29256487][-0.32459939 -0.43574965 -0.28201245][-0.51435304 -0.09181746  0.85161481]]
S:
[20.4062792   6.38531027  2.1935318 ]
VT:
[[-0.50876229 -0.39526907 -0.76480278][-0.39634588  0.89616339 -0.19950219][ 0.7642453   0.20162724 -0.61259741]]

这里要提一下的是full_matrices这个参数，默认值是1，即如果输入矩阵是mxn的，则返回的左奇异矩阵是mxm的，右奇异矩阵是nxn的。
如果full_matrices=0，则返回的左奇异矩阵是mxk的，右奇异矩阵是kxn的，而k是m和n中较小的那个。画个表格来看就比较明白了。

	左奇异矩阵	右奇异矩阵
m>n	mxn	nxn
m<n	mxm	mxn

当m>n时，左奇异矩阵是mxn的，右奇异矩阵是xn的。
那如果我不想用这个现成的方法，能不能从求特征值得方法来求分解呢，也是可以的，而且计算结果是一样的。下面直接上代码：

import numpy as np
A=np.asarray([[6,1,9],[4,9,9],[4,0,6],[7,4,7]],dtype=np.float32)
#先求左奇异矩阵和对应的奇异值
#np.linalg.eigh是求特征值和特征向量的方法
#第一个返回值是特征值列表
#第二个返回值是特征向量组成的矩阵
EIGENV_U,U = np.linalg.eigh(np.dot(A,A.T))
#np.linalg.eigh的返回值不会像svd方法一样按特征值大小排好序，所以必须自己排序
#获得对奇异值的排序的位置
EIGENV_U_sort_idx = np.argsort(EIGENV_U)[::-1]
#重新对奇异值按大小排序
EIGENV_U = EIGENV_U[EIGENV_U_sort_idx]
#对奇异矩阵里的特征向量也进行一样的顺序调整
U=U[:,EIGENV_U_sort_idx]
#对拍好序的特征值求平方差，得到奇异值
S = np.sqrt(EIGENV_U)
#这里有个细节处理需要注意
#因为AxAT是mxm的，我们会得到m个特征值，如果m>n,我们实际上只需要取前n个
#因为求逆阵必须是满秩的，因此先对对角方阵求逆，
#再将其填入一个nxm的矩阵中，A是mxn的，此处逆阵的行列数记得反一下
S= np.diag(S[:min(A.shape[0],A.shape[1])])
S= np.linalg.inv(S)
S_diag = np.zeros((A.shape[1],A.shape[0]))
S_diag[:min(A.shape[0],A.shape[1]),:min(A.shape[0],A.shape[1])] = S
VT=np.dot(np.dot(S_diag,np.linalg.inv(U)),A)
print('A:')
print(A)
print('U:')
print(U)
print('S:')
print(S)
print('VT:')
print(VT)
output=========
A:
[[6. 1. 9.][4. 9. 9.][4. 0. 6.][7. 4. 7.]]
U:
[[-0.50626904  0.5132768  -0.3310997  -0.6087788 ][-0.61136556 -0.7336475  -0.29256487  0.0489822 ][-0.3245994   0.43574965 -0.28201243  0.79071265][-0.51435304  0.09181746  0.85161483  0.04198474]]
S:
[20.406279  6.38531   2.193532  0.      ]
VT：
[[-0.5087623  -0.39526909 -0.76480279][ 0.39634593 -0.89616334  0.19950226][ 0.76424524  0.20162717 -0.61259742]]

从求得结果来看和直接用svd是一样的，但是VT中第二行的值是正负相反的。但是如果我们验证一下结果，算一下 U S V T USV^T USVT的结果正好也等于 A A A。这其实也算是分解成功的，只是有些情况下，向量值会正负相反。

sklearn里的TruncatedSVD方法

最后，我想说一下sklearn里的TruncatedSVD方法，这个方法创建是设置参数n_components可以直接实现降维，他会舍弃你设置的维度以后的所有数值。
比如说n_components=2，那么他返回给你的就是 U S US US 得到的矩阵的前两列。
来看个代码例子就更清楚了。

from sklearn.decomposition import TruncatedSVD
svd = TruncatedSVD(n_components=2)
M_reduced=svd.fit_transform(A)
print('sklearn TruncatedSVD')
print(M_reduced)
#U是之前用svd直接求得的左奇异矩阵
#S_diag是由奇异值构成的mxn对角矩阵
print('np svd U x S_diag, get left 2 cols')
print(np.dot(U,S_diag))[:,:2])
output=========
sklearn TruncatedSVD
[[10.331068  -3.27743  ][12.475696   4.6845646][ 6.6238675 -2.7823956][10.496031  -0.5862823]]np svd U x S_diag, get left 2 cols
[[-10.33106704  -3.27743167][-12.47569588   4.68456701][ -6.62386549  -2.78239667][-10.49603137  -0.58628297]]

你看，值是不是一模一样。使用TruncatedSVD可以直接提取到我们想保留的最重要的特征维度，舍弃不需要的数据，提升计算性能，在处理大量数据是非常的有用。

参考文章
[1] 奇异值分解原理和应用（SVD和TruncatedSVD）
[2]使用Python求解特征值、特征向量及奇异值分解（SVD）