KPCA原理及演示

1、主成份（Principal Component Analysis）分析是降维（Dimension Reduction）的重要手段。每一个主成分都是数据在某一个方向上的投影，在不同的方向上这些数据方差Variance的大小由其特征值（eigenvalue）决定。一般我们会选取最大的几个特征值所在的特征向量（eigenvector），这些方向上的信息丰富，一般认为包含了更多我们所感兴趣的信息。当然，这里面有较强的假设：（1）特征根的大小决定了我们感兴趣信息的多少。即小特征根往往代表了噪声，但实际上，向小一点的特征根方向投影也有可能包括我们感兴趣的数据； （2）特征向量的方向是互相正交（orthog

2、onal）的，这种正交性使得PCA容易受到Outlier的影响，例如在【1】中提到的例子（3）难于解释结果。例如在建立线性回归模型（Linear Regression Model）分析因变量（response）和第一个主成份的关系时，我们得到的回归系数（Coefficiency）不是某一个自变量（covariate）的贡献，而是对所有自变量的某个线性组合（Linear Combination）的贡献。在Kernel PCA分析之中，我们同样需要这些假设，但不同的地方是我们认为原有数据有更高的维数，我们可以在更高维的空间（Hilbert Space）中做PCA分析（即在更高维空间里，把原始数据向

3、不同的方向投影）。这样做的优点有：对于在通常线性空间难于线性分类的数据点，我们有可能再更高维度上找到合适的高维线性分类平面。我们第二部分的例子就说明了这一点。本文写作的动机是因为作者没有找到一篇好的文章（看了wikipedia和若干google结果后）深层次介绍PCA和Kernel PCA之间的联系，以及如何以公式形式来解释如何利用Kernel PCA来做投影，特别有些图片的例子只是展示了结果和一些公式，这里面具体的过程并没有涉及。希望这篇文章能做出较好的解答。1. Kernel Principal Component Analysis 的矩阵基础我们从解决这几个问题入手：传统的PCA如何做？

4、在高维空间里的PCA应该如何做？如何用Kernel Trick在高维空间做PCA？如何在主成分方向上投影？如何Centering 高维空间的数据？1.1 传统的PCA如何做？让我先定义如下变量：X=x1,x2,xN是一个dN矩阵，代表输入的数据有N个，每个sample的维数是d。我们做降维，就是想用k维的数据来表示原始的d维数据（kd）。当我们使用centered的数据（即ixi=0）时，可定义协方差矩阵C为：C=1NxixTi=1NXXT做特征值分解，我们可以得到：CU=UC=UUT=aauauTa注意这里的C,U,的维数都是dd, 且U=u1,u2,ud,=diag(1,2,d)。当我们做

5、降维时，可以利用前k个特征向量Uk=u1,u2,uk。则将一个d维的xi向k维的主成分的方向投影后的yi=UTkxi（这里的每一个ui都是d维的，代表是一个投影方向，且uTiui=1，表示这是一个旋转变量）1.2 在高维空间里的PCA应该如何做？高维空间中，我们定义一个映射:XdF，这里F表示Hilbert泛函空间。现在我们的输入数据是(xi),i=1,2,n， 他们的维数可以说是无穷维的（泛函空间）。在这个新的空间中，假设协方差矩阵同样是centered，我们的协方差矩阵为：C=1N(xi)(xi)T=1N(X)(X)T这里有一个陷阱，我跳进去过：在对Kernel trick一知半解的时候，

6、我们常常从形式上认为C可以用Ki,j=K(xi,xj)来代替，因此对K=(Kij)做特征值分解，然后得到K=UUT，并且对原有数据降维的时候，定义Yi=UTkXi。但这个错误的方法有两个问题：一是我们不知道矩阵C的维数；二是UTkXi从形式上看不出是从高维空间的(Xi)投影，并且当有新的数据时，我们无法从理论上理解UTkXnew是从高维空间的投影。如果应用这种错误的方法，我们有可能得到看起来差不多正确的结果，但本质上这是错误的。正确的方法是通过Kernel trick将PCA投影的过程通过内积的形式表达出来，详细见1.31.3 如何用Kernel Trick在高维空间做PCA？在1.1节中，通

7、过PCA，我们得到了U矩阵。这里将介绍如何仅利用内积的概念来计算传统的PCA。首先我们证明U可以由x1,x2,xN展开（span)：Cua=auaua=1aCu=1a(ixixTi)u=1aixi(xTiu)=1ai(xTiu)xi=ixTiuaxi=iaixi这里定义ai=xTiua。因为xTiu是一个标量（scala），所以ai也是一个标量，因此ui是可以由xi张成。进而我们显示PCA投影可以用内积运算表示，例如我们把xi向任意一个主成分分量ua进行投影，得到的是uTaxi，也就是xTiua。作者猜测写成这种形式是为了能抽出xTixj=的内积形式。xTiCuaxTi1NjxjxTjkakx

8、kjakk(xTixj)(xTjxk)=axTiua=axTikakxk=Nakak(xTixk)当我们定义Kij=xTixj时，上式可以写为K2=NaKa（这里a定义为a1,a2,aNT.）进一步，我们得到解为：K=awitha=NaK矩阵包含特征值和a，我们可以通过可以计算得到ua，注意特征值分解时Eigendecomposition，a只代表一个方向，它的长度一般为1，但在此处不为1。这里计算出a的长度（下面将要用到）：因为ua的长度是1，我们有：1=uTaua=(iaixi)T(jajxj)=ijaiajxTixTj=(a)TKa=(a)T(Naa)=Na(aTa)a=1/Na=1/a

9、在上面的分析过程中，我们只使用了内积。因此当我们把Kij=xTixj推广为Kij=(xi)T(xj)时，上面的分析结果并不会改变。1.4 如何在主成分方向上投影？投影时，只需要使用U矩阵，假设我们得到的新数据为t，那么t在ua方向的投影是：uTat=iaixTit=iai(xTit)对于高维空间的数据(xi),(t)，我们可以用Kernel trick，用K(xi,t)来带入上式：uTat=iaiK(xi,t)1.5 如何Centering 高维空间的数据？在我们的分析中，协方差矩阵的定义需要centered data。在高维空间中，显式的将(xi)居中并不简单,因为我们并不知道的显示表达。但

10、从上面两节可以看出，所有的计算只和K矩阵有关。具体计算如下：令i=(xi)，居中Ci=i1NkkKCij=(i1Nkk)T(j1Nll）=Tij1NlTil1NkTkj+1N2klTkl=Kij1NlKil1NkKkj+1N2klKkl不难看出，KC=K1NKK1N+1NK1N其中1N为NN的矩阵，其中每一个元素都是1/N对于新的数据，我们同样可以K(xi,t)C=(i1Nkk)T(t1Nll）=Tit1NlTil1NkTkt+1N2klTkl=K(xi,t)1NlKil1NkK(xk,t)+1N2klKkl2. 演示 (R code)首先我们应该注意输入数据的格式，一般在统计中，我们要求X矩

11、阵是Nd的，但在我们的推导中，X矩阵是dN。这与统计上的概念并不矛盾：在前面的定义下协方差矩阵为XTX，而在后面的定义中是XXT。另外这里的协方差矩阵是样本（Sample）的协方差矩阵，我们的认为大写的X代表矩阵，而不是代表一个随机变量。另外，在下面的结果中，Gaussian 核函数（kernel function）的标准差（sd）为2。在其他取值条件下，所得到的图像是不同的。KPCA图片：R 源代码（Source Code）：链接到完整的代码KernelPCAKernel PCA部分代码：12345678910111213141516171819202122232425# Kernel PC

12、A# Polynomial Kernel# k(x,y) = t(x) %*% y + 1k1 = function (x,y) (x1 * y1 + x2 * y2 + 1)2 K = matrix(0, ncol = N_total, nrow = N_total)for (i in 1:N_total) for (j in 1:N_total) Ki,j = k1(Xi, Xj,)ones = 1/N_total* matrix(1, N_total, N_total)K_norm = K - ones %*% K - K %*% ones + ones %*% K %*% onesre

13、s = eigen(K_norm)V = res$vectorsD = diag(res$values)rank = 0for (i in 1:N_total) if (Di,i 1e-6) break V,i = V,i / sqrt (Di,i)rank = rank + 1Y = K_norm %*% V,1:rankplot(Y,1, Y,2, col = rainbow(3)label, main = Kernel PCA (Poly), xlab=First component, ylab=Second component)3. 主要参考资料【1】A Tutorial on Principal Component Analysis ,Jonathon Shlens,Shlens03【2】Wikipedia： /wiki/Kernel_principal_component_analysis【3】 Original KPCA Paper：Kernel principal component analysis，Bernhard Schlkopf, Alexander Smola and Klaus-Robert Mllerhttp:/www.sp