lda.m

It is for Face Recognition

%function [EigenVal, LDAMat]=lda(Sample, nVecNum)
function [EigenVal, LDAMat]=lda(Sample, nVecNum)
% LDA Linear Discriminant Analysis.

if nargin < 1,
    disp('Invalid parameter number');
    return;
end
if nargin == 1,
    nVecNum = size(Sample,2)-1;
end

% ====== Initialization
data_n = size(Sample, 1);
feature_n = size(Sample,2)-1;
featureMatrix = Sample(:, 1:end-1); 
classLabel = Sample(:, end);
clear Sample;
[diffClassLabel, classSize] = countele(classLabel);
class_n = length(diffClassLabel);
sampleMean = mean(featureMatrix);

% ====== Compute B and W
% ====== B: between-class scatter matrix
% ====== W:  within-class scatter matrix
B = zeros(feature_n,feature_n);
W = zeros(feature_n,feature_n);
for i = 1:class_n,
	index = find(classLabel==diffClassLabel(i));
    X_i = featureMatrix(index, :);
	classMean = mean(X_i,1);
    % B = B + length(index) * (classMean - sampleMean)' * (classMean - sampleMean);
    % candidate: B = B + (classMean - sampleMean)' * (classMean - sampleMean);
    B = B + (classMean - sampleMean)' * (classMean - sampleMean);
    % same as W = W + (length(index)-1) * cov(X_i);
    for j =1 : feature_n
        X_i(:,j) = X_i(:,j) - classMean(1,j);
    end
    W = W + X_i' * X_i;
end

% compute the generalized eigenvectors and eigenvalues
[VMat, DValMat] = eig(B, W);
DValMat = diag(DValMat); % transform eigenvalues to one vector
%ordering
[Temp, nIndex] = sort(abs(DValMat));
clear Temp;
nIndex = flipud(nIndex);
EigenVal = DValMat(nIndex);
LDAMat = VMat(:,nIndex);

% ====== Find the best discriminant vectors
%invW = inv(W);
%Q = invW*B;
%D = [];
%EigenVal = [];
%for i = 1:nVecNum,
%	[eigVec, eigVal] = eig(Q);
%	[nMaxEigVal, nIndex] = max(abs(diag(eigVal)));  
%	D = [D, eigVec(:, nIndex)];	% Each col of D is a eigenvector
%    EigenVal = [EigenVal, nMaxEigVal];
%	Q = (eye(feature_n)-invW*D*inv(D'*invW*D)*D')*invW*B;
%end
% charged: newSample = [featureMatrix*D(:,1:nVecNum) classLabel];
%discrim_vec = D;

% compute the eigenvectors and eigenvalues
%[LDAMat, EigenVal] = eig( inv(W) * B);
%EigenVal = diag(EigenVal); % transform eigenvalues to one vector
% the same as: [LDAMat, SCORE, EigenVal, tsquare] = princomp(fTrainMatrix)

% Compute projection matrix
%[U, EigenVal, V] = svd( inv(W) * B );
%LDAMat = V(:,1:nVecNum);
%EigenVal = diag(EigenVal);

% get the first nVecNum vectors
size(EigenVal)
LDAMat = LDAMat(:,1:nVecNum);

% translation
% model.b = - LDAMat * sampleMean;

return;