Julia言語 Kernel PCA実装

カーネルPCAを実装した．参考にしたのは 赤穂先生のカーネル本とこの修士論文．

「RKHS内でのΦ(x)の平均を簡単のためにゼロにする」という巷の説明をあちこちで見かけるが，無限次元空間でΦ(x)の平均をゼロにする方法は自明じゃない．この修論のappendixがとても参考になった．

https://amzn.to/3zFF2wV

まず，カーネル行列を計算するための関数を書く． （この辺，後々，testデータとtrainデータ間のカーネルも必要になったりするので，地味に書き方が難しい．カーネル法のライブラリってこの辺，どうしてるんだろう）

using LinearAlgebra

function rbf_inners(X,Y,c)
    rbf_inner(x,y,c) = exp( -c*norm(x-y)^2 )
    d, n = size(X)
    d, m = size(Y)
    K = zeros(m,n)
    for j = 1:n
        for i = 1:m
           K[i,j] = rbf_inner(Y[:,i],X[:,j],c) 
        end
    end
    return K
end

function poly_inners(X,Y,c,deg)
    poly_inner(x,y,c,deg) = (c+x'*y)^deg
    d, n = size(X)
    d, m = size(Y)
    K = zeros(m,n)
    for j = 1:n
        for i = 1:m
           K[i,j] = poly_inner(Y[:,i],X[:,j],c,deg) 
        end
    end
    return K
end

function lapras_inners(X,Y,c)
    lapras_inner(x,y,c) = exp(-c*sum( abs.(x .- y)))
    d, n = size(X)
    d, m = size(Y)
    K = zeros(m,n)
    for j = 1:n
        for i = 1:m
           K[i,j] = lapras_inner(Y[:,i],X[:,j],c) 
        end
    end
    return K 
end

function Kernel(X, Y; c=2.0, deg=3.0, kernel="rbf")
    if kernel == "rbf"
        K = rbf_inners(X,Y,c) 
    elseif kernel == "poly" 
        K = poly_inners(X,Y,c,deg)
    elseif kernel == "linear"
        K = poly_inners(X,Y,0,1)
    elseif kernel == "lapras"
        K = lapras_inners(X,Y,c)
    else
        error("Kernel $kernel is not defined")
    end
    return K
end

次に，カーネルPCAの本体を書いていく．KからKbarを得る更新式は，‪Bernhard Schölkopfさんの修論のappendixに書いてある．良い性質のカーネルを選ばないと，Φ空間での特徴量ベクトルの平均が0になっていることはない(と思う)ので，centered_processedは基本的にはfalseにするのが良いと思う．

rankがターゲットランクで，無限次元の特徴量空間につくる有限次元の線形部分空間の次元に対応する．訓練に使ったデータを，Φ空間のrank次元部分空間に射影した先のrank次元の座標がprojected_lowdim_pointsになる．射影した点に興味はなくて，基底だけほしい場合は，projection=falseにする．

function Kernel_PCA(X ; rank=2, kernel="rbf", c=2.0, deg=3.0, centered_processed=false, projection=true)
    K = Kernel(X, X, c=c, deg=deg, kernel=kernel)
    n, n = size(K)
    if centered_processed == true 
        Kbar = K
    else
        one = ones(n,n)
        Kbar = K - 1/n * one * K - 1/n * K * one + 1/(n^2) * one * K * one
    end
    Kbar = Symmetric(Kbar)
    
    f = eigen(Kbar, sortby=-)
    lambda = f.values ./ n
    alpha = f.vectors
    
    projected_lowdim_points = zeros(rank,n)
    if projection
        projected_lowdim_points = n * Diagonal(lambda[1:rank]) * alpha[:,1:rank]' 
    end
    return K, lambda, alpha, projected_lowdim_points
end

途中のSymmetricはなくてもよいはずだが，いれないと，微妙に固有値が複素数になることがあるので，念のため．基本的にはこれで実装完了で，ちょっと実験してみる．とりあえず適当に線形分離できない2次元データを作ってプロットしてみる．

using Plots
X = zeros(2,200)
for j = 1:100
    X[:,j] = [cos(2*pi*j/100.0)+0.1*rand() sin(2*pi*j/100.0)+0.1*rand()]
end
for j = 101:200
    X[:,j] = [0.25*cos(2*pi*j/100.0)+0.1*rand() 0.25*sin(2*pi*j/100.0)+0.1*rand()]
end

plt=scatter(X[1,1:100],X[2,1:100], c="blue", legend=false)
scatter!(plt,X[1,101:200], X[2,101:200], c="red")

もちろん，こんなのは，古典的なPCAじゃ手に負えない．これをKernel-PCAして，Φ空間での2次元空間にプロットしてみる．

K, lambda, alpha, lowdims=Kernel_PCA(X, rank=2, kernel="rbf", c=2.0);
plt=scatter(lowdims[1,1:100],lowdims[2,1:100], c="blue", legend=false)
scatter!(plt,lowdims[1,101:200], lowdims[2,101:200], c="red")

素晴らしい！！！ きちんと線形分離できとる！！

データXに加えて，データXtestが観測された場合を考える． 素朴に考えれば，XとXtestをconcateして，同じことを繰り返せばよいが，それは手間である．固有値を求めるのは計算コストが高いので同じことを繰り返したくない．．．そこで，やはり彼の修論のappendixを参照にして，fitさせる関数をつくる．

function fit_Kernel_PCA(K, Xtrain, Xtest, alpha; rank, kernel="rbf", c=2.0, deg=3.0, centered_processed=false)
    Ktest = Kernel(Xtrain, Xtest, c=c, deg=deg, kernel=kernel)
    m, l = size(Ktest)
    if centered_processed 
        Ktest_bar = Ktest
    else
        one_ml = ones(m,l)
        one_ll = ones(l,l)
        Ktest_bar = Ktest - 1/l * Ktest*one_ll - 1/l * one_ml * K + 1/(l^2) * one_ml * K * one_ll
    end
    projected_lowdim_points = (Ktest_bar * alpha[:,1:rank])' 
    return projected_lowdim_points
end

これを使うと，固有値計算を繰り返す必要はない（テストサンプルとトレーニングサンプル間のカーネルを計算する必要はあるが，これは，固有値計算に比べればずっと軽い）．実験してみよう．×印が，新たに観測されたデータで，〇印が，Φでrank次元線形部分空間をつくるのに使ったデータ．

Xtest = zeros(2,100)
for j = 1:50
    Xtest[:,j] = [cos(2*pi*j/50.0)+0.1*rand() sin(2*pi*j/50.0)+0.1*rand()]
end
for j = 51:100
    Xtest[:,j] = [0.25*cos(2*pi*j/50.0)+0.1*rand() 0.25*sin(2*pi*j/50.0)+0.1*rand()]
end

plt=scatter(X[1,1:100],X[2,1:100], c="blue", legend=false)
scatter!(plt,X[1,101:200], X[2,101:200], c="red")
scatter!(plt,Xtest[1,1:50], Xtest[2,1:50], c="blue", markershape=:cross)
scatter!(plt,Xtest[1,51:100], Xtest[2,51:100], c="red", markershape=:cross)

lowdims_test =fit_Kernel_PCA(K, X, Xtest, alpha; rank=2, kernel="rbf", c=2.0, deg=3.0, centered_processed=false)
plt=scatter(lowdims[1,1:100],lowdims[2,1:100], c="blue", legend=false)
scatter!(plt,lowdims[1,101:200], lowdims[2,101:200], c="red")
scatter!(plt,lowdims_test[1,1:50],lowdims_test[2,1:50], c="blue", markershape=:cross)
scatter!(plt,lowdims_test[1,51:100],lowdims_test[2,51:100], c="red", markershape=:cross)

正しく実装できてる．