1.前言
第一次接觸奇異值分解還是在本科期間,那個時候要用到點對點的剛體配准,這是查文獻剛好找到了四元數理論用於配准方法(點對點配准可以利用四元數方法,如果點數不一致更建議應用ICP算法)。一直想找個時間把奇異值分解理清楚、弄明白,直到今天才系統地來進行總結。
上一次學習過關於PCA的文章,PCA的實現一般有兩種,一種是用特征值分解去實現的,一種是用奇異值分解去實現的。特征值和奇異值在大部分人的印象中,往往是停留在純粹的數學計算中。而且線性代數或者矩陣論里面,也很少講任何跟特征值與奇異值有關的應用背景。奇異值分解是一個有着很明顯的物理意義的一種方法,它可以將一個比較復雜的矩陣用更小更簡單的幾個子矩陣的相乘來表示,這些小矩陣描述的是矩陣的重要的特性。就像是描述一個人一樣,給別人描述說這個人長得濃眉大眼,方臉,絡腮胡,而且帶個黑框的眼鏡,這樣寥寥的幾個特征,就讓別人腦海里面就有一個較為清楚的認識,實際上,人臉上的特征是有着無數種的,之所以能這么描述,是因為人天生就有着非常好的抽取重要特征的能力,讓機器學會抽取重要的特征,SVD是也一個重要的方法。在機器學習領域,有相當多的應用與奇異值都可以扯上關系,比如做feature reduction的PCA,做數據壓縮(以圖像壓縮為代表)的算法,還有做搜索引擎語義層次檢索的LSI(Latent Semantic Indexing)。
本文主要關注奇異值的一些特性,還會稍稍提及奇異值的計算。另外,本文里面有部分不算太深的線性代數的知識,如果完全忘記了線性代數,看文可能會有些困難。
2.奇異值分解詳解
特征值分解和奇異值分解兩者有着很緊密的關系,特征值分解和奇異值分解的目的都是一樣,就是提取出一個矩陣最重要的特征。先談談特征值分解吧:
1 特征值:
如果說一個向量v是方陣A的特征向量,將一定可以表示成下面的形式:
這時候λ就被稱為特征向量v對應的特征值,一個矩陣的一組特征向量是一組正交向量。特征值分解是將一個矩陣分解成下面的形式:
其中Q是這個矩陣A的特征向量組成的矩陣,Σ是一個對角陣,每一個對角線上的元素就是一個特征值。首先,要明確的是,一個矩陣其實就是一個線性變換,因為一個矩陣乘以一個向量后得到的向量,其實就相當於將這個向量進行了線性變換。比如說下面的一個矩陣:
它其實對應的線性變換是下面的形式:
因為這個矩陣M乘以一個向量(x,y)的結果是:
上面的矩陣是對稱的,所以這個變換是一個對x,y軸的方向一個拉伸變換(每一個對角線上的元素將會對一個維度進行拉伸變換,當值>1時拉長,當值<1時縮短),當矩陣不是對稱的時候,假如說矩陣是下面的樣子:
它所描述的變換是下面的樣子:
這其實是在平面上對一個軸進行的拉伸變換(如藍色的箭頭所示),在圖中,藍色的箭頭是一個最主要的變化方向(變化方向可能有不止一個),如果我們想要描述好一個變換,那我們就描述好這個變換主要的變化方向就好了。反過頭來看看之前特征值分解的式子,分解得到的Σ矩陣是一個對角陣,里面的特征值是由大到小排列的,這些特征值所對應的特征向量就是描述這個矩陣變化方向(從主要的變化到次要的變化排列)
當矩陣是高維的情況下,那么這個矩陣就是高維空間下的一個線性變換,這個線性變化可能沒法通過圖片來表示,但是可以想象,這個變換也同樣有很多的變換方向,我們通過特征值分解得到的前N個特征向量,那么就對應了這個矩陣最主要的N個變化方向。我們利用這前N個變化方向,就可以近似這個矩陣(變換)。也就是之前說的:提取這個矩陣最重要的特征。總結一下,特征值分解可以得到特征值與特征向量,特征值表示的是這個特征到底有多重要,而特征向量表示這個特征是什么,可以將每一個特征向量理解為一個線性的子空間,我們可以利用這些線性的子空間干很多的事情。不過,特征值分解也有很多的局限,比如說變換的矩陣必須是方陣。
2 奇異值:
下面重點談談奇異值分解。特征值分解是一個提取矩陣特征很不錯的方法,但是它只是對方陣而言的,在現實的世界中,我們看到的大部分矩陣都不是方陣,比如說有N個學生,每個學生有M科成績,這樣形成的一個N * M的矩陣就不可能是方陣,我們怎樣才能描述這樣普通的矩陣呢的重要特征呢?奇異值分解可以用來干這個事情,奇異值分解是一個能適用於任意的矩陣的一種分解的方法:
假設A是一個N * M的矩陣,那么得到的U是一個N * N的方陣(里面的向量是正交的,U里面的向量稱為左奇異向量),Σ是一個N * M的矩陣(除了對角線的元素都是0,對角線上的元素稱為奇異值),V’(V的轉置)是一個N * N的矩陣,里面的向量也是正交的,V里面的向量稱為右奇異向量),如下圖所示:
那么奇異值和特征值是怎么對應起來的呢?首先,我們將一個矩陣A的轉置 * A,將會得到一個方陣,我們用這個方陣求特征值可以得到:
這里得到的v,就是我們上面的右奇異向量。此外我們還可以得到:
這里的σ就是上面說的奇異值,u就是上面說的左奇異向量。奇異值σ跟特征值類似,在矩陣Σ中也是從大到小排列,而且σ的減少特別的快,在很多情況下,前10%甚至1%的奇異值的和就占了全部的奇異值之和的99%以上了。也就是說,我們也可以用前r大的奇異值來近似描述矩陣,這里定義一下部分奇異值分解:
r是一個遠小於m、n的數,這樣矩陣的乘法看起來像是下面的樣子:
右邊的三個矩陣相乘的結果將會是一個接近於A的矩陣,在這兒,r越接近於n,則相乘的結果越接近於A。而這三個矩陣的面積之和(在存儲觀點來說,矩陣面積越小,存儲量就越小)要遠遠小於原始的矩陣A,我們如果想要壓縮空間來表示原矩陣A,我們存下這里的三個矩陣:U、Σ、V就好了。
3.如何計算奇異值
奇異值的計算是一個難題,是一個O(N^3)的算法。在單機的情況下當然是沒問題的,matlab在一秒鍾內就可以算出1000 * 1000的矩陣的所有奇異值,但是當矩陣的規模增長的時候,計算的復雜度呈3次方增長,就需要並行計算參與了。
其實SVD還是可以用並行的方式去實現的,在解大規模的矩陣的時候,一般使用迭代的方法,當矩陣的規模很大(比如說上億)的時候,迭代的次數也可能會上億次,如果使用Map-Reduce框架去解,則每次Map-Reduce完成的時候,都會涉及到寫文件、讀文件的操作。個人猜測Google雲計算體系中除了Map-Reduce以外應該還有類似於MPI的計算模型,也就是節點之間是保持通信,數據是常駐在內存中的,這種計算模型比Map-Reduce在解決迭代次數非常多的時候,要快了很多倍。
Lanczos迭代就是一種解對稱方陣部分特征值的方法(之前談到了,解A’* A得到的對稱方陣的特征值就是解A的右奇異向量),是將一個對稱的方程化為一個三對角矩陣再進行求解。
由於奇異值的計算是一個很枯燥,純數學的過程,而且前人的研究成果(論文中)幾乎已經把整個程序的流程圖給出來了。更多的關於奇異值計算的部分,將在后面的參考文獻中給出,這里不再深入,我還是focus在奇異值的應用中去。
4.奇異值分解應用
奇異值與主成分分析(PCA)
這里主要談談如何用SVD去解PCA的問題。PCA的問題其實是一個基的變換,使得變換后的數據有着最大的方差。方差的大小描述的是一個變量的信息量,我們在講一個東西的穩定性的時候,往往說要減小方差,如果一個模型的方差很大,那就說明模型不穩定了。但是對於我們用於機器學習的數據(主要是訓練數據),方差大才有意義,不然輸入的數據都是同一個點,那方差就為0了,這樣輸入的多個數據就等同於一個數據了。以下面這張圖為例子:
這個假設是一個攝像機采集一個物體運動得到的圖片,上面的點表示物體運動的位置,假如我們想要用一條直線去擬合這些點,那我們會選擇什么方向的線呢?當然是圖上標有signal的那條線。如果我們把這些點單純的投影到x軸或者y軸上,最后在x軸與y軸上得到的方差是相似的(因為這些點的趨勢是在45度左右的方向,所以投影到x軸或者y軸上都是類似的),如果我們使用原來的xy坐標系去看這些點,容易看不出來這些點真正的方向是什么。但是如果我們進行坐標系的變化,橫軸變成了signal的方向,縱軸變成了noise的方向,則就很容易發現什么方向的方差大,什么方向的方差小了。
一般來說,方差大的方向是信號的方向,方差小的方向是噪聲的方向,我們在數據挖掘中或者數字信號處理中,往往要提高信號與噪聲的比例,也就是信噪比。對上圖來說,如果我們只保留signal方向的數據,也可以對原數據進行不錯的近似了。
PCA的全部工作簡單點說,就是對原始的空間中順序地找一組相互正交的坐標軸,第一個軸是使得方差最大的,第二個軸是在與第一個軸正交的平面中使得方差最大的,第三個軸是在與第1、2個軸正交的平面中方差最大的,這樣假設在N維空間中,我們可以找到N個這樣的坐標軸,我們取前r個去近似這個空間,這樣就從一個N維的空間壓縮到r維的空間了,但是我們選擇的r個坐標軸能夠使得空間的壓縮使得數據的損失最小。
還是假設我們矩陣每一行表示一個樣本,每一列表示一個feature,用矩陣的語言來表示,將一個m * n的矩陣A的進行坐標軸的變化,P就是一個變換的矩陣從一個N維的空間變換到另一個N維的空間,在空間中就會進行一些類似於旋轉、拉伸的變化。
而將一個m * n的矩陣A變換成一個m * r的矩陣,這樣就會使得本來有n個feature的,變成了有r個feature了(r < n),這r個其實就是對n個feature的一種提煉,我們就把這個稱為feature的壓縮。用數學語言表示就是:
但是這個怎么和SVD扯上關系呢?之前談到,SVD得出的奇異向量也是從奇異值由大到小排列的,按PCA的觀點來看,就是方差最大的坐標軸就是第一個奇異向量,方差次大的坐標軸就是第二個奇異向量…我們回憶一下之前得到的SVD式子:
在矩陣的兩邊同時乘上一個矩陣V,由於V是一個正交的矩陣,所以V轉置乘以V得到單位陣I,所以可以化成后面的式子:
將后面的式子與A * P那個m * n的矩陣變換為m * r的矩陣的式子對照看看,在這里,其實V就是P,也就是一個變化的向量。這里是將一個m * n 的矩陣壓縮到一個m * r的矩陣,也就是對列進行壓縮,如果我們想對行進行壓縮(在PCA的觀點下,對行進行壓縮可以理解為,將一些相似的sample合並在一起,或者將一些沒有太大價值的sample去掉)怎么辦呢?同樣我們寫出一個通用的行壓縮例子:
這樣就從一個m行的矩陣壓縮到一個r行的矩陣了,對SVD來說也是一樣的,我們對SVD分解的式子兩邊乘以U的轉置U':
這樣我們就得到了對行進行壓縮的式子。可以看出,其實PCA幾乎可以說是對SVD的一個包裝,如果我們實現了SVD,那也就實現了PCA了,而且更好的地方是,有了SVD,我們就可以得到兩個方向的PCA,如果我們對A’A進行特征值的分解,只能得到一個方向的PCA。