核PCA——从理论到实现

核PCA——从理论到实现

1. PCA方法：

设 X=[x1,x2,…xN],x∈Rd ，为d维样本空间的N个向量 h=i 。

协方差矩阵： C=1NXXT

可以做特征值分解 C=UΛUT ，其中U为特征矩阵， Λ 是特征值矩阵。
若只保留k个维度，只需要以k个特征向量为基矢，将所有的向量向他们投影，就可以了。
对于二维的情况如下图（来自wikipedia，[2]）

图中画出了两个特征向量，将所有点向特征向量上投影，就得到中心在0的，1维的点。

概念：
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方法是一个统计学方法，它使用正交变换将可能存在相性相关关系的变量转换为线性无关的变量，这些线性无关的向量称为principal components。也就是说，在上图中的两个向量(也就是wikipedia所说的变量)为坐标的基矢，他们在样本空间中是线性相关的，通过PCA方法，他们转变为样本空间中的两个特征向量，他们在该样本空间中是线性无关的（这是显然的）。一种比较直观的理解是，PCA方法找到特征值最大的基矢，使得样本点在这些基矢上的投影的方差尽可能地大。

2. KPCA

主成份（Principal Component Analysis）分析是降维（Dimension Reduction）的重要手段。每一个主成分都是数据在某一个方向上的投影，在不同的方向上这些数据方差Variance的大小由其特征值（eigenvalue）决定。一般我们会选取最大的几个特征值所在的特征向量（eigenvector），这些方向上的信息丰富，一般认为包含了更多我们所感兴趣的信息。当然，这里面有较强的假设：（1）特征根的大小决定了我们感兴趣信息的多少。即小特征根往往代表了噪声，但实际上，向小一点的特征根方向投影也有可能包括我们感兴趣的数据； （2）特征向量的方向是互相正交（orthogonal）的，这种正交性使得PCA容易受到Outlier的影响，例如在【1】中提到的例子（3）难于解释结果。例如在建立线性回归模型（Linear Regression Model）分析因变量（response）和第一个主成份的关系时，我们得到的回归系数（Coefficiency）不是某一个自变量（covariate）的贡献，而是对所有自变量的某个线性组合（Linear Combination）的贡献。

在Kernel PCA分析之中，我们同样需要这些假设，但不同的地方是我们认为原有数据有更高的维数，我们可以在更高维的空间（Hilbert Space）中做PCA分析（即在更高维空间里，把原始数据向不同的方向投影）。这样做的优点有：对于在通常线性空间难于线性分类的数据点，我们有可能再更高维度上找到合适的高维线性分类平面。我们第二部分的例子就说明了这一点。

之前转载了一篇博客，介绍了KPCA的推导，作者主要是参考、翻译了另一篇客http://zhanxw.com/blog/2011/02/kernel-pca-原理和演示/，这里讲KPCA的整个过程很详细，还有R实现的代码，很不错！建议保存下来看，这个网站有时候上不去。

    理论明白了，如果想在matlab里实现KPCA或者真正去用它，可以参考如下两个代码：

    博客http://blog.sina.com.cn/s/blog_7671b3eb01012d9s.html中实现了KPCA的整个过程，但把特征向量的去单位化注释掉了，实际上这段代码是必要的（因为eig解得的特征向量是单位化的，但实际要求的向量长度不是1，详见公式推导）；代码的最后还有新数据在特征空间的投影。

    还有一个实现的代码是matlab exchange上的http://www.mathworks.com/matlabcentral/fileexchange/27319-kernel-pca，代码的每一步都有英文说明，代码挺规范的，只是最后没有新数据向特征空间的投影代码。

    上面的有助于加深对KPCA的理解，如果真正要用，可以用stprtool模式识别工具箱的kpca函数，使用很方便，如下例子，这段代码用KPCA实现对一个二维圆的降维处理，可以看到效果很好：

% 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');

 

    另尝试编写的一个代码如下，实现对3类数据特征空间的聚类，如下图，红、绿、蓝三种颜色分别表示3类数据，经过rbf核映射到新的空间后，分别聚成了3类：

    这里使用的是rbf核，还可以使用多项式核或sigmoid核，要想得到好的结果就要不断改变核参数，rbf效果好调一些，sigmoid不易找到好的参数，多项式核在一些范围效果不错，但两个参数调整需要好好找找，如下为一个多项式核函数的结果

    注意，如果你测试代码，结果不一样很正常，因为gencircledata函数产生的数据有随机成分在。

% 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);

核函数为高斯径向基函数RBF的KPCA

function [train_kpca,test_kpca] = kpcaFordata(train,test,threshold,rbf_var)
%% Data kpca processing
%% 源地址：http://blog.sina.com.cn/lssvm
%% 函数默认设置
if nargin <4
rbf_var=10000;%？
end
if nargin <3
threshold = 90;
end
%% 数据处理
patterns=zscore(train); %训练数据标准化
test_patterns=zscore(test); %测试数据标准化
train_num=size(patterns,1); %train_num是训练样本的个数
test_num=size(test_patterns,1);%test_num是测试样本的个数
cov_size = train_num; %cov_size是训练样本的个数
%% 计算核矩阵
for i=1:cov_size,
for j=i:cov_size,
K(i,j) = exp(-norm(patterns(i,:)-patterns(j,:))^2/rbf_var); %核函数 rbf_var ？？
K(j,i) = K(i,j);
end
end
unit = ones(cov_size, cov_size)/cov_size;%cov_size是样本的个数
%% 中心化核矩阵
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为排序后结果，index为索引
evals=flipud(x) ;% flipud函数实现矩阵的上下翻转
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是训练样本的个数
%% evecs单位化后的特征矩阵,K_n训练数据的中心化核矩阵
train_kpca=[K_n * evectors(:,1:index(1))];
%% 重建测试数据
unit_test = ones(test_num,cov_size)/cov_size;%cov_size是训练样本的个数
K_test = zeros(test_num,cov_size); %test_num是测试样本的个数， cov_size是训练样本的个数
for i=1:test_num, %test_num是测试样本的个数
for j=1:cov_size,%cov_size是训练样本的个数
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是测试样本的个数
test_kpca = [K_test_n * evectors(:,1:index(1))];