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 個特徵值  λ1 <= λ2 <=  λn  ,以及這 n 個特徵值所對應的特徵向量  {W1,  W2, ...Wn},如果這 n 個特徵向量線性無關,那麼矩陣 A 就可以用下式的特徵分解表示:

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

  一般我們會把 W 的這 n 個特徵向量標準化,即滿足 ||wi||2 = 1,或者說 wiTwi=1,此時 W 的 n 個特向量為標準正交基,滿足WTW = I,即 WT = W-1

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

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

2,特徵分解和奇異值分解(SVD)的區別

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

   所有的矩陣都可以進行奇異值分解,但只有方陣才可以進行特徵值分解。當所給的矩陣是對稱的方陣,AT=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 的一個方陣 ATA。既然 ATA 是方陣,那麼我們就可以進行特徵分解,得到的特徵值和特徵向量滿足下式:

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

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

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

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

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

  我們注意到:

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

  上面還有一個問題就是我們說 ATA 的特徵向量組成的就是我們SVD中的 V矩陣,而 AAT的特徵向量組成的就是我們 SVD中的 U矩陣,這有什麼根據嗎?這個其實很容易證明,我們以V矩陣的證明為例:

  上式證明使用了:UTU = I, ΣTΣ = Σ2。可以看出 ATA 的特徵向量組成的的確就是我們SVD中的 V 矩陣。類似的方法可以得到 AAT 的特徵向量組成的就是我們 SVD中的 U 矩陣。

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

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

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

注意1:通過ATA的特徵值取平方根求解奇異值的注意點

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

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

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

注意2:酋矩陣的定義

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

  其中 In 為 n 階單位矩陣,UT 為 U的共軛轉置,則U為酋矩陣(即矩陣U為酋矩陣,當且僅當其共軛轉置UT 為其逆矩陣:U-1 = UT

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降維中,需要找到樣本協方差矩陣 XTX 的最大的 d 個特徵向量,然後用這最大的 d 個特徵向量張成的矩陣來做低維投影降維。可以看出,在這個過程中需要先求出協方差矩陣 XTX,當樣本數多樣本特徵數也多的時候,這個計算量是很大的。

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

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

  假設我們那的樣本是 m*n 的矩陣 X,如果我們通過 SVD找到了 矩陣 XXT 最大的 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如下:

  我們首先求出 ATA和 AAT

  進行求出 ATA 的特徵值和特徵向量(此處的特徵向量取的是單位向量):

  則 V 為:

  接著求 AAT 的特徵值和特徵向量(此處的特徵向量取的是單位向量):

  則 U 為:

  利用 Aviiui,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的步驟推導等等

相關文章