转自:http://zhanxw.com/blog/2011/02/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,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,我们的协方差矩阵为:
这里有一个陷阱,我跳进去过:
在对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):
这里定义 αai=xTiuλa 。
因为 xTiu 是一个标量(scala),所以 αai 也是一个标量,因此 ui 是可以由 xi 张成。
进而我们显示PCA投影可以用内积运算表示,例如我们把 xi 向任意一个主成分分量 ua 进行投影,得到的是 uTaxi ,也就是 xTiua 。作者猜测写成这种形式是为了能抽出 xTixj=<xi,xj> 的内积形式。
当我们定义 Kij=xTixj 时,上式可以写为 K2α=NλaKαa
(这里 αa 定义为 [αa1,αa2,…,αaN]T .)
进一步,我们得到解为:
K 矩阵包含特征值 λ~ 和 αa ,我们可以通过 α 可以计算得到 ua ,
注意特征值分解时Eigendecomposition, αa 只代表一个方向,它的长度一般为1,但在此处不为1。
这里计算出 αa 的长度(下面将要用到):
因为 ua 的长度是1,我们有:
在上面的分析过程中,我们只使用了内积。因此当我们把 Kij=xTixj 推广为 Kij=<Φ(xi),Φ(xj>=Φ(xi)TΦ(xj) 时,上面的分析结果并不会改变。
1.4 如何在主成分方向上投影?
投影时,只需要使用
U
矩阵,假设我们得到的新数据为
t
,那么
t
在
ua
方向的投影是:
对于高维空间的数据 Φ(xi),Φ(t) ,我们可以用Kernel trick,用 K(xi,t) 来带入上式:
1.5 如何Centering 高维空间的数据?
在我们的分析中,协方差矩阵的定义需要centered data。在高维空间中,显式的将
Φ(xi)
居中并不简单,
因为我们并不知道
Φ
的显示表达。但从上面两节可以看出,所有的计算只和
K
矩阵有关。具体计算如下:
令
Φi=Φ(xi)
,居中
ΦCi=Φi–1N∑kΦk
不难看出,
其中 1N 为 N×N 的矩阵,其中每一个元素都是 1/N
对于新的数据,我们同样可以
2. 演示 (R code)
首先我们应该注意输入数据的格式,一般在统计中,我们要求
X
矩阵是
N×d
的,但在我们的推导中,
X
矩阵是
d×N
。
这与统计上的概念并不矛盾:在前面的定义下协方差矩阵为
XTX
,而在后面的定义中是
XXT
。另外这里的协方差矩阵是样本(Sample)的协方差矩阵,我们的认为大写的X代表矩阵,而不是代表一个随机变量。
另外,在下面的结果中,Gaussian 核函数(kernel function)的标准差(sd)为2。在其他取值条件下,所得到的图像是不同的。
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