请问谁有smo算法程序
SVM中的序列最小最优化算法(Sequetial Minimal Optimization)，我看了Platt的文章，
上面有个伪代码，
感觉不是很清楚。
在网上搜了一下，也没有搜到相关代码。
请大虾指教。
谢谢。
就帖在后面吧，
我想应该不只我一个人需要才对。。
C的或者matlab的都可以。 SMO算法的Matlab实现
http://www.blog.edu.cn/user2/huangbo929/archives/2007/1713347.shtml
基于Gauss核函数的SVM原码（使用SMO算法作为优化算法）
程序由研学的luke提供
#include "iostream.h"
#include "fstream.h"
#include "ctype.h"
#include "math.h"
#include "stdlib.h"
#include "time.h"
#define N 14104//样本数总数
#define end_support_i 12090
#define first_test_i 12090
#define d 340//维数
//////////global variables
//int N= //样本数
//int d= //维数
float C=100;//惩罚参数
float tolerance=0.001;//ui与边界0,C之间的公差范围
float eps=1e-3;//一个近似0的小数
float two_sigma_squared=10;//RBF核函数中的参数
float alph;//Lagrange multipiers
float b;//threshold
float error_cache;//存放non-bound样本误差
int target;//训练与测试样本的目标值
float precomputed_self_dot_product;//预存dot_product_func(i,i)的值，以减少计算量
int dense_points;//存放训练与测试样本，0-end_support_i-1训练;first_test_i-N测试
//函数的申明
int takeStep(int,int);
int examineNonBound(int);
int examineBound(int);
int examineFirstChoice(int,float);
int examineExample(int);
float kernel_func(int,int);
float learned_func(int);
float dot_product_func(int,int);
float error_rate();
void setX();
void setT();
void initialize();
///////////////////////////////////////////////////////
void main()
{
ofstream os("d:\\smo1.txt");//存放训练的后的参数
int numChanged=0,examineAll=1;
//srand((unsigned int)time(NULL));
initialize();
//以下两层循环，开始时检查所有样本，选择不符合KKT条件的两个乘子进行优化，选择成功，返回1，否则，返回0
//所以成功了，numChanged必然>0，从第二遍循环时，就不从整个样本中去寻找不符合KKT条件的两个乘子进行优化，
//而是从边界的乘子中去寻找，因为非边界样本需要调整的可能性更大，边界样本往往不被调整而始终停留在边界上。
//如果，没找到，再从整个样本中去找，直到整个样本中再也找不到需要改变的乘子为止，此时，算法结束。
while(numChanged>0||examineAll)
{
numChanged=0;
if(examineAll)
{
for(int k=0;k<end_support_i;k++)
numChanged+=examineExample(k);//examin all exmples
}
else{
for(int k=0;k<end_support_i;k++)
if(alph!=0&&alph!=C)
numChanged+=examineExample(k);//loop k over all non-bound lagrange multipliers
}
if(examineAll==1)
examineAll=0;
else if(numChanged==0)
examineAll=1;
}
//存放训练后的参数
{
os<<"d="<<d<<endl;
os<<"b="<<b<<endl;
os<<"two_sigma_squared="<<two_sigma_squared<<endl;
int n_support_vectors=0;
for(int i=0;i<end_support_i;i++)
if(alph>0)//此地方是否可以改为alph>tolerance?????????????????
n_support_vectors++;
os<<"n_support_vectors="<<n_support_vectors<<"rate="<<(float)n_support_vectors/first_test_i<<endl;
for(i=0;i<end_support_i;i++)
if(alph>0)
os<<"alph["<<i<<"]="<<alph<<" ";
os<<endl;
}
//输出，以及测试
cout<<error_rate()<<endl;
}
/////////examineExample程序
//假定第一个乘子ai(位置为i1)，examineExample(i1)首先检查，如果它超出tolerance而违背KKT条件，那么它就成为第一个乘子；
//然后，寻找第二个乘子(位置为i2),通过调用takeStep(i1,i2)来优化这两个乘子。
int examineExample(int i1)
{
float y1,alph1,E1,r1;
y1=target;
alph1=alph;
if(alph1>0&&alph1<C)
E1=error_cache;
else
E1=learned_func(i1)-y1;//learned_func为计算输出函数
r1=y1*E1;
//////违反KKT条件的判断
if((r1>tolerance&&alph1>0)||(r1<-tolerance&&alph1<C))
{
/////////////使用三种方法选择第二个乘子
//1：在non-bound乘子中寻找maximum fabs(E1-E2)的样本
//2：如果上面没取得进展,那么从随机位置查找non-boundary 样本
//3：如果上面也失败，则从随机位置查找整个样本,改为bound样本
if (examineFirstChoice(i1,E1))//1
{
return 1;
}
if (examineNonBound(i1))//2
{
return 1;
}
if (examineBound(i1))//3
{
return 1;
}
}
///没有进展
return 0;
}
////////////1：在non-bound乘子中寻找maximum fabs(E1-E2)的样本
int examineFirstChoice(int i1,float E1)
{
int k,i2;
float tmax;
for(i2=-1,tmax=0.0,k=0;k<end_support_i;k++)//*******************************end_support_i
if(alph>0&&alph<C)
{
float E2,temp;
E2=error_cache;
temp=fabs(E1-E2);
if(temp>tmax)
{
tmax=temp;
i2=k;
}
}
if(i2>=0)
{
if(takeStep(i1,i2))
return 1;
}
return 0;
}
/////////////// 2：如果上面没取得进展,那么从随机位置查找non-boundary 样本
int examineNonBound(int i1)
{
int k,k0 = rand()%end_support_i;
int i2;
for (k = 0; k < end_support_i; k++)
{
i2 = (k + k0) % end_support_i;//从随机位开始
if (alph > 0.0 && alph < C)
{
if (takeStep(i1, i2))
{
return 1;
}
}
}
return 0;
}
////////////3：如果上面也失败，则从随机位置查找整个样本,(改为bound样本)
int examineBound(int i1)
{
int k,k0 = rand()%end_support_i;
int i2;
for (k = 0; k < end_support_i; k++)
{
i2 = (k + k0) % end_support_i;//从随机位开始
//if (alph= 0.0 || alph=C)//修改******************
{
if (takeStep(i1, i2))
{
return 1;
}
}
}
return 0;
}
///////////takeStep()
//用于优化两个乘子，成功，返回1，否则，返回0
int takeStep(int i1,int i2)
{
int y1,y2,s;
float alph1,alph2;//两个乘子的旧值
float a1,a2;//两个乘子的新值
float E1,E2,L,H,k11,k22,k12,eta,Lobj,Hobj,delta_b;
if(i1==i2) return 0;//不会优化两个同一样本
//给变量赋值
alph1=alph;
alph2=alph;
y1=target;
y2=target;
if(alph1>0&&alph1<C)
E1=error_cache;
else
E1=learned_func(i1)-y1;//learned_func(int)为非线性的评价函数，即输出函数
if(alph2>0&&alph2<C)
E2=error_cache;
else
E2=learned_func(i2)-y2;
s=y1*y2;
//计算乘子的上下限
if(y1==y2)
{
float gamma=alph1+alph2;
if(gamma>C)
{
L=gamma-C;
H=C;
}
else
{
L=0;
H=gamma;
}
}
else
{
float gamma=alph1-alph2;
if(gamma>0)
{
L=0;
H=C-gamma;
}
else
{
L=-gamma;
H=C;
}
}//计算乘子的上下限
if(L==H) return 0;
//计算eta
k11=kernel_func(i1,i1);//kernel_func(int,int)为核函数
k22=kernel_func(i2,i2);
k12=kernel_func(i1,i2);
eta=2*k12-k11-k22;
if(eta<0)
{
a2=alph2+y2*(E2-E1)/eta;//计算新的alph2
//调整a2，使其处于可行域
if(a2<L) a2=L;
if(a2>H) a2=H;
}
else//此时，得分别从端点H,L求目标函数值Lobj,Hobj，然后设a2为求得最大目标函数值的端点值
{
float c1=eta/2;
float c2=y2*(E2-E1)-eta*alph2;
Lobj=c1*L*L+c2*L;
Hobj=c1*H*H+c2*H;
if(Lobj>Hobj+eps)//eps****************
a2=L;
else if(Lobj<Hobj-eps)
a2=H;
else
a2=alph2;//加eps的目的在于，使得Lobj与Hobj尽量分开，如果，很接近，就认为没有改进(make progress)
}
if(fabs(a2-alph2)<eps*(a2+alph2+eps))
return 0;
a1=alph1-s*(a2-alph2);//计算新的a1
if(a1<0)//调整a1,使其符合条件*****??????????????????????????????????????????
{
a2+=s*a1;
a1=0;
}
else if(a1>C)
{
float t=a1-C;
a2+=s*t;
a1=C;
}
//更新阀值b
{
float b1,b2,bnew;
if(a1>0&&a1<C)
bnew=b+E1+y1*(a1-alph1)*k11+y2*(a2-alph2)*k12;
else{
if(a2>0&&a2<C)
bnew=b+E2+y1*(a1-alph1)*k12+y2*(a2-alph2)*k22;
else{
b1=b+E1+y1*(a1-alph1)*k11+y2*(a2-alph2)*k12;
b2=b+E2+y1*(a1-alph1)*k12+y2*(a2-alph2)*k22;
bnew=(b1+b2)/2;
}
}
delta_b=bnew-b;
b=bnew;
}
//对于线性情况，要更新权向量，这里不用了
//更新error_cache，对取得进展的a1,a2,所对应的i1,i2的error_cache=error_cache=0
{
float t1=y1*(a1-alph1);
float t2=y2*(a2-alph2);
for(int i=0;i<end_support_i;i++)
if(0<alph&&alph<C)
error_cache+=t1*kernel_func(i1,i)+t2*(kernel_func(i2,i))-delta_b;
error_cache=0.0;
error_cache=0.0;
}
alph=a1;//store a1,a2 in the alpha array
alph=a2;
return 1;//说明已经取得进展
}
//////////learned_func(int)
//评价分类学习函数
float learned_func(int k)
{
float s=0.0;
for(int i=0;i<end_support_i;i++)
if(alph>0)
s+=alph*target*kernel_func(i,k);
s-=b;
return s;
}
//计算点积函数dot_product_func(int,int)
float dot_product_func(int i1,int i2)
{
float dot=0;
for(int i=0;i<d;i++)
dot+=dense_points*dense_points;
return dot;
}
//径向机核函数RBF:kernel_func(int,int)
float kernel_func(int i1,int i2)
{
float s=dot_product_func(i1,i2);
s*=-2;
s+=precomputed_self_dot_product+precomputed_self_dot_product;
return exp(-s/two_sigma_squared);
}
//初始化initialize()
void initialize()
{
//初始化阀值b为0
b=0.0;
//初始化alph[]为0
for(int i=0;i<end_support_i;i++)
alph=0.0;
//设置样本值矩阵
setX();
//设置目标值向量
setT();
//设置预计算点积
{for(int i=0;i<N;i++)
precomputed_self_dot_product=dot_product_func(i,i);//注意，这是对训练样本的设置，对于测试样本也应考虑?????????????????
}
}
//计算误差率error_rate()
float error_rate()
{ ofstream to("d:\\smo_test.txt");
int tp=0,tn=0,fp=0,fn=0;
float ming=0,te=0,total_q=0,temp=0;
for(int i=first_test_i;i<N;i++)
{ temp=learned_func(i);
if(temp>0&&target>0)
tp++;
else if(temp>0&&target<0)
fp++;
else if(temp<0&&target>0)
fn++;
else if(temp<0&&target<0)
tn++;
to<<i<<"实际输出"<<temp<<endl;
}
total_q=(float)(tp+tn)/(float)(tp+tn+fp+fn);//总精度
ming=(float)tp/(float)(tp+fn);
te=(float)tp/(float)(tp+fp);
to<<"---------------测试结果-----------------"<<endl;
to<<"tp="<<tp<<" tn="<<tn<<"fp="<<fp<<"fn="<<fn<<endl;
to<<"ming="<<ming<<"te="<<te<<"total_q="<<total_q<<endl;
return (1-total_q);
}
//计算样本X[]
void setX(){
ifstream ff("D:\\17_smo.txt");////////////////////////
if(!ff)
{
cerr<<"error!!"<<endl;
exit(1);
}
int i=0,j=0;
char ch;
while(ff.get(ch))
{
if(isspace(ch)) continue;
if(isdigit(ch))
{
dense_points=(int)ch-48;
j++;
if(j==d)
{
j=0;
i++;
}
}
}
}
//set targetT[]
void setT()
{
for(int i=0;i<6045;i++)
target=1;
for(i=6045;i<12090;i++)
target=-1;
for(i=12090;i<13097;i++)
target=1;
for(i=13097;i<14104;i++)
target=-1;
} 另外有人将其改为IDL语言，http://hi.baidu.com/hilbertspace/blog/item/febd1dd5a8a738c550da4b70.html
代码如下：
实现了硬最大线性间隔分界(最简单的情况)
;;;;;SVM的SMO优化实现
;;;;借用研学论坛上luke的代码
;;Introduction: x: 二维样本数据.
;;;;;;;;;;;;;;;;n: 数据规模
;;;;;;;;;;;;;;;;y; 数据标志
;;;;;;;;;;;;;;;;w; 边界线斜率
;;;;;;;;;;;;;;;;b: 边界线截距
;;;;;;;;;;;;;;;;所以 y=wx+b 是边界线
PRO svmlineartest
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;generate the training set;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
x=fltarr(2,40)
y=intarr(40)
x=randomu(1.5,40)
x=randomu(2.5,40)+0.5
y=-1
y=1
n=40
plot,x,x,psym=1,xrange=,yrange=
oplot,x,x,psym=1,color=255
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;finished generate the training set;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;minimise w^2, subject to: yi*(wxi+b) >= 1, then it can be converted to
;;;;;;sum wi*alphi=0, & , w- sum yi*alphi*xi =0 , then get
;;;;;;maxmize L= sum alphi-0.5* sum (yi*yj*alphi*alphj*<xi.*xj> subject to: sum yi*alphi =0, ai>=0
;;;;;;sumj alphi*yi*yj*xi,*xj=2
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
w=fltarr(2)
w=svmw(x,y,n)
END
Function svmw,x,y,n
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;intialization;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;使用指针方便全局操作;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
alph=fltarr(n)
C=100
tol=0.001
eps=1e-3
error_cache=fltarr(n)
alphptr=ptr_new(alph)
b=0.0
bptr=ptr_new(b)
errorptr=ptr_new(error_cache)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;Finish intialization;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
numChange=0
examineAll=1
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;以下两层循环，开始时检查所有样本，选择不符合KKT条件的两个乘子进行优化
;;;;;;;选择成功，返回1，否则，返回
;;;;;;;成功了，numChanged必然>0，从第二遍循环时，就不从整个样本中去寻找不符合KKT条件的两个乘子进行优化，
;;;;;;;而是从边界的乘子中去寻找，因为非边界样本需要调整的可能性更大，边界样本往往不被调整而始终停留在边界上。
;;;;;;;如果，没找到，再从整个样本中去找，直到整个样本中再也找不到需要改变的乘子为止，此时，算法结束。
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
while ((numChange gt 0) or examineAll) do begin
numChanged =0
if (examineAll) then begin
for i=0,n-1 do begin
numChanged+=examineExample(i,x,y,alphptr,n,bptr,errorptr)
endfor
endif
if (examineAll eq 1) then examineAll=0 else if (numChanged eq 0) then examineAll=1
endwhile
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;计算w,b并画出图;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
w=fltarr(2)
w= total((*alphptr)*x*y)
w= total((*alphptr)*x*y)
b=-(max(w*x+w*x)+min(w*x+w*x))/2
print,w
print,b
x1new=findgen(100)*0.01*2
ynew=-(w*x1new+b)/w
oplot, x1new,ynew
End
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;takestep用于优化两个乘子，成功，返回1，否则，返回0;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
Function takestep,i1,i2,x,y,alphptr,n,bptr,errorptr
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;取出乘子和x,y,定义C值;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
alph=*alphptr
b=*bptr
C=100
tol=0.001
eps=1e-3
error_cache=*errorptr
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;完成取出乘子和x,y;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if (i1 eq i2) then return, 0
alph1 = alph
y1 = y
alph2 = alph
y2 = y
s=y1*y2
if ((alph1 gt 0) and (alph1 lt C)) then E1=error_cache else E1=learned_func(i1,x,y,alph,n)-y1;//learned_func(int)为非线性的评价函数，即输出函数
if ((alph2 gt 0) and (alph2 lt C)) then E2=error_cache else E2=learned_func(i2,x,y,alph,n)-y2;
if y1 ne y2 then begin
L= 0 > (alph-alph)
H= C < (C+alph-alph)
endif
if y1 eq y2 then begin
L= 0 > (alph+alph-C)
H= C < (alph+alph)
endif
if (L eq H) then return, 0
k11=kernel_func(i1,i1,x,y,n);//kernel_func(int,int)为核函数
k22=kernel_func(i2,i2,x,y,n);
k12=kernel_func(i1,i2,x,y,n)
eta = 2*k12-k11-k22
if (eta lt 0) then begin
a2=alph2 - y2*(E1-E2)/eta
if (a2 lt L) then a2=L else if (a2 gt H) then a2=H
endif else begin
c1=eta/2;
c2=y2*(E2-E1)-eta*alph2;
Lobj=c1*L*L+c2*L
Hobj=c1*H*H+c2*H
if (Lobj gt (Hobj+eps)) then begin
a2=L
endif else if (Lobj lt Hobj-eps) then begin
a2=H
endif else begin
a2=alph2
endelse
endelse
if (abs(a2-alph2) lt eps*(a2+alph2+eps)) then return, 0
a1=alph1+s*(alph2-a2)
if (a1 lt 0) then begin
a2+=s*a1
a1=0
endif else if (a1 gt C) then begin
t=a1-C
a2+=s*t
a1=C
endif
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;更新阈值b;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if (a1 gt 0) and (a1 lt C) then begin
bnew=b+E1+y1*(a1-alph1)*k11+y2*(a2-alph2)*k12;
endif else if (a2 gt 0) and (a2 lt C) then begin
bnew=b+E2+y1*(a1-alph1)*k12+y2*(a2-alph2)*k22
endif else begin
b1=b+E1+y1*(a1-alph1)*k11+y2*(a2-alph2)*k12
b2=b+E2+y1*(a1-alph1)*k12+y2*(a2-alph2)*k22
bnew=(b1+b2)/2
endelse
delta_b=bnew-b
b=bnew
*bptr=b
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;完成更新阈值b;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Update weight vector to reflect change in a1&a2,if linear SVM;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;更新error cache;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
t1=y1*(a1-alph1);
t2=y2*(a2-alph2);
for i=0,n-1 do begin
if (0 lt alph and alph lt C) then begin
error_cache+=t1*kernel_func(i1,i,x,y,n)+t2*(kernel_func(i2,i,x,y,n))-delta_b;
error_cache=0.0;
error_cache=0.0;
endif
endfor
*errorptr=error_cache
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;完成更新error cache;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;更新lagrange乘子;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
alph=a1;Store a1 in the alpha array
alph=a2;Store a2 in the alpha array
*alphptr=alph
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;完成更新lagrange乘子;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
return, 1
End
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;examineExample程序;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;假定第一个乘子ai(位置为i1)，examineExample(i1)首先检查;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;如果它超出tolerance而违背KKT条件，那么它就成为第一个乘子;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;然后，寻找第二个乘子(位置为i2),通过调用takeStep(i1,i2)来优化这两个乘子;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
Function examineExample,i1,x,y,alphptr,n,bptr,errorptr
alph=(*alphptr)
tol=0.001
C=100
eps=1e-3
error_cache=*errorptr
y1=y
alph1=alph
if ((alph1 gt 0) and (alph1 lt C)) then E1=error_cache else E1=learned_func(i1,x,y,alph,n)-y1;//learned_func为计算输出函数
r1=y1*E1
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;违反KKT条件的判断;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if (((r1 lt -tol) and (alph1 lt C)) or (( r1 gt tol) and (alph1 gt 0))) then begin
;if(size(where (alph ne 0)) gt 1 ) then begin
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;使用三种方法选择第二个乘子;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;1：在non-bound乘子中寻找maximum fabs(E1-E2)的样本;;;;;;;;;;;;;;
;;;;;;;;;;2：如果上面没取得进展,那么从随机位置查找non-boundary 样本;;;;;;
;;;;;;;;;;3：如果上面也失败，则从随机位置查找整个样本,改为bound样本;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if (examineFirstChoice(i1,E1,x,y,alphptr,n,bptr,errorptr)) then return, 1
if (examineNonBound(i1,x,y,alphptr,n,bptr,errorptr)) then return, 1
if (examineBound(i1,x,y,alphptr,n,bptr,errorptr)) then return, 1
endif
;endif
return, 0
End
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;1：在non-bound乘子中寻找maximum fabs(E1-E2)的样本;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
Function examineFirstChoice,i1,E1,x,y,alphptr,n,bptr,errorptr
C=100
alph=(*alphptr)
error_cache=*errorptr
for k=0,n-1 do begin
i2=-1
tmax=0.0
if(alph gt 0) and (alph lt C) then begin
E2=error_cache;
temp=abs(E1-E2);
if (temp gt tmax) then begin
tmax=temp;
i2=k;
endif
endif
endfor
if (i2 ge 0) then begin
if (takeStep(i1,i2,x,y,alphptr,n,bptr,errorptr)) then return, 1
endif
return, 0
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;2：如果上面没取得进展,那么从随机位置查找non-boundary 样本;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function examineNonBound,i1,x,y,alphptr,n,bptr,errorptr
C=100
alph=(*alphptr)
error_cache=*errorptr
k0 = randomu(seed,1) mod n-1
for k = 0, n-1 do begin
i2 = (k + k0) mod n-1;//从随机位开始
if (alph gt 0.0) and (alph lt C) then begin
if (takeStep(i1, i2,x,y,alphptr,n,bptr,errorptr)) then begin
return, 1;
endif
endif
endfor
return,0
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;如果上面也失败，则从随机位置查找整个样本,(改为bound样本);;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function examineBound,i1,x,y,alphptr,n,bptr,errorptr
C=100
alph=(*alphptr)
error_cache=*errorptr
k0 = randomu(seed,1) mod (n-1);
for k = 0, n-1 do begin
i2 = (k + k0) mod (n-1);//从随机位开始
if (alph eq 0.0 )or (alph eq C) then begin
if (takeStep(i1, i2,x,y,alphptr,n,bptr,errorptr)) then return, 1
endif
endfor
return,0
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;compute the learn function(y=wx+b) //test;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
Function learned_func ,i,x,y,alphlinshi,n
s=fltarr(n)
for j=0,n-1 do begin
s=alphlinshi*alphlinshi*y*y*(x[*,i] ## transpose(x[*,j]))
endfor
learned_func=total(s)
return, learned_func
End
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;核函数kernel_func(int,int);;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function kernel_func,i1,i2,x,y,n
func=(x[*,i1] ## transpose(x[*,i2]))
return,func
;float s=dot_product_func(i1,i2);
;s*=-2;
;s+=precomputed_self_dot_product+precomputed_self_dot_product;
;return exp(-s/two_sigma_squared);
end 太感谢了。 ((em:08))非常感谢 很有借鉴性 什么都不说了，感谢的一塌糊涂 谢谢楼主无私分享，很好很强大！！
页:
[1]