Kernel PCA 原理和演示
转载自：https://zhuanlan.zhihu.com/p/21583787，版权归原作者所有。
主成份（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分析（即在更高维空间里，把原始数据向不同的方向投影）。这样做的优点有：对于在通常线性空间难于线性分类的数据点，我们有可能再更高维度上找到合适的高维线性分类平面。我们第二部分的例子就说明了这一点。
本文写作的动机是因为作者没有找到一篇好的文章（看了wikipedia和若干google结果后）深层次介绍PCA和Kernel PCA之间的联系，以及如何以公式形式来解释如何利用Kernel PCA来做投影，特别有些图片的例子只是展示了结果和一些公式，这里面具体的过程并没有涉及。希望这篇文章能做出较好的解答。
1. Kernel Principal Component Analysis 的矩阵基础
我们从解决这几个问题入手：传统的PCA如何做？在高维空间里的PCA应该如何做？如何用Kernel Trick在高维空间做PCA？如何在主成分方向上投影？如何Centering 高维空间的数据？
1.1 传统的PCA如何做？
让我先定义如下变量： 是一个
矩阵，代表输入的数据有
个，每个sample的维数是
。我们做降维，就是想用
维的数据来表示原始的
维数据(
)。
当我们使用centered的数据（即）时，可定义协方差矩阵
为：
做特征值分解，我们可以得到：
注意这里的
当我们做降维时，可以利用前
1.2 在高维空间里的PCA应该如何做？
高维空间中，我们定义一个映射，这里
表示Hilbert泛函空间。
现在我们的输入数据是，他们的维数可以是无穷维的（泛函空间），为了方便讨论，设其空间的维数为
(
)。
在这个新的空间(特征空间)中，将中心化的数据记为，即
其中
将中心化的核矩阵记为(具体求法见*****)，且定义
中心化的数据（
）的协方差矩阵为
：
这里有一个陷阱：
在对Kernel trick一知半解的时候，我们常常从形式上认为
因此对
但这个错误的方法有两个问题：一是我们不知道矩阵
如果应用这种错误的方法，我们有可能得到看起来差不多正确的结果，但本质上这是错误的。
正确的方法是通过Kernel trick将PCA投影的过程通过内积的形式表达出来，详细见1.3
1.3 如何用Kernel Trick在高维空间做PCA？
在1.1节中，通过PCA，我们得到了矩阵。这里将介绍如何仅利用内积的概念来计算传统的PCA。
首先我们证明U可以由展开（span)：
这里定义。
因为 是一个标量(scala)，所以
也是一个标量，因此
是可以由
张成。
进而我们显示PCA投影可以用内积运算表示，例如我们把向任意一个主成分分量
进行投影，得到的是
，也就是
。作者猜测写成这种形式是为了能抽出
的内积形式。
也就是说，对于，
将其写成矩阵的形式为：
当我们定义时，上式可以写为
（这里
定义为
.）
进一步，我们得到解为：
注意特征值分解(Eigen decomposition)时，
因为的长度是1，我们有：
所以
在上面的分析过程中，我们使用了内积的思想完成了传统PCA的过程，即可以通过核矩阵的特征值
和特征向量
计算协方差矩阵的特征值
和特征向量
（PS：不要问为什么不直接用协方差矩阵进行特征值分解(SVD)，得到想要的主方向。请记住在特征空间的协方差矩阵式是不能得到的，所以不能进行大侠想的特征分解）。
1.4 KPCA实现
下面我们重复1.3过程，看看KPCA中的有关推导。
设的特征值和对应的特征向量分别为
和
，即
由1.3知，，即存在
，使得：
用同时作用于
得
上式左边可以化为：
其中
,
右边可以化为：
让依次取遍
的所有数，得到上面推到的“左边=右边”的矩阵表示：
化简得到
即之前假设的为中心化的核矩阵
的第
个特征向量（对应的特征值为
）。现在我们可以用所求的
和
去表示协方差矩阵
的特征值
和向量
：
，
，其中
现在还有个小问题没有解决，在传统PCA方法中，协方差矩阵得到的特征向量是单位正交的。同样，我们希望也是单位向量，即：
因此，将的特征向量
的长度归一化到
，便能保证
的长度为1。
下面证明正交性：
得证！
最后总结下之前的结论，特征空间中协方差矩阵的特征值为
，对应的特征向量为
,
，其中
和
为
的
特征值和特征向量。
由于映射未知，直接在特征空间内进行数据中心化不可行，直接将输入数据在特征空间中心化的过程见1.6。
1.5 如何在主成分方向上投影？
传统PCA投影时，只需要使用矩阵，假设我们得到的新数据为
，那么
在
方向的投影是：
对于高维空间的数据,
，我们可以用Kernel
trick，用
来带入上式：
一般不能显式得到，从而需要计算输入数据
的中心化核矩阵，记为
。
1.5 如何Centering 高维空间的数据？
在我们的分析中，协方差矩阵的定义需要centered data。在高维空间中，显式的将中心化并不简单, 因为我们并不知道
的显示表达。但从上面两节可以看出，所有的计算只和
矩阵有关。具体计算如下：
令，居中
：
不难看出，
其中
对于新的数据
所以：
如果数据，则上式中有
，且
。
2. 演示 (R code)
首先我们应该注意输入数据的格式，一般在统计中，我们要求矩阵是
的，但在我们的推导中，
矩阵是
。
这与统计上的概念并不矛盾：在前面的定义下协方差矩阵为，而在后面的定义中是
。另外这里的协方差矩阵是样本（Sample）的协方差矩阵，我们的认为大写的
代表矩阵，而不是代表一个随机变量。
另外，在下面的结果中，Gaussian 核函数（kernel function）的标准差（sd）为2。在其他取值条件下，所得到的图像是不同的。
KPCA图片：
R 源代码（Source Code）：链接到完整的代码 KernelPCA
Kernel PCA部分代码：
# Kernel PCA
# Polynomial Kernel
# k(x,y) = t(x) %*% y + 1
k1 = function (x,y) { (x[1] * y[1] + x[2] * y[2] + 1)^2 }
K = matrix(0, ncol = N_total, nrow = N_total)
for (i in 1:N_total) {
for (j in 1:N_total) {
K[i,j] = k1(X[i,], X[j,])
}}
ones = 1/N_total* matrix(1, N_total, N_total)
K_norm = K - ones %*% K - K %*% ones + ones %*% K %*% ones
res = eigen(K_norm)
V = res$vectors
D = diag(res$values)
rank = 0
for (i in 1:N_total) {
if (D[i,i] < 1e-6) { break }
V[,i] = V[,i] / sqrt (D[i,i])
rank = rank + 1
}
Y = K_norm %*% V[,1:rank]
plot(Y[,1], Y[,2], col = rainbow(3)[label], main = "Kernel PCA (Poly)"
, xlab="x", ylab="y")
function demo_kpca_ushape
%% PARAMETERS
N_split = [30, 30, 60]; % number of data points in 2 straight parts and curve of "U"
corners = [3 4 2 4]; % corners: x1 x2 y1 y2
nvar = 0.05; % noise variance
kernel.type = 'gauss';
kernel.par = .5;
numeig = 4; % number of eigenvalues to plot
Ntest = [50,50]; % number of test grid divisions, in each dimension
border = [0 6 0 6]; % test grid border
%% PROGRAM
tic
%% generate data
N = sum(N_split);
N1 = N_split(1); N2 = N_split(2); N3 = N_split(3);
X1 = [linspace(corners(1),corners(2),N1)' repmat(corners(3),N1,1)]; % lower straight part
X2 = [linspace(corners(1),corners(2),N2)' repmat(corners(4),N2,1)]; % upper straight part
angles = linspace(pi/2,3*pi/2,N3)';
rad = (corners(4)-corners(3))/2;
m3 = [corners(1) corners(3)+rad];
X3 = repmat(m3,N3,1) + rad*[cos(angles) sin(angles)];
n = nvar*randn(N,2);
X = [X1; X2; X3] + n;
%% generate test grid data
Nt1 = Ntest(1); Nt2 = Ntest(2);
Xtest = zeros(Nt1*Nt2,2);
absc = linspace(border(1),border(2),Nt1);
ordi = linspace(border(3),border(4),Nt2);
for i=1:Nt1,
for j=1:Nt2,
Xtest((i-1)*Nt2+j,:) = [absc(i) ordi(j)];
end
end
%% calculate kernel principal components and projections of Xtest
[E,v] = km_kpca(X,numeig,kernel.type,kernel.par);
Kt = km_kernel(X,Xtest,kernel.type,kernel.par); % kernels of test data set
Xtestp = E'*Kt; % projections of test data set on the principal directions
Y = cell(numeig,1);
for i=1:numeig,
Y{i} = reshape(Xtestp(i,:),Nt2,Nt1); % shape into 2D grid
end
toc
%% OUTPUT
subplot(3,2,1)
plot(X(:,1),X(:,2),'.')
% axis(border)
title('original data of shape U')
% fireworks!
for i=1:numeig,
subplot(2,3,i+2); hold on
[C,h] = contourf(-interp2(Y{i},2),25);
plot(X(:,1)/border(2)*(Nt1*4-1)+1,X(:,2)/border(4)*(Nt2*4-1)+1,'o',...
'MarkerFaceColor','White','MarkerEdgeColor','Black')
title(['projections to the ',num2str(i), 'th principal directions'])
set(gca, 'XTick', [])
set(gca, 'YTick', [])
end
%Additional function
function [E,v] = km_kpca(X,m,ktype,kpar)
n = size(X,1);
K = km_kernel(X,X,ktype,kpar);
[E,V] = eig(K);
v = diag(V); % eigenvalues
[v,ind] = sort(v,'descend');
v = v(1:m);
E = E(:,ind(1:m)); % principal components
for i=1:m
E(:,i) = E(:,i)/sqrt(n*v(i)); % normalization
end
function K = km_kernel(X1,X2,ktype,kpar)
switch ktype
case 'gauss' % Gaussian kernel
sgm = kpar; % kernel width
dim1 = size(X1,1);
dim2 = size(X2,1);
norms1 = sum(X1.^2,2);
norms2 = sum(X2.^2,2);
mat1 = repmat(norms1,1,dim2);
mat2 = repmat(norms2',dim1,1);
distmat = mat1 + mat2 - 2*X1*X2'; % full distance matrix
K = exp(-distmat/(2*sgm^2));
case 'gauss-diag' % only diagonal of Gaussian kernel
sgm = kpar; % kernel width
K = exp(-sum((X1-X2).^2,2)/(2*sgm^2));
case 'poly' % polynomial kernel
p = kpar(1); % polynome order
c = kpar(2); % additive constant
K = (X1*X2' + c).^p;
case 'linear' % linear kernel
K = X1*X2';
otherwise % default case
error ('unknown kernel type')
end
结果
3. 主要参考资料
[1] A Tutorial on Principal Component Analysis ,Jonathon Shlens, Shlens03
[2] Wikipedia： http://en.wikipedia.org/wiki/Kernel_principal_component_analysis
[3] Original KPCA Paper：Kernel principal component analysis，Bernhard Schölkopf, Alexander Smola and Klaus-Robert Müller http://www.springerlink.com/content/w0t1756772h41872/fulltext.pdf
[4] Max Wellings’s classes notes for machine learning Kernel Principal Component Analaysis http://www.ics.uci.edu/~welling/classnotes/papers_class/Kernel-PCA.pdf
转载（纠正部分错误）http://zhanxw.com/blog/2011/02/kernel-pca-原理和演示/