function [eigvector, eigvalue] = SLDA(gnd,data)
% SLDA: Subspace Linear Discriminant Analysis
% the generalized optimal discriminant vectors should be selected from 
% the complementary subspace of the null space of St

% PCA transformation is employed to discard the null space of the 
% within-class scatter matrix Sw or Sb, following which LDA is applied.

% ====== Initialization
[nSmp,nFea] = size(data);
classLabel = unique(gnd);
nClass = length(classLabel);


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

%======================================
if nSmp > nFea
    ddata = data'*data;
    ddata = max(ddata,ddata');

    [eigvector_PCA, eigvalue_PCA] = eig(ddata);
    eigvalue_PCA = diag(eigvalue_PCA);
    clear ddata;

    maxEigValue = max(abs(eigvalue_PCA));
    eigIdx = find(eigvalue_PCA/maxEigValue < 1e-12);
    eigvalue_PCA(eigIdx) = [];
    eigvector_PCA(:,eigIdx) = [];

    [junk, index] = sort(-eigvalue_PCA);
    eigvalue_PCA = eigvalue_PCA(index);
    eigvector_PCA = eigvector_PCA(:, index);

    eigvalue_PCA = eigvalue_PCA.^-.5;
    data = (data*eigvector_PCA).*repmat(eigvalue_PCA',nSmp,1);
else
    ddata = data*data';
    ddata = max(ddata,ddata');

    [eigvector, eigvalue_PCA] = eig(ddata);
    eigvalue_PCA = diag(eigvalue_PCA);
    clear ddata;

    maxEigValue = max(eigvalue_PCA);
    eigIdx = find(eigvalue_PCA/maxEigValue < 1e-12);
    eigvalue_PCA(eigIdx) = [];
    eigvector(:,eigIdx) = [];

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

    eigvalue_PCA = eigvalue_PCA.^-.5;
    eigvector_PCA = (data'*eigvector).*repmat(eigvalue_PCA',nFea,1);

    data = eigvector;
    clear eigvector;
end

[nSmp,nFea] = size(data);

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

[dumpVec,eigvalue,eigvector] = svd(Hb,'econ');
eigvalue = diag(eigvalue);
eigIdx = find(eigvalue < 1e-3);
eigvalue(eigIdx) = [];
eigvector(:,eigIdx) = [];

eigvalue = eigvalue.^2;
eigvector = eigvector_PCA*(repmat(eigvalue_PCA,1,length(eigvalue)).*eigvector);

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

