%PCAKLM Principal Component Analysis/Karhunen-Loeve Mapping
% (PCA or MCA of overall/mean covariance matrix)
%
% [W,FRAC] = PCAKLM(TYPE,A,N)
% [W,N] = PCAKLM(TYPE,A,FRAC)
%
% INPUT
% A Dataset
% TYPE Type of mapping: 'pca' or 'klm'. Default: 'pca'.
% N or FRAC Number of dimensions (>= 1) or fraction of variance (< 1)
% to retain; if > 0, perform PCA; otherwise MCA.
% Default: N = inf.
%
% OUTPUT
% W Affine Karhunen-Loeve mapping
% FRAC or N Fraction of variance or number of dimensions retained.
%
% DESCRIPTION
% Performs a principal component analysis (PCA) or minor component analysis
% (MCA) on the overall or mean class covariance matrix (weighted by the
% class prior probabilities). It finds a rotation of the dataset A to an
% N-dimensional linear subspace such that at least (for PCA) or at most (for
% MCA) a fraction FRAC of the total variance is preserved.
%
% PCA is applied when N (or FRAC) >= 0; MCA when N (or FRAC) < 0. If N is
% given (abs(N) >= 1), FRAC is optimised. If FRAC is given (abs(FRAC) < 1),
% N is optimised.
%
% Objects in a new dataset B can be mapped by B*W, W*B or by A*KLM([],N)*B.
% Default (N = inf): the features are decorrelated and ordered, but no
% feature reduction is performed.
%
% ALTERNATIVE
%
% V = PCAKLM(A,TYPE,0)
%
% Returns the cumulative fraction of the explained variance. V(N) is the
% cumulative fraction of the explained variance by using N eigenvectors.
%
% This function should not be called directly, only trough PCA or KLM.
% Use FISHERM for optimizing the linear class separability (LDA).
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% MAPPINGS, DATASETS, PCLDC, KLLDC, PCAM, KLM, FISHERM
% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands
% $Id: pcaklm.m,v 1.15 2010/02/08 15:29:48 duin Exp $
function [w,truefrac] = pcaklm (type,a,frac)
truefrac = [];
% Default: preserve all dimensions (identity mapping).
if (nargin < 3) | (isempty(frac))
frac = inf;
prwarning (3,'no dimensionality given, only decorrelating and ordering dimensions');
end
% Default: perform PCAM.
if (nargin < 1) | (isempty(type))
type = 'pcam';
prwarning (3,'no type given, assuming PCA');
end
if (strcmp(type,'pcam'))
mapname = 'PCA';
elseif (strcmp(type,'klm'))
mapname = 'Karhunen-Loeve Mapping';
else
error('Unknown type specified');
end
%DXD Make the name a bit more informative:
if isfinite(frac)
if (frac<1)
mapname = [mapname sprintf(' ret. %4.1f%% var',100*frac)];
else
mapname = [mapname sprintf(' to %dD',frac)];
end
end
% Empty mapping: return straightaway.
if (nargin < 2) | (isempty(a))
w = prmapping(type,frac);
w = setname(w,mapname);
return
end
%nodatafile(a);
if ~isdataset(a)
if isa(a,'double')
a = prdataset(a,1); % make sure we have a dataset
else
error('nodatafile','Data should be given in a dataset or as doubles')
end
end
islabtype(a,'crisp','soft');
isvaldfile(a,1); % at least 1 object per class
a = setfeatdom(a,[]); % get rid of domain testing
[m,k,c] = getsize(a);
p = getprior(a);
a = setprior(a,p); % make class frequencies our prior
% If FRAC < 0, perform minor component analysis (MCA) instead of
% principal component analysis.
mca = (frac < 0); frac = abs(frac);
% Shift mean of data to origin.
b = a*scalem(a);
% If there are less samples M than features K, first perform a lossless
% projection to the (M-1) dimensional space spanned by the samples.
if (m <= k)
testdatasize(b,'objects');
u = reducm(b); b = b*u;
korg = k; [m,k] = size(b);
frac = min(frac,k);
else
testdatasize(b,'features');
u = [];
end
% Calculate overall or average class prior-weighted covariance matrix and
% find eigenvectors F.
if (strcmp(type,'pcam'))
if (c==0) | ~islabtype(a,'crisp') % we have unlabeled data!
G = prcov(+b); % use all
else
bb = [];
classsiz = classsizes(b);
for j = 1:c
%bb = [bb; seldat(b,j)*filtm([],'double')*p(j)/classsiz(j)];
bb = [bb; double(+seldat(b,j))*p(j)/classsiz(j)];
end
%[U,G] = meancov(remclass(setnlab(bb*m,1)));
[U,G] = meancov(bb);
end
else
%DXD For high dimensional dataset with many classes, we cannot
%store all individual cov. matrices in memory (like in the next
%line), but we have to compute them one by one:
%[U,GG] = meancov(b,1);
G = zeros(k,k);
for i = 1:c
%G = G + p(i)*GG(:,:,i);
[U,GG] = meancov(seldat(b,i),1);
G = G + p(i)*GG;
end
end
[F,V] = preig(G); % overdone if reducm has been called
% v = V(I) contains the sorted eigenvalues:
% descending for PCAM, ascending for MCA.
if (mca)
[v,I] = sort(diag(V));
else
[v,I] = sort(-diag(V)); v = -v;
end
if (frac == inf) % Return all dimensions, decorrelated and ordered.
n = k; truefrac = k;
elseif (frac == 0) % Just return cumulative retained variance.
w = cumsum(v)/sum(v);
return
elseif (frac >= 1) % Return FRAC dimensions.
n = abs(frac);
if (n > k),
error('illegal dimensionality requested');
end
I = I(1:n);
sv = sum(v);
if (sv ~= 0),
truefrac = cumsum(v(1:n))/sv;
else,
truefrac = 0;
end;
elseif (frac > 0) % Return the N dimensions that retain at least (PCA)
% or at most (MCA) FRAC variance.
J = find(cumsum(v)/sum(v) > frac);
if (mca), n = J(1)-1; else, n = J(1); end;
truefrac = n; I = I(1:n);
end
% If needed, apply pre-calculated projection to (M-1) dimensional subspace.
if (~isempty(u))
rot = u.data.rot*F(:,I);
off = u.data.offset*F(:,I);
else
rot = F(:,I);
off = -mean(a*F(:,I));
end
% Construct affine mapping.
w = affine(rot,off,a);
w = setdata(w,v,'eigenvalues');
w = setname(w,mapname);
return