一步步教你輕鬆學奇異值分解SVD降維演算法

阿新 • • 發佈：2018-12-17

摘要：奇異值分解（singular value decomposition）是線性代數中一種重要的矩陣分解，在生物資訊學、訊號處理、金融學、統計學等領域有重要應用，SVD都是提取資訊的強度工具。在機器學習領域，很多應用與奇異值都有關係，比如推薦系統、資料壓縮（以影象壓縮為代表）、搜尋引擎語義層次檢索的LSI等等。（本文原創，轉載必須註明出處.）

奇異值分解原理

什麼是奇異值分解(SVD）

奇異值分解

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

$$ M_{m×n}=U_{m×m} \Sigma_{m×n} V^T_{n×n} $$

其中U是m×m階酉矩陣；Σ是m×n階非負實數對角矩陣；而$V^T$，即V的共軛轉置，是n×n階酉矩陣。這樣的分解就稱作M的奇異值分解。Σ對角線上的元素$Σ_i$,i即為M的奇異值。常見的做法是將奇異值由大而小排列。如此Σ便能由M唯一確定了。（雖然U和V仍然不能確定。）

V的列組成一套對$M$的正交"輸入"或"分析"的基向量。這些向量是$M^*M$的特徵向量。
U的列組成一套對$M$的正交"輸出"的基向量。這些向量是$MM^*$的特徵向量。
Σ對角線上的元素是奇異值，可視為是在輸入與輸出間進行的標量的"膨脹控制"。這些是$ MM^* $及 $ M^* M $的特徵值的非負平方根，並與U和V的行向量相對應。

SVD 的計算方法

SVD 與特徵值

現在，假設矩陣 $ \mathbf M_{m\times n} $ 的 SVD 分解是 $ \mathbf M = \mathbf U\mathbf\Sigma\mathbf V^{\mathsf H}; $那麼，我們有

$$
\begin{aligned}
\mathbf M\mathbf M^{\mathsf H} &{}= \mathbf U\mathbf\Sigma\mathbf V^{\mathsf H}\mathbf V\mathbf\Sigma^{\mathsf H}\mathbf U^{\mathsf H} = \mathbf U(\mathbf\Sigma\mathbf\Sigma^{\mathsf H})\mathbf U^{\mathsf H} \\
\mathbf M^{\mathsf H}\mathbf M &{}= \mathbf V\mathbf\Sigma^{\mathsf H}\mathbf U^{\mathsf H}\mathbf U\mathbf\Sigma\mathbf V^{\mathsf H} = \mathbf V(\mathbf\Sigma^{\mathsf H}\mathbf\Sigma)\mathbf V^{\mathsf H}\
\end{aligned}
$$

這也就是說，$ \mathbf U $ 的列向量（左奇異向量），是$ \mathbf M\mathbf M^{\mathsf H} $ 的特徵向量；同時，$ \mathbf V $ 的列向量（右奇異向量），是 $ \mathbf M^{\mathsf H}\mathbf M $ 的特徵向量；另一方面，$ \mathbf M $ 的奇異值（$ \mathbf\Sigma $ 的非零對角元素）則是 $ \mathbf M\mathbf M^{\mathsf H} $ 或者 $ \mathbf M^{\mathsf H}\mathbf M $ 的非零特徵值的平方根。

如何計算 SVD

有了這些知識，我們就能手工計算出任意矩陣的 SVD 分解了；具體來說，演算法如下

計算 $ \mathbf M\mathbf M^{\mathsf H} $ 和 $ \mathbf M^{\mathsf H}\mathbf M $；
分別計算 $ \mathbf M\mathbf M^{\mathsf H} $ 和 $ \mathbf M^{\mathsf H}\mathbf M $ 的特徵向量及其特徵值；
$ \mathbf M\mathbf M^{\mathsf H} $ 的特徵向量組成 $ \mathbf U $；而 $ \mathbf M^{\mathsf H}\mathbf M $ 的特徵向量組成 $ \mathbf V $；
對 $ \mathbf M\mathbf M^{\mathsf H} $ 和 $ \mathbf M^{\mathsf H}\mathbf M $ 的非零特徵值求平方根，對應上述特徵向量的位置，填入 $ \mathbf\Sigma $ 的對角元。

實際計算看看

現在，我們來試著計算 $ \mathbf M = \begin{bmatrix}2 & 4 \\ 1 & 3 \\ 0 & 0 \\ 0 & 0\end{bmatrix} $ 的奇異值分解。計算奇異值分解，需要計算 $ \mathbf M $ 與其共軛轉置的左右積；這裡主要以 $ \mathbf M\mathbf M^{\mathsf H} $ 為例。
首先，我們需要計算 $ \mathbf M\mathbf M^{\mathsf H} $，

$$
\mathbf W = \mathbf M\mathbf M^{\mathsf H} = \begin{bmatrix}2 & 4 \\ 1 & 3 \ 0 & 0 \\ 0 & 0\end{bmatrix}\begin{bmatrix}2 & 1 & 0 & 0 \\ 4 & 3 & 0 & 0\end{bmatrix} = \begin{bmatrix}20 & 14 & 0 & 0 \\ 14 & 10 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0\end{bmatrix}.
$$

現在，我們要求 $ \mathbf W $ 的特徵值與特徵向量。根據定義$ \mathbf W\vec x = \lambda \vec x $；因此 $ (\mathbf W - \lambda\mathbf I)\vec x = \vec 0 $。這也就是說

$$
\begin{bmatrix}
20 - \lambda & 14 & 0 & 0 \\
14 & 10 - \lambda & 0 & 0 \\
0 & 0 & -\lambda & 0 \\
0 & 0 & 0 & -\lambda
\end{bmatrix}\vec x = \vec 0.
$$

根據線性方程組的理論，若要該關於$ \vec x $ 的方程有非零解，則要求係數矩陣的行列式為 0；也就是

$$
\begin{vmatrix}
20 - \lambda & 14 & 0 & 0 \\
14 & 10 - \lambda & 0 & 0 \\
0 & 0 & -\lambda & 0 \\
0 & 0 & 0 & -\lambda
\end{vmatrix} =
\begin{vmatrix}
20 - \lambda & 14 \\
14 & 10 - \lambda \\
\end{vmatrix}\begin{vmatrix}
-\lambda & 0 \\
0 & -\lambda \\
\end{vmatrix}
= 0,
$$

這也就是 $ \bigl((20 - \lambda)(10 - \lambda) - 196\bigr)\lambda^2 = 0$；解得 $\lambda_1 = \lambda_2 = 0 $, $ \lambda_{3} = 15 + \sqrt{221} \approx 29.866 $, $ \lambda_{4} = 15 - \sqrt{221} \approx 0.134 $。將特徵值代入原方程，可解得對應的特徵向量；這些特徵向量即作為列向量，形成矩陣

$$\mathbf U = \begin{bmatrix}-0.82 & -0.58 & 0 & 0 \\ -0.58 & 0.82 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\end{bmatrix}.$$

同理可解得（注意，$ \mathbf M\mathbf M^{\mathsf T} $ 和 $ \mathbf M^{\mathsf T}\mathbf M $ 的特徵值相同）

$$\mathbf V = \begin{bmatrix}-0.40 & -0.91 \\ -0.91 & 0.40\end{bmatrix}.$$

以及 $ \mathbf\Sigma $ 上的對角線元素由 $ \mathbf W $ 的特徵值的算術平方根組成；因此有

$$\mathbf\Sigma = \begin{bmatrix}5.46 & 0 \\ 0 & 0.37 \\ 0 & 0 \\ 0 & 0\end{bmatrix}.$$

因此我們得到矩陣 $ \mathbf M $ 的 SVD 分解（數值上做了近似）：

$$\begin{bmatrix}2 & 4 \\ 1 & 3 \\ 0 & 0 \\ 0 & 0\end{bmatrix} \approx \begin{bmatrix}-0.82 & -0.58 & 0 & 0 \\ -0.58 & 0.82 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\end{bmatrix}\begin{bmatrix}5.46 & 0 \\ 0 & 0.37 \\ 0 & 0 \\ 0 & 0\end{bmatrix}\begin{bmatrix}-0.40 & -0.91 \\ -0.91 & 0.40\end{bmatrix}$$

幾何上的直觀解釋

我們先來看一個例子。假設 $ \mathbf M $ 是一個 $ m\times n $ 的矩陣，而 $ \mathbf x $ 是線性空間 $ \mathbb K^n $ 中的向量，則 $ \mathbf y = \mathbf M\cdot\mathbf x $ 是線性空間 $ \mathbb K^m $ 中的向量。這樣一來，矩陣 $ \mathbb A $ 就對應了一個從$ \mathbb K^n $ 到 $ \mathbb K^m $ 的變換$ T: \mathbb K^n \to \mathbb K^m $，具體來說既是 $ \mathbf x\mapsto \mathbf M\cdot\mathbf x $。這也就是說，線上性代數中，任意矩陣都能看做是一種變換。這樣一來，我們就統一了矩陣和變換。

SVD 場景

隱性語義檢索

資訊檢索-隱性語義檢索（Lstent Semantic Indexing, LSI）或隱形語義分析（Latent Semantic Analysis, LSA）
隱性語義索引：矩陣 = 文件 + 詞語
最早的 SVD 應用之一，我們稱利用 SVD 的方法為隱性語義索引（LSI）或隱性語義分析（LSA）。

推薦系統

利用 SVD 從資料中構建一個主題空間。
再在該空間下計算其相似度。(從高維-低維空間的轉化，在低維空間來計算相似度，SVD 提升了推薦系統的效率。)

影象壓縮

例如：32*32=1024 => 32*2+2*1+32*2=130(2*1表示去掉了除對角線的0), 幾乎獲得了10倍的壓縮比。

SVD 工作原理

矩陣分解

矩陣分解是將資料矩陣分解為多個獨立部分的過程。
矩陣分解可以將原始矩陣表示成新的易於處理的形式，這種新形式是兩個或多個矩陣的乘積。（類似代數中的因數分解）
舉例：如何將12分解成兩個數的乘積？（1，12）、（2，6）、（3，4）都是合理的答案。

SVD 是矩陣分解的一種型別，也是矩陣分解最常見的技術

SVD 將原始的資料集矩陣 Data 分解成三個矩陣 U、∑、V
舉例：如果原始矩陣 $Data_{m*n} $ 是m行n列，
- $U_{m * k}$ 表示m行k列
- $∑_{k * k}$ 表示k行k列
- $V_{k * n}$ 表示k行n列。

$$ Data_{m×n} = U_{m×k} * ∑_{k×k} * V_{k×n} $$

具體的案例：

$$
\begin{vmatrix}
0 & -1.6 & 0.6 \\
0 & 1.2 & 0.8 \\
0 & 0 & 0 \\
0 & 0 & 0 \\
\end{vmatrix} =
\begin{vmatrix}
0.8 & 0.6 & 0 & 0 \\
-0.6 & 0.8 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{vmatrix} *
\begin{vmatrix}
2 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 0 \\
\end{vmatrix} *
\begin{vmatrix}
0 & 0 & 1 \\
-1 & 0 & 0 \\
0 & 1 & 0 \\
\end{vmatrix}
$$

上述分解中會構建出一個矩陣∑，該矩陣只有對角元素，其他元素均為0(近似於0)。另一個慣例就是，∑的對角元素是從大到小排列的。這些對角元素稱為奇異值。
奇異值與特徵值(PCA 資料中重要特徵)是有關係的。這裡的奇異值就是矩陣 $Data * Data^T$ 特徵值的平方根。
普遍的事實：在某個奇異值的數目(r 個=>奇異值的平方和累加到總值的90%以上)之後，其他的奇異值都置為0(近似於0)。這意味著資料集中僅有 r 個重要特徵，而其餘特徵則都是噪聲或冗餘特徵。

SVD 演算法特點

優點：簡化資料，去除噪聲，優化演算法的結果
缺點：資料的轉換可能難以理解
使用的資料型別：數值型資料

基於SVD的程式碼實現

'''分析 Sigma 的長度取值
根據自己的業務情況，就行處理，設定對應的 Singma 次數
通常保留矩陣 80% ～ 90% 的能量，就可以得到重要的特徵並取出噪聲。
'''
def analyse_data(Sigma, loopNum=20):
    # 總方差的集合（總能量值）
    Sig2 = Sigma**2
    SigmaSum = sum(Sig2)
    for i in range(loopNum):
        SigmaI = sum(Sig2[:i+1])
        print('主成分：%s, 方差佔比：%s%%' % (format(i+1, '2.0f'), format(SigmaI/SigmaSum*100, '.2f')))


'''基於SVD的評分估計
Args:
    dataMat         訓練資料集
    user            使用者編號
    simMeas         相似度計算方法
    item            未評分的物品編號
Returns:
    ratSimTotal/simTotal     評分（0～5之間的值）
'''
def svdEst(dataMat, user, simMeas, item):
    # 物品數目
    n = shape(dataMat)[1]
    # 對資料集進行SVD分解
    simTotal = 0.0 ;  ratSimTotal = 0.0
    # 奇異值分解,只利用90%能量值的奇異值，奇異值以NumPy陣列形式儲存
    U, Sigma, VT = la.svd(dataMat)
    # 分析 Sigma 的長度取值
    # analyse_data(Sigma, 20)

    # 如果要進行矩陣運算，就必須要用這些奇異值構建出一個對角矩陣
    Sig4 = mat(eye(4) * Sigma[: 4]) # eye對角矩陣
    # 利用U矩陣將物品轉換到低維空間中，構建轉換後的物品(物品+4個主要的特徵)
    xformedItems = dataMat.T * U[:, :4] * Sig4.I # I 逆矩陣
    # print('dataMat', shape(dataMat))
    # print('U[:, :4]', shape(U[:, :4]))
    # print('Sig4.I', shape(Sig4.I))
    # print('VT[:4, :]', shape(VT[:4, :]))
    # print('xformedItems', shape(xformedItems))

    # 對於給定的使用者，for迴圈在使用者對應行的元素上進行遍歷
    # 和standEst()函式的for迴圈一樣，這裡相似度計算在低維空間下進行的。
    for j in range(n):
        userRating = dataMat[user, j]
        if userRating == 0 or j == item:
            continue
        # 相似度的計算方法也會作為一個引數傳遞給該函式
        similarity = simMeas(xformedItems[item, :].T, xformedItems[j, :].T)
        # for 迴圈中加入了一條print語句，以便了解相似度計算的進展情況。如果覺得累贅，可以去掉
        # print('the %d and %d similarity is: %f' % (item, j, similarity))
        # 對相似度不斷累加求和
        simTotal += similarity
        # 對相似度及對應評分值的乘積求和
        ratSimTotal += similarity * userRating
    if simTotal == 0:
        return 0
    else:
        # 計算估計評分
        return ratSimTotal/simTotal

排序獲取最後的推薦結果

'''recommend函式推薦引擎，預設呼叫standEst函式，產生最高的N個推薦結果
Args:
    dataMat         訓練資料集
    user            使用者編號
    simMeas         相似度計算方法
    estMethod       使用的推薦演算法
Returns:  返回最終 N 個推薦結果
'''
def recommend(dataMat, user, N=3, simMeas=cosSim, estMethod=standEst):
    # 尋找未評級的物品,對給定的使用者建立一個未評分的物品列表
    unratedItems = nonzero(dataMat[user, :].A == 0)[1] # .A: 矩陣轉陣列
    # 如果不存在未評分物品，那麼就退出函式
    if len(unratedItems) == 0:
        return 'you rated everything'
    # 物品的編號和評分值
    itemScores = []
    # 在未評分物品上進行迴圈
    for item in unratedItems:
        # 獲取 item 該物品的評分
        estimatedScore = estMethod(dataMat, user, simMeas, item)
        itemScores.append((item, estimatedScore))
    # 按照評分得分 進行逆排序，獲取前N個未評級物品進行推薦
    return sorted(itemScores, key=lambda jj: jj[1], reverse=True)[: N]

測試和專案呼叫

測試程式碼

# 計算相似度的方法
myMat = mat(loadExData3())
# 計算相似度的第一種方式
# print(recommend(myMat, 1, estMethod=svdEst))
# 計算相似度的第二種方式
# print(recommend(myMat, 1, estMethod=svdEst, simMeas=pearsSim))

# 預設推薦（菜館菜餚推薦示例）
print(recommend(myMat, 2))

執行結果

菜館菜餚推薦結果： [(3, 4.0), (5, 4.0), (6, 4.0)]

***Repl Closed***

分析結果，我們不難發現，分別對3烤牛肉，5魯賓三明治、6印度烤雞給我4星好評，推薦給我們的使用者。

要點補充

基於內容(content-based)的推薦

通過各種標籤來標記菜餚
將這些屬性作為相似度計算所需要的資料
這就是：基於內容的推薦。

構建推薦引擎面臨的挑戰

問題

1）在大規模的資料集上，SVD分解會降低程式的速度
2）存在其他很多規模擴充套件性的挑戰性問題，比如矩陣的表示方法和計算相似度得分消耗資源。
3）如何在缺乏資料時給出好的推薦-稱為冷啟動【簡單說：使用者不會喜歡一個無效的物品，而使用者不喜歡的物品又無效】

建議

1）在大型系統中，SVD分解(可以在程式調入時執行一次)每天執行一次或者其頻率更低，並且還要離線執行。
2）在實際中，另一個普遍的做法就是離線計算並儲存相似度得分。(物品相似度可能被使用者重複的呼叫)
3）冷啟動問題，解決方案就是將推薦看成是搜尋問題，通過各種標籤／屬性特徵進行基於內容的推薦。

專案案例: 基於SVD的影象壓縮

收集並準備資料

將文字資料轉化為矩陣

'''影象壓縮函式'''
def imgLoadData(filename):
    myl = []
    for line in open(filename).readlines():
        newRow = []
        for i in range(32):
            newRow.append(int(line[i]))
        myl.append(newRow)
    # 矩陣調入後，就可以在螢幕上輸出該矩陣
    myMat = mat(myl)
    return myMat

分析資料: 分析Sigma的長度個數

通常保留矩陣 80% ～ 90% 的能量，就可以得到重要的特徵並去除噪聲。

'''分析 Sigma 的長度取值
根據自己的業務情況，就行處理，設定對應的 Singma 次數
通常保留矩陣 80% ～ 90% 的能量，就可以得到重要的特徵並取出噪聲。
'''
def analyse_data(Sigma, loopNum=20):
    # 總方差的集合（總能量值）
    Sig2 = Sigma**2
    SigmaSum = sum(Sig2)
    for i in range(loopNum):
        SigmaI = sum(Sig2[:i+1])
        print('主成分：%s, 方差佔比：%s%%' % (format(i+1, '2.0f'), format(SigmaI/SigmaSum*100, '.2f')))

使用演算法: 對比使用 SVD 前後的資料差異對比，對於儲存大家可以試著寫寫

例如：32*32=1024 => 32*2+2*1+32*2=130(2*1表示去掉了除對角線的0), 幾乎獲得了10倍的壓縮比。

'''列印矩陣
由於矩陣保護了浮點數，因此定義淺色和深色，遍歷所有矩陣元素，當元素大於閥值時列印1，否則列印0
'''
def printMat(inMat, thresh=0.8):
    for i in range(32):
        for k in range(32):
            if float(inMat[i, k]) > thresh:
                print(1)
            else:
                print(0)
        print('')


'''實現影象壓縮，允許基於任意給定的奇異值數目來重構影象
Args:
    numSV       Sigma長度
    thresh      判斷的閾值
'''
def imgCompress(numSV=3, thresh=0.8):
    # 構建一個列表
    myMat = imgLoadData('./0_5.txt')

    print("****original matrix****")
    # 對原始影象進行SVD分解並重構影象e
    printMat(myMat, thresh)

    # 通過Sigma 重新構成SigRecom來實現
    # Sigma是一個對角矩陣，因此需要建立一個全0矩陣，然後將前面的那些奇異值填充到對角線上。
    U, Sigma, VT = la.svd(myMat)
    # SigRecon = mat(zeros((numSV, numSV)))
    # for k in range(numSV):
    #     SigRecon[k, k] = Sigma[k]

    # 分析插入的 Sigma 長度
    # analyse_data(Sigma, 20)

    SigRecon = mat(eye(numSV) * Sigma[: numSV])
    reconMat = U[:, :numSV] * SigRecon * VT[:numSV, :]
    print("****reconstructed matrix using %d singular values *****" % numSV)
    printMat(reconMat, thresh)

參考文獻

完整程式碼下載

原始碼請進【機器學習和自然語言QQ群：436303759】檔案下載：

一步步教你輕鬆學奇異值分解SVD降維演算法

目錄

奇異值分解原理

什麼是奇異值分解(SVD）

SVD 的計算方法

SVD 與特徵值

如何計算 SVD

實際計算看看

幾何上的直觀解釋

SVD 場景

SVD 工作原理

矩陣分解

SVD 是矩陣分解的一種型別，也是矩陣分解最常見的技術

SVD 演算法特點

推薦系統

推薦系統場景

推薦系統要點

相似度計算

推薦系統的評價

推薦系統原理

專案實戰: 餐館菜餚推薦系統

收集並準備資料

分析資料

訓練演算法: 通過呼叫 recommend() 函式進行推薦

排序獲取最後的推薦結果

測試和專案呼叫

要點補充

專案案例: 基於SVD的影象壓縮

收集並準備資料

分析資料: 分析Sigma的長度個數

使用演算法: 對比使用 SVD 前後的資料差異對比，對於儲存大家可以試著寫寫

參考文獻

完整程式碼下載

相關推薦