Kernel PCA 原理和演示

Kernel PCA 原理和演示

主成份（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如何作？

让我先定义以下变量： X=[x1,x2,…,xN] 是一个d×N矩阵，表明输入的数据有N 个，每一个sample的维数是d。咱们作降维，就是想用k维的数据来表示原始的d维数据（k≤d）。
当咱们使用centered的数据（即∑ixi=0）时，可定义协方差矩阵C为：

C=1NxixTi=1NXXT

作特征值分解，咱们能够获得：

CU=UΛ⇒C=UΛUT=∑aλauauTa

注意这里的C,U,Λ的维数都是d×d, 且U=[u1,u2,…,ud], Λ=diag(λ1,λ2,…,λd)。
当咱们作降维时，能够利用前k个特征向量Uk=[u1,u2,…,uk]。则将一个d维的xi向k维的主成分的方向投影后的yi=UTkxi （这里的每个ui都是d维的，表明是一个投影方向，且uTiui=1，表示这是一个旋转变量）

 

1.2 在高维空间里的PCA应该如何作？

高维空间中，咱们定义一个映射Φ:Xd→F，这里F表示Hilbert泛函空间。
如今咱们的输入数据是Φ(xi),i=1,2,…n， 他们的维数能够说是无穷维的（泛函空间）。
在这个新的空间中，假设协方差矩阵一样是centered，咱们的协方差矩阵为：

C¯=1NΦ(xi)Φ(xi)T=1NΦ(X)Φ(X)T

这里有一个陷阱，我跳进去过：
在对Kernel trick只知其一;不知其二的时候，咱们经常从形式上认为C¯能够用Ki,j=K(xi,xj)来代替，
所以对K=(Kij)作特征值分解，而后获得K=UΛUT，而且对原有数据降维的时候，定义Yi=UTkXi。
但这个错误的方法有两个问题：一是咱们不知道矩阵C¯的维数；二是UTkXi从形式上看不出是从高维空间的Φ(Xi)投影，而且当有新的数据时，咱们没法从理论上理解UTkXnew是从高维空间的投影。
若是应用这种错误的方法，咱们有可能获得看起来差很少正确的结果，但本质上这是错误的。
正确的方法是经过Kernel trick将PCA投影的过程经过内积的形式表达出来，详细见1.3

 

1.3 如何用Kernel Trick在高维空间作PCA？

 

在1.1节中，经过PCA，咱们获得了U矩阵。这里将介绍如何仅利用内积的概念来计算传统的PCA。
首先咱们证实U能够由x1,x2,…,xN展开（span)：

Cua=λaua

ua=1λaCu=1λa(∑ixixTi)u=1λa∑ixi(xTiu)=1λa∑i(xTiu)xi=∑ixTiuλaxi=∑iαaixi

这里定义αai=xTiuλa。
由于xTiu 是一个标量（scala），因此αai也是一个标量，所以ui 是能够由xi张成。

 

进而咱们显示PCA投影能够用内积运算表示，例如咱们把xi向任意一个主成分份量ua进行投影，获得的是uTaxi，也就是xTiua 。做者猜想写成这种形式是为了能抽出xTixj=<xi,xj>的内积形式。

 

xTiCuaxTi1N∑jxjxTj∑kαakxk∑jαak∑k(xTixj)(xTjxk)=λaxTiua=λaxTi∑kαakxk=Nλa∑kαak(xTixk)

当咱们定义Kij=xTixj时，上式能够写为K2α=NλaKαa
（这里αa定义为[αa1,αa2,…,αaN]T.）
进一步，咱们获得解为：

Kα=λ~aαwithλ~a=Nλa

K矩阵包含特征值λ~和αa，咱们能够经过α能够计算获得ua，
注意特征值分解时Eigendecomposition，αa只表明一个方向，它的长度通常为1，但在此处不为1。
这里计算出αa的长度（下面将要用到）：
由于ua的长度是1，咱们有：

1=uTaua=(∑iαaixi)T(∑jαajxj)=∑i∑jαaiαajxTixTj=(αa)TKαa=(αa)T(Nλaαa)=Nλa(αaTαa)⇒∥αa∥=1/Nλa−−−−√=1/λ~a−−√

 

在上面的分析过程当中，咱们只使用了内积。所以当咱们把Kij=xTixj推广为Kij=<Φ(xi),Φ(xj>=Φ(xi)TΦ(xj)时，上面的分析结果并不会改变。

1.4 如何在主成分方向上投影？

投影时，只须要使用U矩阵，假设咱们获得的新数据为t，那么t在ua方向的投影是：

uTat=∑iαaixTit=∑iαai(xTit)

对于高维空间的数据Φ(xi),Φ(t)，咱们能够用Kernel trick，用K(xi,t)来带入上式：

uTat=∑iαaiK(xi,t)

 

1.5 如何Centering 高维空间的数据？

在咱们的分析中，协方差矩阵的定义须要centered data。在高维空间中，显式的将Φ(xi)居中并不简单,
由于咱们并不知道Φ的显示表达。但从上面两节能够看出，全部的计算只和K矩阵有关。具体计算以下：
令Φi=Φ(xi)，居中ΦCi=Φi–1N∑kΦk

KCij=<ΦCiΦCj>=(Φi–1N∑kΦk)T(Φj–1N∑lΦl）=ΦTiΦj–1N∑lΦTiΦl–1N∑kΦTkΦj+1N2∑k∑lΦTkΦl=Kij–1N∑lKil–1N∑kKkj+1N2∑k∑lKkl

不难看出，

KC=K–1NK–K1N+1NK1N

其中1N 为N×N的矩阵，其中每个元素都是1/N
对于新的数据，咱们一样能够

K(xi,t)C=<ΦCiΦCt>=(Φi–1N∑kΦk)T(Φt–1N∑lΦl）=ΦTiΦt–1N∑lΦTiΦl–1N∑kΦTkΦt+1N2∑k∑lΦTkΦl=K(xi,t)–1N∑lKil–1N∑kK(xk,t)+1N2∑k∑lKkl

 

2. 演示 (R code)

首先咱们应该注意输入数据的格式，通常在统计中，咱们要求X矩阵是N×d的，但在咱们的推导中，X矩阵是d×N。
这与统计上的概念并不矛盾：在前面的定义下协方差矩阵为XTX，而在后面的定义中是XXT。另外这里的协方差矩阵是样本（Sample）的协方差矩阵，咱们的认为大写的X表明矩阵，而不是表明一个随机变量。
另外，在下面的结果中，Gaussian 核函数（kernel function）的标准差（sd）为2。在其余取值条件下，所获得的图像是不一样的。

KPCA图片：

R 源代码（Source Code）：连接到完整的代码 KernelPCA

Kernel PCA部分代码：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

# 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= "First component" , ylab= "Second component" )

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

No related posts.