greedyappx.m

用matlab实现的统计模式识别工具箱

function [sel_inx,Alpha,Z,kercnt,MsErr,MaxErr]=...
  greedyappx(X,ker,arg,m,p,mserr,maxerr,verb) 
% GREEDYAPPX Kernel greedy data approximation.
%
% Synopsis:
%  [inx,Alpha,kercnt,mserr,maxerr]=...
%    greedyappx(X,ker,arg,m,m2,mserr,maxerr,verb)
%
% Description:
%  This function aims to select a subset S of input data X such
%  that the feature space representation of X can be well 
%  approximated by feature space representation of S.
%  The feature represenation of data is by the use of
%  specified kernel function.
%
%  The greedy algortihm is used to seletect the subset S. 
%  The algorithm iterates until on of the following stopping 
%  conditions is achieved:
%    - number of vectors of S achieves m 
%    - maximal reconstruction error is less than maxerr 
%    - mean squared sum of reconstruction errors less than mserr. 
%  
% Input:
%  X [dim x num_data] Input data.
%  ker [string] Kernel identifier. See 'help kernel' for more info.
%  arg [...] Argument of selected kernel.
%  m [1x1] Maximal number of vector used for approximation.
%  p [1x1] Depth of search for the best basis vector.
%  mserr [1x1] Desired mean sum of squared reconstruction errors.
%  maxerr [1x1] Desired maximal reconstruction error.
%  verb [1x1] If 1 then infor about process is displayed.
%
% Output:
%  inx [1 x n] Indices of selected vector, i.e., S = X(:,inx).
%  Alpha [m x m] Koefficient of the kernel projection of data on the
%    found base vectors, i.e., z = Alpha*kernel(S,x,ker,arg).
%  Z [m x num_data] Training data projected on the found base vectors.
%  kercnt [1x1] Number of used kernel evaluations.
%  MsErr [1 x m] Sum of squared reconstruction errors for corresponding
%   number of base vectors.
%  MaxErr [1 x m] Maximal squared reconstruction error for crresponding 
%
% See also 
%  GREEDYKPCA.
%

% About: Statistical Pattern Recognition Toolbox
% (C) 1999-2003, Written by Vojtech Franc and Vaclav Hlavac
% <a href="http://www.cvut.cz">Czech Technical University Prague</a>
% <a href="http://www.feld.cvut.cz">Faculty of Electrical Engineering</a>
% <a href="http://cmp.felk.cvut.cz">Center for Machine Perception</a>

% Modifications:
% 10-dec-2004, VF, tmp(find(Errors<=0)) = -inf; added to evoid num errors.
% 5-may-2004, VF
% 13-mar-2004, VF
% 10-mar-2004, VF
% 9-mar-2004, addopted from greedyappx

if nargin < 5, mserr=1e-6; end
if nargin < 6, maxerr=1e-6; end
if nargin < 7, verb=0; end

[dim,num_data]=size(X);

if verb, 
  fprintf('Greedy data approximation.\n');
  fprintf('Setting: ker=%s, arg=%f, m=%d, eps=%f\n', ker,arg,m,maxerr); 
end

kercnt=0;
Errors = diagker(X,ker,arg)';  kercnt = kercnt+num_data;
Z = zeros(m,num_data);

MsErr = [];
MaxErr = [];
Alpha=zeros(m,m);
SV = zeros(dim,m);
sel_inx=[];
work_inx = [1:num_data];

for i=1:m,

  % call greedyappx2
  if i == 1,
    [tmp_sel_inx,tmp_Alpha,tmp_Z,tmp_kercnt,tmp_MsErr,tmp_MaxErr]=...
       ordinary_greedyappx(X,ker,arg,p,0,1e-12,verb);
%    [tmp_sel_inx,tmp_Alpha,tmp_Z,tmp_kercnt,tmp_MsErr,tmp_MaxErr]=...
%       ordinary_greedyappx(X,ker,arg,p,mserr,maxerr,verb);
    kercnt = kercnt+tmp_kercnt;
  else
    init_model.Alpha = Alpha(1:i-1,1:i-1);
    init_model.Z = Z(1:i-1,:);
    init_model.Errors = Errors;

    [tmp_sel_inx,tmp_Alpha,tmp_Z,tmp_kercnt,tmp_MsErr,tmp_MaxErr]=...
       ordinary_greedyappx(X,ker,arg,p,0,1e-12,verb,init_model);
%    [tmp_sel_inx,tmp_Alpha,tmp_Z,tmp_kercnt,tmp_MsErr,tmp_MaxErr]=...
%       ordinary_greedyappx(X,ker,arg,p,mserr,maxerr,verb,init_model);
    kercnt = kercnt+tmp_kercnt;
  end

  tmp_Z = tmp_Z(i:size(tmp_Z,1),:);
  A = tmp_Z*tmp_Z';

  tmp1 = sum(tmp_Z.^2,1);
  tmp1(find(tmp1==0))=inf;
  tmp = sum((A*tmp_Z).*tmp_Z,1)./tmp1;
  tmp(sel_inx) = -inf; tmp(find(Errors<=0)) = -inf;
  [dummy,new_inx]=max(tmp);
 
%  [V,D] = eig(A);
%  D=diag(D);
%  [dummy,inx]=max(D);
%  z = V(:,inx);
%  dummy = z'*A*z/(z'*z)  
   
  % orthonormalization
  sel_inx = [sel_inx new_inx];
  tmp = kernel( X(:,new_inx), X(:,work_inx), ker, arg );
  kercnt=kercnt+num_data-i;

  if i > 1,
     Z(i,work_inx) = ...
       (tmp - Z(1:i-1,new_inx)'*Z(1:i-1,work_inx))/sqrt(Errors(new_inx));

     Alpha(i,:) = - Z(1:i-1,new_inx)'*Alpha(1:i-1,:);
     Alpha(i,i) = 1;
     Alpha(i,:) = Alpha(i,:)/sqrt(Errors(new_inx));
  else
     Z(1,:) = tmp/sqrt(Errors(new_inx));
     Alpha(1,1) = 1/sqrt(Errors(new_inx));
  end

%  SV(:,i) = x;
  % Error(i) = k(i,i)-k'(i,i)
  Errors(work_inx) = Errors(work_inx) - Z(i,work_inx).^2;
  Errors(find(Errors<0)) = 0;
%  Errors(sel_inx)=zeros(1,length(sel_inx));
  work_inx(find(work_inx==new_inx)) = [];
  
  % store errors
  MsErr = [MsErr, sum(Errors)/num_data];
  MaxErr = [MaxErr, max(Errors)];
  
  if verb,
    fprintf('%d: maxerr=%f, mserr=%f, inx=%d\n', ...
        i,MaxErr(end), MsErr(end), new_inx);
  end
  
  % evaluate stopping conditions:
  if maxerr >= MaxErr(end) | mserr >= MsErr(end),
    break;
  end
end

%if MsErr(end) < 0, 
%  i = i-1
%  MaxErr = MaxErr(1:end-1);
%  MsErr = MsErr(1:end-1);
%  sel_inx = sel_inx(1:end-1);
%end % Patch to avaid nummerical errors

% cut off non-used memory if number of used base vector is less than m
Alpha=Alpha(1:i,1:i);
Z = Z(1:i,:);

return;


%=================================================
function [sel_inx,Alpha,Z,kercnt,MsErr,MaxErr]=...
  ordinary_greedyappx(X,ker,arg,m,mserr,maxerr,verb,init_model) 

[dim,num_data]=size(X);

kercnt=0;
sel_inx=[];              % indices of seleted base vectors
work_inx = [1:num_data]; % indices of the rest 
MsErr = [];
MaxErr = [];

if nargin < 8,
  Errors = diagker(X,ker,arg)';
  Z = zeros(m,num_data);
  Alpha=zeros(m,m);
  curr_m = 0;
else
  Errors = init_model.Errors;
  curr_m = size(init_model.Z,1);
  m = m + curr_m;
  Z = zeros(m,num_data);
  Alpha=zeros(m,m);

  Z(1:curr_m,:) = init_model.Z;
  Alpha(1:curr_m,1:curr_m) = init_model.Alpha;
end

if verb, fprintf('('); end

for i=curr_m+1:m,
  
  % find vector with highest reconstruction error
  [curr_maxerr,new_inx] = max( Errors );
  sel_inx = [sel_inx,new_inx];

  % orthonormalization
  tmp = kernel( X(:,new_inx), X(:,work_inx), ker, arg ); 
  kercnt = kercnt + num_data - i;

  if i > 1,
     Z(i,work_inx) = ...
       (tmp - Z(1:i-1,new_inx)'*Z(1:i-1,work_inx))/sqrt(Errors(new_inx)); 

     Alpha(i,:) = - Z(1:i-1,new_inx)'*Alpha(1:i-1,:);
     Alpha(i,i) = 1;
     Alpha(i,:) = Alpha(i,:)/sqrt(Errors(new_inx));
  else
     Z(1,:) = tmp/sqrt(Errors(new_inx));
     Alpha(1,1) = 1/sqrt(Errors(new_inx));
  end

  % Error(i) = k(i,i)-k'(i,i)
  Errors(work_inx) = Errors(work_inx) - Z(i,work_inx).^2;  
%  Errors(sel_inx)=zeros(1,length(sel_inx));
  work_inx(find(work_inx==new_inx)) = [];
    
  % store errors
  MsErr = [MsErr, sum(Errors)/num_data];
  MaxErr = [MaxErr, max(Errors)];
  
  if verb,
    fprintf('.', i, MsErr(end) );
  end
  
  % evaluate stopping conditions:
  if maxerr >= MaxErr(end) | mserr >= MsErr(end),
    break;
  end
end

if verb, fprintf(')\n'); end

% cut off non-used memory if number of used base vector is less than m
Alpha=Alpha(1:i,1:i);
Z = Z(1:i,:);

return;