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

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

线性代数里的SVD

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

$A{A}^T=US{V}^T(US{V}^T{)}^T=US{V}^{T}V{S}^T{U}^T=US{S}^T{U}^T$
同样的：
$A^{T}A=(US{V}^T{)}^{T}US{V}^T=V{S}^T{U}^{T}US{V}^T=V{S}^{T}S{V}^T$
说明：
$S{S}^T和S^{T}S$都是对角阵，其对角线上的元素是奇异值的平方，但$S{S}^T和S^{T}S$的维度是不同的，$S{S}^T$是mxm的，而$S^{T}S$是nxn的。这是，我们对于A的奇异值个数应该是取m和n中较小的一个，也就是说多出来的值需要被舍弃，因为我们最终需要的$S$是一个mxn的矩阵，主对角线上的奇异值的个数只能是min(m,n)。

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

但是，要注意的是当我们已经求得左奇异矩阵U后，再求右奇异矩阵V时不能再使用$A^{T}A$了，因为对于同一个特征值其特征向量是不唯一的，如果$u_i$是$A^{T}A$对应于${\lambda }_i$特征向量，那么$-u_i$也是其特征向量。因为：
$A(u_i)={\lambda }_i(u_i),A(-u_i)={\lambda }_i(-u_i)$
因此，如果我们使用$A^{T}A$来求右奇异矩阵V，你很可能会遇到一组或几组特征向量的符号是反的，从而导致你求得的$US{V}^T$不等于$A$。
那么为了确保得到正确的右奇异矩阵V，我们在求得左奇异矩阵U后，直接利用$A=US{V}^T,V^T=S^{-1}{U}^{-1}A$来得到对应的$V^T$。此处，$S^{-1}{U}^{-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中第二行的值是正负相反的。但是如果我们验证一下结果，算一下$US{V}^T$的结果正好也等于$A$。这其实也算是分解成功的，只是有些情况下，向量值会正负相反。

sklearn里的TruncatedSVD方法

最后，我想说一下sklearn里的TruncatedSVD方法，这个方法创建是设置参数n_components可以直接实现降维，他会舍弃你设置的维度以后的所有数值。
比如说n_components=2，那么他返回给你的就是 $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）