kerskfmat.m

该工具箱是用于统计模式识别的

function [Alpha,bias,sol,t,flps,margin,up,lo]=...
    kerskfmat(X,I,epsilon,ker,arg,tmax,C)
% KERSKFMAT fast version of KERNELSK. 
% [Alpha,bias,sol,t,flps,margin,up,lo]=kerskfmat(X,I,epsilon,ker,arg,tmax,C)
%
% Faster version (with respect to number of floating point operations) 
% of KERNELSK algorithm. It does not compute the whole kernel matrix.
% 
% Inputs:
%   X [NxL] training patterns, N is dimension and L number of patterns.
%   I [1xL] labels, 1 for 1st class and 2 for 2nd class.
%   epsilon [1x1] precision of found solution. The margin of found 
%       hyperplane is less than the optimal margin at most by epsilon. 
%   ker [string] kernel, see 'help kernel'.
%   arg [...] argument of given kernel, see 'help kernel'.
%   tmax [int] maximal number of iterations.
%   C [real] trade-off between margin and training error.
%  
% Outputs:
%   Alpha [1xL] Weights (Lagrangians) of patterns.
%   bias [real] bias (threshold) of found decision rule.
%   sol [int] 1 solution is found
%             0 algorithm stoped (t == tmax) before converged.
%            -1 hyperplane with margin greater then epsilon 
%               does not exist.
%   t [int] number of iterations.
%   margin [real] margin between classes.
%   flps [int] number of used floating point operations.
%   up [1,t] evolution of the upper bound on the optimal margin.
%   lo [1,t] evolution of the lower bound on the optimal margin.
%   CE [real] classification error on training patterns.
%
% See also KERNELSK, SVM.
%

% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz
% Written Vojtech Franc (diploma thesis) 02.11.1999, 13.4.2000
% Modifications
% 19-September-2001, V.Franc, comments changed.
% 18-August-2001, V.Franc, up and lo bounds added
% 9-August-2001, V.Franc, version without computing kernel matrices
% 13-July-2001, V.Franc, comments
% 12-July-2001, V.Franc, C, bias and normal vect. normalized.
% 11-July-2001, V.Franc, Rosta Horcik proved that the computation 
%                        of threshold is OK.
% 10-July-2001, V.Franc, derived from kekozinec2

flops(0);

% set default values of the input argiments
if nargin < 7,
  C = inf;
end

% maximal number of iteraions
if nargin < 6,
  tmax = inf;
end

% indexes of pattens in the 1st and 2nd class
xinx1 = find(I == 1);
xinx2 = find(I == 2);

X1=X(:,xinx1);      % patters from 1st class
X2=X(:,xinx2);      % patters from 1st class

l1 = size(X1,2);    % number os patterns
l2 = size(X2,2);

% make problem lin-separable in high dimensional space
if C ~= 0,
  kadd = 1/(2*C);
else
  kadd = 0;
end

% convex coeficients defining normal of the decision hyperplane
% (they correspond to the Lagrangian multiplyers).
s1 = zeros(l1, 1);
s2 = zeros(l2, 1);

% initial sol
s1(1)=1;       % take the 1st pattern from the 1st class
s2(1)=1;       % take the 2nd pattern from the 2st class

sol=0;
t = 0;

minXDA1=-inf;
minXDA2=-inf;

% -- Inicalization of temp. variables ---------------------------------
Di1 = zeros(l1,1);
Fi1 = zeros(l1,1);
for i=1:l1,  
  Di1(i) = kernel(X1(:,1),X1(:,i), ker,arg);
  Fi1(i) = kernel(X2(:,1),X1(:,i), ker,arg);
end
Di1(1) = Di1(1) + kadd;

Di2 = zeros(l2,1);
Fi2 = zeros(l2,1);
for i=1:l2,  
  Di2(i) = kernel(X1(:,1),X2(:,i), ker,arg);
  Fi2(i) = kernel(X2(:,1),X2(:,i), ker,arg);
end
Fi2(1) = Fi2(1)+kadd;

a = kernel(X1(:,1),X1(:,1),ker,arg) + kadd;
b = kernel(X2(:,1),X2(:,1),ker,arg) + kadd;
c = kernel(X1(:,1),X2(:,1),ker,arg );

% upper and lower bounds on the optimal margin
up=[];
lo=[];

% main cycle
while sol == 0 & tmax > t,

   t = t + 1;
   sol = 1;
   
   % -- compute auxciliary variables --   
 
   if sqrt( a -2*c +b) <= 0,    
     % algorithm has converged to the zero vector --> classes overlap
     sol = -1;
     break;
   end

   [emin,inx1] = min(Di1-Fi1);
   [gmin,inx2] = min(Fi2-Di2);
           
   % projection x \in X_1 on (w_1 - w_2)
   proj1 = (emin + b -c )/sqrt(a-2*c+b);
   
   % projection x \in X_2 on (w_2 - w_1)
   proj2 = (gmin + a - c)/sqrt(a-2*c+b);
     
  
   % --- compute stop condition for the alpha1 (1st class) ------

   % (proj1 < proj2) ~ the worst point will be used for update 
   if (proj1 < proj2) & (proj1 <= (sqrt(a-2*c+b) - epsilon)),
    
     % -- Adaptation phase of vector alpha1 ----------------------------
     
     k = (a - emin - c)/...
       (a+kernel(X1(:,inx1),X1(:,inx1),ker,arg)+kadd-2*(Di1(inx1)-Fi1(inx1)) );
  
     k = min( 1, k );
    
     s1 = s1 * (1-k);
     s1(inx1) = s1(inx1) + k;
    
     sol = 0; 
     
     % -------------------------------------------------------------
     a = a*(1-k)^2 + 2*(1-k)*k*Di1(inx1) + ...
          k^2 * (kernel(X1(:,inx1),X1(:,inx1),ker,arg)+kadd );
     
     c = c*(1-k) + k*Fi1(inx1);
     
     for i=1:l1,
       Di1(i) = Di1(i)*(1-k) + k*kernel(X1(:,i),X1(:,inx1),ker,arg);
     end
     Di1(inx1) = Di1(inx1) + k*kadd;
     for i=1:l2,
       Di2(i) = Di2(i)*(1-k) + k*kernel(X2(:,i),X1(:,inx1),ker,arg);    
     end
     
   else        
      % --- compute stop condition for the alpha2 (2st class) ------   

      if proj2 <= (sqrt(a-2*c+b) - epsilon ),

         % -- Adaptation phase ----------------------------------

        k = (b - gmin -c)/...
        (b+kernel(X2(:,inx2),X2(:,inx2),ker,arg)+kadd-2*(Fi2(inx2)-Di2(inx2)));

         k = min( 1, k );
    
         s2 = s2 * (1-k);
         s2(inx2) = s2(inx2) + k;

         sol = 0;

         % ------------------------------------------------------
         b = b*(1-k)^2 + 2*(1-k)*k*Fi2(inx2) + ...
          k^2 * (kernel(X2(:,inx2),X2(:,inx2),ker,arg)+kadd );

         c = c*(1-k) + k*Di2(inx2);
         
         for i=1:l1,
           Fi1(i) = (1-k)*Fi1(i) + k*kernel(X2(:,inx2),X1(:,i),ker,arg);
         end
         for i=1:l2,
           Fi2(i) = (1-k)*Fi2(i) + k*kernel(X2(:,inx2),X2(:,i),ker,arg);
         end
         Fi2(inx2) = Fi2(inx2) + k*kadd;
         
      end  
   end
   
   % store up=||w||/2 and current margin m(w1-w2,theta) = min( m1, m2) 
   m = min([proj1,proj2]) - 0.5*sqrt(a-2*c+b);
   up = [up,sqrt(a-2*c+b)/2];       
   lo = [lo,m ];   
   
%%   disp(sprintf('step = %d', t ));
   
end

if sol == 1 & (proj1 < 0 | proj2 < 0),
   sol = -2;
end

% -- compute threshold -----------------------

% sqared margin in transfromed space
margin2 = a - 2*c + b;

% threshold after normalization
theta = (a-b)/margin2;

% solution (normal vect. in the transformed space) after normalization
s1=2*s1/margin2;
s2=2*s2/margin2;

% -- make SVM-like output --------------------

Alpha=zeros(1,l1+l2);
Alpha(xinx1)=s1;
Alpha(xinx2)=s2;

bias = -theta;

% -- compute margin -----------------------------------
if nargout >= 6,
  margin = 0;
  for i=find(Alpha ~= 0),
    for j=find(Alpha ~= 0 ),
      margin = margin + Alpha(i)*Alpha(j)*itosgn(I(i))*...
               itosgn(I(j))* kernel(X(:,i),X(:,j),ker,arg);
    end
  end
  margin = 2/sqrt(margin);
end

% overall number of used float point operations
flps=flops;

return;