Kernel PCA 原理和演示

来源：互联网发布：ubuntu cuda9.0 caffe 编辑：程序博客网时间：2024/04/30 09:02

主成份（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 = 1 N x i x T i = 1 N X X T

做特征值分解，我们可以得到：

C U = U Λ \Rightarrow C = U Λ U T = \sum a λ a u a u T a

注意这里的

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 ¯ = 1 N Φ (x i) Φ (x i) T = 1 N Φ (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)：

C u a = λ a u a

u a = 1 λ a C u = 1 λ a (\sum i x i x T i) u = 1 λ a \sum i x i (x T i u) = 1 λ a \sum i (x T i u) x i = \sum i x T i u λ a x i = \sum i α a i x i

这里定义

αai=xTiuλa。
因为

xTiu 是一个标量（scala），所以

αai也是一个标量，因此

ui 是可以由

xi张成。

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

x T i C u a x T i 1 N \sum j x j x T j \sum k α a k x k \sum j α a k \sum k (x T i x j) (x T j x k) = λ a x T i u a = λ a x T i \sum k α a k x k = N λ a \sum k α a k (x T i x k)

当我们定义

Kij=xTixj时，上式可以写为

K2α=NλaKαa
（这里

αa定义为

[αa1,αa2,…,αaN]T.）
进一步，我们得到解为：

K α = λ ~ a α w i t h λ ~ a = N λ a

K矩阵包含特征值

λ~和

αa，我们可以通过

α可以计算得到

ua，
注意特征值分解时Eigendecomposition，

αa只代表一个方向，它的长度一般为1，但在此处不为1。
这里计算出

αa的长度（下面将要用到）：
因为

ua的长度是1，我们有：

1 = u T a u a = (\sum i α a i x i) T (\sum j α a j x j) = \sum i \sum j α a i α a j x T i x T j = (α a) T K α a = (α a) T (N λ a α a) = N λ a (α a T α a) \Rightarrow ∥ α a ∥ = 1 / N λ a - - - - \sqrt = 1 / λ ~ a - - \sqrt

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

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

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

u T a t = \sum i α a i x T i t = \sum i α a i (x T i t)

对于高维空间的数据

Φ(xi),Φ(t)，我们可以用Kernel trick，用

K(xi,t)来带入上式：

u T a t = \sum i α a i K (x i, t)

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

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

K C i j = < Φ C i Φ C j > = (Φ i - 1 N \sum k Φ k) T (Φ j - 1 N \sum l Φ l ） = Φ T i Φ j - 1 N \sum l Φ T i Φ l - 1 N \sum k Φ T k Φ j + 1 N 2 \sum k \sum l Φ T k Φ l = K i j - 1 N \sum l K i l - 1 N \sum k K k j + 1 N 2 \sum k \sum l K k l

不难看出，

K C = K - 1 N K - K 1 N + 1 N K 1 N

其中

1N 为

N×N的矩阵，其中每一个元素都是

1/N
对于新的数据，我们同样可以

K (x i, t) C = < Φ C i Φ C t > = (Φ i - 1 N \sum k Φ k) T (Φ t - 1 N \sum l Φ l ） = Φ T i Φ t - 1 N \sum l Φ T i Φ l - 1 N \sum k Φ T k Φ t + 1 N 2 \sum k \sum l Φ T k Φ l = K (x i, t) - 1 N \sum l K i l - 1 N \sum k K (x k, t) + 1 N 2 \sum k \sum l K k l

2. 演示 (R code)

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

KPCA图片：

源地址：http://zhanxw.com/blog/2011/02/kernel-pca-原理和演示/

R 源代码（Source Code）：

# 3 groups of sample, each size is 20N_in_group = 40N_group = 3N_total = N_in_group * N_groupx=matrix(0, ncol = 2, nrow = N_total)label = rep(c(1,2,3), rep(N_in_group,N_group))# 4 graphslayout(matrix(seq(4), 2,2))# use polar coordinate to generate sample# theta ~ UNIF(0, 2pi)r = rep(c(1, 3, 6), rep(N_in_group,N_group))theta = runif(N_total)*2*pix[,1] = (r  ) * cos(theta) + rnorm(N_total, sd = .2)x[,2] = (r  ) * sin(theta) + rnorm(N_total, sd = .2)plot(x[,1], x[,2], col = rainbow(3)[label], main = "Origin", xlab="First dimension", ylab="Second dimension")X = xXtX = t(X) %*% Xres = eigen(XtX)V = res$vectorsD = diag(res$values)# verify eigen decop# sum(abs(XtX %*% V - V %*% (D)))Y = X%*% Vplot(Y[,1], Y[,2], col = rainbow(3)[label], main = "Traditional PCA" , xlab="First component", ylab="Second component")# Kernel PCA# Polynomial Kernel# k(x,y) = t(x) %*% y + 1k1 = 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 %*% onesres = eigen(K_norm)V = res$vectorsD = diag(res$values)Y = K %*% Vplot(Y[,1], Y[,2], col = rainbow(3)[label], main = "Kernel PCA (Poly)", xlab="First component", ylab="Second component")# Gaussian Kernel# k(x,y) = exp(-sum((x-y)^2)))k2 = function (x,y) { dnorm(norm(matrix(x-y), type="F"))}K = matrix(0, ncol = N_total, nrow = N_total)for (i in 1:N_total) {  for (j in 1:N_total) {     K[i,j] = k2(X[i,], X[j,])}}ones = 1/N_total* matrix(1, N_total, N_total)K_norm = K - ones %*% K - K %*% ones + ones %*% K %*% onesres = eigen(K_norm)V = res$vectorsD = diag(res$values)Y = K %*% Vplot(Y[,1], Y[,2], col = rainbow(3)[label], main = "Kernel PCA (Gaussian)", xlab="First component", ylab="Second component")