function [eigvector, Sb,St,eigvalue, elapse] = RLDA(gnd, data, lam)
% LDA: Linear Discriminant Analysis
%
%             Input:
%               data  - Data matrix. nSam x nfea
%               gnd   - Label. nSam x 1
%             Output:
%               eigvector - Each column is an embedding function, for a new
%                           data point (row vector) x,  y = x*eigvector
%                           will be the embedding result of x.
%               eigvalue  - The sorted eigvalue of LDA eigen-problem.
%               elapse    - Time spent on different steps
%

tmp_T = cputime;

% ====== Initialization===========================
[nSmp, nFea] = size(data);
if length(gnd) ~= nSmp
    error('gnd and data mismatch!');
end

classLabel = unique(gnd);
nClass = length(classLabel);

sampleMean = mean(data,1);
data = (data - repmat(sampleMean,nSmp,1));
%%=================================================

St = data'*data;
St = max(St,St');  % total scatter matrix St

Hb = zeros(nClass,nFea);
for i = 1 : nClass
    index = find(gnd==classLabel(i));
    classMean = mean(data(index,:),1);
    Hb(i,:) = sqrt(length(index))*classMean;
end
Sb = Hb'*Hb;
Sb = max(Sb,Sb');  % between-class scatter matrix Sb

%-------------------------------------------------------------
%----- St = Sb + Sw
%----- argmax (a'*Sb*a)/(a'*Sw*a) <----> argmax (a'*Sb*a)/(a'*St*a)
%----- St^{-1}*Sb eigvector -----Sb*eigvector=St*eigvector*eigvalue

St = St - Sb;  % Sw
St = St + lam*eye(size(St));


[eigvector, eigvalue] = eig(Sb, St);
eigvalue = diag(eigvalue);
eigIdx = find(eigvalue < 1e-3);
eigvalue(eigIdx) = [];
eigvector(:,eigIdx) = [];

[junk, index] = sort(-eigvalue);
eigvalue = eigvalue(index);
eigvector = eigvector(:,index);

for i = 1:size(eigvector,2)
    eigvector(:,i) = eigvector(:,i)./norm(eigvector(:,i));
end

elapse = cputime - tmp_T;


