multisvc.m

支持向量机 这是我的老师编写的matlab源文件 希望对大家有用

function  [nsv, alpha, b0,t] = multisvc (X,Y,ker,C)
%MULTISVC Support Vector Classification
%
%  Usage: [nsv alpha bias] = multisvc(X,Y,ker,C)
%
%  Parameters: X      - Training inputs
%              Y      - Training targets
%              ker    - kernel function
%              C      - upper bound (non-separable case)
%              nsv    - number of support vectors
%              alpha  - Lagrange Multipliers
%              b0     - bias term
%
%  Author: Ren weihua
 if (nargin <2 | nargin>5) % check correct number of arguments
    help multisvc
 else
      fprintf('Support Vector Multiclassification\n')
      fprintf('_____________________________\n')
      if (nargin<4) C=Inf;, end
      if (nargin<3) ker='linear';, end
       % tolerance for Support Vector Detection
       epsilon = svtol(C(end,end));
       % Construct the Kernel matrix
       fprintf('Constructing ...\n')
       t=0;
       svi=[];
       for i=0:9
           Yn=((Y==i)-(Y~=i));
           n = size(X,1);
           H = zeros(n,n);
           for i=1:n
               for j=1:n
                   H1(i,j)=svkernel(ker,X(i,:),X(j,:));
               end 
           end
           Dist = zeros(n,n); 
           for i=1:n
               for j=1:n
                   Dist=H1(i,i)+H(j,j)-2*H(i,j);
               end 
           end
           if size(Y,2)~=1
               Y=Y'
           end
           H = (Y'*Y).*H1;%*svkernel(ker,X(i,:),X(j,:));
           c = -ones(n,1); 
           Rmax=max(max(Dist));%2-2*sum(sum(H1))/(n*n);
           clear H1;
           % Add small amount of zero order regularisation to 
           % avoid problems when Hessian is badly conditioned. 
           H = H+1e-10*eye(size(H));
           % Set up the parameters for the Optimisation problem
           vlb = zeros(n,1);      % Set the bounds: alphas >= 0
           if size(C,1)==1 & size(C,2)==1     %  alphas <= C
               vub = C*ones(n,1);
           elseif size(C,1)==1
               vub=C';
           elseif size(C,2)==1
               vub=C;
           end
           x0 = zeros(n,1);       % The starting point is [0 0 0   0]
           neqcstr = nobias(ker); % Set the number of equality constraints (1 or 0)
           if neqcstr
               A = Y';, b = 0;     % Set the constraint Ax = b
           else
               A = [];, b = [];
           end
           % Solve the Optimisation Problem
           fprintf('Optimising ...\n');
           st = cputime;
           %[alpha lambda how] = qp(H, c, A, b, vlb, vub, x0, neqcstr);
           options=optimset('MaxIter',5000);
           [alpha,fval,exitflag]=quadprog(H,c,[],[],A,b,vlb,vub,x0,options);
           t=t+cputime - st;
           fprintf('Status : %d\n',exitflag);
           w2 = alpha'*H*alpha;
           fprintf('Margin    : %f\n',2/sqrt(w2));
           fprintf('Sum alpha : %f\n',sum(alpha));
           fprintf('the max R : %f,     %f\n',Rmax,Rmax*w2);
           nsv = length(svi);
           svi = [svi,find( alpha > epsilon)];
           b0 = 0;
           % Explicit bias, b0 
           if nobias(ker) ~= 0
               % find b0 from average of support vectors on margin
               % SVs on margin have alphas: 0 < alpha < C
               svii = find( alpha > epsilon & alpha < (C - epsilon));
               if length(svii) > 0
                   b0 =  (1/length(svii))*sum(Y(svii) - H(svii,svi)*alpha(svi).*Y(svii));
               else 
                   fprintf('No support vectors on margin - cannot compute bias.\n');
               end
           end
       end
       fprintf('Execution time: %4.1f seconds\n',t);
       % Compute the number of Support Vectors
       nsv = length(svi);
       fprintf('Support Vectors : %d (%3.1f%%)\n',nsv,100*nsv/n);
 end