Kernel PCA 原理和演示

04-26

主成份（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 = [{x_1},{x_2}, cdots ,{x_N}]{ m{ }}$ 是一個 $d imes N$ 矩陣，代表輸入的數據有 $N$ 個，每個sample的維數是 $d$ 。我們做降維，就是想用 $k$ 維的數據來表示原始的 $d$ 維數據( $kleq d$ )。

當我們使用centered的數據（即 $sum olimits_i {{x_i} } = 0$ ）時，可定義協方差矩陣 $C$ 為：

$[C = frac{1}{N}sumlimits_i {{x_i}x_i^T} = frac{1}{N}X{X^T}]$

做特徵值分解，我們可以得到：

$[CU = ULambda Rightarrow C = ULambda {U^T} = sumlimits_a {{lambda _a}{u_a}u_a^T} ]$

注意這裡的 $C$ , $U$ , $Lambda$ 的維數都是 $d imes d$ , 且 $[U = [{u_1},{u_2}, cdots ,{u_d}]]$ , $Lambda = diag({lambda _1},{lambda _2}, cdots ,{lambda _d})$ 。

當我們做降維時，可以利用前 $k$ 個特徵向量 $[U_k = [{u_1},{u_2}, cdots ,{u_k}]]$ 。則將一個 $d$ 維的 $x_i$ 向 $k$ 維的主成分的方向投影后的數據為 $y_i=U_k^{T}x_i$ （這裡的每一個 $u_i$ 都是 $d$ 維的，代表是一個投影方向，且 $u_i^{T}u_i=1$ ，表示這是一個旋轉變數）

1.2 在高維空間里的PCA應該如何做？

高維空間中，我們定義一個映射 $Phi :{X^d} o mathcal{F}$ ，這裡 $mathcal{F}$ 表示Hilbert泛函空間。

現在我們的輸入數據是 $Phi(x_i), i=1,2,cdots,n$ ，他們的維數可以是無窮維的（泛函空間），為了方便討論，設其空間的維數為 $D$ ( $Dgg d$ )。

在這個新的空間(特徵空間)中，將中心化的數據記為 $ilde Phi(x_i)$ ，即

$ilde Phi(x_i)=Phi(x_i)-ar Phi$

其中

$ar Phi = frac{1}{N}sum_i^{N}{Phi(x_i)}$

將中心化的核矩陣記為 $ilde K$ (具體求法見*****)，且定義 $ilde K_{ij}=< ilde Phi(x_i), ilde Phi(x_j)>= ilde Phi(x_i)^T ilde Phi(x_j)$

$ilde Phi(x_i)$ 中心化的數據（ $sum_i^N{ ildePhi(x_i)}=0$ ）的協方差矩陣為 $ar C$ ：

$ar C = frac{1}{N}sumlimits_i { ilde Phi ({x_i}) ilde Phi {{({x_i})}^T}} = frac{1}{N} ilde Phi (X) ilde Phi {(X)^T}$

這裡有一個陷阱：

在對Kernel trick一知半解的時候，我們常常從形式上認為 $ar{C}$ 可以用 $ilde K_{i,j}= ilde K(x_i,x_j)$ 來代替，

因此對 $ilde K=( ilde K_{ij})$ 做特徵值分解，然後得到 $ilde K = U Lambda U^T$ ，並且對原有數據降維的時候，定義 $Y_i = U_k^{T}X_i$ 。

但這個錯誤的方法有兩個問題：一是我們不知道矩陣 $ar{C}$ 的維數，也即不知道 $ar C$ 的具體形式；二是 $U_k^{T}X_i$ 從形式上看不出是從高維空間的 $Phi(x_i)$ 的投影，並且當有新的數據時 $X_{new}$ ，我們無法從理論上理解 $U_k^T X_{new}$ 是從高維空間的投影，後面（11111）發現，投影方向不是這裡的 $U$ 。

如果應用這種錯誤的方法，我們有可能得到看起來差不多正確的結果，但本質上這是錯誤的。

正確的方法是通過Kernel trick將PCA投影的過程通過內積的形式表達出來，詳細見1.3

1.3 如何用Kernel Trick在高維空間做PCA？

在1.1節中，通過PCA，我們得到了 $U$ 矩陣。這裡將介紹如何僅利用內積的概念來計算傳統的PCA。

首先我們證明U可以由 $x_1,x_2,cdots,x_N$ 展開（span)：

由於 $Cu_a=lambda_a u_a$ 得

$[{u_a} = frac{1}{{{lambda _aN}}}C{u_a}{kern 1pt} = frac{1}{{{lambda _aN}}}(sumlimits_i {{x_i}x_i^T} ){u_a} = frac{1}{{{lambda _aN}}}sumlimits_i {{x_i}(x_i^T{u_a})} {kern 1pt} = frac{1}{{{lambda _aN}}}sumlimits_i {(x_i^T{u_a}){x_i}} = sumlimits_i {frac{{x_i^T{u_a}}}{{{Nlambda _a}}}{x_i}} {kern 1pt} = sumlimits_i {alpha _i^a{x_i}} ]$

這裡定義 ${alpha _i^a}={frac{{x_i^T{u_a}}}{{{Nlambda _a}}}}$ 。

因為 $x_i^{T}u$ 是一個標量(scala)，所以 $alpha_i^a$ 也是一個標量，因此 $u_a$ 是可以由 $x_i$ 張成。

進而我們顯示PCA投影可以用內積運算表示，例如我們把 $x_i$ 向任意一個主成分分量 $u_a$ 進行投影，得到的是 $u_a^Tx_i$ ，也就是 $x^Tu_a$ 。作者猜測寫成這種形式是為了能抽出 $x_i^Tx_j=<x_i,x_j>$ 的內積形式。

$[x_i^TC{u_a} = {lambda _a}x_i^T{u_a}]$

$x_i^Tfrac{1}{N}sumlimits_j {{x_j}} x_j^Tsumlimits_k {alpha _k^a{x_k}} = {lambda _a}x_i^Tsumlimits_k {alpha _k^a{x_k}}$

$sumlimits_k {alpha _k^asumlimits_j {(x_i^T{x_j})(x_j^T{x_k})} } = N{lambda _a}sumlimits_k {alpha _k^a(x_i^T{x_k})}$

也就是說，對於 $i$ ，

$egin{split}&alpha_1^a(x_i^Tx_1x_1^Tx_1+x_i^Tx_2x_2^Tx_1+cdots+x_i^Tx_Nx_N^Tx_1)+\&alpha_2^a(x_i^Tx_1x_1^Tx_2+x_i^Tx_2x_2^Tx_2+cdots+x_i^Tx_Nx_N^Tx_2)+cdots\&alpha_N^a(x_i^Tx_1x_1^Tx_N+x_i^Tx_2x_2^Tx_N+cdots+x_i^Tx_Nx_N^Tx_N)+\&=\&Nlambda_a(alpha_1^ax_i^Tx_1+x_i^Tx_2+cdots+x_i^Tx_N)end{split}$

將其寫成矩陣的形式為：

$left(egin{array}{cccc}x_1^Tx_1 & x_1^Tx_2 & cdots & x_1^Tx_N\x_2^Tx_1 & x_2^Tx_2 & cdots & x_2^Tx_N \vdots & vdots & ddots & vdots \x_N^Tx_1 & x_N^Tx_2& cdots & x_N^Tx_Nend{array} ight)left(egin{array}{cccc}x_1^Tx_1 & x_1^Tx_2 & cdots & x_1^Tx_N\x_2^Tx_1 & x_2^Tx_2 & cdots & x_2^Tx_N \vdots & vdots & ddots & vdots \x_N^Tx_1 & x_N^Tx_2& cdots & x_N^Tx_Nend{array} ight)left(egin{array}{c}alpha_1^a \alpha_2^2 \vdots \alpha_N^a end{array} ight)=Nlambda_aleft(egin{array}{cccc}x_1^Tx_1 & x_1^Tx_2 & cdots & x_1^Tx_N\x_2^Tx_1 & x_2^Tx_2 & cdots & x_2^Tx_N \vdots & vdots & ddots & vdots \x_N^Tx_1 & x_N^Tx_2& cdots & x_N^Tx_Nend{array} ight)left(egin{array}{c}alpha_1^a \alpha_2^2 \vdots \alpha_N^a end{array} ight)$

當我們定義 $K_{ij}=x_i^Tx_j$ 時，上式可以寫為 $K^2alpha^a=Nlambda_aKalpha^a$ （這裡 $alpha^a$ 定義為 $[alpha_1^a,alpha_2^a,cdots,alpha_N^a]^T$ .）

進一步，我們得到解為：

$Kalpha^a= ilde lambda_aalpha^a$ ，其中 $ilde lambda_a=Nlambda_a$

$K$ 矩陣包含特徵值 $ilde lambda_a$ 和特徵向量 $alpha^a$ ，我們可以通過 $alpha^a$ 可以計算得到特徵向量 $u_a$ 。

注意特徵值分解(Eigen decomposition)時， $alpha^a$ 只代表一個方向，它的長度(模長)一般為1，但在此處不為1。下面計算出 $alpha^a$ 的長度（下面將要用到）：

因為 $u_a$ 的長度是1，我們有：

$egin{split}1&=u_a^tu_a\&=(sum_{i}{alpha_i^ax_i})^T(sum_{j}{alpha_j^ax_j}) \&=sum_isum_j{(alpha_i^a)^Talpha_j^ax_i^Tx_j}\&=(alpha^a)^TKalpha^a\&=(alpha^a)^T(Nlambda_aalpha^a)\&=Nlambda_a(alpha^a)^Talpha^aend{split}$

所以

$left|| alpha^a ight|| = frac{1}{sqrt{Nlambda_a} } = frac{1}{sqrt{ ilde{lambda}_a} }$

在上面的分析過程中，我們使用了內積的思想完成了傳統PCA的過程，即可以通過核矩陣 $K$ 的特徵值 $ilde lambda_a$ 和特徵向量 $alpha^a$ 計算協方差矩陣的特徵值 $lambda_a=frac{ ilde lambda_a}{N}$ 和特徵向量 $u_a=sum_i^N{frac{alpha _i^a}{sqrt{ ilde lambda_a}}{x_i}}$

（PS：不要問為什麼不直接用協方差矩陣進行特徵值分解(SVD)，得到想要的主方向。請記住在特徵空間的協方差矩陣式是不能得到的，所以不能進行大俠想的特徵分解）。

1.4 KPCA實現

下面我們重複1.3過程，看看KPCA中的有關推導。

設 $ar C$ 的特徵值和對應的特徵向量分別為 $lambda_k$ 和 $v_k$ ，即

$ar Cv_k=lambda_k v_k(k=1,2,cdots,D)(*)$

由1.3知， $v_kin spanleft{ ilde Phi(x_1),cdots, ilde Phi(x_N) ight}$ ，即存在 $alpha_{k1},cdots,alpha_{kN}$ ，使得：

$v_k=sum_i^N{alpha_{ki} ilde Phi(x_i)}$

用 $ilde Phi(x_j)^T$ 同時作用於 $(*)$ 得

$lambda_k ilde Phi(x_j)^T v_k = ilde Phi(x_j)^T ar C v_k$

上式左邊可以化為：

$egin{split}lambda_k ilde Phi(x_j)^T v_k&=lambda_kleft( sum_i^N{alpha_{ki} ilde Phi(x_j)^T ilde Phi(x_i)} ight)\&=lambda_kleft( sum_i^N{alpha_{ki} ilde K_{ji}} ight) \&=lambda_kleft( ilde K alpha_k ight) _jend{split}$

其中

$ilde K = left( ilde K_{ij} ight)_{N imes N}=left( ilde Phi(x_i)^T ilde Phi(x_i) ight) _{N imes N}$ , $alpha_k=left( alpha_{k1},cdots,alpha_{kN} ight)^T$

右邊可以化為：

$egin{split} ilde Phi(x_j)^T ar C v_k&= ilde Phi(x_j)^T frac{1}{N}sum_l^N{ ilde Phi(x_l) ilde Phi(x_l)^T}sum_i^N{alpha_{ki} ilde Phi(x_i)}\&=frac{1}{N}sum_l^N{sum_i^N{alpha_{ki}left( ilde Phi(x_j)^T ilde Phi(x_l) ilde Phi(x_l)^T ilde Phi(x_i) ight) }}\&=frac{1}{N}sum_l^N{sum_i^N{alpha_{ki} ilde K_{jl} ilde K_{li}}}\&=frac{1}{N}left( ilde K^2alpha_k ight)_j end{split}$

讓 $j$ 依次取遍 $1,cdots,N$ 的所有數，得到上面推到的「左邊=右邊」的矩陣表示：

$lambda_k ilde K alpha_k=frac{1}{N} ilde K ^2 alpha_k$

化簡得到

$ilde K alpha_k = N lambda_k alpha_k = ilde lambda_k alpha_k$

即之前假設的 $alpha_k=left( alpha_{k1},cdots,alpha_{kN} ight)^T$ 為中心化的核矩陣 $ilde K$ 的第 $k$ 個特徵向量（對應的特徵值為 $ilde lambda_k$ ）。現在我們可以用所求的 $ilde lambda_k$ 和 $alpha_k$ 去表示協方差矩陣 $ar C$ 的特徵值 $lambda_k$ 和向量 $v_k$ ：

$lambda_k=frac{ ilde lambda_k}{N}$ ， $v_k=sum_i^N{alpha_{ki} ilde Phi(x_i)}= ilde Phi alpha_k$ ，其中 $ilde Phi=left( ilde Phi(x_1),cdots, ilde Phi(x_N) ight)$

現在還有個小問題沒有解決，在傳統PCA方法中，協方差矩陣得到的特徵向量是單位正交的。同樣，我們希望 $v_k$ 也是單位向量，即：

$egin{split}1&=v_k^Tv_k\&=left(sum_{i}^N{alpha_{ki} ilde Phi(x_i)} ight)^Tleft(sum_{j}^N{alpha_{kj}} ilde Phi(x_j) ight) \&=alpha_k^T ilde Phi^T ilde Phi alpha_k\&=alpha_k^T ilde K alpha_k\&=alpha_k^T ilde lambda_k alpha_k\&= ilde lambda_k alpha_k^T alpha_k\&= ilde lambda_k left| {alpha_k} ight|^2\end{split}$

因此，將 $ilde K$ 的特徵向量 $alpha_k$ 的長度歸一化到 $left| {alpha_k} ight| = frac{1}{sqrt{ ilde lambda_k}}$ ，便能保證 $v_k$ 的長度為1。

下面證明正交性：

$egin{split}v_i^Tv_j &=alpha_i^T ilde Phi^T ilde Phi alpha_j\&=alpha_i^T ilde K alpha_j\&=alpha_i^T ilde lambda_j alpha_j\&= ilde lambda_j alpha_i^T alpha_j\&=0end{split}$

得證！

最後總結下之前的結論，特徵空間中協方差矩陣 $ar C$ 的特徵值為 $lambda_k=frac{ ilde lambda_k}{N}$ ，對應的特徵向量為 $v_k = frac{1}{sqrt{ ilde lambda_k}} ilde Phi alpha_k$ , $(k=1,2,cdots,D)$ ，其中 $ilde lambda_k$ 和 $alpha_k$ 為 $ilde K$ 的特徵值和特徵向量。

由於映射 $Phi$ 未知，直接在特徵空間內進行數據中心化不可行，直接將輸入數據在特徵空間中心化的過程見1.6。

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

傳統PCA投影時，只需要使用 $U$ 矩陣，假設我們得到的新數據為 $T=left( t_1,t_2,cdots,t_M ight)$ ，那麼 $t_k$ 在 $u_a$ 方向的投影是：

$u_a^Tt_k=sum_i{alpha_i^ax_i^Tt_k}=sum_i{alpha_i^a(x_i^Tt_k)}$

對於高維空間的數據 $ilde Phi(x_i)$ , $ilde phi(t_k)$ ，我們可以用Kernel trick，用 $ilde Psi (x_i,t_k)= ilde Phi(x_i)^T ilde phi(t_k)$ 來帶入上式：

$egin{split}v_k^T ilde phi(t_k)&=frac{1}{sqrt{ ilde lambda_k}}left( sum_i^N{alpha_{ki} ilde Phi(x_i)} ight)^T ilde phi (t_k)\&=frac{1}{sqrt{ ilde lambda_k}}sum_i^N{alpha_{ki}left( ilde Phi(x_i)^T ilde phi(t_k) ight)}\&=sum_i^N{frac{alpha_{ki}}{sqrt{ ilde lambda_k}} ilde Psi(x_i,t_k)}\end{split}$

一般 $ilde phi(t_k)$ 不能顯式得到，從而需要計算輸入數據 $T$ 的中心化核矩陣，記為 $ilde Psi$ 。

1.5 如何Centering 高維空間的數據？

在我們的分析中，協方差矩陣的定義需要centered data。在高維空間中，顯式的將 $Phi(x_i)$ 中心化並不簡單, 因為我們並不知道 $Phi$ 的顯示錶達。但從上面兩節可以看出，所有的計算只和 $K$ 矩陣有關。具體計算如下：

令 $Phi_i=Phi(x_i)$ ，居中 $ilde Phi_i=Phi_i-frac{1}{N}sum_kPhi_k$ ：

$egin{split} ilde K_{ij}&=< ilde Phi_i, ilde Phi_j>\&=(Phi_i-frac{1}{N}sum_k{Phi_k})^T(Phi_j-frac{1}{N}sum_l{Phi_l})\&=Phi_i^TPhi_j-frac{1}{N}sum_l{Phi_i^TPhi_l}-frac{1}{N}sum_k{Phi_k^TPhi_j}+frac{1}{N^2}sum_ksum_l{Phi_k^TPhi_l}\&=K_{ij}-frac{1}{N}sum_l{K_{il}}-frac{1}{N}sum_k{K_{kj}}+frac{1}{N^2}sum_ksum_l{K_{kl}}end{split}$

不難看出，

$ilde K=K-1_NK-K1_N+1_NK1_N$

其中 $1_N$ 為 $N imes N$ 的矩陣，其中每一個元素都是 $frac{1}{N}$ 。

對於新的數據 $t_k$ ，如果數據 $t_k otin X$ ： $egin{split} ilde Psi_{i,k}&=< ilde Phi_i, ilde phi_k>\&=(Phi_i-frac{1}{N}sum_k{Phi_k})^T(phi_k-frac{1}{M}sum_l{phi_l})\&=Phi_i^Tphi_t-frac{1}{M}sum_l{Phi_i^Tphi_l}-frac{1}{N}sum_k{Phi_k^TPhi_t}+frac{1}{MN}sum_ksum_l{Phi_k^Tphi_l}\&= Psi(x_i,t_k)-frac{1}{N}sum_l{Psi_{il}}-frac{1}{N}sum_k{Psi(x_k,t)}+frac{1}{MN}sum_ksum_l{Psi_{kl}}end{split}$

所以：

$ilde Psi=Psi-1_NK-K1_M+1_NK1_M$

如果數據 $t_k in X$ ，則上式中有 $M=N$ ，且 $ilde K = ilde Psi$ 。

2. 演示 (R code)

首先我們應該注意輸入數據的格式，一般在統計中，我們要求 $X$ 矩陣是 $N imes d$ 的，但在我們的推導中， $X$ 矩陣是 $d imes N$ 。

這與統計上的概念並不矛盾：在前面的定義下協方差矩陣為 $X^TX$ ，而在後面的定義中是 $XX^T$ 。另外這裡的協方差矩陣是樣本（Sample）的協方差矩陣，我們的認為大寫的 $X$ 代表矩陣，而不是代表一個隨機變數。

另外，在下面的結果中，Gaussian 核函數（kernel function）的標準差（sd）為2。在其他取值條件下，所得到的圖像是不同的。

KPCA圖片：

R 源代碼（Source Code）：鏈接到完整的代碼 KernelPCA

Kernel PCA部分代碼：

# 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) rank = 0for (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="x", ylab="y")

實例：

function demo_kpca_ushape%% PARAMETERSN_split = [30, 30, 60]; % number of data points in 2 straight parts and curve of "U"corners = [3 4 2 4]; % corners: x1 x2 y1 y2nvar = 0.05; % noise variancekernel.type = gauss;kernel.par = .5;numeig = 4; % number of eigenvalues to plotNtest = [50,50]; % number of test grid divisions, in each dimensionborder = [0 6 0 6]; % test grid border %% PROGRAMtic%% generate dataN = sum(N_split);N1 = N_split(1); N2 = N_split(2); N3 = N_split(3);X1 = [linspace(corners(1),corners(2),N1) repmat(corners(3),N1,1)]; % lower straight partX2 = [linspace(corners(1),corners(2),N2) repmat(corners(4),N2,1)]; % upper straight partangles = linspace(pi/2,3*pi/2,N3);rad = (corners(4)-corners(3))/2;m3 = [corners(1) corners(3)+rad];X3 = repmat(m3,N3,1) + rad*[cos(angles) sin(angles)];n = nvar*randn(N,2);X = [X1; X2; X3] + n;%% generate test grid dataNt1 = Ntest(1); Nt2 = Ntest(2);Xtest = zeros(Nt1*Nt2,2);absc = linspace(border(1),border(2),Nt1);ordi = linspace(border(3),border(4),Nt2);for i=1:Nt1, for j=1:Nt2, Xtest((i-1)*Nt2+j,:) = [absc(i) ordi(j)]; endend%% calculate kernel principal components and projections of Xtest[E,v] = km_kpca(X,numeig,kernel.type,kernel.par);Kt = km_kernel(X,Xtest,kernel.type,kernel.par); % kernels of test data setXtestp = E*Kt; % projections of test data set on the principal directionsY = cell(numeig,1);for i=1:numeig, Y{i} = reshape(Xtestp(i,:),Nt2,Nt1); % shape into 2D gridendtoc%% OUTPUTsubplot(3,2,1)plot(X(:,1),X(:,2),.)% axis(border)title(original data of shape U)% fireworks!for i=1:numeig, subplot(2,3,i+2); hold on [C,h] = contourf(-interp2(Y{i},2),25); plot(X(:,1)/border(2)*(Nt1*4-1)+1,X(:,2)/border(4)*(Nt2*4-1)+1,o,... MarkerFaceColor,White,MarkerEdgeColor,Black) title([projections to the ,num2str(i), th principal directions]) set(gca, XTick, []) set(gca, YTick, [])end%Additional functionfunction [E,v] = km_kpca(X,m,ktype,kpar)n = size(X,1);K = km_kernel(X,X,ktype,kpar);[E,V] = eig(K);v = diag(V); % eigenvalues[v,ind] = sort(v,descend);v = v(1:m);E = E(:,ind(1:m)); % principal componentsfor i=1:m E(:,i) = E(:,i)/sqrt(n*v(i)); % normalizationendfunction K = km_kernel(X1,X2,ktype,kpar)switch ktype case gauss % Gaussian kernel sgm = kpar; % kernel width dim1 = size(X1,1); dim2 = size(X2,1); norms1 = sum(X1.^2,2); norms2 = sum(X2.^2,2); mat1 = repmat(norms1,1,dim2); mat2 = repmat(norms2,dim1,1); distmat = mat1 + mat2 - 2*X1*X2; % full distance matrix K = exp(-distmat/(2*sgm^2)); case gauss-diag % only diagonal of Gaussian kernel sgm = kpar; % kernel width K = exp(-sum((X1-X2).^2,2)/(2*sgm^2)); case poly % polynomial kernel p = kpar(1); % polynome order c = kpar(2); % additive constant K = (X1*X2 + c).^p; case linear % linear kernel K = X1*X2; otherwise % default case error (unknown kernel type)end

結果

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

轉載（糾正部分錯誤）http://zhanxw.com/blog/2011/02/kernel-pca-原理和演示/