conrvm.m

关于人脸识别的一个VC程序

% CONRVM	The Relevance Vector Machine - constructive version
% 
%	[WEIGHTS,NZ] = CONRVM(X,T,ITS,ALPHA,KERNEL,ETA)
% 
% 
% Output:	WEIGHTS	The weights corresponding to the RV's
%		NZ	A binary 'indicator' (N+1)-vector,
%			specifying the relevance vectors.
%			NZ(1) corresponds to the bias (which may not be
%			present, unlike the SVM) while a 1 in the
%			remaining N locations signifies
%			that the corresponding kernel is utilised. 
%			
%			So LENGTH(WEIGHTS) = SUM(NZ).
% 
% Input:	X	The data, in row vector form.
%
%		T	Target N-vector, containing 1's and 0's.
%
%		ITS	Three element vector containing iteration parameters
%			ITS(1)	The number of 'outer' iterations, for
%				each of which a new kernel is added.
%				The function may return early if there
%				are no 'good' kernels remaining.
%			ITS(2)	The number of iterations of optimisation
%				of the whole model after adding a kernel 
%				'Suggested' good value = 5.
%			ITS(3)	The number of iterations in the IRLS loop.
%				'Suggested' good value = 25.
%			ITS(4)	Number of iterations at the finish.
%				'Suggested' value = minimum 250.
% 
%			Entering '0' (zero) for a value will result in
%			the choice of a 'sensible' default.
% 
%		ALPHA	Initial value of alpha, 'suggested' value = 1e-3.
% 
%		KERNEL	Kernel type:	e.g. 'gauss','poly2','poly3', etc.
%			(See 'kernelFunction.m')
%		ETA	Kernel inverse squared - 'width' parameter. 
%			e.g. for Gaussian, K = ETA * SQUARED_DISTANCE.
% 
% Note: you should set the parameters KERNEL_MEMORY and CACHE_SIZE,
% (probably to the same value) defined at the beginning of the function,
% according to the memory available on your machine.
% 

function [w_nz, nonZero] = ...
    conrvm(X,T,iterations,alphaInit,kernel_,eta)

KERNEL_MEMORY	= 128;	% in MB
CACHE_SIZE	= 128;	% again, in MB

iterations	= [iterations(:) ; zeros(4-length(iterations),1)];
defaultIts	= [100 5 25 250];
i		= iterations==0;
iterations(i)	= defaultIts(i);

outerIts	= iterations(1);
its1		= iterations(2);
its2		= iterations(3);
finalIts	= iterations(4);

FAST_METHOD	= 1;
CHANGE_SLOW	= 0;


W_STOP_CRITERION = 0;

EM_update	= 0;
alpha_factor	= 1e13;
 
if length(its1)==2
  mon	= abs(its1(2));
  if its1(2)<0
    plotmon = 1;
  else
    plotmon = 0;
  end
  its1	= its1(1);
else
  mon		= 0;
  plotmon	= 0;
end
[N d]	= size(X);

%
% Work out if we can afford to compute the matrix of basis outputs
% 
memoryRequirement	= N*(N+1)*8/1e6; % megabytes
if memoryRequirement<=KERNEL_MEMORY
  KERNEL_MATRIX	= 1;
  fprintf('Computing full kernel matrix (%.0fMB)\n', memoryRequirement)
else
  KERNEL_MATRIX = 0;
  fprintf(['Full kernel matrix (%.1fMB) too large for '...
	  'allowed limit (%.1fMB)\n'], memoryRequirement, KERNEL_MEMORY)
  fprintf('\t-> computing kernels on-line with cache size of %.0fMB.\n',...
	  CACHE_SIZE)
end

if KERNEL_MATRIX
  H	= kernelFunction(X,X,kernel_,eta);
else
  n_cache	= CACHE_SIZE * 1e6 / (8*N);
  in_cache	= 1:n_cache;
  CACHE		= kernelFunction(X,X(in_cache,:),kernel_,eta,'NOBIAS');
end


M	= N+1;

nonZero		= logical([1 ; zeros(M-1,1)]);

alpha		= 1e16*ones(M,1);
w		= zeros(M,1);
alpha(nonZero)	= alphaInit;

options		= optimset('MaxIter',its2,'Display',0);
if KERNEL_MATRIX
  H_nz		= H(:,nonZero);
else
  H_nz		= ones(N,1);
end
w(nonZero)	= (H_nz'*H_nz + diag(alpha(nonZero))) \ (H_nz'*(2*T-1));

FINISH		= 0;
tryAlpha	= 1e12;
for o=0:outerIts
  if o
    f	= find(nonZero==0);
    ra	= randperm(length(f));
    for j	= 1:length(f)
      newh	= f(ra(j));
      if KERNEL_MATRIX
	H_new	= H(:,newh); 
      else
	if newh==1
	  H_new	= ones(N,1);
	else
	  %
	  % Check if its cache-d
	  % 
	  ci	= find(in_cache==newh);
	  if isempty(ci)
	    H_new	= kernelFunction(X,X(newh-1,:),kernel_,eta,'NOBIAS');
	  else
	    H_new	= CACHE(:,newh-1);
	  end
	end
      end
      
      if FAST_METHOD
	%
	% Some quick calc code
	% 
	betaH_new	= (vary.*H_new);
	hp	= H_nz'*betaH_new;
	hpp	= (H_new'* betaH_new) + tryAlpha;
	tmp	= Li'*hp;
	Spp	= 1 / (hpp - tmp'*tmp);
	%
	% now for the mean
	% 
	% mup	= Spp * (betaH_new'*T - hp'*w(nonZero));
	y	= H_nz*w_nz;
	betaTlin	= vary.*y + e;
	wTry	= Spp * (H_new'*betaTlin - hp'*w(nonZero));
	grad	= 1/tryAlpha - wTry^2 - Spp;
      else
	[wTry errs diagC] = ...
	    regirls([H_nz H_new],T,[alpha_nz ; tryAlpha],[w_nz ; 0],options);
	grad	= 1/tryAlpha - wTry(end)^2 - diagC(end);
      end

      if grad<0
	fprintf(['Cycle %3d: %2d kernels -> ' ...
		 'introducing kernel %4d \t(%d/%d)\n'],...
		o,sum(nonZero),newh,j,length(f))
	split		= max(find(newh>find(nonZero)));
	if isempty(split)
	  H_nz		= [H_new H_nz];
	else
	  H_nz		= [H_nz(:,1:split) H_new H_nz(:,split+1:end)];
	end
	
	nonZero(newh)	= logical(1);
	alpha(newh)	= alphaInit;
	w(newh)		= wTry(end);
	break
      end
    end
    if grad>=0 | o==outerIts
      
      if grad>=0 & FAST_METHOD & CHANGE_SLOW
	fprintf('Cycle %3d: switching to slow method\n',o)
	FAST_METHOD = 0;
      else
	FINISH	= 1;
	its1	= finalIts;
	mon	= 25;
	plotmon	= mon;
	if grad>=0
	  s_ = 'no good functions remain';
	else
	  s_ = 'natural termination';
	end
	fprintf('Cycle %3d: %s - final optimisation\n',o,s_)
      end
    end
  end
    
  for i=1:its1
    wOld		= w;
    badHess		= 1;
    alpha_cutoff	= min(alpha)*alpha_factor;
    while badHess
      alpha_cutoff	= alpha_cutoff/10;
      oldZero		= nonZero;
      % nonZero		= (alpha<alpha_cutoff) & (alpha>=0);
      nonZero		= nonZero & (alpha<alpha_cutoff) & (alpha>=0);
      numNZ		= sum(nonZero);
      alpha_nz		= alpha(nonZero);
      w_nz		= w(nonZero);

      if KERNEL_MATRIX
	H_nz	= H(:,nonZero);
      else
	%H_nz		= kernelFunction(X,X(nonZero(2:end),:),kernel_,eta);
	%if ~nonZero(1)
	%  H_nz(:,1)	= [];
	%end
	H_nz	= H_nz(:,nonZero(oldZero));
      end
      [w_nz errs diagC logDetH LD badHess Li e vary] = ...
	  regirls(H_nz,T,alpha_nz,w_nz,options);
    end
    
    evidenceMeasure(i)	= 0.5 * (sum(log(alpha_nz)) - ...
				 (w_nz.^2)'*alpha_nz - logDetH) + LD;
    
    % 
    % alpha re-estimation
    % 

    if numNZ>1
      gamma		= zeros(size(alpha_nz));
      gamma		= 1 - alpha_nz.*diagC;
      if ~EM_update
	alpha_nz	= gamma ./ (w_nz.^2);
      else
	alpha_nz	= 1./(diagC + w_nz.^2);
      end
    end
    
    if mon & ~rem(i,mon)
      fprintf('%d\tGamma = %.2f (nz = %d)\n', i, sum(gamma), numNZ)
      if plotmon
	subplot(1,3,1)
	plot(evidenceMeasure)
	subplot(1,3,2)
	stem(w,'.')
	hold off
	set(gca,'Xlim',[-1 M+1])
	subplot(1,3,3)
	hist(log10(alpha_nz),[-12:12])
	set(gca,'Xlim',[-13 13])
      end
      drawnow
    end
  
    alpha(nonZero)	= alpha_nz;
    w(nonZero)		= w_nz;
    delta_w		= max(abs(w-wOld));
    if delta_w<W_STOP_CRITERION
      fprintf('** Stopping at iteration %d: delta_w = %g\n', i, delta_w)
      break
    end
  end
  if FINISH
    break;
  end
end