smo_algorithm.m

核方法＆svm是模式识别是很重要的方法

function outPutPar = SMO_algorithm(class_one,class_two,kPar,C)
% the function for generating the two_class based classification
% InputPar:
% class_one: the first class data;
% class_two: the second class data;
% sizeTrainClassOne: the size of the training number of the first class
% data;
% sizeTrainClassTwo: the size of the training number of the second class
% data;
% OutputPar:
% outPutPar.w:
% outPutPar.b;
% outPutPar.alpha;
% Written by WangZhe on 2005-04-05;
size_class_one = size(class_one,1);
size_class_two = size(class_two,1);
Y=[ones(size_class_one,1);-1*ones(size_class_two,1)];
% input data completed

X=[class_one;class_two];

mat_kernel=(X*X'+1).^kPar;
% calculate the kernel use a kernel

trainSampleNum = size_class_one+size_class_two;
alpha=zeros(trainSampleNum,1);
b=0;
% set the value of alpha ,b and C

littleValue=10^(-7);
% the little value control the precison of this programj

% the allowed times that is against KKT conditions
loop_times=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

while(loop_times<5000)
    sup_index = find(alpha>littleValue & alpha<C-littleValue);
    E = (mat_kernel(:,sup_index)*(alpha(sup_index,1).*Y(sup_index,:))+b)-Y;
    E(abs(E)<littleValue)=0;
    % determine alpha1 and alpha2
    not_alpha1_index = [];
    for i = 1:trainSampleNum
        [alpha1,alpha1_index,success] = determine_alpha1(Y,mat_kernel,alpha,not_alpha1_index,b,C,littleValue);
        if success==1
            break
        end
        [alpha2,alpha2_index,alpha2_new] = determine_alpha2(mat_kernel,E,alpha,Y,littleValue,C,alpha1,alpha1_index);
        
        %can't find alpha2,to exchange alpha1
        if abs(alpha2_new-alpha2) < littleValue
            not_alpha1_index = [not_alpha1_index;alpha1_index];
        else
            break
        end
    end
    if success==1
        break
    end
    
    %optimization
    alpha1_new = alpha1+Y(alpha1_index)*Y(alpha2_index)*(alpha2-alpha2_new);
    alpha(alpha1_index) = alpha1_new;
    alpha(alpha2_index) = alpha2_new;
    b1_new = -E(alpha1_index)-Y(alpha1_index)*(alpha1_new-alpha1)*mat_kernel(alpha1_index,alpha1_index)-Y(alpha2_index)*...
        (alpha2_new-alpha2)*mat_kernel(alpha2_index,alpha1_index)+b;
    b2_new = -E(alpha2_index)-Y(alpha1_index)*(alpha1_new-alpha1)*mat_kernel(alpha1_index,alpha2_index)-Y(alpha2_index)*...
        (alpha2_new-alpha2)*mat_kernel(alpha2_index,alpha2_index)+b;
    b = (b1_new+b2_new)/2;
    if abs(b) < littleValue
        b = 0;
    end
    
    loop_times=loop_times+1;
end

outPutPar.numberSV=length(find(abs(alpha)>littleValue));
outPutPar.w=((alpha.*Y)'*X)';
outPutPar.b=b;
outPutPar.alpha=alpha;