subsc.m

The Matlab Toolbox for Pattern Recognition

%SUBSC Subspace Classifier
%
%   W = SUBSC(A,N)
%   W = SUBSC(A,FRAC)
%
% INPUT
%   A          Dataset
%   N or FRAC  Desired model dimensionality or fraction of retained 
%              variance per class
%
% OUTPUT
%   W          Subspace classifier  
%
% DESCRIPTION
% Each class in the trainingset A is described by linear subspace of
% dimensionality N, or such that at least a fraction FRAC of its variance
% is retained. This is realised by calling PCA(AI,N) or PCA(AI,FRAC) for
% each subset AI of A (objects of class I). For each class a model is
% built that assumes that the distances of the objects to the class
% subspaces follow a one-dimensional distribution. 
%
% New objects are assigned to the class of the nearest subspace.
% Classification by D = B*W, in which W is a trained subspace classifier
% and B is a testset, returns a dataset D with one-dimensional densities
% for each of the classes in its columns.
%
% REFERENCE
% E. Oja, The Subspace Methods of Pattern Recognition, Wiley, New York, 1984.
%
% See DATASETS, MAPPINGS, PCA, FISHERC, FISHERM, GAUSSM

% Copyright: R.P.W. Duin, [e-mail]
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

function W = subsc(A,N)

  name = 'Subspace classf.';
  
  % handle default
  if nargin < 2, N = 1; end
  
  % handle untrained calls like subsc([],3);
  if nargin < 1 | isempty(A)
    W = mapping(mfilename,{N});
    W = setname(W,name);
    return
  end
  
  if isa(N,'double')
    
  % handle training like A*subsc, A*subsc([],3), subsc(A)
  % PRTools takes care that they are all converted to subsc(A,N)
      
    islabtype(A,'crisp');     % allow crisp labels only
    isvaldset(A,1,2);         % at least one object per class, two objects
    [m,k,c] = getsize(A);     % size of the training set
    for j = 1:c               % run over all classes
      B = seldat(A,j);        % get the objects of a single class only
      u = mean(B);            % compute its mean
      B = B - repmat(u,size(B,1),1); % subtract mean
      v = pca(B,N);           % compute PCA for this class
      v = v*v';               % trick: affine mappings in the original space
      B = B - B*v;            % differences of objects and their mappings
      s = mean(sum(B.*B,2));  % mean square error w.r.t. the subspace
      data(j).u = u;          % store mean
      data(j).w = v;          % store mapping
      data(j).s = s;          % store mean square distance
    end
                              % define trained mapping,
                              % store class labels and size
    W = mapping(mfilename,'trained',data,getlablist(A),k,c);
    W = setname(W,name);
    
  elseif isa(N,'mapping')
    
  % handle evaluation of a trained subspace classifier W for a dataset A.
  % The command D = A*W is by PRTools translated into D = subsc(A,W)
  % Such a call is detected here by the fact that N appears to be a mapping.
    
    W = N;                    % avoid confusion: call the mapping W
    m = size(A,1);            % number of test objects
    [k,c] = size(W);          % mapping size: from K features to C classes
    d = zeros(m,c);           % output: C class densities for M objects
    
    for j=1:c                 % run over all classes
      u = W.data(j).u;        % class mean in training set
      v = W.data(j).w;        % mapping to subspace in original space
      s = W.data(j).s;        % mean square distance
      B = A - repmat(u,m,1);  % substract mean from test set
      B = B - B*v;            % differences objects and their mappings
      d(:,j) = sum(B.*B,2)/s; % convert to distance and normalise 
    end
    d = exp(-d/2)/sqrt(2*pi); % convert to normal density
    
    A = dataset(A);           % make sure A is a dataset 
    d = setdata(A,d,getlabels(W)); % take data from D and use 
                              % class labels as given in W
                              % other information in A is preserved   
    W = d;                    % return result in output variable W
    
  else     
    
    error('Illegal call')     % this should not happen
    
  end
  
return