KPCA原理及演示汇编

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

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

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

4、. Kernel Principal Component Analysis的矩阵基础我们从解决这几个问题入手：传统的 PCA 如何做？在高维空间里的 PCA 应该如何做？如 何用 Kernel Trick 在高维空间做 PCA ？如何在主成分方向上投影？如何 Centering 高维空 间的数据？1.1 传统的 PCA 如何做？让我先定义如下变量： X=x1,x2,xN 是一个 dN 矩阵，代表输入的数据有 N 个， 每个 sample 的维数是 d 。我们做降维，就是想用 k 维的数据来表示原始的 d 维数据 (kd)。当我们使用 centered 的数据(即 ixi=0 )时，可定义协方差

5、矩阵 C 为：C=1 Nx ixTi =1 NXX T做特征值分解，我们可以得到：CU=U?C=UUT=aauauTa注意这里的 C,U, 的维数都是 dd, 且U= u1,u2, u,d , =diag ( 1 , 2 , ,d)。 当我们做降维时，可以利用前 k 个特征向量 Uk=u1,u2, ,u k 。则将一个 d 维的 xi 向k 维的主成分的方向投影后的 yi=UTkxi (这里的每一个 ui都是 d 维的，代表是一个投影方向，且 uTiui=1 ，表示这是一个旋转变量)1.2 在高维空间里的 PCA 应该如何做？高维空间中，我们定义一个映射 :Xd F，这里 F 表示 Hilbe

6、rt 泛函空间。 现在我们的输入数据是 (xi),i=1,2, n， 他们的维数可以说是无穷维的 (泛函空间) 。 在这个新的空间中，假设协方差矩阵同样是 centered ，我们的协方差矩阵为：C =N1(xi) x(i )T =1 N(X) X()T这里有一个陷阱，我跳进去过：在对 Kernel trick 一知半解的时候，我们常常从形式上认为 C可以用Ki,j=K(xi,xj)来代替，因此对 K=(Kij )做特征值分解，然后得到 K=UUT，并且对原有数据降维的 时候，定义 Yi=UTkXi。但这个错误的方法有两个问题：一是我们不知道矩阵 C的维数；二是 UTkXi 从形式上看不出是从

7、高维空间的 (Xi) 投影，并且当有新的数据时，我们无法 从理论上理解 U Tk X new 是从高维空间的投影。如果应用这种错误的方法， 我们有可能得到看起来差不多正确的结果， 但本质上 这是错误的。正确的方法是通过 Kernel trick 将 PCA投影的过程通过内积的形式表达出来， 详细见 1.31.3 如何用 Kernel Trick 在高维空间做 PCA ？ 在 1.1 节中，通过 PCA ，我们得到了 U 矩阵。这里将介绍如何仅利用内积的概念来计算传 统的 PCA 。首先我们证明 U 可以由 x1,x2, ,xN 展开( span) ：Cua=auaua=1 aCu=1 a(ix

8、ixTi)u=1 aixi(xTiu)=1 ai(xTiu)xi=ixTiu axi=iaixi这里定义 ai = xTi ua。因为 xTiu 是一个标量( scala )，所以ai 也是一个标量，因此 ui 是可以由 xi 张成。进而我们显示 PCA 投影可以用内积运算表示， 例如我们把 xi 向任意一个主成分分量 ua 进 行投影，得到的是 uTaxi，也就是 xTiua 。作者猜测写成这种形式是为了能抽出 xTixj=的内积形式。xTiCuaxTi1NjxjxTjkakxkjakk(xTixj)(xTjxk)=axTiua=axTikakxk=Nakak(xTixk)当我们定义 Kij

9、=xTixj 时，上式可以写为 K2 =NaKa(这里 a定义为 a 1 , a2 , ,aN T. ) 进一步，我们得到解为：K = awith a=N aK 矩阵包含特征值 和 a，我们可以通过 可以计算得到 ua ， 注意特征值分解时 Eigendecomposition ，a 只代表一个方向，它的长度一般为 1，但在此处不为 1。这里计算出 a 的长度(下面将要用到)：因为 ua 的长度是 1，我们有：1= uTaua=(iaixi)T(jajxj)=ijaiajxTixTj=( a)TKa=(a)T(N aa)=Na(aTa)? a=1/ N a =1/ a- 在上面的分析过程中，我

10、们只使用了内积。因此当我们把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) 居中并不简单 ,因为我们

11、并不知道 的显示表达。 但从上面两节可以看出， 所有的计算只和 K 矩阵有关。 具体计算如下：令 i =(xi )，居中 Ci= i1NkkKCij =( i1 Nkk)T( j1 Nll)=Tij1NlTil1NkTkj+1N2klTkl=Kij1NlKil1 N kK kj +1 N2 k lK kl不难看出，KC=K1NKK1N+1NK1 N其中 1N 为 NN 的矩阵，其中每一个元素都是 1/ N 对于新的数据，我们同样可以K(xi,t)C=( i1Nkk)T( t 1 Nll)=Tit1NlTil1NkTkt+1N2klTkl=K(xi,t)1NlKil1NkK(xk,t)+1 N

12、2 k lK kl2. 演示 (R code)首先我们应该注意输入数据的格式，一般在统计中，我们要求 X 矩阵是 N d 的，但在 我们的推导中， X 矩阵是 d N。这与统计上的概念并不矛盾：在前面的定义下协方差矩阵为X TX ，而在后面的定义中是XXT 。另外这里的协方差矩阵是样本( Sample )的协方差矩阵，我们的认为大写的 X 代 表矩阵，而不是代表一个随机变量。2。在其另外，在下面的结果中， Gaussian 核函数( kernel function )的标准差( sd )为 他取值条件下，所得到的图像是不同的。KPCA 图片：R 源代码( Source Code )：链接到完整

13、的代码 KernelPCAKernel PCA 部分代码：1 # Kernel PCA2 # Polynomial Kernel3 # k(x,y) = t(x) %*% y + 14 k1 = function (x,y) (x1 * y1 + x2 * y2 + 1)2 5 K = matrix (0, ncol = N_total, nrow = N_total)6 for (i in 1:N_total) 7 for (j in 1:N_total) 8 Ki,j = k1 (Xi, Xj,)9 10 ones = 1/N_total*matrix (1, N_total, N_tot

14、al)11 K_norm = K - ones %*% K - K %*% ones + ones %*% K %*% ones12 res = eigen (K_norm)1314 V = res$vectors15 D = diag (res$values)1617 rank = 018 for (i in 1:N_total) 19 if (Di,i 1e-6) break 20 V,i = V,i /sqrt (Di,i)21 rank = rank + 122 23 Y = K_norm %*% V,1:rank24 plot (Y,1, Y,2, col =rainbow (3)l

15、abel, main =Kernel PCA (Poly)25 , 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 analysi