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 greedyappxif nargin < 5, mserr=1e-6; endif nargin < 6, maxerr=1e-6; endif 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); endkercnt=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;  endend%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 mAlpha=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 vectorswork_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;endif verb, fprintf('('); endfor 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;  endendif verb, fprintf(')\n'); end% cut off non-used memory if number of used base vector is less than mAlpha=Alpha(1:i,1:i);Z = Z(1:i,:);return;