矩陣的奇異值分解(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，也是需要考慮的，有機會的話，可以交流哦~~

　　本次分享到此結束，歡迎大家批評與交流~~

參考文獻：

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