PCA 演算法也叫主成分分析(principal components analysis),主要是用於資料降維的。
為什麼要進行資料降維?因為實際情況中我們的訓練資料會存在特徵過多或者是特徵累贅的問題,比如:
- 一個關於汽車的樣本資料,一個特徵是”km/h的最大速度特徵“,另一個是”英里每小時“的最大速度特徵,很顯然這兩個特徵具有很強的相關性
- 拿到一個樣本,特徵非常多,樣本缺很少,這樣的資料用迴歸去你和將非常困難,很容易導致過度擬合
PCA演算法就是用來解決這種問題的,其核心思想就是將 n 維特徵對映到 k 維上(k < n),這 k 維是全新的正交特徵。我們將這 k 維成為主元,是重新構造出來的 k 維特徵,而不是簡單地從 n 維特徵中取出其餘 n-k 維特徵。
PCA 的計算過程
假設我們得到 2 維資料如下:
其中行代表樣例,列代表特徵,這裡有10個樣例,每個樣例有2個特徵,我們假設這兩個特徵是具有較強的相關性,需要我們對其進行降維的。
第一步:分別求 x 和 y 的平均值,然後對所有的樣例都減去對應的均值
這裡求得 x 的均值為 1.81 , y 的均值為 1.91,減去均值後得到資料如下:
注意,此時我們一般應該在對特徵進行方差歸一化,目的是讓每個特徵的權重都一樣,但是由於我們的資料的值都比較接近,所以歸一化這步可以忽略不做
第一步的演算法步驟如下:
本例中步驟3、4沒有做。
第二步:求特徵協方差矩陣
公式如下:
第三步:求解協方差矩陣的特徵值和特徵向量
第四步:將特徵值從大到小進行排序,選擇其中最大的 k 個,然後將其對應的 k 個特徵向量分別作為列向量組成特徵矩陣
這裡的特徵值只有兩個,我們選擇最大的那個,為: 1.28402771 ,其對應的特徵向量為:
注意:matlab 的 eig 函式求解協方差矩陣的時候,返回的特徵值是一個特徵值分佈在對角線的對角矩陣,第 i 個特徵值對應於第 i 列的特徵向量
第五步: 將樣本點投影到選取的特徵向量上
假設樣本列數為 m ,特徵數為 n ,減去均值後的樣本矩陣為 DataAdjust(m*n),協方差矩陣為 n*n ,選取 k 個特徵向量組成後的矩陣為 EigenVectors(n*k),則投影后的資料 FinalData 為:
FinalData (m*k) = DataAdjust(m*n) X EigenVectors(n*k)
得到的結果是:
這樣,我們就將 n 維特徵降成了 k 維,這 k 維就是原始特徵在 k 維上的投影。
整個PCA的過程貌似很簡單,就是求協方差的特徵值和特徵向量,然後做資料轉換。但為什麼協方差的特徵向量就是最理想的 k 維向量?這個問題由PCA的理論基礎來解釋。
PCA 的理論基礎
關於為什麼協方差的特徵向量就是 k 維理想特徵,有3個理論,分別是:
- 最大方差理論
- 最小錯誤理論
- 座標軸相關度理論
這裡簡單描述下最大方差理論:
最大方差理論
訊號處理中認為訊號具有較大的方差,噪聲有較小的方差,訊雜比就是訊號與噪聲的方差比,越大越好。因此我們認為,最好的 k 為特徵既是將 n 維樣本點轉換為 k 維後,每一維上的樣本方差都很大
PCA 處理圖解如下:
降維轉換後:
上圖中的直線就是我們選取的特徵向量,上面例項中PCA的過程就是將空間的2維的點投影到直線上。
那麼問題來了,兩幅圖都是PCA的結果,哪一幅圖比較好呢?
根據最大方差理論,答案是左邊的圖,其實也就是樣本投影后間隔較大,容易區分。
其實從另一個角度看,左邊的圖每個點直線上的距離絕對值之和比右邊的每個點到直線距離絕對值之和小,是不是有點曲線迴歸的感覺?其實從這個角度看,這就是最小誤差理論:選擇投影后誤差最小的直線。
再回到上面的左圖,也就是我們要求的最佳的 u ,前面說了,最佳的 u 也就是最佳的曲線,它能夠使投影后的樣本方差最大或者是誤差最小。
另外,由於我們前面PCA演算法第一步的時候已經執行對樣本資料的每一維求均值,並讓每個資料減去均值的預處理了,所以每個特徵現在的均值都為0,投影到特徵向量上後,均值也為0.因此方差為:
最後的等式中中間的那部分其實就是樣本方差的協方差矩陣(xi 的均值為 0)
由於 u 是單位向量,得到
上式兩邊痛乘以 u,得到:
於是我們得到
最佳投影直線就是特徵值 λ 最大是對應的特徵向量,其次是 λ 第二大對應的特徵向量(求解的到的特徵向量都是正交的)。其中 λ 就是我們的方差,也對應了我們前面的最大方差理論,也就是找到能夠使投影后方差最大的直線。
Matlab 實現
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
function [lowData,reconMat] = PCA(data,K) [row , col] = size(data); meanValue = mean(data); %varData = var(data,1,1); normData = data - repmat(meanValue,[row,1]); covMat = cov(normData(:,1),normData(:,2));%求取協方差矩陣 [eigVect,eigVal] = eig(covMat);%求取特徵值和特徵向量 [sortMat, sortIX] = sort(eigVal,'descend'); [B,IX] = sort(sortMat(1,:),'descend'); len = min(K,length(IX)); eigVect(:,IX(1:1:len)); lowData = normData * eigVect(:,IX(1:1:len)); reconMat = (lowData * eigVect(:,IX(1:1:len))') + repmat(meanValue,[row,1]); % 將降維後的資料轉換到新空間 end |
呼叫方式
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
function testPCA %% clc clear close all %% filename = 'testSet.txt'; K = 1; data = load(filename); [lowData,reconMat] = PCA(data,K); figure scatter(data(:,1),data(:,2),5,'r') hold on scatter(reconMat(:,1),reconMat(:,2),5) hold off end |
效果圖