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.