SVD（奇異值分解）Python實現

摘要： 注：在 《SVD（異值分解）小結 》 中分享了SVD原理，但其中只是利用了 numpy.linalg.svd 函式應用了它，並沒有提到如何自己編寫程式碼實現它，在這裡，我再分享一下如何自已寫一個 SVD函式 。但是這裡會利用到SVD的原理，如何大家還不明白它的原理，可以去看看 《SVD（異值...

注：在 ofollow,noindex" target="_blank">《SVD（異值分解）小結 》 中分享了SVD原理，但其中只是利用了 numpy.linalg.svd 函式應用了它，並沒有提到如何自己編寫程式碼實現它，在這裡，我再分享一下如何自已寫一個 SVD函式 。但是這裡會利用到SVD的原理，如何大家還不明白它的原理，可以去看看 《SVD（異值分解）小結 》 ，或者自行百度/google。

1、SVD演算法實現

1.1 SVD原理簡單回顧

有一個 \(m \times n\) 的實數矩陣 \(A\) ，我們可以將它分解成如下的形式

\[ A = U\Sigma V^T \tag{1-1} \]

其中 \(U\) 和 \(V\) 均為單位正交陣，即有 \(UU^T=I\) 和 \(VV^T=I\) ， \(U\) 稱為 左奇異矩陣 ， \(V\) 稱為 右奇異矩陣 ， \(\Sigma\) 僅在主對角線上有值，我們稱它為 奇異值 ，其它元素均為0。上面矩陣的維度分別為 \(U \in \mathbf{R}^{m\times m},\ \Sigma \in \mathbf{R}^{m\times n}\) , \(\ V \in \mathbf{R}^{n\times n}\) 。

正常求上面的 \(U,V,\Sigma\) 不便於求，我們可以利用如下性質

\[ AA^T=U\Sigma V^TV\Sigma^TU^T=U\Sigma \Sigma^TU^T \tag{1-2} \]

\[ A^TA=V\Sigma^TU^TU\Sigma V^T=V\Sigma^T\Sigma V^T \tag{1-3} \]

1.2 SVD演算法

據 1.1小節 ，對式（1-3）和式（1-4）做特徵值分解，即可得到奇異值分解的結果。但是樣分開求存在一定的問題，由於做特徵值分解的時候，特徵向量的正負號並不影響結果，比如，我們利用式（1-3）和（1-4）做特徵值分解

\[ AA^T\mathbf{u}_i = \sigma_i \mathbf{u}_i\quad \text{or} \quad AA^T(-\mathbf{u}_i) = \sigma_i (-\mathbf{u}_i)\\ A^TA\mathbf{v}_i = \sigma_i \mathbf{v}_i\quad \text{or} \quad A^TA(-\mathbf{v}_i) = \sigma_i (-\mathbf{v}_i) \]

如果在計算過程取，取上面的 \(\mathbf{u}_i\) 組成左奇異矩陣 \(U\) ，取 \(-\mathbf{v}_i\) 組成右奇異矩陣 \(V\) ，此時 \(A\ne U\Sigma V^T\) 。因此求 \(\mathbf{v}_i\) 時，要根據 \(\mathbf{u}_i\) 來求，這樣才能保證 \(A= U\Sigma V^T\) 。因此，我們可以得出如下1.1計算SVD的演算法。它主要是先做特性值分解，再根據特徵值分解得到的左奇異矩陣 \(U\) 間接地求出部分的右奇異矩陣 \(V'\in \mathbf{R}^{m\times n}\) 。

1. 計算特徵值：特徵值分解 \(AA^T\) ，其中 \(A \in \mathbf{R}^{m\times n}\) 為原始樣本資料

   \[ AA^T=U\Sigma \Sigma^TU^T \]

   得到左奇異矩陣 \(U \in \mathbf{R}^{m \times m}\) 和奇異值矩陣 \(\Sigma' \in \mathbf{R}^{m \times m}\)

2. 間接求部分右奇異矩陣：求 \(V' \in \mathbf{R}^{m \times n}\)

   利用 \(A=U\Sigma'V'\) 可得

   \[ V' = (U\Sigma')^{-1}A = (\Sigma')^{-1}U^TA \tag{1-4} \]

3. 返回 \(U,\ \Sigma',\ V'\) ，分別為 左奇異矩陣，奇異值矩陣，右奇異矩陣。

注：這裡得到的 \(\Sigma'\) 和 \(V'\) 與式（1-2）所得到的 \(\Sigma,\ V\) 有區別，它們的維度不一樣。 \(\Sigma'\) 是隻取了前 \(m\) 個奇異值形成的對角方陣，即 \(\Sigma' \in \mathbf{R}^{m \times m}\) ； \(V'\) 不是一個方陣，它只取了 \(V \in \mathbf{R}^{m \times n}\) 的前 \(m\) 行（假設 \(m < n\) ），即有 \(V' = V(:m,\cdot)\) 。這樣一來，我們同樣有類似式（1-1）的數學關係成立，即

\[ A = U\Sigma' (V')^T\tag{1-5} \]

我們可以利用此關係重建原始資料。

2、SVD的Python實現

以下程式碼的執行環境為 python3.6 + jupyter5.4 。

2.1 SVD實現過程

讀取資料

這裡面的資料集大家隨便找一個數據就好，如果有需要我的資料集，可以下在面留言。

import numpy as np
import pandas as pd
from scipy.io import loadmat

# 讀取資料，使用自己資料集的路徑。
train_data_mat = loadmat("../data/train_data2.mat")
train_data = train_data_mat["Data"]
print(train_data.shape)

特徵值分解

# 資料必需先轉為浮點型，否則在計算的過程中會溢位，導致結果不準確
train_dataFloat = train_data / 255.0
# 計算特徵值和特徵向量
eval_sigma1,evec_u = np.linalg.eigh(train_dataFloat.dot(train_dataFloat.T))

計算右奇異矩陣

#降序排列後，逆序輸出
eval1_sort_idx = np.argsort(eval_sigma1)[::-1]
# 將特徵值對應的特徵向量也對應排好序
eval_sigma1 = np.sort(eval_sigma1)[::-1]
evec_u = evec_u[:,eval1_sort_idx]
# 計算奇異值矩陣的逆
eval_sigma1 = np.sqrt(eval_sigma1)
eval_sigma1_inv = np.linalg.inv(np.diag(eval_sigma1))
# 計算右奇異矩陣
evec_part_v = eval_sigma1_inv.dot((evec_u.T).dot(train_dataFloat))

上面的計算出的 evec_u, eval_sigma1, evec_part_v 分別為左奇異矩陣，所有奇異值，右奇異矩陣。

2.2 SVD降維後重建資料

取不同個數的奇異值，重建圖片，計算出均方誤差，如圖2-1所示。從圖中可以看出，隨著奇異值的增加，均方誤差（MSE）在減小，且奇異值和的比率正快速上升，在100維時，奇異值佔總和的53%。

圖2-1 奇值分解維度和均方誤差變化圖

注：均方誤差MSE有如下計算公式

\[ \text{MSE} = \frac{1}{n}\left((y_1-y_1')^2+(y_2-y_2')^2+\cdots+(y_n-y_n')^2\right) \]

\(\text{RMSE}=\sqrt{\text{MSE}}\)

。

將圖和10、50、100維的圖進行比較，如圖3-2所示。在直觀上，100維時，能保留較多的資訊，此時能從圖片中看出車輛形狀。

總結

SVD 與特徵值分解（EVD）非常類似，應該說EVD只是SVD的一種特殊懷況。我們可以通這它們在實際的應用中返過來理解 特徵值/奇異值 的含義：它代表著資料的資訊量，它的值越大，資訊越多。

最近作業是真的多呀，冒著生命危險來分享，希望能給大家帶來幫助:smile: