Decentralized Orthogonal Iteration (DOI) + Spectral Clustering

Decentralized Orthogonal Iteration (DOI) + Spectral clustering 

  

On spectral clustering analysis and an algorithm.pdf <= Spectral clustering에 사용된 논문 

A Decentralized Algorithm for Spectral Analysis.pdf   <= Decentralized Orthogonal Iteration에 사용된 논문 

 DOI 는 Distributed algorithm 중에 하나이다. 이 알고리즘을 통해서 전체 Graph G의 Adjacency matrix A의 top k eigenvector들을 구할 수 있다.  이는 Orthogonal Iteration (OI)를 Push-sum (PS) algorithm이라는 consensus 알고리즘을 통해서 distributed 하게 구현한 것이다. 

 이렇게 구한 k 개의 eigenvector들은 전체 Graph의 spectral clustering에 사용될 수 있다. Spectral clustering이란 노드와 엣지로 이뤄진 (원과 연결선) 그래프 구조를 적절히 나누는 알고리즘을 의미한다. 여기선 node와 node 사이의 거리에 반비례하는 weight를 엣지에 부여했지만 이를 페이스북 같은 것에 고려할 때는 두 사용자 간의  친밀도를 나타낼 수 있는 주고 받은 메세지의 수나 태그된 사진의 수 등으로 표현될 수 있을 것이다. 

 

 DOI 와 같은 분산 처리 기법을 사용할 경우 얻을 수 있는 이점은 크게 두 가지이다. 하나는 하나의 성능 좋은 서버가 모든 연산을 할 필요가 없다는 것이다. 즉 페이스북과 같이 사용자가 수천만에 달하는 시스템의 경우 수천만*수천만 의 엄청난 크기의 matrix의 eigenvector를 구해야하는데 이는 실질적으로 불가능하다. 두 번째 이점은 privacy이다. 최근에 프리즘과 같은 CIA의 감시 도청 시스템에 대한 말이 많이 대두되는 시기에 분산 처리를 할 경우 각 노드에 해당하는 사용자의 PC에서 자신의 neighbor (이 경우는 친구)의 정보만을 사용하므로 중앙 서버에서는 따로 개인 정보를 볼 필요가 없다. 

1. Orthogonal Iteration 

매트랩 코드 main_OI.m

%% Orthogonal Iteration Algorithm - Computing top k eigenvectors

clc;

clear all;

close all;

%% Set matrix A

n   = 10;

A   = randn(n, n);

A   = A + A';  % make symmetrc: make evs real

% A   = A*A';  % make psd: make evs positive

%% Calculate eigenvalue and vectors from A 

tic;

[e_vec_A e_val_A] = eig(A); 

fprintf('\n Eigenvalue decomposition ');toc;

e_vec_A    = fliplr(e_vec_A);   % flip eval and evec to make ascending order 

e_val_A    = rot90(e_val_A, 2); 

diag_sv    = abs(diag(e_val_A));

[~, idx]   = sort(diag_sv, 'ascend');

idx(idx')  = linspace(n, 1, n);

%% Orthogonal Iteration 

k    = 4; % the number of eigenvectors 

Q    = randn(n, k);

maxiter = 1E3;

tic;

for iter = 1:maxiter

    V = A*Q;

    Q = qr_GS(V);

end 

fprintf('\n Orthogonal ITeration '); toc;

%% 결과를 확인한다. 

Qsorted(:, idx) = e_vec_A;

Qref            = Qsorted(:, 1:4);

QorthogonalIter = Q;

fprintf('\n Qref: eigenvectors are sorted by its singular value \n'); 

Qref

QorthogonalIter

error_norm = norm(abs(Qref)-abs(Q));

error_max  = max(max(abs(Qref)-abs(Q)));

fprintf('Error norm: %.5f max: %.5f \n', error_norm, error_max); 

 

2. Push Sum 

매트랩 코드 main_PS.m

%% Push-Sum Algorithm computing Sum of all nodes

clc;

clear all;

close all;

%% 먼저 stage에 agent를 뿌리고 neighor를 설정한다.

make_agent_formation; % <= 여기서 stage와 agent변수를 설정

%% 위에서 설정한 agent들로 Adjacency matrix를 만든다.

sig2d = 3; % 거리에 따른 weight의 decay를 결정하는 값.

A     = make_adjacency_matrix(agents, stage, sig2d);

%% Push-Sum Algorithm

% 1. Stochstic matrix B를 만든다.

B = A > 0;

% B = diag(sum(B)) \ B; % <== Make Row-Stochastic

B = B / diag(sum(B));   % <== Make Column-Stochastic

% 2. 각 node마다 임의의 값을 설정한다. 

ref_val = randi([1 100], stage.nr_agent, 5);

for ai = 1:stage.nr_agent

    agents{ai}.S = ref_val(ai, :);

end

% 3. Push-Sum Initialization

for ai =1:stage.nr_agent

    agents{ai}.w = 0;

end

% 4. Push-Sum Implementation

agents{1}.w = 1; % <== 하나만 1로 바꾼다. 

ps.Iter     = 1;

ps.maxIter  = 1E2; 

tic;

while 1

    % fprintf('[%d/%d] \n', ps.Iter, ps.maxIter);

    

    % (4th, 5th line) Push-Sum Algorithm

    for ai = 1:stage.nr_agent

        agents{ai}.Stemp = zeros(size(agents{ai}.S));

        agents{ai}.wtemp = 0; 

        sumB = 0;

        for nj = agents{ai}.nei

            agents{ai}.Stemp = agents{ai}.Stemp ...

                + B(ai, nj)*agents{nj}.S;

            agents{ai}.wtemp = agents{ai}.wtemp ...

                + B(ai, nj)*agents{nj}.w; 

            sumB = sumB + B(ai, nj);

        end 

        sumB

    end

    

    % temp 를 ref로 덮어쓴다.

    for ai = 1:stage.nr_agent 

        agents{ai}.S = agents{ai}.Stemp;

        agents{ai}.w = agents{ai}.wtemp;

        % fprintf(' [%d] w:%.2f s:%.2f', ai, agents{ai}.w, agents{ai}.S);

    end

    % fprintf('\n');

    

    ps.Iter = ps.Iter + 1;

    if ps.Iter >= ps.maxIter

        break;

    end

end

toc;

% (7th Line) After Iteration

for ai = 1:stage.nr_agent

    agents{ai}.S = agents{ai}.S / agents{ai}.w; 

end

%% Push-Sum의 결과를 분석한다. 

for ai = 1:stage.nr_agent

%     fprintf(' [%02d] Psuh-Sum: %.1f Reference: %.1f (%.1f %%) \n' ...

%         , ai, agents{ai}.S, sum(ref_val) ...

%         , 100*(agents{ai}.S - sum(ref_val))/sum(ref_val));

    fprintf(' [%02d] Diff btw/ Psuh-Sum and Reference: %.1f %% \n' ...

        , ai, 100*(agents{ai}.S - sum(ref_val))/sum(ref_val));

end

fprintf('\n');

fprintf('Reference result: '); disp(sum(ref_val));

fprintf('Push-Sum result: '); disp(agents{1}.S );

 

3. Decentralized Orthogonal Iteration + Spectral Clustering 

*Spectral Clustering Algorithm 

매트랩 코드 main_DOI.m

%% Decentralized Orthogonal Iteration: Calculating k eigenvectors

clc;

clear all;

close all;

%% 먼저 stage에 agent를 뿌리고 neighor를 설정한다.

make_agent_formation; % <= 여기서 stage와 agent변수를 설정

%% 위에서 설정한 agent들로 Adjacency matrix를 만든다.

sig2d = 3; % 거리에 따른 weight의 decay를 결정하는 값.

A     = make_adjacency_matrix(agents, stage, sig2d);

%% 1. Calculate eigenvalue and vectors from A

tic;

[e_vec_A e_val_A] = eig(A);

fprintf('\n Eigenvalue decomposition ');toc;

e_vec_A    = fliplr(e_vec_A);   % flip eval and evec to make ascending order

e_val_A    = rot90(e_val_A, 2);

diag_sv    = abs(diag(e_val_A));

[~, idx]   = sort(diag_sv, 'ascend');

idx(idx')  = linspace(stage.nr_agent, 1, stage.nr_agent);

%% 2. Orthogonal Iteration

k    = 10; % the number of eigenvectors

Q    = randn(stage.nr_agent, k);

maxiter = 1E3;

tic;

for iter = 1:maxiter

    V = A*Q;

    Q = qr_GS(V);

end

fprintf('\n Orthogonal ITeration '); toc;

% Orthogonal Iteration의 결과를 확인한다.

Qsorted(:, idx) = e_vec_A;

Qref            = Qsorted(:, 1:k);

QorthogonalIter = Q;

fprintf('\n Qref: eigenvectors are sorted by its singular value \n');

Qref

QorthogonalIter

error_norm = norm(abs(Qref)-abs(QorthogonalIter));

error_max  = max(max(abs(Qref)-abs(QorthogonalIter)));

fprintf('Error norm: %.5f max: %.5f \n', error_norm, error_max);

%% 3. Decentralized Orthogonal Iteration

% 1. 각 agent에 random k-dimensional vector를 설정한다.

for ai = 1:stage.nr_agent

    agents{ai}.Q = randn(1, k);

end

% 2. Push-Sum을 위한 col-stochastic mtx B를 만든다.

B = A > 0;

B = B / diag(sum(B));   % <== Make Column-Stochastic

% 3. DOI를 수행한다.

doi.Iter    = 1;

doi.maxIter = 5E2;

while 1

    fprintf('[%d/%d] \n', doi.Iter, doi.maxIter);

    

    % (DOI-3rd, 4th line)

    for ai = 1:stage.nr_agent

        % (3rd)

        agents{ai}.V = zeros(1, k);

        for nj = agents{ai}.nei

            agents{ai}.V = agents{ai}.V + A(ai, nj)*agents{nj}.Q;

        end

        

        % (4th)

        agents{ai}.K = agents{ai}.V'*agents{ai}.V;

    end

    

    % (DOI-5th line) Push-Sum

    for ai = 1:stage.nr_agent

        agents{ai}.S = agents{ai}.K;

        agents{ai}.w = 0;

    end

    agents{1}.w = 1; % <== 하나만 1로 바꾼다.

    ps.Iter     = 1;

    ps.maxIter  = 5E2;

    while 1

        % fprintf('[%d/%d] \n', ps.Iter, ps.maxIter);

        

        % (PS-4th, 5th line) Push-Sum Algorithm

        for ai = 1:stage.nr_agent

            agents{ai}.Stemp = zeros(size(agents{ai}.S));

            agents{ai}.wtemp = 0;

            sumB = 0;

            for nj = agents{ai}.nei

                agents{ai}.Stemp = agents{ai}.Stemp ...

                    + B(ai, nj)*agents{nj}.S;

                agents{ai}.wtemp = agents{ai}.wtemp ...

                    + B(ai, nj)*agents{nj}.w;

                sumB = sumB + B(ai, nj);

            end

        end

        

        % temp 를 ref로 덮어쓴다.

        for ai = 1:stage.nr_agent

            agents{ai}.S = agents{ai}.Stemp;

            agents{ai}.w = agents{ai}.wtemp;

            % fprintf(' [%d] w:%.2f s:%.2f', ai, agents{ai}.w, agents{ai}.S);

        end

        % fprintf('\n');

        

        % Push-Sum 에 대한 종료 조건

        ps.Iter = ps.Iter + 1;

        if ps.Iter >= ps.maxIter

            break;

        end

    end % 여기서 Push-Sum이 끝났다.

    for ai = 1:stage.nr_agent

        agents{ai}.S = agents{ai}.S / agents{ai}.w;

    end

    

    % (DOI-6th and 7th line)

    for ai = 1:stage.nr_agent

        agents{ai}.R = chol(agents{ai}.S);

        agents{ai}.Q = agents{ai}.V / agents{ai}.R;

    end

    

    % DOI 에 대한 종료 조건

    doi.Iter = doi.Iter + 1;

    if doi.Iter > doi.maxIter

        break;

    end

end

%% DOI 결과를 저장한다. 

Qdoi = [];

for ai = 1:stage.nr_agent

    Qdoi = [Qdoi ; agents{ai}.Q];

end 

save('Qresults.mat', 'Qref', 'Qdoi', 'QorthogonalIter');

%% 최종 결과를 plot한다.

clc;

clear all;

test_index = 1;

load(['MultiAgentFormation' num2str(test_index) '.mat']);

load(['Qresults' num2str(test_index) '.mat']); 

[Qref Qdoi QorthogonalIter];

error_norm = norm(abs(Qref)-abs(QorthogonalIter));

error_max  = max(max(abs(Qref)-abs(QorthogonalIter)));

fprintf('Error norm: %.3f max: %.3f \n', error_norm, error_max);

error_norm = norm(abs(Qref)-abs(Qdoi));

error_max  = max(max(abs(Qref)-abs(Qdoi)));

fprintf('Error norm: %.3f max: %.3f rate: %.2f %% \n' ...

    , error_norm, error_max, 100*error_norm/norm(abs(Qref)) );

figure();

hold on;

for ai = 1:stage.nr_agent

    plot(agents{ai}.pos(1), agents{ai}.pos(2)...

        , 's' ...

        , 'LineWidth', 2 ...

        , 'MarkerEdgeColor','k' ...

        , 'MarkerFaceColor','g' ...

        , 'MarkerSize',10)

    str = sprintf('  [%d]', ai);

    text(agents{ai}.pos(1), agents{ai}.pos(2), str);

    for nj = agents{ai}.nei

        plot([agents{ai}.pos(1) agents{nj}.pos(1)], ...

            [agents{ai}.pos(2) agents{nj}.pos(2)], 'b-');

    end

end

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

grid on; hold off;

title('Multi-Agent Coordination');

%% Spectral Clustering을 한 결과를 보여준다. 

nr_k  = 10;

colors = rand(nr_k, 3);

Qtest = Qref; % Qdoi, Qref

[kmeans_label kmeans_center] ...

    = kmeans(Qtest, nr_k ...

    , 'Distance','cosine' ... % sqEuclidean , cosine

    , 'Replicates', 5); 

fig_kmeans = figure(3);

hold on;

for ai = 1:stage.nr_agent

    for nj = agents{ai}.nei

        plot([agents{ai}.pos(1) agents{nj}.pos(1)], ...

            [agents{ai}.pos(2) agents{nj}.pos(2)], '--' ...

            , 'Color', [0.5 0.6 0.8] ...

            , 'LineWidth', 1 );

    end

end

for ai = 1:stage.nr_agent

    plot(agents{ai}.pos(1), agents{ai}.pos(2)...

        , 's' ...

        , 'LineWidth', 2 ...

        , 'MarkerEdgeColor','k' ...

        , 'MarkerFaceColor', colors(kmeans_label(ai), :) ...

        , 'MarkerSize',10)

    % str = sprintf('  [%d] %d', ai, kmeans_label(ai));

    str = sprintf('   <%d>', kmeans_label(ai));

    text(agents{ai}.pos(1), agents{ai}.pos(2), str);

    

end

hold off;

title('k means result');

 

추가적으로 필요한 코드들 

1. make_agent_formation.m 

%% Multi-agent formation을 만들자. 

make_new_formation = 0;

plot_formation     = 0;

%% Multi-agent formation을 만든다. (혹은 불러온다.)

if make_new_formation

    stage.nr_agent  = 100;    % Agent의 수

    stage.xmin = 0; stage.xmax = 10; stage.w = 10;   

    stage.ymin = 0; stage.ymax = 10; stage.h = 10;

    stage.nei_radius  = 2;

    agents = cell(stage.nr_agent, 1);

    for ai = 1:stage.nr_agent

        agents{ai}.pos = [stage.xmin+stage.w*rand() stage.ymin+stage.h*rand()];

    end

    for ai = 1:stage.nr_agent 

        agents{ai}.nei = [];

        for aj = 1:stage.nr_agent

            if ai == aj, continue; end;

            dist_ij = norm(agents{ai}.pos-agents{aj}.pos);

            if dist_ij <= stage.nei_radius

                agents{ai}.nei = [agents{ai}.nei aj ];

            end

        end

    end

    save('MultiAgentFormation.mat', 'stage', 'agents');

else

    load('MultiAgentFormation.mat');

end

%% agent들과 neighbor들을 plot한다.

if plot_formation

    figure();

    hold on;

    for ai = 1:stage.nr_agent

        for nj = agents{ai}.nei

            plot([agents{ai}.pos(1) agents{nj}.pos(1)], ...

                [agents{ai}.pos(2) agents{nj}.pos(2)], 'b-');

        end

    end

    

    for ai = 1:stage.nr_agent

        plot(agents{ai}.pos(1), agents{ai}.pos(2)...

            , 's' ...

            , 'LineWidth', 2 ...

            , 'MarkerEdgeColor','k' ...

            , 'MarkerFaceColor','g' ...

            , 'MarkerSize',10)

        str = sprintf('  [%d]', ai);

        text(agents{ai}.pos(1), agents{ai}.pos(2), str);    

    end

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

    grid on; hold off;

    title('Multi-Agent Coordination');

end 

2. make_adjacency_matrix.m 

function A = make_adjacency_matrix(agents, stage, sigma_2)

%% 주어진 agent와 stage에 대해서 adjacency matrix를 만든다. 

A       = zeros(stage.nr_agent, stage.nr_agent);

for ai = 1:stage.nr_agent

    for aj = 1:stage.nr_agent

        if ai == aj, continue; end;

        dist_ij = norm(agents{ai}.pos-agents{aj}.pos);

        if dist_ij <= stage.nei_radius

            % 거리가 nei_radius보다 작으면 neighbor이다.

            A(ai, aj) = stage.nei_radius*exp(-dist_ij^2/sigma_2);

        end

    end

end 

3. qr_GS.m 

function [Q,R] = qr_GS(A)

[m,n] = size(A);

% compute QR using Gram-Schmidt

for j = 1:n

   v = A(:,j);

   for i=1:j-1

        R(i,j) = Q(:,i)'*A(:,j);

        v = v - R(i,j)*Q(:,i);

   end

   R(j,j) = norm(v);

   Q(:,j) = v/R(j,j);

end

 DOI + Spectral clustering 결과 

위의 그림에서 같은 색과 숫자가 부여된 노드 (네모) 들이 같은 cluster로 분류된 애들이다. 대충 잘 된다. kmeans cluster의 distance measure는 cosine distance를 사용하였다. 이상하게 유클리디안 distance를 사용하면 결과가 잘 나오지 않는다..