smo算法matlab实现

Posted 2020-08-19 holy_black_cat

看完CSDN上结构之法，算法之道的支持向量机通俗导论（理解SVM的三层境界）

http://blog.csdn.net/v_july_v/article/details/7624837

    参考了台湾的林智仁教授写了一个封装SVM算法的libsvm库，下载地址：

http://www.csie.ntu.edu.tw/~cjlin/libsvm/，此外下载了一份libsvm的注释文档，下载地址：

http://www.pami.sjtu.edu.cn/people/gpliu/document/libsvm_src.pdf

 

    SVM的原理在三层境界里已经讲的很清楚了，然而SMO算法的实现却仍然困惑着很多人。机器学习课程中有个作业：

“1. 需要Matlab2010b以上版本的运行环境；

2. my_svm.m 与my_svmtrain.m 是两个与此作业相关的文件；

3. 请尝试设计一个序列最优化函数，替代原函数自带的seqminopt函数。my_svmtrain.m中对应的代码部分如下：

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

[alpha bias] = seqminopt(training, groupIndex, ...

        boxconstraint, tmp_kfun, smo_opts);

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

”

 这个函数就是SMO算法源码。函数第一个输入参数的意义依次是分类数据，分类label，参数b（bias），核函数和smo参数（包括tol，MaxIter）。

   要SMO算法，必须知道算法的流程，简单来说就是个二次规划问题。

   for i=1:iter

a. 根据预先设定的规则，从所有样本中选出两个

b. 保持其他拉格朗日乘子不变，更新所选样本对应的拉格朗日乘子

end

下面解这个只有两个变量的二次规划问题：

步骤二计算二阶导数也可以写成

    按照算法的步骤，写下SMO算法的实现代码。之前有看到一份python实现的代码

http://www.tuicool.com/articles/RRZvYb

 

    下面给出matlab代码，有错误请指正：

function [alphas offset] = my_seqminopt(data, targetLabels, boxConstraints, ...

    kernelFunc, smoOptions)

%kernelFunc is not used because

%initialization

tol=smoOptions.TolKKT;

maxIter=smoOptions.MaxIter;

m=size(data,1);

alphas = zeros(size(data,1), 1);

itCount=1;

offset = 0;

%iteration

while itCount

    for i=1:m

        num_alpha_change=0;

        gI=(alphas.*targetLabels)‘ *(data*data(i,:)‘)+offset;

        %violate the KKT condition

        eI=gI-targetLabels(i);

        if ((targetLabels(i)*eI<-tol)&& (alphas(i)tol) && (alphas(i)>0) )

            j=randi([i,m]);%select randomly

            gJ=(alphas.*targetLabels)‘ *(data*data(j,:)‘)+offset;

            eJ=gJ-targetLabels(j);

            

            %oldAlphaI,oldAlphaJ

            oldAlphaI=alphas(i);

            oldAlphaJ=alphas(j);

            

            if targetLabels(i)~=targetLabels(j)

                L=max(0,alphas(j)-alphas(i));

                H=min(boxConstraints(i),boxConstraints(i)+alphas(j)-alphas(i));

            else

                L=max(0,alphas(j)+alphas(i)-boxConstraints(i));

                H=min(boxConstraints(i),alphas(j)+alphas(i));

            end

            

            eta=2.0*data(i,:)*data(j,:)‘-data(i,:)*data(i,:)‘-data(j,:)*data(j,:)‘;

            if eta>=0

                continue;

            end

            alphas(j)=alphas(j)-targetLabels(j)*(eI-eJ)/eta;

            if alphas(j)

                alphas(j)=L;

            elseif alphas(j)>H

                    alphas(j)=H;

            end

             

            alphaChange=alphas(j)-oldAlphaJ;

            if abs(alphaChange)<1e-5

                continue;

            end

            alphas(i)=alphas(i)+targetLabels(j)*targetLabels(i)*(oldAlphaJ-alphas(j));

            b1=boxConstraints(i)-eI-targetLabels(i)*(alphas(i)-oldAlphaI)*data(i,:)*data(i,:)‘-targetLabels(j)*(alphas(j)-oldAlphaJ)*data(i,:)*data(j,:)‘;

            b2=boxConstraints(i)-eJ-targetLabels(i)*(alphas(i)-oldAlphaI)*data(i,:)*data(j,:)‘-targetLabels(j)*(alphas(j)-oldAlphaJ)*data(j,:)*data(j,:)‘;

            if ((alphas(i)>0) && (alphas(i)

                offset=b1;

            elseif ((alphas(j)>0) && (alphas(j)

                offset=b2;

            else

                offset=(b1+b2)/2.0;    

            end

            num_alpha_change=num_alpha_change+1;

            

        end

    end

        if num_alpha_change == 0

            itCount=itCount+1;

        else

            itCount=0;

        end

end