前面我們講了 QR 分解有一些優良的特性，但是 QR 分解僅僅是對矩陣的行進行操作（左乘一個酉矩陣），可以得到列空間。這一小節的 SVD 分解則是將行與列同等看待，既左乘酉矩陣，又右乘酉矩陣，可以得出更有意思的信息。奇異值分解( SVD, Singular Value Decomposition ) 在計算矩陣的偽逆( pseudoinverse )，最小二乘法最優解，矩陣近似，確定矩陣的列向量空間，秩以及線性系統的解集空間都有應用。
1. SVD 的形式
對於一個任意的 m×n 的矩陣 A，SVD 將一個矩陣分解為三個特殊矩陣的乘積, 國外還做了一個視頻——SVD 之歌：
其中， U 和
是酉矩陣，
是對角線矩陣。注意酉矩陣只是坐標的轉換，數據本身分布的形狀並沒有改變，而對角矩陣，則是對數據進行了拉伸或者壓縮。由於 m >= n，又可以寫成下面這種 thin SVD 形式：
2. SVD 的幾何解釋
考慮 A 是 2×2 的簡單情況。我們知道，一個幾何形狀左乘一個矩陣 A 實際上就是將該形狀進行旋轉、對稱、拉伸變換等錯切變換, A 就是所謂的 shear 矩陣。比如可以將一個圓通過左乘 A 得到一個旋轉后的橢圓。
舉個例子，假設 A 對平面上的一個圓進行變換:
如果只看矩陣 A， 我們幾乎無法直觀看到圓是如何變換的。變換前的圓是：
C, M, Y, K 分別表示第一、二、三、四象限。左乘 A 矩陣后得到：
我們的問題是，旋轉了多少度？伸縮的方向是多少？最大伸縮比例是多少？
利用這個在線工具對 A 進行 SVD 分解:
明顯，我們看到 A 變換實際上是先對圓順時針旋轉 45°(可以看做是坐標軸逆時針宣戰了 45°，主成分方向)，再關於 x 軸對稱(第一行乘以 -1）, 即左乘 V^T：
然后在 x 方向拉伸 3 倍(左乘S)：
最后再順時針旋轉 45°，再關於 x 軸對稱(第一行乘以 -1, 兩次對稱操作抵消了）, 即左乘 U：
上述繪圖過程的 python 代碼如下：
1 # -*- coding=utf-8 2 #!/usr/bin/python2.7 3 4 from pylab import * 5 6 7 def plotCircle(before, M=matrix([[1, 0], [0, 1]])): 8 '''before: 變換前的矩陣 9 M： 變換矩陣,默認為單位矩陣 10 返回變換之后的矩陣 11 ''' 12 eclMat = M * before # M 變換 13 eclX = array(eclMat[0]).reshape(-1) 14 eclY = array(eclMat[1]).reshape(-1) 15 axis('equal') 16 axis([-3, 3, -3, 3]) 17 grid(True) 18 plot(eclX[:25], eclY[:25], 'c', linewidth=3) 19 plot(eclX[25:50], eclY[25:50], 'm', linewidth=3) 20 plot(eclX[50:75], eclY[50:75], 'y', linewidth=3) 21 plot(eclX[75:100], eclY[75:100], 'k', linewidth=3) 22 show() 23 return eclMat 24 25 26 ang = linspace(0, 2*pi, 100) 27 28 x = cos(ang) 29 y = sin(ang) 30 cirMat = matrix([x, y]) # 2 × 100 的圓圈矩陣 31 32 # 畫最開始的圖形——圓 33 plotCircle(cirMat) 34 35 # 畫變換之后的橢圓 36 M = matrix([[2, 1], [1, 2]]) # 2 × 2 的變換矩陣 37 clMat = plotCircle(cirMat, M) 38 39 # 將 M 矩陣進行 svd 分解 40 U, s, V = np.linalg.svd(M, full_matrices=True) 41 S = np.diag(s) 42 43 # SVD 矩陣對圓的變換 44 plotCircle(cirMat) 45 Tran1 = plotCircle(cirMat, V) 46 Tran2 = plotCircle(Tran1, S) 47 Tran3 = plotCircle(Tran2, U)
3. SVD 向量空間
假設矩陣 A 的秩是 r， 那么對角線矩陣的秩也是 r (乘以酉矩陣不會改變矩陣的秩）， 我們假設：
-
那么， Ax = 0 的 null space 也就是解空間是什么呢？答案如下：
證明非常簡單，直接帶入 SVD 分解三個矩陣的右邊，計算的時候正交的向量相乘都等於0, 不為 0 的恰好都被對角線矩陣的 0 元素歸零，等號成立。如果 A 是列滿秩的，顯然符合條件的只剩下零向量了。同樣的，你可以知道 A' 的 null space。 -
矩陣 A 的線性子空間是啥？設想有一個很高的矩陣 A， 它的列向量很有可能不是正交的。對於任意一個坐標 x，矩陣 A 的線性子空間可以定義為：
R(A) = {y | y = Ax， x 是任意的坐標}
那么， u1, u2,..., ur 是 R(A) 的正交基。
這個想法也非常直觀，無論來了一個什么向量（或者叫坐標）， 經過 V^T 和 對角線矩陣變換之后還是一個坐標，這個坐標就是 U 矩陣列向量線性組合的系數。
4. SVD 的計算
這就是正定矩陣的對角化。計算過程如下：
- 計算 A' (A 的轉置)和 A'A
- 計算 A'A 的特征值，將特征值按照遞減的順序排列, 求均方根，得到 A 的奇異值
- 由奇異值可以構建出對角線矩陣 S， 同時求出 S 逆，以備后面的計算
- 有上述排好序的特征值可以求出對應的特征向量，以特征向量為列得到矩陣 V， 轉置后得到 V'
- U = AVS逆，求出 U，完畢。
Lanczos algorithm 是一種迭代的計算方法，沒來得及細看。
感覺 SVD 計算水很深，要用到的時候再看，現在暫不深入了。
5. PCA
主成分分析經常用於減少數據集的維數，同時保持數據集中的對方差貢獻最大的特征。假設上面的橢圓中（二維空間，兩個坐標值）均勻分布了非常多的點，如何在一維空間（一個坐標值）里面就能最大程度將這些點區分開來呢？這時候就要用到 PCA：
可以看到，上圖中將所有的點投影到 +45° 直線上，將二維空間映射到一維空間，可以最大程度地區分開這些點，即投影后的樣本分布方差最大，這個方向就是 v1 向量的方向。就上圖來說，假設 100×2 的矩陣 X 表示樣本點在平面上的坐標，將這些點投影到 v1 上保留最大的樣本方差：
這樣，就將 100 個 2 維空間的樣本點壓縮為 100 個 1 維空間的樣本點，這里的列壓縮實際上是對特征進行了壓縮, 而不是簡單地丟棄。PCA 是不是只能在樣本點個數大於特征個數的時候才能壓縮呢？例如，現在 3 個 100 維特征的樣本用 100×3 的矩陣表示，用 SVD 分解后可以得到三個矩陣的乘積，做如下變換：
同樣也最大限度保留了主成分。
總結起來，一個矩陣 A， 如果想對行進行壓縮並保留主成分，那么左乘 u1'， 如果相對列進行壓縮並保留主成分(讓我聯想起了稀疏表示)，那么右乘 v1。
當然，上面是簡單的保留第一個主成分，PCA的全部工作簡單點說，就是對原始的空間中順序地找一組相互正交的坐標軸，第一個軸是使得方差最大的，第二個軸是在與第一個軸正交的平面中使得方差最大的，第三個軸是在與第1、2個軸正交的平面中方差最大的，這樣假設在N維空間中，我們可以找到N個這樣的坐標軸，我們取前r個去近似這個空間，這樣就從一個N維的空間壓縮到r維的空間了，但是我們選擇的r個坐標軸能夠使得空間的壓縮使得數據的損失最小。
6. 最小二乘問題
想法和 QR 分解的辦法類似，主要是利用酉矩陣變換的長度不變性：
得到最優解是：
