核PCA——從理論到實(shí)現(xiàn)

1. PCA方法：

設(shè)X=[x1,x2,…xN],x∈Rd，為d維樣本空間的N個(gè)向量h=i。

協(xié)方差矩陣：C=1NXXT

可以做特征值分解C=UΛUT，其中U為特征矩陣，Λ是特征值矩陣。
若只保留k個(gè)維度，只需要以k個(gè)特征向量為基矢，將所有的向量向他們投影，就可以了。
對(duì)于二維的情況如下圖（來(lái)自wikipedia，[2]）

圖中畫出了兩個(gè)特征向量，將所有點(diǎn)向特征向量上投影，就得到中心在0的，1維的點(diǎn)。

概念：
Principal component:主成分

From wikipedia: Principal component analysis?(PCA) is a statistical procedure that uses?orthogonal transformation?to convert a set of observations of possibly correlated variables into a set of values of?linearly uncorrelated?variables called?principal components.

PCA方法是一個(gè)統(tǒng)計(jì)學(xué)方法，它使用正交變換將可能存在相性相關(guān)關(guān)系的變量轉(zhuǎn)換為線性無(wú)關(guān)的變量，這些線性無(wú)關(guān)的向量稱為principal components。也就是說(shuō)，在上圖中的兩個(gè)向量(也就是wikipedia所說(shuō)的變量)為坐標(biāo)的基矢，他們?cè)跇颖究臻g中是線性相關(guān)的，通過(guò)PCA方法，他們轉(zhuǎn)變?yōu)闃颖究臻g中的兩個(gè)特征向量，他們?cè)谠摌颖究臻g中是線性無(wú)關(guān)的（這是顯然的）。一種比較直觀的理解是，PCA方法找到特征值最大的基矢，使得樣本點(diǎn)在這些基矢上的投影的方差盡可能地大。

2. KPCA

主成份（Principal Component Analysis）分析是降維（Dimension Reduction）的重要手段。每一個(gè)主成分都是數(shù)據(jù)在某一個(gè)方向上的投影，在不同的方向上這些數(shù)據(jù)方差Variance的大小由其特征值（eigenvalue）決定。一般我們會(huì)選取最大的幾個(gè)特征值所在的特征向量（eigenvector），這些方向上的信息豐富，一般認(rèn)為包含了更多我們所感興趣的信息。當(dāng)然，這里面有較強(qiáng)的假設(shè)：（1）特征根的大小決定了我們感興趣信息的多少。即小特征根往往代表了噪聲，但實(shí)際上，向小一點(diǎn)的特征根方向投影也有可能包括我們感興趣的數(shù)據(jù)； （2）特征向量的方向是互相正交（orthogonal）的，這種正交性使得PCA容易受到Outlier的影響，例如在【1】中提到的例子（3）難于解釋結(jié)果。例如在建立線性回歸模型（Linear Regression Model）分析因變量（response）和第一個(gè)主成份的關(guān)系時(shí)，我們得到的回歸系數(shù)（Coefficiency）不是某一個(gè)自變量（covariate）的貢獻(xiàn)，而是對(duì)所有自變量的某個(gè)線性組合（Linear Combination）的貢獻(xiàn)。

在Kernel PCA分析之中，我們同樣需要這些假設(shè)，但不同的地方是我們認(rèn)為原有數(shù)據(jù)有更高的維數(shù)，我們可以在更高維的空間（Hilbert Space）中做PCA分析（即在更高維空間里，把原始數(shù)據(jù)向不同的方向投影）。這樣做的優(yōu)點(diǎn)有：對(duì)于在通常線性空間難于線性分類的數(shù)據(jù)點(diǎn)，我們有可能再更高維度上找到合適的高維線性分類平面。我們第二部分的例子就說(shuō)明了這一點(diǎn)。

之前轉(zhuǎn)載了一篇博客，介紹了KPCA的推導(dǎo)，作者主要是參考、翻譯了另一篇客http://zhanxw.com/blog/2011/02/kernel-pca-原理和演示/，這里講KPCA的整個(gè)過(guò)程很詳細(xì)，還有R實(shí)現(xiàn)的代碼，很不錯(cuò)！建議保存下來(lái)看，這個(gè)網(wǎng)站有時(shí)候上不去。

????理論明白了，如果想在matlab里實(shí)現(xiàn)KPCA或者真正去用它，可以參考如下兩個(gè)代碼：

????博客http://blog.sina.com.cn/s/blog_7671b3eb01012d9s.html中實(shí)現(xiàn)了KPCA的整個(gè)過(guò)程，但把特征向量的去單位化注釋掉了，實(shí)際上這段代碼是必要的（因?yàn)閑ig解得的特征向量是單位化的，但實(shí)際要求的向量長(zhǎng)度不是1，詳見(jiàn)公式推導(dǎo)）；代碼的最后還有新數(shù)據(jù)在特征空間的投影。

????還有一個(gè)實(shí)現(xiàn)的代碼是matlab exchange上的http://www.mathworks.com/matlabcentral/fileexchange/27319-kernel-pca，代碼的每一步都有英文說(shuō)明，代碼挺規(guī)范的，只是最后沒(méi)有新數(shù)據(jù)向特征空間的投影代碼。

????上面的有助于加深對(duì)KPCA的理解，如果真正要用，可以用stprtool模式識(shí)別工具箱的kpca函數(shù)，使用很方便，如下例子，這段代碼用KPCA實(shí)現(xiàn)對(duì)一個(gè)二維圓的降維處理，可以看到效果很好：

% 3.3 Kernel Principal Component Analysis
clc
clear
close all

% generate circle data
X = gencircledata([1;1],5,250,1);

% compute kernel PCA
options.ker = 'rbf'; % use RBF kernel
options.arg = 4; % kernel argument
options.new_dim = 2; % output dimension
model = kpca(X,options);

XR = kpcarec(X,model); % compute reconstruced data

% Visualization
figure;
h1 = ppatterns(X);
h2 = ppatterns(XR, '+r');
legend([h1 h2],'Input vectors','Reconstructed');

?

????另嘗試編寫的一個(gè)代碼如下，實(shí)現(xiàn)對(duì)3類數(shù)據(jù)特征空間的聚類，如下圖，紅、綠、藍(lán)三種顏色分別表示3類數(shù)據(jù)，經(jīng)過(guò)rbf核映射到新的空間后，分別聚成了3類：

????這里使用的是rbf核，還可以使用多項(xiàng)式核或sigmoid核，要想得到好的結(jié)果就要不斷改變核參數(shù)，rbf效果好調(diào)一些，sigmoid不易找到好的參數(shù)，多項(xiàng)式核在一些范圍效果不錯(cuò)，但兩個(gè)參數(shù)調(diào)整需要好好找找，如下為一個(gè)多項(xiàng)式核函數(shù)的結(jié)果

????注意，如果你測(cè)試代碼，結(jié)果不一樣很正常，因?yàn)間encircledata函數(shù)產(chǎn)生的數(shù)據(jù)有隨機(jī)成分在。

% 3.3 Kernel Principal Component Analysis
clc
clear
close all

% generate circle data
X0 = gencircledata([1;1],1,250,0.1);
X1 = gencircledata([1;1],3,250,0.1);
X2 = gencircledata([1;1],6,250,0.1);
X0 = X0 - repmat(mean(X0, 2), 1, 250);
X1 = X1 - repmat(mean(X1, 2), 1, 250);
X2 = X2 - repmat(mean(X2, 2), 1, 250);

X = [X0 X1 X2];
y = [ones(1, size(X0, 2)) 2*ones(1, size(X1, 2)) 3*ones(1, size(X2, 2))];
data.X = X;
data.y = y;

figure
ppatterns(data);

% compute kernel PCA
kernelflag = 1;
if kernelflag == 1
????options.ker = 'rbf';
????options.arg = 2; % kernel argument
elseif kernelflag == 2
????options.ker = 'sigmoid';
????options.arg = [3 10]; % kernel argument
else
????options.ker = 'poly';
????options.arg = [2 0]; % kernel argument
end

options.new_dim = 2; % output dimension
model = kpca(data.X, options);
kpca_data = kernelproj(data, model);

figure
ppatterns(kpca_data);

核函數(shù)為高斯徑向基函數(shù)RBF的KPCA

function [train_kpca,test_kpca] = kpcaFordata(train,test,threshold,rbf_var)
%% Data kpca processing
%% 源地址：http://blog.sina.com.cn/lssvm
%% 函數(shù)默認(rèn)設(shè)置
if nargin <4
rbf_var=10000;%？
end
if nargin <3
threshold = 90;
end
%%?數(shù)據(jù)處理
patterns=zscore(train); %訓(xùn)練數(shù)據(jù)標(biāo)準(zhǔn)化
test_patterns=zscore(test); %測(cè)試數(shù)據(jù)標(biāo)準(zhǔn)化
train_num=size(patterns,1); %train_num是訓(xùn)練樣本的個(gè)數(shù)
test_num=size(test_patterns,1);%test_num是測(cè)試樣本的個(gè)數(shù)
cov_size = train_num; %cov_size是訓(xùn)練樣本的個(gè)數(shù)
%%?計(jì)算核矩陣
for i=1:cov_size,
for j=i:cov_size,
K(i,j) = exp(-norm(patterns(i,:)-patterns(j,:))^2/rbf_var); %核函數(shù) rbf_var ？？
K(j,i) = K(i,j);
end
end
unit = ones(cov_size, cov_size)/cov_size;%cov_size是樣本的個(gè)數(shù)
%%?中心化核矩陣
K_n = K - unit*K - K*unit + unit*K*unit;% 中心化核矩陣
%%?特征值分解
[evectors_1,evaltures_1] = eig(K_n/cov_size);
[x,index]=sort(real(diag(evaltures_1))); %sort每行按從小到大排序，x為排序后結(jié)果，index為索引
evals=flipud(x) ;% flipud函數(shù)實(shí)現(xiàn)矩陣的上下翻轉(zhuǎn)
index=flipud(index);
%%?將特征向量按特征值的大小順序排序
evectors=evectors_1(:,index);
%%?單位化特征向量
% for i=1:cov_size
% evecs(:,i) = evectors(:,i)/(sqrt(evectors(:,i)));
% end
train_eigval = 100*cumsum(evals)./sum(evals);
index = find(train_eigval >threshold);

train_kpca = zeros(train_num, index(1)); %train_num是訓(xùn)練樣本的個(gè)數(shù)
%% evecs單位化后的特征矩陣,K_n訓(xùn)練數(shù)據(jù)的中心化核矩陣
train_kpca=[K_n * evectors(:,1:index(1))];
%% 重建測(cè)試數(shù)據(jù)
unit_test = ones(test_num,cov_size)/cov_size;%cov_size是訓(xùn)練樣本的個(gè)數(shù)
K_test = zeros(test_num,cov_size); %test_num是測(cè)試樣本的個(gè)數(shù)， cov_size是訓(xùn)練樣本的個(gè)數(shù)
for i=1:test_num, %test_num是測(cè)試樣本的個(gè)數(shù)
for j=1:cov_size,%cov_size是訓(xùn)練樣本的個(gè)數(shù)
K_test(i,j) = exp(-norm(test_patterns(i,:)-patterns(j,:))^2/rbf_var);
end
end
K_test_n = K_test - unit_test*K - K_test*unit + unit_test*K*unit;
test_kpca = zeros(test_num, index(1));%test_num是測(cè)試樣本的個(gè)數(shù)
test_kpca = [K_test_n * evectors(:,1:index(1))];

總結(jié)

以上是生活随笔為你收集整理的核PCA——从理论到实现的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。