矩阵的奇异值分解（SVD）（理论）

　　矩阵的奇异值分解(Singular Value Decomposition,SVD)是数值计算中的精彩之处，在其它数学领域和机器学习领域得到了广泛的应用，如矩阵的广义逆，主分成分析(PCA),自然语言处理(NLP)中的潜在语义索引（Latent Semantic Indexing）,推荐算法等。

　　鉴于实际应用，本次分享中的数域为实数域，即我们只在实数范围内讨论。我们假定读者具有大学线性代数的水平。那么，矩阵的奇异值分解定理如下：

（定理）（奇异值分解定理）任意一个$m imes n$矩阵A可分解为

$$A=PDQ$$

其中P是$m imes m$正交矩阵，D是$m imes n$对角阵，Q是$n imes n$正交矩阵。

证明：矩阵$A^{T}A$是$n imes 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}.$当$ige r+1$时，$sigma_{i}=0$,从而$Au_{i}=0$.

       用${u_{1}^{T},u_{2}^{T},...,u_{n}^{T}}$作为行构成一个$n imes n$矩阵$Q$.接着，定义

$$v_{i}=sigma_{i}^{-1}Au_{i} (1le i le r).$$

当$1le 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 eq j$时，$delta=0$.

       我们选择额外的向量$v_{i}$使得${v_{1},v_{2},...,v_{m}}$为$mathbb{R}^{m}$的标准正交基。设P是$m imes m$矩阵，其列是$v_{1},v_{2},...,v_{m}$.设D是$m imes n$对角阵，$sigma_{1},sigma_{2},...sigma_{r}$在其对角线上，其余地方均为0.于是有

$$A=PDQ.$$

这是因为$(P^{T}AQ^{T})_{ij}=v_{i}^{T}Au_{j}$,当$jge r+1$时，该式为0，当$jle 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模块中，已经实现了矩阵的奇异值分解。以下为示例的应用代码：

import numpy as np
#generate a random 3*4 matrix 
A =  np.random.randint(5, size=(3, 4))
#parameter full_matrices: control the size of P and Q
#d returns as numpy.ndarray, not matrix 
P,d,Q = np.linalg.svd(A, full_matrices=True)
print('A:',A)
print('P:',P)
#D return as diagonal 3*4 matrix
D = np.zeros(12).reshape(3,4)
for i in range(len(d)):
    D[i][i] = d[i]
print('D:',D)
print('Q:',Q)
#check if P*D*Q == A
print('P*D*Q:',np.dot(P,np.dot(D,Q)))

输入结果如下：

　　至于如何用原始算法来实现矩阵的SVD，也是需要考虑的，有机会的话，可以交流哦~~

　　本次分享到此结束，欢迎大家批评与交流~~

---

参考文献：

1. SVD 维基百科：https://en.wikipedia.org/wiki/SVD
2. 数值分析  机械工业出版社 作者：萨奥尔(Timothy Sauer)  译者：裴玉茹 
3. numpy的svd实现函数： https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.svd.html
4.  奇异值分解(SVD)原理与在降维中的应用：https://www.cnblogs.com/pinard/p/6251584.html
5.  奇异值分解SVD应用——LSI：http://blog.csdn.net/abcjennifer/article/details/8131087
6. 论文：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