Gaussian Process Approximations (PITC, D2FAS)

 Gaussian process 는 non-parametric Bayesian regression의 대표적인 알고리즘으로 주어진 데이터를 가지고 새로운 위치에서 값을 regression할 때, 사용될 수 있는 알고리즘이다. 

 수학적인 얘기는 집어 치우고, 결국 이 알고리즘이 할 수 있는 것은 sparse한 (띄엄띄엄있는 의 pedantic 한 단어) 기상청에서 온도를 쟀을 때, 서울 전체에 대한 온도 지도를 만들 수 있다. 이 때 기상청의 위치들의 모음을 X라 하고, 해당하는 온도를 Y라 할 때, 새로운 위치 Z에서 온도를 구할 수 있는 알고리즘이 Gaussian process regression이다. 

 수학적으로 조금 들어가보면, Y = f(X)의 관계를 만족시키는 f를 구하는 것인데, f(Z) ~ N(muZ, sig2Z)에서 muZ와 sig2Z를 구하는데, 이를 위해서 (X, Z)에 대한 joint Gaussian distribution을 정의한다. 이렇게 하면, Gauss느님에 정의에 따라 f(Z)|Y에 대한 conditional Gaussian distribution을 자동으로 구할 수 있다. 

 '$f(Z)\sim N(\mu_Z, var_Z)$'

이고, 

 '$mu_Z = k(Z, X)(k(X, X)+\sigma^2_wI)^{-1}Y$'

 '$var_Z = z(Z, Z)+k(Z, X)(k(X, X)+\sigma^2_wI)^{-1}k(X, Z)$'

로 정의된다. 

여기서 문제가 '$(k(X, X)+\sigma^2_wI)^{-1}$'를 구하기 위해서 들어가는 computational cost가 '$O(N^3)$'이라는 점이다. '$N$'은 데이터의 수 (기상청의 수) 이다. 

 이러한 문제를 해결하기 위해서 여러 Approximation 기법들이 고안되었다. 이들 대부분은 SoD (Subset of Data) 기법으로 불릴 수 있는데, support set '$U\in R^m$'를 사용하는 알고리즘이다. ('$m < N$') 

 이번 포스팅에서 설명할 내용은 Partially Independent Training Conditional (PITC) 과 이 알고리즘을 분산화 시킨 Decentralized Data Fusion and Active Sensing (D2FAS) 이다. 

 자세한 내용은 논문을 참조 

"Decentralized Data Fusion and Active Sensing with Mobile Sensors for Modeling and Predicting Spatiotemporal Traffic Phenomena"

Abstract

 The problem of modeling and predicting spatiotemporal traffic phenomena over an urban road network is important to many traffic applications such as detecting and forecasting congestion hotspots. This paper presents a decentralized data fusion and active sensing (D2FAS) algorithm for mobile sensors to actively explore the road network to gather and assimilate the most informative data for predicting the traffic phenomenon. We analyze the time and communication complexity of D2FAS and demonstrate that it can scale well with a large number of observations and sensors. We provide a theoretical guarantee on its predictive performance to be equivalent to that of a sophisticated centralized sparse approximation for the Gaussian process (GP) model: The computation of such a sparse approximate GP model can thus be parallelized and distributed among the mobile sensors (in a Google-like MapReduce paradigm), thereby achieving efficient and scalable prediction. We also theoretically guarantee its active sensing performance that improves under various practical environmental conditions. Empirical evaluation on real-world urban road network data shows that our D2FAS algorithm is significantly more time-efficient and scalable than state-of-the-art centralized algorithms while achieving comparable predictive performance.

 

알고리즘 설명은 귀찮고, 그냥 논문에 있는 중요한 Definition과 Theorem을 올리자. 

일반적인 GPR, PITC, 그리고 D2FAS를 적용한 결과가 다음과 같다. 

결과

[GPR with 1000 input data 500 groups 10 support set] 

GP takes 2.361sec / GP-PITC takes 0.106sec (22.221 times faster) 

GP takes 2.361sec / GP-D2FAS takes 0.089sec (26.592 times faster) 

GP error: 0.009 / SGP error: 0.231 

코드

1. 메인 코드 (demo_GPR_PITC_D2FAS.m)

%% 1D GP와 SGP를 한다. GP의 경우 hyperparameter opt도 한다. 

% clc;

% clear all;

close all;

plot_ref = 0;      % Reference를 그린다. 

plot_learning = 1; % Hyperparameter Learning curve를 그린다. 

%% 중요 파라미터들 

nrData  = 1E3;    % GP 데이터의 수 

nrGroup = 5E0;    % 몇 개의 group으로 나눌지 

nrU     = 1E1;    % Suppot set의 크기 

nrZ     = 3E1;    % Output의 크기 

%% Reference function을 만든다. 

xmin = 1;

xmax = 7;

% f = @(x)(3+cos(0.02*exp(x)));

f = @(x)(3+cos(3*x));

Xref = linspace(xmin, xmax, 500)';

Yref = f(Xref); 

if plot_ref == 1

    figure(1);

    plot(Xref, Yref, 'r-' ...

        ,'LineWidth', 2);

    title('Reference graph');

end

%% GP Data로 사용될 데이터를 모은다.  

Xgp = linspace(xmin, xmax, nrData)';

Ygp = f(Xgp);

% Regression을 할 위치들 

Z = linspace(xmin, xmax, nrZ)';

%% Xgp와 Ygp를 nrGroup 만큼 나눈다. (PITC, D2FAS) 

% Xgp 를 k개로 나눈다. (PITC, D2FAS) 

X_PITC.k    = nrGroup; % Number of Groups 

X_PITC.N    = nrData;  % Total Number of Data

X_PITC.all  = Xgp;     % Total Data 

[X_PITC.data, X_PITC.idxs] = split_vec(Xgp, X_PITC.k); 

X_D2FAS = X_PITC;

% Ygp 를 k개로 나눈다. (D2FAS)

Y_D2FAS.k    = nrGroup; % Number of Groups 

Y_D2FAS.N    = nrData;  % Total Number of Data

Y_D2FAS.all  = Ygp;     % Total Data 

[Y_D2FAS.data, Y_D2FAS.idxs] = split_vec(Ygp, Y_D2FAS.k); 

%% Support set를 구한다. 

U   = linspace(xmin, xmax, nrU)'; % Total support set

%% GP의 hyp의 초기값을 설정한다. 

hyp.sig2f = 1E+0;

hyp.sig2w = 1E-2;

hyp.sigx  = 5E-1; 

%% hyperparameter를 구한다. 

% maxIter = 2E2;

maxIter = 0; % hyperparameter optimization은 하지 않는다. (일단)

if maxIter ~= 0 

    % Hyperparameter estimation을 위한 step size를 정한다. 

    stepSize.sig2f = 1E-3;

    stepSize.sig2w = 1E-3;

    stepSize.sigx  = 1E-3;

    

    logLiks = zeros(maxIter, 1); % Learning curve를 위한 변수 

    fprintf('Start optimizing the hyperparameters. \n');

    for i = 1:maxIter 

        if rem(i, maxIter/10) == 0, 

            fprintf('%d%% ', floor(i/maxIter*100)); 

        end;

        [hyp_der, logLik] = se_kernel_logderivative(Xgp, Ygp, hyp);

        hyp.sig2f = hyp.sig2f + stepSize.sig2f * hyp_der.sig2f;

        hyp.sig2w = hyp.sig2w + stepSize.sig2w * hyp_der.sig2w;

        hyp.sigx  = hyp.sigx  + stepSize.sigx  * hyp_der.sigx;

        % log likelihood 

        logLiks(i) = logLik;

    end

    fprintf(' Done \n');

    % Learning curve를 그려본다. 

    if plot_learning

        figure(2);

        plot(logLiks, 'b-');

        xlabel('Iterations');

        ylabel('log marginal likelihood');

        title('Learning curve of hyperparameters of GP');

    end 

    save('GPlearning.mat', 'hyp');

else

    % load('GPlearning.mat');

end

%% GPR, GPR_PITC, GPR_D2FAS를 수행한다. 

tic;

[mu_gp, var_gp] = gpr(Z, Xgp, Ygp, hyp);

gp_toc = toc;

%% PITC

tic;

[mu_pitc, var_pitc] = gpr_pitc(Z, X_PITC, Ygp, U, hyp);   

gp_pitc_toc = toc;

%% D2FAS

tic;

[mu_d2fas, var_d2fas] = gpr_d2fas(Z, X_D2FAS, Y_D2FAS, U, hyp);

gp_d2fas_toc = toc;

% mu_d2fas = mu_pitc; var_d2fas = var_pitc;

%% 그림을 그린다. 

fig_gp_all = figure('Position', [100 100 1400 800]);

hold on;

hdata  = plot(Xgp, Ygp, 'ro');

href   = plot(Xref, Yref, 'r-', 'LineWidth', 2);

hgpr   = plot(Z, mu_gp, 'g*-', 'LineWidth', 2);

plot(Z, mu_gp+var_gp, 'g--');

plot(Z, mu_gp-var_gp, 'g--');

hpitc  = plot(Z, mu_pitc, 'bo-', 'LineWidth', 2);

plot(Z, mu_pitc+var_pitc, 'b--');

plot(Z, mu_pitc-var_pitc, 'b--');

hd2fas  = plot(Z, mu_d2fas, 'c+-', 'LineWidth', 2);

plot(Z, mu_d2fas+abs(var_d2fas), 'c--');

plot(Z, mu_d2fas-abs(var_d2fas), 'c--');

hu = plot(U, zeros(length(U), 1) ...

    , 'kx', 'MarkerSize', 15);

hold off;

title(sprintf('GPR with %d input data %d groups %d support set'...

    , nrData, nrGroup, nrU) ...

    , 'FontSize', 20);

legend([href hgpr hpitc hd2fas hu hdata] ...

    , 'Reference', 'GPR'...

    , 'GPR-PITC', 'GPR-D2FAS' ...

    , 'Support Set', 'GP Data'...

    , -1);

axis([xmin xmax min(Yref)-2 max(Yref)+2]);

% 그림을 파일로 저장을 해보자. 

str_name = sprintf('results/%d_input_%d_groups_%d_supports' ...

    , nrData, nrGroup, nrU);

set(fig_gp_all,'PaperPositionMode','auto')

print (fig_gp_all , '-dpng', [str_name, '.png']) ;

%%

gpErr = norm(mu_gp - f(Z));

sgpErr = norm(mu_pitc - f(Z));

str1 = sprintf('[GPR with %d input data %d groups %d support set] \n'...

    , nrData, nrGroup, nrU);

str2 = sprintf('GP takes %.3fsec / GP-PITC takes %.3fsec (%.3f times faster) \n' ...

    , gp_toc, gp_pitc_toc, gp_toc/gp_pitc_toc );

str3 = sprintf('GP takes %.3fsec / GP-D2FAS takes %.3fsec (%.3f times faster) \n' ...

    , gp_toc, gp_d2fas_toc, gp_toc/gp_d2fas_toc);

str4 = sprintf('GP error: %.3f / SGP error: %.3f \n' ...

    , gpErr, sgpErr);

strAll = [str1, str2, str3, str4];

fileID = fopen([ str_name , '.txt'], 'w');

fprintf(fileID, '%s \r\n', strAll);

fclose(fileID);

fprintf('%s \r\n', strAll);

2. PITC 코드 (gpr_pitc.m)

function [muhatZ, varhatZ] = gpr_pitc(Z, Xpitc, Y, U, hyp)

% 

% partially independent training conditional (PITC) Algorithm

% 

% Z: unobserved points

% Xpitc: Input data structure

% Y: Output data

% U: Support set

% hyp: Hyperparameters

%

nrInput = Xpitc.N;

KZZ = se_kernel(Z, Z, hyp);

KZU = se_kernel(Z, U, hyp);

KUZ = KZU';

KUU = se_kernel(U, U, hyp);

KXU = se_kernel(Xpitc.all, U, hyp);

KUX = KXU';

PHI = KUU \ KUX; 

LAMBDA = zeros(nrInput, nrInput);

for i = 1:Xpitc.k

    if i == 1, frIdx = 1;

    else frIdx = Xpitc.idxs(i-1)+1; 

    end

    toIdx = Xpitc.idxs(i);

    LAMBDA(frIdx:toIdx, frIdx:toIdx) ...

        = se_kernel(Xpitc.data{i}, Xpitc.data{i}, hyp);

end

PSI = KZU*PHI/(KXU*PHI + LAMBDA);

% 평균을 미리 구한다. 

muY = repmat(mean(Y), size(Y, 1), 1);

muZ = repmat(mean(Y), size(Z, 1), 1);

% muhat과 varhat을 구한다. 

muhatZ  = muZ + PSI*(Y-muY);

varhatZ = KZZ - PSI*KXU/KUU*KUZ;

varhatZ = diag(varhatZ); 

3. D2FAS 코드 (gpr_d2fas.m)

function [muhatZ, varhatZ] = gpr_d2fas(Z, X_D2FAS, Y_D2FAS, U, hyp)

% 

% Distributed Data Fusion Algorithm (D2FAS) Algorithm

% 

% Z: unobserved points

% Xd2fas: Input data structure

% Y: Output data

% U: Support set

% hyp: Hyperparameters

%

% 0. Precompute the inverse of KUU

KUU = se_kernel(U, U, hyp);

invKUU = inv(KUU);

%% 1. Make Local Summary 

localSummary.zu = cell(X_D2FAS.k, 1);

localSummary.sigu = cell(X_D2FAS.k, 1);

for i = 1:X_D2FAS.k

    % For each agent in a group

    Xk = X_D2FAS.data{i}; 

    Yk = Y_D2FAS.data{i};

    

    % muYk = mean(Yk)*ones(size(Yk));

    muYk = mean(Y_D2FAS.all)*ones(size(Yk));

    KUXk = se_kernel(U, Xk, hyp);

    KXkU = KUXk';

    KXkXk = se_kernel(Xk, Xk, hyp);

    % SDkDk_U = (KXkXk - KXkU/KUU*KUXk);

    % Gamma = KUXk/SDkDk_U;

    Gamma = KUXk/KXkXk;

    

    % Local summary of z

    localSummary.zu{i} = Gamma*(Yk-muYk);

    localSummary.sigu{i} = Gamma*KXkU;

end

%% 2. Make Global Summary 

nrSupportSet = size(U, 1);

zu_bar = zeros(nrSupportSet, 1);

sigu_bar = zeros(nrSupportSet, nrSupportSet);

for i = 1:X_D2FAS.k

    zu_bar = zu_bar + localSummary.zu{i};

    sigu_bar = sigu_bar + localSummary.sigu{i};

end

sigu_bar = sigu_bar + KUU;

%% 3. Make regression using D2FAS 

KZU = se_kernel(Z, U, hyp);

KUZ = KZU';

KZZ = se_kernel(Z, Z, hyp);

muZ = mean(Y_D2FAS.all)*ones(length(Z), 1);

% estimated mean

muhatZ  = muZ + KZU/sigu_bar*zu_bar;

% estimated variance

varhatZ = KZZ - KZU*(invKUU + inv(sigu_bar))*KUZ;

varhatZ = diag(varhatZ);

4. GPR 코드 (gpr.m)

function [muhatZ, varhatZ] = gpr(Z, X, Y, hyp)

%

% Gaussian Process Regression 

% 

% Z: unobserved points

% X: Input data

% Y: Output data

% hyp: Hyperparameters

%

KXX = se_kernel(X, X, hyp);

KZX = se_kernel(Z, X, hyp);

KZZ = se_kernel(Z, Z, hyp);

muY = repmat(mean(Y), size(Y, 1), 1);

muZ = repmat(mean(Y), size(Z, 1), 1);

muhatZ = muZ + KZX/KXX*(Y-muY);

varhatZ = KZZ - KZX/KXX*KZX';

varhatZ = diag(varhatZ); 

5. 커널 관련 코드

a. se_kernel.m

function K = se_kernel(X1, X2, hyp)

sig2f = hyp.sig2f;

sig2w = hyp.sig2w;

sigx = hyp.sigx;

nr_X1 = size(X1, 1);

nr_X2 = size(X2, 1);

K = zeros(nr_X1, nr_X2);

for i = 1:nr_X1

    for j = 1:nr_X2

        curX1 = X1(i, :);

        curX2 = X2(j, :);

        dist = norm(curX1-curX2);

        

        val = sig2f*exp(-1/2*(dist/sigx)^2 );

        

        % if isequal(curX1, curX2)

        if nr_X1 == nr_X2 && i == j

            val = val + sig2w;

        end 

        

        K(i,j) = val;

    end

end

b. se_kernel_logderivative.m

function [hyp__derivative, logLik] = se_kernel_logderivative(X, Y, hyp)

%% 1. 먼저 K를 구한다. 

K = se_kernel(X, X, hyp);

%% 2. K\y를 계산한다. 

invK = inv(K);

invKY = K \ Y;

invKYsq = invKY*invKY';

%% 3. 각 log partial derivative를 구한다. 

sig2f = hyp.sig2f;

sig2w = hyp.sig2w;

sigx  = hyp.sigx;

nr_X = size(X, 1);

dK_dSig2f = zeros(nr_X, nr_X);

dK_dSig2w = zeros(nr_X, nr_X);

dK_dSigx  = zeros(nr_X, nr_X);

for i = 1:nr_X

    for j = 1:nr_X

        curXi = X(i, :);

        curXj = X(j, :);

        dij = norm(curXi - curXj);

        

        % 3.1. sig2f에 대해서 구하자. 

        dK_dSig2f(i, j) = 2*sqrt(sig2f)*exp(-0.5*(dij/sigx)^2);

        

        % 3.2. sigx에 대해서 구하자. 

        dK_dSigx(i, j)  = sig2f*exp(-0.5*(dij/sigx)^2)*dij^2/sigx^3;

        

        % 3.3. sig2w에 대해서 구하자. 

        if i == j

            dK_dSig2w(i, j) = 2*sqrt(sig2w);

        else

            dK_dSig2w(i, j) = 0;

        end

    end

end 

%% 4. log derivative를 구한다. 

hyp__derivative.sig2f = 0.5*trace( (invKYsq - invK) * dK_dSig2f);

hyp__derivative.sig2w = 0.5*trace( (invKYsq - invK) * dK_dSig2w);

hyp__derivative.sigx  = 0.5*trace( (invKYsq - invK) * dK_dSigx);

%% 5. Log likelihood를 구한다. 

logLik = -0.5*Y'*invKY - 0.5*log(det(K)) - nr_X/2*log(2*pi);

6. 추가 유틸 코드

a. split_vec.m

function [divided_vec, divided_idxs] = split_vec(ref, k)

% ref의 vector를 k 개로 나눠서 ret에 matrix로 저장한다. (데이터는 row)

nr_data  = length(ref);

unit_n   = (nr_data / k);

% Assert conditions

assert(unit_n > 1, '[split_vec] Too many groups!');

% assert(mod(unit_n, 1) == 0, '[split_vec] Groups cannot be divided ');

% Start algorithm  

unit_n = floor(unit_n);

divided_vec = cell(k, 1);

divided_idxs = zeros(k, 1);

curr_idx = 1;

for i = 1:k-1

    next_idx = curr_idx + unit_n - 1;

    divided_vec{i}   = ref(curr_idx : next_idx, :);

    curr_idx = next_idx + 1;

    

    divided_idxs(i) = next_idx;

end

% 마지막을 채워넣는다. 

divided_vec{k} = ref(curr_idx : end, :);

divided_idxs(k) = nr_data;

b. varycolor.m

function ColorSet=varycolor(NumberOfPlots)

error(nargchk(1,1,nargin))%correct number of input arguements??

error(nargoutchk(0, 1, nargout))%correct number of output arguements??

%Take care of the anomolies

if NumberOfPlots<1

    ColorSet=[];

elseif NumberOfPlots==1

    ColorSet=[0 1 0];

elseif NumberOfPlots==2

    ColorSet=[0 1 0; 0 1 1];

elseif NumberOfPlots==3

    ColorSet=[0 1 0; 0 1 1; 0 0 1];

elseif NumberOfPlots==4

    ColorSet=[0 1 0; 0 1 1; 0 0 1; 1 0 1];

elseif NumberOfPlots==5

    ColorSet=[0 1 0; 0 1 1; 0 0 1; 1 0 1; 1 0 0];

elseif NumberOfPlots==6

    ColorSet=[0 1 0; 0 1 1; 0 0 1; 1 0 1; 1 0 0; 0 0 0];

else %default and where this function has an actual advantage

    %we have 5 segments to distribute the plots

    EachSec=floor(NumberOfPlots/5); 

    

    %how many extra lines are there? 

    ExtraPlots=mod(NumberOfPlots,5); 

    

    %initialize our vector

    ColorSet=zeros(NumberOfPlots,3);

    

    %This is to deal with the extra plots that don't fit nicely into the

    %segments

    Adjust=zeros(1,5);

    for m=1:ExtraPlots

        Adjust(m)=1;

    end

    

    SecOne   =EachSec+Adjust(1);

    SecTwo   =EachSec+Adjust(2);

    SecThree =EachSec+Adjust(3);

    SecFour  =EachSec+Adjust(4);

    SecFive  =EachSec;

    for m=1:SecOne

        ColorSet(m,:)=[0 1 (m-1)/(SecOne-1)];

    end

    for m=1:SecTwo

        ColorSet(m+SecOne,:)=[0 (SecTwo-m)/(SecTwo) 1];

    end

    

    for m=1:SecThree

        ColorSet(m+SecOne+SecTwo,:)=[(m)/(SecThree) 0 1];

    end

    

    for m=1:SecFour

        ColorSet(m+SecOne+SecTwo+SecThree,:)=[1 0 (SecFour-m)/(SecFour)];

    end

    for m=1:SecFive

        ColorSet(m+SecOne+SecTwo+SecThree+SecFour,:)=[(SecFive-m)/(SecFive) 0 0];

    end

    

end