矩阵的奇异值分解（SVD）（理论）
矩阵的奇异值分解(Singular Value Decomposition,SVD)是数值计算中的精彩之处，在其它数学领域和机器学习领域得到了广泛的应用，如矩阵的广义逆，主分成分析(PCA),自然语言处理(NLP)中的潜在语义索引（Latent Semantic Indexing）,推荐算法等。
鉴于实际应用，本次分享中的数域为实数域，即我们只在实数范围内讨论。我们假定读者具有大学线性代数的水平。那么，矩阵的奇异值分解定理如下：
（定理）（奇异值分解定理）任意一个$m \times n$矩阵A可分解为
$$A=PDQ$$
其中P是$m \times m$正交矩阵，D是$m \times n$对角阵，Q是$n \times n$正交矩阵。
证明：矩阵$A^{T}A$是$n \times n$对称矩阵，因为$(A^{T}A)^{T}=A^{T}(A^{T})^{T}=A^{T}A$.又因为
$$x^{T}（A^{T}A）x=(Ax)^{T}(Ax)\ge0,$$
所以$A^{T}A$是半正定矩阵，从而，$A^{T}A$的特征值为非负数。
假设$A^{T}A$的特征值为$\sigma_{1}^{2},\sigma_{2}^{2},...,\sigma_{n}^{2}$,其中，$\sigma_{1}^{2},\sigma_{2}^{2},...,\sigma_{r}^{2}$都是正的，$\sigma_{r+1}^{2},\sigma_{r+2}^{2},...,\sigma_{n}^{2}$都是0，$r$为$A^{T}A$的秩。设$\{u_{1},u_{2},...,u_{n}\}$为$A^{T}A$的标准正交特征向量集，则
$$A^{T}Au_{i}=\sigma_{i}^{2}u_{i} (i=1,2,...,n)$$
于是$(Au_{i})^{T}(Au_{i})=u_{i}^{T}(A^{T}A)u_{i}=u_{i}^{T}\sigma_{i}^{2}u_{i}=\sigma_{i}^{2}.$当$i\ge r+1$时，$\sigma_{i}=0$,从而$Au_{i}=0$.
用$\{u_{1}^{T},u_{2}^{T},...,u_{n}^{T}\}$作为行构成一个$n\times n$矩阵$Q$.接着，定义
$$v_{i}=\sigma_{i}^{-1}Au_{i} (1\le i \le r).$$
当$1\le i,j \le r$时，$v_{i}$构成一个标准正交系，这是因为
$$v_{i}^{T}v_{j}=\sigma_{i}^{-1}(Au_{i})^{T}\sigma_{j}^{-1}(Au_{j})=(\sigma_{i}\sigma_{j})^{-1}(u_{i}^{T}A^{T}Au_{j})=(\sigma_{i}\sigma_{j})^{-1}(u_{i}^{T}\sigma_{j}^{2}u_{j})=\delta_{ij},$$
其中$\delta_{ij}$为Kronecker符号，即当$i=j$时，$\delta=1$,当$i\neq j$时，$\delta=0$.
我们选择额外的向量$v_{i}$使得$\{v_{1},v_{2},...,v_{m}\}$为$\mathbb{R}^{m}$的标准正交基。设P是$m\times m$矩阵，其列是$v_{1},v_{2},...,v_{m}$.设D是$m\times n$对角阵，$\sigma_{1},\sigma_{2},...\sigma_{r}$在其对角线上，其余地方均为0.于是有
$$A=PDQ.$$
这是因为$(P^{T}AQ^{T})_{ij}=v_{i}^{T}Au_{j}$,当$j\ge r+1$时，该式为0，当$j\le r$时，该式为$v_{i}^{T}\sigma v_{j}=\sigma_{j}\delta_{ij}$,从而$P^{T}AQ^{T}=D$.又因$P,Q$为正交矩阵，因此$$A=PDQ.$$
证毕。
在上面证明中，我们称实数$\sigma_{1},\sigma_{2},...,\sigma_{n}$(取非负数)为矩阵A的奇异值，它们是$A^{T}A$的特征值的非负平方根。定理中的分解$A=PDQ$就是一个奇异值分解。由上面的证明，我们可以知道：矩阵的奇异值分解并不唯一，因为$\sigma_{1},\sigma_{2},...,\sigma_{n}$的次序及$v_{r+1},v_{r+2},...,v_{n}$的选择并不唯一。
在Python中的Numpy模块中，已经实现了矩阵的奇异值分解。以下为示例的应用代码：
1 import numpy as np 2 #generate a random 3*4 matrix 3 A = np.random.randint(5, size=(3, 4)) 4 #parameter full_matrices: control the size of P and Q 5 #d returns as numpy.ndarray, not matrix 6 P,d,Q = np.linalg.svd(A, full_matrices=True) 7 print('A:',A) 8 print('P:',P) 9 #D return as diagonal 3*4 matrix 10 D = np.zeros(12).reshape(3,4) 11 for i in range(len(d)): 12 D[i][i] = d[i] 13 print('D:',D) 14 print('Q:',Q) 15 #check if P*D*Q == A 16 print('P*D*Q:',np.dot(P,np.dot(D,Q)))
输入结果如下：
至于如何用原始算法来实现矩阵的SVD，也是需要考虑的，有机会的话，可以交流哦~~
本次分享到此结束，欢迎大家批评与交流~~
参考文献：
- SVD 维基百科：https://en.wikipedia.org/wiki/SVD
- 数值分析 机械工业出版社 作者：萨奥尔(Timothy Sauer) 译者：裴玉茹
- numpy的svd实现函数： https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.svd.html
- 奇异值分解(SVD)原理与在降维中的应用：https://www.cnblogs.com/pinard/p/6251584.html
- 奇异值分解SVD应用——LSI：http://blog.csdn.net/abcjennifer/article/details/8131087
- 论文：CALCULATING THE SINGULAR VALUES AND PSEUDO-INVERSE OF A MATRIX, G. GOLUB AND W. KAHAN, J. SIAM llrM,B. AfeArd.Ser. B, Vol. 2, No. 2, 1965
