function [A, elapse] = KDA_EVD (K, label) 
% KDA_EVD method 
% written by Guang Dai
% Jan 29th, 2009
% email: daiguang116@gmail.com

% input
% K:        kernel matrix for training data set, each row vector corresponds to a data point.
% label:    class label, a row vector with the same columns as K

% output
% A:        transformation matrix
% elapse:   computational time

% inside
% Sb:       between-class scatter matrix
% Sw:       within-class scatter matrix

% ====== Initialization ======

% size checking
K = K';
clear data;
[nFea, nSmp] = size(K);
if length(label) ~= nSmp
    error('label and data mismatch!');
end

tmp_T = cputime;

% obtain num of labes
classLabel = unique(label);
nClass = length(classLabel);
Dim = nClass - 1;

% check if it is a spars matrix 
if issparse(K)
    K = full(K);
end

% center the data
sampleMean = mean(K, 2);

P = [ones(1, nSmp); K];
P = P';
PP = pinv(P);
clear P;

Y = zeros(nSmp, nClass);
if 1==1
    beta = ones(nSmp,nClass);
else
    beta = ones(nSmp,nClass);
    for i=1:nClass
        label_j = classLabel(i);
        beta(:,label_j) = nSmp/length(find(label==label_j));
    end
end
for i = 1:nClass,
    label_j = classLabel(i);
    index = find(label==label_j);
    Y(index, label_j)= beta(index, label_j);
end
W = PP * Y;
W1= W(2:end,:);
W2= W1';
% Configure Hb
Hb = zeros(nFea, nClass);
for i = 1:nClass,
    label_j = classLabel(i);
    index = find(label==label_j);
    classMean = mean(K(:,index),2);
    sr = sqrt(length(index));
    Hb (:,label_j) = sr*(classMean-sampleMean);
    W2 (label_j,:) = W2(label_j,:)/(sr*beta(1,label_j));
end
B = W2 * Hb;
% EVD ???
[U, S, V] = svd(B);
Q1 = U(:,1:end-1);

A = inv(sqrt(S(1:end-1,1:end-1))) * Q1' * W2;
A = A';
elapse = cputime - tmp_T;