function [eigvector, eigvalue, elapse] =
LDA(gnd,options,data)
% LDA: Linear Discriminant Analysis
%
% [eigvector, eigvalue] = LDA(data, gnd, options)
%
% Input:
% data - Data matrix. Each row
vector of fea is a data point.
% gnd - Colunm vector of the
label information for each
% data
point.
% options - Struct value in Matlab. The fields in
options
% that can
be set:
%
% Regu -
1: regularized solution,
% a* = argmax
(a'data'WXa)/(a'data'DXa+ReguAlpha*I)
% 0: solve
the sinularity problem by SVD
% Default:
1
%
% ReguAlpha
- The regularization parameter. Valid
% when Regu==1. Default value is
0.1.
%
% ReguType -
'Ridge': Tikhonov regularization
% 'Custom': User provided
% regularization matrix
% Default:
'Ridge'
% regularizerR -
(nFea x nFea) regularization
% matrix
which should be provided
% if ReguType
is 'Custom'. nFea is
% the feature
number of data
% matrix
% Fisherface
- 1:
Fisherface approach
% PCARatio = nSmp - nClass
% Default:
0
%
% PCARatio - The percentage of
principal
% component kept in the PCA
% step. The percentage is
% calculated based on the
% eigenvalue. Default is 1
% (100%, all the non-zero
% eigenvalues will be kept.
% If PCARatio > 1, the PCA step
% will keep exactly PCARatio principle
% components (does not exceed the
% exact number of non-zero components).
%
%
% 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
%
% Examples:
%
% fea
= rand(50,70);
% gnd
= [ones(10,1);ones(15,1)*2;ones(10,1)*3;ones(15,1)*4];
% options = [];
% options.Fisherface = 1;
% [eigvector, eigvalue] = LDA(gnd, options, fea);
% Y =
fea*eigvector;
%
%
% See also LPP, constructW, LGE
if ~exist('data','var')
global data;
end
if (~exist('options','var'))
options = [];
end
if ~isfield(options,'Regu') | ~options.Regu
bPCA = 1;
if
~isfield(options,'PCARatio')
options.PCARatio = 1;
end
else
bPCA = 0;
if
~isfield(options,'ReguType')
options.ReguType = 'Ridge';
end
if
~isfield(options,'ReguAlpha')
options.ReguAlpha = 0.1;
end
end
tmp_T = cputime;
% ====== Initialization
[nSmp,nFea] = size(data);
if length(gnd) ~= nSmp
error('gnd and data
mismatch!');
end
classLabel = unique(gnd);
nClass = length(classLabel);
Dim = options.dim;
% Dim = nClass - 1;
if Dim > nClass - 1
error('the reduced
dimension is large!');
end
if bPCA & isfield(options,'Fisherface') &
options.Fisherface
options.PCARatio = nSmp
- nClass;
end
if issparse(data)
data = full(data);
end
sampleMean = mean(data,1);
data = (data - repmat(sampleMean,nSmp,1));
bChol = 0;
if bPCA & (nSmp > nFea+1) & (options.PCARatio >=
1)
DPrime =
data'*data;
DPrime =
max(DPrime,DPrime');
[R,p] =
chol(DPrime);
if p == 0
bPCA = 0;
bChol = 1;
end
end
%======================================
% SVD
%======================================
if bPCA
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);
%=======================================
if options.PCARatio > 1
idx =
options.PCARatio;
if idx
< length(eigvalue_PCA)
eigvalue_PCA =
eigvalue_PCA(1:idx);
eigvector_PCA =
eigvector_PCA(:,1:idx);
end
elseif options.PCARatio < 1
sumEig =
sum(eigvalue_PCA);
sumEig =
sumEig*options.PCARatio;
sumNow =
0;
for idx =
1:length(eigvalue_PCA)
sumNow = sumNow +
eigvalue_PCA(idx);
if sumNow >= sumEig
break;
end
end
eigvalue_PCA = eigvalue_PCA(1:idx);
eigvector_PCA = eigvector_PCA(:,1:idx);
end
%=======================================
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);
%=======================================
if options.PCARatio > 1
idx =
options.PCARatio;
if idx
< length(eigvalue_PCA)
eigvalue_PCA =
eigvalue_PCA(1:idx);
eigvector =
eigvector(:,1:idx);
end
elseif options.PCARatio < 1
sumEig =
sum(eigvalue_PCA);
sumEig =
sumEig*options.PCARatio;
sumNow =
0;
for idx =
1:length(eigvalue_PCA)
sumNow = sumNow +
eigvalue_PCA(idx);
if sumNow >= sumEig
break;
end
end
eigvalue_PCA = eigvalue_PCA(1:idx);
eigvector
= eigvector(:,1:idx);
end
%=======================================
eigvalue_PCA = eigvalue_PCA.^-.5;
eigvector_PCA =
(data'*eigvector).*repmat(eigvalue_PCA',nFea,1);
data = eigvector;
clear eigvector;
end
else
if ~bChol
DPrime = data'*data;
% options.ReguAlpha =
nSmp*options.ReguAlpha;
switch lower(options.ReguType)
case
{lower('Ridge')}
for i=1:size(DPrime,1)
DPrime(i,i) = DPrime(i,i) +
options.ReguAlpha;
end
case
{lower('Tensor')}
DPrime = DPrime +
options.ReguAlpha*options.regularizerR;
case
{lower('Custom')}
DPrime = DPrime +
options.ReguAlpha*options.regularizerR;
otherwise
error('ReguType does not
exist!');
end
DPrime = max(DPrime,DPrime');
end
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
elapse.timeW = 0;
elapse.timePCA = cputime - tmp_T;
tmp_T = cputime;
if bPCA
[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);
else
WPrime = Hb'*Hb;
WPrime =
max(WPrime,WPrime');
dimMatrix =
size(WPrime,2);
if Dim >
dimMatrix
Dim = dimMatrix;
end
if
isfield(options,'bEigs')
if options.bEigs
bEigs =
1;
else
bEigs =
0;
end
else
if (dimMatrix > 1000 & Dim <
dimMatrix/10) | (dimMatrix > 500 & Dim < dimMatrix/20) |
(dimMatrix > 250 & Dim <
dimMatrix/30)
bEigs =
1;
else
bEigs =
0;
end
end
if bEigs
%disp('use eigs to speed up!');
option = struct('disp',0);
if bChol
option.cholB = 1;
[eigvector, eigvalue] = eigs(WPrime,R,Dim,'la',option);
else
[eigvector, eigvalue] = eigs(WPrime,DPrime,Dim,'la',option);
end
eigvalue = diag(eigvalue);
else
[eigvector, eigvalue] =
eig(WPrime,DPrime);
eigvalue = diag(eigvalue);
[junk, index] = sort(-eigvalue);
eigvalue = eigvalue(index);
eigvector = eigvector(:,index);
if Dim < size(eigvector,2)
eigvector
= eigvector(:, 1:Dim);
eigvalue =
eigvalue(1:Dim);
end
end
end
for i = 1:size(eigvector,2)
eigvector(:,i) =
eigvector(:,i)./norm(eigvector(:,i));
end
elapse.timeMethod = cputime - tmp_T;
elapse.timeAll = elapse.timePCA + elapse.timeMethod;