Python機器學習筆記：奇異值分解（SVD）演算法

阿新 • • 發佈：2021-01-19

完整程式碼及其資料，請移步小編的GitHub

　　傳送門：請點選我

　　如果點選有誤：https://github.com/LeBron-Jian/MachineLearningNote

　　奇異值分解（Singular Value Decomposition，後面簡稱 SVD）是線上性代數中一種重要的矩陣分解，它不光可用在降維演算法中（例如PCA演算法）的特徵分解，還可以用於推薦系統，以及自然語言處理等領域，在機器學習，訊號處理，統計學等領域中有重要應用。

　　比如之前的學習的PCA，掌握了SVD原理後再去看PCA是非常簡單的，因為我最近在整理學習線性代數基礎，並溫習了一遍特徵值與特徵向量，所以這裡同時也溫習一下SVD演算法，如果想看我之前的PCA降維的文章，可以參考下面：

Python機器學習筆記：主成分分析（PCA）演算法

Python機器學習筆記：使用scikit-learn工具進行PCA降維

　　好了，話不多說，這裡再對SVD演算法做一個總結（這裡總結主線是參考劉建平大佬和老師的網課視訊學習，首先在這裡提出感謝）

1，基變換與特徵值分解

1.1 向量的表示與基變換

　　首先看看向量，向量可以表示為（3, 2），實際上表示線性組合：

　　基表示：（1, 0）和（0, 1）叫做二維空間中的一組基。

　　一組基就確定了對向量進行表示的空間，而向量的表示就是基於這組基進行座標表示。選定的基不同，空間自然不同，那麼向量表示的座標自然也不同。一般情況下，我們會預設傳統的正交座標系為標準的基選擇，但其實對傳統正交座標系進行選擇，放縮或剪下變換後得到的依然是一個座標系。

　　基變換：資料與一個基做內積運算，結果作為第一個新的座標分量，然後與第二個基做內積運算，結果作為第二個新座標的分量。

　　將一組基A下的座標轉換成另一組集B下的座標——核心思想就是將A中的基用B對應的座標系進行表示，將該向量在A下的座標表示變換成關於基A的線性組合，進行代換即可。

　　資料（3， 2）對映到基中座標為：

1.2 特徵值分解

　　首先回歸下特徵值和特徵向量的定義，其定義如下：

　　其中A是一個 n * n 的矩陣，x 是一個 n 維向量，則我們說 λ 是矩陣A的一個特徵值，而 x 是矩陣A的特徵值 λ 所對應的特徵向量。

　　求出特徵值和特徵向量有什麼好處呢？就是我們可以將矩陣 A 特徵分解。如果我們求出了矩陣 A 的 n 個特徵值 λ₁

<= λ₂ <= λ_n ，以及這 n 個特徵值所對應的特徵向量 {W₁, W₂, ...W_n}，如果這 n 個特徵向量線性無關，那麼矩陣 A 就可以用下式的特徵分解表示：

　　其中 W 是這 n 個特徵向量所長成的 n * n 維矩陣，而 Σ 為這個 n 個特徵值為主對角線的 n * n 維矩陣。

　　一般我們會把 W 的這 n 個特徵向量標準化，即滿足 ||w_i||₂ = 1，或者說 w_i^Tw_i=1，此時 W 的 n 個特向量為標準正交基，滿足W^TW = I，即 W^T = W^-1。

　　這樣我們特徵分解表示式可以寫成：

　　注意到要進行特徵分解，矩陣A必須為方陣。那麼如果A不是方陣，即行和列不相同時，我們還可以對矩陣進行分解嗎？答案是可以，就是下面我們要學習的SVD演算法。

2，特徵分解和奇異值分解（SVD）的區別

　　特徵值分解和奇異值分解的目的都是一樣，就是提取出一個矩陣最重要的特徵。

　　所有的矩陣都可以進行奇異值分解，但只有方陣才可以進行特徵值分解。當所給的矩陣是對稱的方陣，A^T=A，二者的結果是相同的。也就是說對稱矩陣的特徵值分解是所有奇異值分解的一個特例。但是二者還是存在一些小的差異，奇異值分解需要對奇異值進行從大到小的排序，而且全部是大於等於零。

　　對於如上的奇異值分解公式，我們可以看到奇異值分解同時包含了旋轉，縮放和投影三種作用，並且U和V都起到了對A進行旋轉的作用，而 Σ 起到了對 A 縮放的作用，特徵值分解只有縮放的效果。

　　在應用上，奇異值分解主要用於資料矩陣，而特徵值分解主要用於方型的相關矩陣。

3，奇異值分解的定義

　　特徵值分解釋一個提取矩陣特徵很不錯的方法，但是它只適用於方陣，而在現實的世界中，我們看到的大部分矩陣都不是方陣，比如說有M個學生，每個學生有N個成績，這樣就形成了一個M*N的矩陣就可能不是方陣了，我們怎麼樣才能像描述特徵值一樣描述這樣一般矩陣的重要特徵呢？奇異值分解就是幹這個事情，奇異值分解是一個能適用於任何的矩陣的一種分解的方法。

3.1 SVD的理論描述

　　假設 A 是一個 m*n 階矩陣，其中的元素全部屬於域 K，也就是實數域或複數域。如此則存在一個分解使得：

　　其中 U 是 m*m 階酋矩陣，Σ 是半正定 m*n 階對角矩陣，而 V* 是V的共軛轉置，是 n*n 階酋矩陣，這樣的分解就稱作 M 的奇異值分解。Σ（是對角矩陣，但不一定是方陣）對角線上的元素 Σi 即為 M 的奇異值。

　　一般的 Σ 有如下形式：

　　上述分解中會構建出一個矩陣 Σ ，該矩陣只有對角元素，其他元素均為0，另一個慣例是，Σ 的對角元素是從大到小排列的。這些對角元素稱為奇異值。在科學和工程中，一直存在這個一個普遍事實：在某個奇異值的數目（r個）之後，其他奇異值都置為零，這就意味著資料集中僅有 r 個重要特徵，其餘的都是噪聲或冗餘資料。

　　從下圖可以形象的看出上面SVD的定義：

　　那麼我們如何求出SVD分解後的 U，Σ，V 這三個矩陣呢？

　　如果我們將 A 的轉置和 A 做矩陣乘法，那麼會得到 n*n 的一個方陣 A^TA。既然 A^TA 是方陣，那麼我們就可以進行特徵分解，得到的特徵值和特徵向量滿足下式：

　　這樣我們就可以得到矩陣 A^TA 的 n 個特徵值和對應的 n 個特徵向量 v 了。將 A^TA 的所有特徵向量張成一個 n*n 的矩陣 V，就是我們SVD公式裡面的V矩陣了。一般我們將 V 中的每個特徵向量叫做 A 的右奇異向量。

　　如果我們將 A 和 A 的轉置做矩陣乘法，那麼會得到 m*m的一個方陣 AA^T。既然 AA^T 是方陣，那麼我們就可以進行特徵分解，得到的特徵值和特徵向量滿足下面公式：

　　這樣我們就可以得到矩陣 AA^T 的 m 個特徵值和對應的 m 個特徵向量 u了。將 AA^T 的所有特徵張成一個 m*m 的矩陣 U，就是我們 SVD 公式裡面的 U矩陣了。一般我們將 U 中的每個特徵向量叫做 A 的左奇異向量。

　　U 和 V 我們都求出來了，現在就剩下奇異值矩陣 Σ 沒有求出了。

　　由於 Σ 除了對角線上是奇異值其他位置都是0，所以我們只需要求出每個奇異值 σ 就可以了。

　　我們注意到：

　　這樣我們可以求出我們的每一個奇異值，進而求出奇異值矩陣 Σ。

　　上面還有一個問題就是我們說 A^TA 的特徵向量組成的就是我們SVD中的 V矩陣，而 AA^T的特徵向量組成的就是我們 SVD中的 U矩陣，這有什麼根據嗎？這個其實很容易證明，我們以V矩陣的證明為例：

　　上式證明使用了：U^TU = I, Σ^TΣ = Σ²。可以看出 A^TA 的特徵向量組成的的確就是我們SVD中的 V 矩陣。類似的方法可以得到 AA^T 的特徵向量組成的就是我們 SVD中的 U 矩陣。

　　進一步我們還可以看出我們的特徵值矩陣等於奇異值矩陣的平方，也就是說特徵值和奇異值滿足如下關係：

　　這樣也就是說，我們可以不用上面推導的式子來計算奇異值，上式如下：

　　我們可以直接通過求 A^TA 的特徵值取平方根來求奇異值。

注意1：通過A^TA的特徵值取平方根求解奇異值的注意點

　　我們知道，正常求 U，V，Σ 不便求，所以，我們利用如下性質（公式上面有提到，這裡再複述一下）：

　　需要注意的是：這裡的 ΣΣ^T 和 Σ^TΣ 在矩陣的角度上來講，他們是不相等的，因為他們的維度不同 ΣΣ^T €R^m*m，而 Σ^TΣ€R^n*n，但是他們在主對角線的奇異值是相等的，即有：

　　所以對於 ΣΣ^T 和 Σ^TΣ 中的特徵值開方，可以得到所有的奇異值。

注意2：酋矩陣的定義

　　若n行n列的複數矩陣 U滿足：

　　其中 In 為 n 階單位矩陣，U^T 為 U的共軛轉置，則U為酋矩陣（即矩陣U為酋矩陣，當且僅當其共軛轉置U^T 為其逆矩陣：U^-1 = U^T）

3.2 SVD的幾何層面理解

　　下面從幾何層面上去理解二維的SVD：對於任意的 2*2 矩陣，通過 SVD 可以將一個互相垂直的網路（orthogonal grid）變換到另外一個互相垂直的網路。

　　我們可以通過向量的方式來描述這個事實：首先，選擇兩個互相正交的單位向量 v1 和 v2，向量 Mv1 和 Mv2 正交。

　　u1 和 u2 分別表示 Mv1 和 Mv2 的單位向量，σ1*u1=Mv1 和 σ2*u2=Mv2，σ1 和 σ2分別表示這不同方向向量上的模，也稱為作為矩陣 M 的奇異值。

　　這樣就有了如下關係式：

　　我們現在可以簡單描述下經過 M 線性變換後的向量 x 的表達形式。由於向量 v1 和 v2 是正交的單位向量，我們可以得到如下式子：

　　這就意味著：

　　向量內積可以用向量的轉置來標色，如下所示：

　　最終的式子為：

　　上述的式子經常表示為：

　　u 矩陣的列向量分別是 u1 和 u2，Σ 是一個對角矩陣，對角元素分別是對應的 σ1 和 σ2，v 矩陣的列向量分別為 v1 和 v2。

　　這就表明了任意的矩陣 M 是可以分解成三個矩陣，v 表示原始域的標準正交基， u 表示經過 M 變換後的 co-domain 的標準正交基，Σ 表示了 v 中的向量與 u 中相對應向量之間的關係。

3.3 SVD的推導

　　這裡直接拿了百度百科的內容。

3.3 SVD的一些性質

　　上面我們對SVD的定義和計算做了詳細的描述，下面對其一些性質進行分析。

　　對於奇異值，它和我們特徵分解中的特徵值類似，在奇異值矩陣中也是按照從大到小排列，而且奇異值的減少特別的塊，在很多情況下，前10%甚至1%的奇異值的和就佔了全部的奇異值之和的99%以上的比例。也就是說，我們可以用最大的k個奇異值和對應的左右奇異值向量來近似描述矩陣。也就是說：

　　其中 k 要比 n 小很多，也就是一個大的矩陣 A 可以用三個小的矩陣來標色，如下圖所示，現在我們的矩陣 A 只需要灰色的部分的三個小矩陣就可以近似描述了。

　　由於奇異值的特點就是衰減特別快，也就是說僅僅有前幾個奇異值特別大，後面的值都很小，這一點可以用於資料壓縮和去噪，PCA降維，也可以用於推薦演算法，將使用者和喜好對應的矩陣分解，進而得到隱含的使用者需求來做推薦。

4，SVD演算法的應用

4.1 矩陣近似值

　　奇異值分解在統計中的主要應用為主成分分析（PCA），一種資料分析方法，用來找出大量資料中所隱含的“模式”，PCA演算法的作用是把資料集對映到低維空間中去。資料集的特徵值（在SVD中用奇異值表徵）按照重要性排列，降維的過程就是捨棄不重要的特徵向量的過程，而剩下的特徵向量組成的空間即為降維後的空間。

　　在PCA降維中，需要找到樣本協方差矩陣 X^TX 的最大的 d 個特徵向量，然後用這最大的 d 個特徵向量張成的矩陣來做低維投影降維。可以看出，在這個過程中需要先求出協方差矩陣 X^TX，當樣本數多樣本特徵數也多的時候，這個計算量是很大的。

　　注意到SVD也可以得到協方差矩陣 X^TX 最大的 d 個特徵向量張成的矩陣，但是SVD有個好處，有一些SVD的實現演算法可以不先求出協方差矩陣 X^TX ，也能求出我們的右奇異矩陣 V。也就是說，我們的 PCA演算法可以不用做特徵分解，而是做 SVD來完成。這個方法在樣本量很大的時候很有效。實際上，scikit-learn 演算法的背後真正的實現就是SVD，而不是我們認為的暴力特徵分解。

　　另一方面，注意到 PCA僅僅使用了我們 SVD的右奇異矩陣，沒有使用左奇異矩陣，那麼左奇異矩陣有什麼用呢？

　　假設我們那的樣本是 m*n 的矩陣 X，如果我們通過 SVD找到了矩陣 XX^T 最大的 d 個特徵向量張成的 m*d 維矩陣 U，則我們如果進行如下處理：

　　可以得到一個 d*n 的矩陣 X'，這個矩陣和我們原來的 m*n 維樣本矩陣 X 相比，行數從 m 減到了 d，可見對行數進行了壓縮，也就是說，左奇異矩陣可以用於行數的壓縮。相對的，右奇異矩陣可以用於列數集特徵維度的壓縮，也就是我們的PCA降維。

　　舉個例子，我們蒐集的資料中總是存在噪聲：無論採用的裝置多精密，方法有多好，總是會存在一些誤差的。大的奇異值對應了矩陣中的主要資訊的話，運用SVD進行資料分析，提取其中的主要部分的話，還是相當合理的。

　　作為例子，假設我們蒐集的資料如下所示：

　　我們將資料用矩陣的形式表示：

　　經過奇異值分解後，得到：

　　由於第一個奇異值遠比第二個要大，資料中有包含一些噪聲，第二個奇異值在原始矩陣分解相對應的部分可以忽略。經過SVD分解後，保留了主要樣本點如圖所示：

　　就保留主要資料來看，該過程與PCA降維有一些聯絡，PCA 也是要了SVD去檢測資料間依賴和冗餘資訊。

4.2 平行奇異值

　　把頻率選擇性衰落通道進行分解。

4.3 求偽逆

　　奇異值分解可以被用來計算矩陣的偽逆，若矩陣 M 的奇異值分解為 M = UΣV*，那麼 M 的偽逆為： M+ = VΣ⁺U*。

　　其中 Σ⁺是 Σ 的偽逆，並將其主對角線上每個非零元素都求導數之後再轉置得到的。求偽逆通常用來求解線性最小平方，最小二乘法問題。

4.4 在資料表達上的應用

　　下面我們來看一個奇異值分解在資料表達上的應用，假設我們有如下一張 15*25的影象資料：

　　如圖所示，該影象主要由下面三部分構成：

　　我們將影象表示為 15*25 的矩陣，矩陣的元素對應著影象不同畫素，如果畫素是白色的話，就取1，黑色的話就取0，我們得到了一個具有 375個元素的矩陣，如下圖所示：

　　如果我們對矩陣 M 進行奇異值分解以後，得到的奇異值分別是：

　　矩陣 M 就可以表示成：

　　vi 具有 15個元素， ui 具有 25 個元素，σi 對應不同的奇異值，如上圖所示，我們就可以用 123個元素來表示具有 375個元素的影象資料了。

4.5 降噪（noise reduction）

　　前面例子的奇異值都不為零，或者都還算比較大，下面我們來探索一下擁有零或者非常小的奇異值的情況。通常來講，大的奇異值對應的部分會包含更多的資訊。比如，我們有一張掃描，帶有噪聲的影象，如下圖所示：

　　我們採用跟例項二相同的處理方式處理該掃描影象，得到影象矩陣的奇異值：

　　很明顯，前面三個奇異值遠遠比後面的奇異值要大，這樣矩陣 M 的分解方式就可以如下：

　　經過奇異值分解後，我們得到了一張降噪後的影象。

5， SVD計算例項

5.1 例項1

　　對M進行奇異值分解，M矩陣如下：

　　我們看一下 M 矩陣變換後的效果，如下圖所示：

　　在這個例子中，第二個奇異值為0，因此經過變換後只有一個方向上有表達。

　　換句話說，如果某些奇異值非常小的話，其相對應的幾項就可以不同出現在矩陣 M 的分解式中。因此，我們可以看到矩陣 M 的秩大小等於非零奇異值的個數。

5.2 例項2

　　這裡我們用一個例子來說明矩陣是如何進行奇異值分解的。

　　我們需要進行奇異值分解的矩陣A如下：

　　我們首先求出 A^TA和 AA^T：

　　進行求出 A^TA 的特徵值和特徵向量（此處的特徵向量取的是單位向量）：

　　則 V 為：

　　接著求 AA^T 的特徵值和特徵向量（此處的特徵向量取的是單位向量）：

　　則 U 為：

　　利用 Av_i=σ_iu_i，i=1,2 求奇異值：

　　當然，我們也可以用 σ_i= √λ_i 直接求出奇異值為 √3 和 1。

　　即通過：

　　我們代數矩陣即可求出奇異值。

　　最終得到 A 的奇異值分解為：

6，利用Python實現SVD降維

6.1 Python中SVD的函式

　　Python中的 Numpy 包內有一個 linalg 的線性工具，可以使用它來直接計算 SVD 的結果，不需要專門研究對應的線性代數的知識，我們拿上面 5.2 的例項來看，矩陣如下：

　　我們使用該函式來計算它的SVD。

　　得到如下結果：

　　可以看到 Σ 是以行向量的形式出現，因為該矩陣除了對角元素以外其他元素均為0，這樣的格式更節省空間。

　　我們和之前計算的結果進行比對：

　　我們發現 Σ 和我們計算出來的是一致的。

　　當然人可以計算2*2, 3*3的，但是當矩陣的維度越來越大呢？

　　接下來再看一個大的資料集，我們建立如下的資料集：

　　我們求出奇異值為：

　　我們發現它後面兩個值接近於0，所以我們就可以將最後兩個值去掉了，接下來我們就可以去掉一部分元素的矩陣近似表示原始資料：

　　這樣我們重構的近似矩陣如下：

[[ 1.00000000e+00  1.00000000e+00  1.00000000e+00  7.75989921e-16
   7.71587483e-16]
 [ 2.00000000e+00  2.00000000e+00  2.00000000e+00  3.00514919e-16
   2.77832253e-16]
 [ 1.00000000e+00  1.00000000e+00  1.00000000e+00  2.18975112e-16
   2.07633779e-16]
 [ 5.00000000e+00  5.00000000e+00  5.00000000e+00  3.00675663e-17
  -1.28697294e-17]
 [ 1.00000000e+00  1.00000000e+00 -5.48397422e-16  2.00000000e+00
   2.00000000e+00]
 [ 3.21319929e-16  4.43562065e-16 -3.48967188e-16  3.00000000e+00
   3.00000000e+00]
 [ 9.71445147e-17  1.45716772e-16 -1.52655666e-16  1.00000000e+00
   1.00000000e+00]]

　　四捨五入之後與原始資料基本保持一致。

　　我們如果知道僅僅需要保留到前3個奇異值呢？其中一個典型的做法就是保留矩陣中的 90%的能量資訊，我們將所有的奇異值取平方和，直到加到總和的90%為止。

6.2 在影象壓縮中的應用

　　一個影象矩陣，我們總可以將它分解為以下形式，通過選取不同個數的 Σ 中的奇異值，就可以實現影象的壓縮：

　　如果只想實現影象壓縮，我們可以使用python numpy 庫中的 linalg.svd 對影象矩陣進行分解，進而提取前K個奇異值便能實現SVD影象壓縮的效果，下面我們看一下程式碼：

#_*_ coding:utf-8_*_
import numpy as np
import cv2

img = cv2.imread('harden.jpg')
print('origin image shape is ', img.shape)
# 表示 RGB 中各有一個矩陣，都為300*532
#  origin image shape is  (300, 532, 3)


def svd_compression(img, k):
    res_image = np.zeros_like(img)
    for i in range(img.shape[2]):
        # 進行奇異值分解, 從svd函式中得到的奇異值sigma 是從大到小排列的
        U, Sigma, VT = np.linalg.svd(img[:,:,i])
        res_image[:, :, i] = U[:,:k].dot(np.diag(Sigma[:k])).dot(VT[:k,:])

    return res_image


# 保留前 k 個奇異值
res1 = svd_compression(img, k=300)
res2 = svd_compression(img, k=200)
res3 = svd_compression(img, k=100)
res4 = svd_compression(img, k=50)

row11 = np.hstack((res1, res2))
row22 = np.hstack((res3, res4))
res = np.vstack((row11, row22))

cv2.imshow('img', res)
cv2.waitKey(0)
cv2.destroyAllWindows()

　　我們這裡分別提取了前300， 200， 100， 50 的奇異值，圖如下：

　　可以看到，當我們取到前面300個奇異值來重構圖片時，基本與原圖看不出來差別，甚至兩百的都是都還OK。

　　所以從上面的壓縮結果來看，奇異值可以被看作是一個矩陣的代表值，或者說奇異值能夠代表這個矩陣的資訊，當奇異值越大時，它代表的資訊越多。因此我們取前面若干個最大的奇異值，就可以基本還原出資料本身。

參考文獻：https://blog.csdn.net/xiahouzuoxin/article/details/41118351

http://www.ams.org/publicoutreach/feature-column/fcarc-svd

https://www.cnblogs.com/pinard/p/6251584.html

https://jingyan.baidu.com/article/9f63fb916ba5e1c8400f0eca.html

http://blog.sciencenet.cn/blog-696950-699432.html

百度百科，SVD的步驟推導