
function [G, L, DD] = iakdaInit(K, gnd, opts)

% function [G, L, DD] = iakdaInit(K, gnd, opts)
%
% Accelerated kernel discriminant analysis (AKDA)
%
% Description:
%  Accelerated computation of expansion coefficient matrix G and
%  eigenvalue diagonal DD for dimensionality reduction
%
% IN: 
%   K: N x N kernel Gram matrix
%   gnd: class labels
%   opts : additional options
%     reg_alpha: regularization parameter for kernel matrix
%
% OUT:
%  G : N x D projection matrix
%  L : N x N Cholesky factor
%  DD: eigenvalue diagonal matrix
%
% Author: Nikolaos Gkalelis - CERTH-ITI
% Email: gkalelis@iti.gr
%
% First creation (Version 0.0) 1 Apr 2017.
%

fprintf('akdaInit>>> Entering\n');

reguAlpha = opts.regAlpha;

Clbl = unique(gnd); % class labels
if ~issorted(Clbl);
    fprintf('>>> Class labels are not sorted, sorting...\n');
    Clbl = sort(Clbl);
end

C = length(Clbl); % number of classes
N = size(K,1); % number of observations

% compute number of observations of i-th class
Ni = zeros(C,1);
for i=1:C
    Ni(i) = sum(gnd == i);
end

% construct elementary matrices
n_sqrt = zeros(C,1); % square root of cardinality vector
for i=1:C
    n_sqrt(i) = sqrt(Ni(i));
end

% construct core matrix
O = eye(C) -  n_sqrt* n_sqrt.' * (1/N); % (n_sqrt.' * n_sqrt) = N

N2 = zeros(C,C); % inverse of square root of cardinality matrix
for i=1:C
    N2(i,i) = 1/sqrt(Ni(i));
end

Rc = zeros(N,C); % indicator matrix
for n=1:N
    for i=1:C % get class ID
        if gnd(n) == Clbl(i)
            Rc(n,i) = 1;
            break;
        end
    end
end

opts_da.disp=0;
[KSI, DD] = eigs(full(O), rank(O), 'LM', opts_da); % rank(O) = C-1

TH = Rc* N2* KSI;

% regularize Gram matrix
K = max(K,K');
K = K + eye(N, N) * reguAlpha;

L = chol(K); % compute Cholesky factor
W = L\(L'\TH);

% normalize
tmpNorm = sqrt(sum((W'*K).*W',2));
G = W./repmat(tmpNorm',size(W,1),1);
DD = diag(DD);

fprintf('akdaInit>>> Exiting\n');

