Kernel Principal Component Analysis (KPCA)

어떤 데이터들이 들어왔을 때, 이들을 좀 더 분류하기 쉽게 변환, 혹은 처리하기 쉽게 차원을 낮춰주는 역활을 하는 대표적인 방법은 Principal Component Analysis (PCA) 이다. 이번 포스팅은 PCA가 아니라 이 PCA에 kernel trick을 적용한 Kernel PCA란 것이 있다. 

자세한 설명은 아래 slide를 통해서 알 수 있다. 

포항대에서 만든 KPCA자료 (잘 설명되어있다. )

 이것이 무엇인고 하니, 기존의 PCA에 kernel trick을 적용한 것이다. 

 그렇다면 Kernel trick은 무엇인고 하니, 조금 어려운 것을 의미한다. 수학적으로 증명을 하기 위해서 Reproducing Kernel Hilbert Space (RKHS) 까지 나오지만 않다면 Kernel Trick 자체는 어려운 것이 아니다. 

 우리가 마음껏 Affine 연산을 하는 유클리디언 공간에서는 어떤 데이터 간의 관계를 내적으로 표현할 수 있다. 이 내적의 의미는 두 데이터 사이의 유사도를 의미한다. 두 벡터의 내적이란 각 벡터의 element들의 곱의 합이다. 이 때 각 element가 비슷할 수록 (특히 부호) 내적의 값은 크게 나온다. 

 이러한 유사도의 개념을 좀 더 범용적으로 적용하기 위해서 kernel function이 등장한다. 이름만 봐서는 매우 어려워 보이지만, kernel function은 단순히 두 개의 벡터를 입력으로 가지는 함수이다. 그리고 일반적으로 두 벡터 사이의 거리 (보통, 두 벡터의 차의 유클리디안 norm) 에 반비례한다. 기존 formulation의 내적 부분을 이 kernel function을 통해 만든 kernel matrix로 대채하는 것이 바로 Kernel trick이다. 

 그 의미를 조금 들여다 보면, 우리의 데이터가 살고 있는 n 차원의 유클리디언 공간을 매우 차원이 큰 공간으로 mapping하는 nonlinear function이 있다고 가정하는 것이다. 그리고, 그 공간에서 내적은 우리가 정의한 kernel function의 출력값이 된다고 가정하는 것이다. 이러한 가정을 증명할 때 필요한 개념이 앞의 RKHS이다. 

  활용적인 면에서 Kernel을 적용할 경우 가장 큰 장점은 주어진 두 데이터의 상관 관계를 우리가 정한 kernel function을 통해 나타낼 수 있다는 점이다. Kernel trick이 가장 활발히 사용되는 곳은 nonlinear SVM으로 해보면 알겠지만 classification 성능이 비약적으로 상승한다. 이유는 간단한데 기존의 linear SVM은 decision boundary를 칼로 자르는 것과 같이 나누지만, kernel SVM은 학습 데이터에 맞는 곡선으로 나누기 때문이다. 더 flexible하다고나 할까. 

 이번 포스팅의 주제인 Kernel PCA 역시 기존의 PCA에 Kernel trick을 적용한 것이다. PCA를 unsupervised learning에 적용할 경우, 주어진 데이터를 우리가 정한 basis에 projection시키고, 각 basis들에 project된 값이 주어진 데이터의 feature가 되는 구조이다. 이것에 kernel trick을 적용하였기 때문에 기존의 linear projection이 nonlinear한 projection으로 바뀌게 된다. 

그리고 그 결과는 다음과 같다. 

 위의 그림은 6개의 feature를 뽑아서 그린 것이고, 아래 그림은 상위 세 개의 feature를 R, G, B에 대응시킨 것이다. 

  재밌는 점은 원래 데이터의 dimension이 2였는데, kernel PCA를 통해선 그 dimension을 데이터의 갯수에 해당하는 N까지 오히려 올릴 수 있다는 것이다. 이러한 경향은 nonlinear SVM에서도 나타난다. 

데이터 (꼭 필요는 없음. 새롭게 생성할 수 있음.)

kpca_stage.mat

main_kpca.m

clc;

clear all;

close all;

%% 

opt.make_new_data = 0;

opt.plot_data     = 0;

if opt.make_new_data

    kpca.d = 2;     % 데이터의 dimension

    kpca.k = 3;     % 클러스터의 수

    kpca.N = 100;   % 한 클러스터의 데이터의 수

    kpca.Total = kpca.k*kpca.N;

    % stage를 만든다.

    kpca.stage.xmin = 0;

    kpca.stage.xmax = 5;

    kpca.stage.xsize = kpca.stage.xmax-kpca.stage.xmin;

    kpca.stage.xres = 20;

    kpca.stage.ymin = 0;

    kpca.stage.ymax = 5;

    kpca.stage.ysize = kpca.stage.ymax-kpca.stage.ymin;

    kpca.stage.yres = 20;

    kpca.stage.xgrid = linspace(kpca.stage.xmin, kpca.stage.xmax, kpca.stage.xres);

    kpca.stage.ygrid = linspace(kpca.stage.ymin, kpca.stage.ymax, kpca.stage.yres);

    [mesh_X, mesh_Y] = meshgrid(kpca.stage.xgrid, kpca.stage.ygrid);

    kpca.stage.points = [mesh_X(:) mesh_Y(:)];

    clear('mesh_X', 'mesh_Y');

    kpca.stage.nr_points = length(kpca.stage.points);

    % stage 안에 데이터를 만든다.

    kpca.data = [];

    for i = 1:kpca.k

        % center를 만들고, 

        margin = 1;

        kpca.cluster(i).center ...

            = [kpca.stage.xmin+margin + (kpca.stage.xsize-2*margin)*rand() ...

            , kpca.stage.ymin+margin + (kpca.stage.ysize-2*margin)*rand() ];

        % center주변으로 데이터를 뿌린다.

        concentrate_factor = 20; % 이 값이 클수록 데이터가 center에 몰린다.

        kpca.cluster(i).data = repmat(kpca.cluster(i).center, kpca.N, 1) ...

            + [kpca.stage.xsize/concentrate_factor*randn(kpca.N, 1) ...

            kpca.stage.ysize/concentrate_factor*randn(kpca.N, 1)];

        kpca.data = [kpca.data ; kpca.cluster(i).data];

        % etc

        kpca.cluster(i).color = rand(1,3);

    end

    save('kpca_stage.mat', 'kpca');

else

    load kpca_stage.mat;

end

% 데이터들을 그려보자.

if opt.plot_data

    hold on;

    for i = 1:kpca.k

        plot(kpca.cluster(i).data(:,1), kpca.cluster(i).data(:,2)...

            ,'o', 'Color', kpca.cluster(i).color);

    end

    hold off;

    axis([kpca.stage.xmin kpca.stage.xmax kpca.stage.ymin kpca.stage.ymax]);

end 

%% kernel pca 를 한다. 

kpca.opt_kernel.useSE_K    = 1;

kpca.opt_kernel.sigma2f    = 1;

kpca.opt_kernel.sigma2w    = 0;

kpca.opt_kernel.sigma2x    = 5E-1;

kpca.K = kpca_kernel(kpca.data, kpca.data, kpca.opt_kernel);

unit_N = 1/kpca.Total*ones(kpca.Total, kpca.Total);

kpca.Ktilt = kpca.K - unit_N*kpca.K - kpca.K*unit_N + unit_N*kpca.K*unit_N;

% Kernel의 eigenvalue를 구하자.

[kpca.evecs, kpca.evals] = eig(kpca.Ktilt);

kpca.evals = real(diag(kpca.evals));

for i = 1:kpca.Total,

    kpca.evecs(:,i) = kpca.evecs(:,i)/(sqrt(kpca.evals(i)));

end

% Option. eig가 아니라 svd로 문제를 푼다.

% [kpca.evecs, kpca.evals] = svd(kpca.Ktilt);

%% grid의 각 점에 대해서 test를 한다. 

kpca.Ktest = kpca_kernel(kpca.stage.points, kpca.data, kpca.opt_kernel);

unit_Ntest = 1/kpca.Total*ones(kpca.stage.nr_points, kpca.Total);

kpca.Ktest_tilt = kpca.Ktest - unit_Ntest*kpca.K - kpca.Ktest*unit_N + unit_Ntest*kpca.K*unit_N;

% feature를 뽑는다.

kpca.nr_features = 6;

kpca.test_features = kpca.Ktest_tilt * kpca.evecs(:,1:kpca.nr_features);

%% plot한다.

% 1. 총 6개의 feature에 대해서 어떻게 분류되나 확인 

figure( 'Position' , [100 200 800 1000]);

for n = 1:kpca.nr_features,

    subplot(3, 2, n);

    

    axis([kpca.stage.xmin kpca.stage.xmax kpca.stage.ymin kpca.stage.ymax]);

    kpca.imag = reshape(kpca.test_features(:,n), kpca.stage.yres, kpca.stage.xres);

    axis('xy');

    colormap(gray);

    hold on;

    pcolor(kpca.stage.xgrid, kpca.stage.ygrid, kpca.imag);

    shading interp

    contour(kpca.stage.xgrid, kpca.stage.ygrid, kpca.imag, 9, 'b');

    plot(kpca.data(:,1), kpca.data(:,2), 'r.')

    %title(sprintf('Eigenvalue=%4.3f', kpca.evals(n))); 

    title(sprintf(' %d th eigenvalue = %4.3f ', n, kpca.evals(n) )); 

    hold off; colorbar; 

    

end

% 2. 세 개의 가장 큰 eigenvector로 분류된 feature의 값을 차례로 RGB에 넣고 plot 

figure( 'Position' , [1000 200 1200 1000]);

rgb_field = zeros(kpca.stage.yres, kpca.stage.xres, 3);

for n = 1:3

    kpca.imag = reshape(kpca.test_features(:,n), kpca.stage.yres, kpca.stage.xres);

    rgb_field(:, :, n) = kpca.imag;

end

min_rgb = min(min(min(rgb_field)));

max_rgb = max(max(max(rgb_field)));

rgb_field = (rgb_field - min_rgb)/(max_rgb-min_rgb);

subplot(221);

imagesc(kpca.stage.xgrid, kpca.stage.ygrid, rgb_field); colorbar;

hold on;plot(kpca.data(:,1), kpca.data(:,2), 'k.');hold off;

set(gca,'YDir','normal'); axis('xy');

title('RGB values by 3 features');

countour_step = 30;

subplot(222);

pcolor(kpca.stage.xgrid, kpca.stage.ygrid, rgb_field(:,:,1));shading interp; colorbar;

hold on;plot(kpca.data(:,1), kpca.data(:,2), 'k.');

contour(kpca.stage.xgrid, kpca.stage.ygrid, rgb_field(:,:,1), countour_step, 'Color', [.5 .5 .5]);hold off;

title('1st eigenvector');

subplot(223);

pcolor(kpca.stage.xgrid, kpca.stage.ygrid, rgb_field(:,:,2));shading interp; colorbar;

hold on;plot(kpca.data(:,1), kpca.data(:,2), 'k.');

contour(kpca.stage.xgrid, kpca.stage.ygrid, rgb_field(:,:,2), countour_step, 'Color', [.5 .5 .5]);hold off;

title('2nd eigenvector');

subplot(224);

pcolor(kpca.stage.xgrid, kpca.stage.ygrid, rgb_field(:,:,3));shading interp; colorbar;

hold on;plot(kpca.data(:,1), kpca.data(:,2), 'k.');

contour(kpca.stage.xgrid, kpca.stage.ygrid, rgb_field(:,:,3), countour_step, 'Color', [.5 .5 .5]);hold off;

title('3rd eigenvector');

kpca_kernel.m

function K = kpca_kernel(X1, X2, option)

%

% X1: 데이터는 col방향으로 들어있다.

% X2: 데이터는 col방향으로 들어있다.

%

% 여기에는 sigma2w는 포함되어있지 않다!! 

%

% <SE kernel 예제 1>

% option_gp.useSE_K    = 1;

% option_gp.sigma2f    = 1;

% option_gp.sigma2w    = 1E-3;

% option_gp.sigma2x    = 5E-1;

% option_gp.Grad = 0; <== 미분을 안한다. 할 경우는 일단 scalar

% 

% <D kernel 예제 2>

% option_gp.useD_K     = 1;

% option_gp.sigma2f    = 1;

% option_gp.sigma2w    = 1E-3;

% option_gp.r          = 0.5;

% option_gp.Grad = 0;

%

% 

%

if ~isfield(option, 'Grad')

    option.Grad = 0;

end

if ~isfield(option, 'useSE_K')

    option.useSE_K = 0;

end

if ~isfield(option, 'useD_K')

    option.useD_K = 0;

end

len_X1  = size(X1, 1);

len_X2  = size(X2, 1);

K       = zeros(len_X1, len_X2);

for i = 1:len_X1

    for j  = 1:len_X2

        

        if option.Grad

            % X1에서 미분한 gradient를 구한다. (이 때는 앞에서 정의한 K를 사용하지 않을 것이다.)

            % 미분을 일반적으로 하나의 X1과 X2를 가정한다. (예: 두 위치, len_X1 = len_X2 = 1)

            if option.useSE_K

                K = se_kernel(X1(i,:), X2(j,:), option)*-(X1(i,:)-X2(i,:))/option.sigma2x;

            elseif option.useD_K

                diff  = X1(i,:) - X2(j,:);

                dist  = norm(diff);

                sig2f = option.sigma2f;

                r     = option.r;

                if dist == 0

                    K   = zeros(size(diff));

                elseif dist > r

                    K   = zeros(size(diff));

                else

                    K   = sig2f*pi*(dist/r-1)*sin(pi*dist/r)*diff/dist/r;                    

                end

            end

        else

            % kernel의 값을 구한다. 

            if option.useSE_K

                K_data = se_kernel(X1(i,:), X2(j,:), option);

            elseif option.useD_K

                K_data = distribute_kernel(X1(i,:), X2(j,:), option);

            end

            K(i, j) = K_data; 

        end

    end

end

end

%%

function output = se_kernel(x1, x2, option)

    output = option.sigma2f*exp(-norm(x1-x2)^2/option.sigma2x/2);

end

%%

function output = distribute_kernel(x1, x2, option)

    output = option.sigma2f*Lambda(norm(x1-x2)/option.r);

end

function output_lambda = Lambda(h)

    if h <= 1 

        output_lambda = (1-h)*cos(pi*h) + 1/pi*sin(pi*h);

    else 

        output_lambda = 0;

    end

end