function predicted = cosmo_classify_lda(samples_train, targets_train, samples_test, opt)
% linear discriminant analysis classifier - without prior
%
% predicted=cosmo_classify_lda(samples_train, targets_train, samples_test[,opt])
%
% Inputs:
% samples_train PxR training data for P samples and R features
% targets_train Px1 training data classes
% samples_test QxR test data
% opt Optional struct with optional field:
% .regularization Used to regularize covariance matrix. Default .01
% .max_feature_count Used to set the maximum number of features,
% defaults to 5000. If R is larger than this
% value, an error is raised. This is a (conservative)
% safety limit to avoid huge memory consumption
% for large values of R, because training the
% classifier typically requires in the order of 8*R^2
% bytes of memory
%
% Output:
% predicted Qx1 predicted data classes for samples_test
%
% Example:
% ds=cosmo_synthetic_dataset('ntargets',5,'nchunks',15);
% test_chunk=1;
% te=cosmo_slice(ds,ds.sa.chunks==test_chunk);
% tr=cosmo_slice(ds,ds.sa.chunks~=test_chunk);
% pred=cosmo_classify_lda(tr.samples,tr.sa.targets,te.samples,struct);
% % show targets and predicted labels
% disp([te.sa.targets pred])
% %|| 1 1
% %|| 2 2
% %|| 3 3
% %|| 4 4
% %|| 5 5
%
% Notes:
% - this classifier does not support a prior, that is it assumes that all
% classes have the same number of samples. If that is not the case an
% error is thrown.
% - a safety limit of opt.max_feature_count is implemented because a large
% number of features can crash Matlab / Octave, and/or make it very slow.
%
% (Contributions by Joern Diedrichsen, Tobias Wiestler, Nikolaas Oosterhof)
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
if nargin<4 || isempty(opt)
opt=struct();
end
if ~isfield(opt, 'regularization'), opt.regularization=.01; end
if ~isfield(opt, 'max_feature_count'), opt.max_feature_count=5000; end
% support repeated testing on different data after training every time
% on the same data. This is achieved by caching the training data
% and associated model
persistent cached_targets_train;
persistent cached_samples_train;
persistent cached_opt;
persistent cached_model;
if isequal(cached_targets_train, targets_train) && ...
isequal(cached_opt, opt) && ...
isequal(cached_samples_train, samples_train)
% use cache
model=cached_model;
else
% train classifier
model=train(samples_train, targets_train, opt);
% store model
cached_targets_train=targets_train;
cached_samples_train=samples_train;
cached_opt=opt;
cached_model=model;
end
% test classifier
predicted=test(model, samples_test);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function model=train(samples_train, targets_train, opt)
% train LDA classifier
[ntrain, nfeatures]=size(samples_train);
ntrain_=numel(targets_train);
if ntrain_~=ntrain
error(['size mismatch: samples_train has %d rows, '...
'targets_train has %d values'],ntrain,ntrain_);
end
if nfeatures>opt.max_feature_count
% compute size requirements for numbers.
% tyically this is 8 bytes per number, unless single precision
% is used.
w=whos('samples_train');
size_per_number=w.bytes/numel(samples_train);
% the covariance matrix is nfeatures x nfeatures
mem_required=nfeatures^2*size_per_number;
mem_required_mb=round(mem_required)/1e6; % in megabytes
error(['A large number of features (%d) was found, '...
'exceeding the safety limit max_feature_count=%d '...
'for the %s function.\n'...
'This limit is imposed because computing the '...
'covariance matrix will require in the order of '...
'%.0f MB of memory, and inverting it may take '...
'a long time.\n'...
'The safety limit can be changed by setting '...
'the ''max_feature_count'' option to another '...
'value, but large values ***may freeze or crash '...
'the machine***.'],...
nfeatures,opt.max_feature_count,...
mfilename(),mem_required_mb);
end
classes=fast_vector_unique(targets_train);
nclasses=numel(classes);
class_mean=zeros(nclasses,nfeatures); % class means
class_cov=zeros(nfeatures); % within-class variability
% compute mean and (co)variance
for k=1:nclasses;
% select data in this class
msk=targets_train==classes(k);
% number of samples in k-th class
n=sum(msk);
if k==1
if n<2
error(['Need at least two samples per class '...
'in training']);
end
nfirst=n; % keep track of number of samples
elseif nfirst~=n
error(['Different number of classes (%d and %d) - this '...
'is not supported. When using partitions, '...
'consider using cosmo_balance_partitions'], ...
n, nfirst);
end
class_samples=samples_train(msk,:);
class_mean(k,:) = sum(class_samples,1)/n; % class mean
res = bsxfun(@minus,class_samples,class_mean(k,:)); % residuals
class_cov = class_cov+res'*res; % estimate common covariance matrix
end;
% apply regularization
regularization=opt.regularization;
class_cov=class_cov/ntrain;
reg=eye(nfeatures)*trace(class_cov)/max(1,nfeatures);
class_cov_reg=class_cov+reg*regularization;
% linear discriminant
class_weight=class_mean/class_cov_reg;
class_offset=sum(class_weight .* class_mean,2);
model=struct();
model.class_offset=class_offset;
model.class_weight=class_weight;
model.classes=classes;
function predicted=test(model, samples_test)
% test LDA classifier
class_offset=model.class_offset;
class_weight=model.class_weight;
classes=model.classes;
if size(samples_test,2)~=size(class_weight,2)
error('test set has %d features, train set has %d',...
size(class_weight,2),size(samples_test,2));
end
class_proj=bsxfun(@plus,-.5*class_offset,class_weight*samples_test');
[unused,class_idxs]=max(class_proj);
predicted=classes(class_idxs);
function unq_xs=fast_vector_unique(xs)
xs_sorted=sort(xs(:));
idxs=([true;diff(xs_sorted)>0]);
unq_xs=xs_sorted(idxs);
