PCA演算法(主成分分析)
寫在前面
Principle Component Analysis 顧名思義,是通過分析向量空間的主成分,將主成分提取出來,不重要的成分略去,從而達到降維壓縮資訊的目的。
那什麼才是主成分呢?大家應該知道,一個空間會有自己的一組基向量,空間中的任何一個向量都可以通過基向量的組合來表示。
舉個栗子,假如我們有一組2維的點,那麼我們可以找到2個基向量,讓這一組點對著這2個基向量進行投影。每個基向量上,投影后的點越分散,這個基向量就越重要。也就是主成分。不明白的看下圖。(原諒我的靈魂畫作)
紅色的點就是已有的一組2維點,黃色和藍色的線代表兩個基向量的方向。如果我們要把紅色的點從2維降到1維,向哪個方向投影才能儲存更多的資訊呢?直觀上看,當然是黃向量比較好,紅點在黃線上的投影很分散,而在藍線上的投影很多會重疊在一起。黃向量就是比藍向量更主要的成分。
好了,我們現在知道什麼是主成分了,那怎麼才能把它提取出來呢?一般是通過對紅點各維的協方差矩陣進行特徵分解,選取最大的幾個特徵值對應的特徵向量,組成低維空間的一組基向量。再將原始的紅點投影到這組低維基向量上,就得到了紅點降維後的座標。這中間涉及到一些數學的計算和解釋,就放在下面了。
投影
向量的內積相當於投影,
A * B = |A||B|cos(a)
|B|=1時,
A * B= |A|cos(a)
此時,A和B的內積相當於A向B方向投影的長度。
另外,當新基的原點與原基的原點不同時,要先進行變換使原點重合。
討論二維空間基變換
在我們常用的笛卡爾座標系中,一點(x,y)可以分解為
x = r * cos(θ) y = r * sin(θ)
相當於把向量(x,y)投影到兩個基向量上,
(1,0) 和 (0,1)
x和y相當於投影后的長度。
一個向量的準確描述需要給出一組基,以及向量在各個基方向上的投影長度。
注意:不管變換什麼樣的基向量,向量相對於原點的位置是不變的。
例:將基本基上的座標(3,2)變換為基(1/√2,1/√2)、(-1/√2,(1/√2))的座標。
(注意,新基的座標是在基本基上確定的)
多維空間基變換
M個N維向量(組成矩陣A),將其變到R個N維基(組成矩陣B)表示的新空間中
矩陣AB中第m列為A中第m列的變換結果。
兩個矩陣相乘的意義是將右邊矩陣中的每一列列向量變換到左邊矩陣中每一行行向量為基所表示的空間中去。
方差
為實現降維,就要使降維後的結果儘量的分散(重疊太多就不能很好地復原回去),也就是降維後的方差儘量大。
協方差
兩個欄位的協方差表示其相關性,
- 如果結果為正值,則說明兩者是正相關的。
- 如果結果為負值,則說明兩者是負相關的。
- 如果為0,也是就是統計上說的“相互獨立”。
理解協方差矩陣的關鍵就在於牢記它計算的是不同維度之間的協方差,而不是不同樣本之間,拿到一個樣本矩陣,我們最先要明確的就是一行是一個樣本還是一個維度。
a,b兩個欄位組成X,
(相當於m個點的x和y值)
對角線上的為每個欄位各自的方差,其它元素為協方差。
我們要使最後各個方向(欄位)的差距大,如果各個方向相似,這個投影是沒有太大作用的。也就是說協方差儘可能的為0。
PCA計算
X:原始矩陣
C : X對應的協方差矩陣
P : 新基的矩陣
Y = PX
Y : X對P做基變換後的矩陣
D : Y的協方差矩陣
只要找到一組基的矩陣P,使得D為對角矩陣(Y的協方差為0),並且對角元素按從大到小依次排列,那麼P的前K行就是要尋找的基(從N維降到K維)。
又因為C是對稱矩陣,滿足,
(E為正交矩陣)
可得,
所以,我們可以通過對X的協方差矩陣C的特徵分解,來得到E。取E中特徵值最大的幾個特徵值對應的基向量組成P,從而得到降維後的矩陣Y。
在網上看到了一個PCA演算法步驟,覺得挺清楚就貼上過來了。
注:k的選擇。一般選擇佔所有能量99%的特徵值。
程式實現
自己寫了一個簡單的python函式實現PCA,具體的解釋可以看函式前面的註釋~
"""
Author: totodum
Program: Linear_PCA.py
Description: Linear Principle Component Analysis
"""
from numpy import *
from numpy.linalg import *
'''
PCA takes an array of points, each row is a single point. and reduce the dimension
of the points to k-D.
return points of k-D.
'''
def pca(x, k):
m = x.shape[0]
c = dot(transpose(x), x) / m # coefficient
[u, s, v] = svd(c)
p = u[:, 0:k]
y = dot(x, p)
return y
這個是做作業的效果圖(3維降到2維)
三維:
二維: