smo1forsvr.m

一种快速的SVM算法。初学者有用。MATLAB程序。

function [beta, bias] = smo1forsvr(X,Y,ker,C,loss,e,tau)
%SVR Support Vector Regression
%
%  Usage: [nsv beta bias] = svr(X,Y,ker,C,loss,e)
%
%  Parameters: X      - Training inputs
%              Y      - Training targets
%              ker    - kernel function
%              C      - upper bound (non-separable case)
%              loss   - loss function
%              e      - insensitivity
%              nsv    - number of support vectors
%              beta   - Difference of Lagrange Multipliers
%              bias   - bias term
%              tau    - violating pair's tolerance parameter
%  Author: Steve Gunn ([e-mail])


  %if (nargin < 3 | nargin > 7) % check correct number of arguments
    %help svr
%  else

    fprintf('Sequence Minimization Optimization for SVR ....\n')
    fprintf('______________________________\n')
    
    %if (nargin<6) e=0.0;, end
    %if (nargin<5) loss='eInsensitive';, end
    %if (nargin<4) C=Inf;, end
    %if (nargin<3) ker='linear';, end  

    % tolerance for Support Vector Detection
    epsilon = svtol(C);

    % Construct the Kernel matrix
    
    
    %fprintf('Constructing ...\n');
    
    n = size(X,1);
    H = zeros(n,n);  
    for i=1:n
       for j=1:n
          H(i,j) = svkernel(ker,X(i,:),X(j,:));
       end
    end

    beta=zeros(n, 1);
    k=0;
    HH=zeros(n, 1);
    
    
    for i=1:n
        HH(i)=-Y(i)+H(i,:)*beta;
    end
    
    hplus=HH+e;
    hminus=HH-e;
    
    if size(find(beta>0 & beta<=C), 1)==0
       Max=max(hminus(find(beta>-C & beta<=0)));
    elseif size(find(beta>-C & beta<=0), 1)==0
       Max=max(hplus(find(beta>0 & beta<=C)));
    elseif size(find(beta>0 & beta<=C), 1)~=0 & size(find(beta>-C & beta<=0), 1)~=0
       Max=max(max(hplus(find(beta>0 & beta<=C))), max(hminus(find(beta>-C & beta<=0))));
    end
    
    if size(find(beta>=0 & beta<C), 1)==0
       Min=min(hminus(find(beta>=-C & beta<0)));
    elseif size(find(beta>=-C & beta<0), 1)==0
       Min=min(hplus(find(beta>=0 & beta<C)));
    elseif size(find(beta>=0 & beta<C), 1)~=0 & size(find(beta>=-C & beta<0), 1)~=0
       Min=min(min(hplus(find(beta>=0 & beta<C))), min(hminus(find(beta>=-C & beta<0))));
    end
    
       
    while Max-tau > Min;
        Max
        Min
        for i=1:n
            for j=[i+1:n 1:i-1] 
                flag=0;
                eta=H(i,i)+H(j,j)-2*H(i,j);
            if beta(i)>-C & beta(i)<=C & beta(j)>=-C & beta(j)<C
                if beta(i)>-C & beta(i)<=0
                    hup=hminus(i);
                else
                    hup=hplus(i);
                end
                if beta(j)>=-C & beta(j)<0
                    hlow=hminus(j);
                else
                    hlow=hplus(j);
                end
                hup
                hlow
                if hup > hlow+tau
                        fprintf('found a violating pair...........\n')
                        [beta(i), beta(j)] = ntonplusone(beta(i),beta(j),HH(i), HH(j), eta,e,tau,C); %beta(k+1)=...beta(k)
                        beta
                        flag=1;
                end
            end
            if flag==0
               if beta(i)>=-C & beta(i)<C & beta(j)>-C & beta(j)<=C       
                  if beta(i)>=-C & beta(i)<0
                     hlow=hminus(i);
                  else
                     hlow=hplus(i);
                  end
                  if beta(j)>-C & beta(j)<=0
                    hup=hminus(j);
                  else
                    hup=hplus(j);
                  end
                  hlow
                  hup
                  if hlow < hup-tau
                          fprintf('found a violating pair...........\n')
                          [beta(i), beta(j)] = ntonplusone(beta(i),beta(j),HH(i), HH(j), eta,e,tau,C); %beta(k+1)=...beta(k)
                          beta
                          flag=1;
                          %break;
                  end
               end
            end
           
              for i=1:n
                 HH(i)=-Y(i)+H(i,:)*beta;
              end
    
              hplus=HH+e;
              hminus=HH-e;
    
              if size(find(beta>0 & beta<=C), 1)==0
                 Max=max(hminus(find(beta>-C & beta<=0)));
              elseif size(find(beta>-C & beta<=0), 1)==0
                 Max=max(hplus(find(beta>0 & beta<=C)));
              elseif size(find(beta>0 & beta<=C), 1)~=0 & size(find(beta>-C & beta<=0), 1)~=0
                 Max=max(max(hplus(find(beta>0 & beta<=C))), max(hminus(find(beta>-C & beta<=0))));
              end
     
              if size(find(beta>=0 & beta<C), 1)==0
                 Min=min(hminus(find(beta>=-C & beta<0)));
              elseif size(find(beta>=-C & beta<0), 1)==0
                 Min=min(hplus(find(beta>=0 & beta<C)));
              elseif size(find(beta>=0 & beta<C), 1)~=0 & size(find(beta>=-C & beta<0), 1)~=0
                 Min=min(min(hplus(find(beta>=0 & beta<C))), min(hminus(find(beta>=-C & beta<0))));
              end
              
              Max
              Min
              
              if Max-tau <= Min
                  break;
              end
              
            end    %j
            
            if Max-tau <= Min
                break;
            end
           % break;                                        
        end    %i    
                       
       %break;   
    end        %while
    
        
    % Implicit bias, b0
    bias = 0;

    % Explicit bias, b0 
    if nobias(ker) ~= 0
      switch lower(loss)
        case 'einsensitive',
          % find bias from average of support vectors with interpolation error e
          % SVs with interpolation error e have alphas: 0 < alpha < C
          svii = find( abs(beta) > epsilon & abs(beta) < (C - epsilon));
          if length(svii) > 0
            bias = (1/length(svii))*sum(Y(svii) - e*sign(beta(svii)) - H(svii,svi)*beta(svi));
          else 
            fprintf('No support vectors with interpolation error e - cannot compute bias.\n');
            bias = (max(Y)+min(Y))/2;
          end
        case 'quadratic',
            bias = mean(Y - H*beta);
      end 
    end

  end