svc.m

SVM 的代码,工具箱的使用内部有详细的说明,请读者看看Matlab代码中参数使用,将训练数据以矩阵形式输入,即可使用

function [nsv, alpha, b0] = svc(X,Y,ker,C)
%SVC Support Vector Classification
%
%  Usage: [nsv alpha bias] = svc(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: Steve Gunn ([e-mail])

  if (nargin <2 | nargin>4) % check correct number of arguments
    help svc
  else

    fprintf('Support Vector Classification\n')
    fprintf('_____________________________\n')
    n = size(X,1); %n是训练样本的数目
    if (nargin<4) C=Inf;, end %如果不给出C的值则取默认值 无穷大
    if (nargin<3) ker='linear';, end

    % tolerance for Support Vector Detection
    epsilon = svtol(C); %epsilon是一个很小的值，与C有关
    
    % Construct the Kernel matrix
    %生成H矩阵
    fprintf('Constructing ...\n');
    H = zeros(n,n);  
    for i=1:n
       for j=1:n
          H(i,j) = Y(i)*Y(j)*svkernel(ker,X(i,:),X(j,:));
       end
    end
    %优化问题中的c*alpha c为1，即第二项
    c = -ones(n,1);  

    % Add small amount of zero order regularisation to 
    % avoid problems when Hessian is badly conditioned. 
    %对H矩阵的对角部分加以扰动，以防止无法解出的情况
    H = H+1e-10*eye(size(H));
    
    % Set up the parameters for the Optimisation problem
    %0<=alpha<=C
    vlb = zeros(n,1);      % Set the bounds: alphas >= 0
    vub = C*ones(n,1);     %                 alphas <= C
    %初始搜索点
    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,此处实际是指alpha'*Y=0,优化中不需要给出alpha,只给出系数
    else
       A = [];, b = [];
    end

    % Solve the Optimisation Problem
    
    fprintf('Optimising ...\n');
    st = cputime;%取得当前cpu时间
    %解这个优化问题，其表达式大概是
    %  min 1/2(alpha'* H alpha)+c'alpha
    %  s.t.  alpha'Y=0
    %        C>=alpha>=0
    %注意：1/2虽然公式中有，但qp中并没有，这不是qp的错误，其内部考虑了这个问题.
    %关于这一问题，我曾问过作者，他说的。
    [alpha lambda how] = qp(H, c, A, b, vlb, vub, x0, neqcstr);
    
    fprintf('Execution time: %4.1f seconds\n',cputime - st);
    fprintf('Status : %s\n',how);
    %计算w^2
    w2 = alpha'*H*alpha;
    fprintf('|w0|^2    : %f\n',w2);
    %margin表示最优分类面最近点之间的距离,当然是越大越好
    fprintf('Margin    : %f\n',2/sqrt(w2));
    fprintf('Sum alpha : %f\n',sum(alpha));
    
        
    % Compute the number of Support Vectors
    %理论上alpha>0的就可以认为是sv,但实际中因为是数值解法，不可能做到理论上的精确
    %要求alpha大于一个很小值就可以认为是sv.
    svi = find( alpha > epsilon);
    %nsv:sv的数目
    nsv = length(svi);
    fprintf('Support Vectors : %d (%3.1f%%)\n',nsv,100*nsv/n);

    % Implicit bias, b0
    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
      %对于不是线性可分的情况，只有0<alpha<C对应的支持向量才真正在边界上，
      %具体可参加我论文第32页,由于在分类边界上的支持向量满足等式：
      % alpha(y[(w,x)+b]-1)=0,把w的表达式代入就可得到b,这里实际是对所有分类边界上的
      %支持向量的结果求平均值。
      %为了领会这些问题，你可以用svcplot对二维的情况画图，以观察各种情况。
      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