SVD（奇異值分解）Python實現

注：在《SVD（奇異值分解）小結 》中分享了SVD原理，但其中只是利用了numpy.linalg.svd函數應用了它，並沒有提到如何自己編寫代碼實現它，在這里，我再分享一下如何自已寫一個SVD函數。但是這里會利用到SVD的原理，如果大家還不明白它的原理，可以去看看《SVD（奇異值分解）小結 》，或者自行百度/google。數據集：https://pan.baidu.com/s/1ZmpUSIscy4VltcimwwIWew。

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.1：SVD

輸入：樣本數據
輸出：左奇異矩陣，奇異值矩陣，右奇異矩陣

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維的圖進行比較，如圖2-2所示。在直觀上，100維時，能保留較多的信息，此時能從圖片中看出車輛形狀。

圖2-2 原圖與降維重建后的圖比較

總結

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

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