请问谁有smo算法程序

jytan

SVM中的序列最小最优化算法(Sequetial Minimal Optimization)，
我看了Platt的文章，
上面有个伪代码，
感觉不是很清楚。
在网上搜了一下，也没有搜到相关代码。
请大虾指教。
谢谢。
就帖在后面吧，
我想应该不只我一个人需要才对。。
C的或者matlab的都可以。

风花雪月

SMO算法的Matlab实现

http://www.blog.edu.cn/user2/huangbo929/archives/2007/1713347.shtml

风花雪月

基于Gauss核函数的SVM原码（使用SMO算法作为优化算法）

程序由研学的luke提供

1. #include "iostream.h"
   
2. #include "fstream.h"
   
3. #include "ctype.h"
   
4. #include "math.h"
   
5. #include "stdlib.h"
   
6. #include "time.h"
   
7. 
   
8. #define N 14104//样本数总数
   
9. #define end_support_i 12090
   
10. #define first_test_i 12090
    
11. #define d 340//维数
    
12. //////////global variables
    
13. //int N=     //样本数
    
14. //int d=     //维数
    
15. float C=100;//惩罚参数
    
16. float tolerance=0.001;//ui与边界0,C之间的公差范围
    
17. float eps=1e-3;//一个近似0的小数
    
18. float two_sigma_squared=10;//RBF核函数中的参数
    
19. float alph[end_support_i];//Lagrange multipiers
    
20. float b;//threshold
    
21. float error_cache[end_support_i];//存放non-bound样本误差
    
22. int target[N];//训练与测试样本的目标值
    
23. float precomputed_self_dot_product[N];//预存dot_product_func(i,i)的值，以减少计算量
    
24. int dense_points[N][d];//存放训练与测试样本，0-end_support_i-1训练;first_test_i-N测试
    
25. 
    
26. //函数的申明
    
27. int takeStep(int,int);
    
28. int examineNonBound(int);
    
29. int examineBound(int);
    
30. int examineFirstChoice(int,float);
    
31. int examineExample(int);
    
32. float kernel_func(int,int);
    
33. float learned_func(int);
    
34. float dot_product_func(int,int);
    
35. float error_rate();
    
36. void setX();
    
37. void setT();
    
38. void initialize();
    
39. ///////////////////////////////////////////////////////
    
40. void main()
    
41. {
    
42.     ofstream os("d:\\smo1.txt");//存放训练的后的参数
    
43.     int numChanged=0,examineAll=1;
    
44.     //srand((unsigned int)time(NULL));
    
45.     initialize();
    
46.     //以下两层循环，开始时检查所有样本，选择不符合KKT条件的两个乘子进行优化，选择成功，返回1，否则，返回0
    
47.     //所以成功了，numChanged必然>0，从第二遍循环时，就不从整个样本中去寻找不符合KKT条件的两个乘子进行优化，
    
48.     //而是从边界的乘子中去寻找，因为非边界样本需要调整的可能性更大，边界样本往往不被调整而始终停留在边界上。
    
49.     //如果，没找到，再从整个样本中去找，直到整个样本中再也找不到需要改变的乘子为止，此时，算法结束。
    
50.     while(numChanged>0||examineAll)
    
51.     {
    
52.         numChanged=0;
    
53.         if(examineAll)
    
54.         {
    
55.             for(int k=0;k<end_support_i;k++)
    
56.                 numChanged+=examineExample(k);//examin all exmples
    
57.         }
    
58.         else{
    
59.             for(int k=0;k<end_support_i;k++)
    
60.                 if(alph[k]!=0&&alph[k]!=C)
    
61.                     numChanged+=examineExample(k);//loop k over all non-bound lagrange multipliers
    
62.             }
    
63.         if(examineAll==1)
    
64.             examineAll=0;
    
65.         else if(numChanged==0)
    
66.             examineAll=1;
    
67.     }
    
68.     //存放训练后的参数
    
69.     {
    
70.         os<<"d="<<d<<endl;
    
71.         os<<"b="<<b<<endl;
    
72.         os<<"two_sigma_squared="<<two_sigma_squared<<endl;
    
73.         int n_support_vectors=0;
    
74.         for(int i=0;i<end_support_i;i++)
    
75.             if(alph[i]>0)//此地方是否可以改为alph[i]>tolerance?????????????????
    
76.                 n_support_vectors++;
    
77.         os<<"n_support_vectors="<<n_support_vectors<<"rate="<<(float)n_support_vectors/first_test_i<<endl;
    
78.         for(i=0;i<end_support_i;i++)
    
79.             if(alph[i]>0)
    
80.             os<<"alph["<<i<<"]="<<alph[i]<<" ";
    
81.         os<<endl;
    
82.     }
    
83.     //输出，以及测试
    
84.     cout<<error_rate()<<endl;
    
85. }
    
86. 
    
87. /////////examineExample程序
    
88. //假定第一个乘子ai(位置为i1)，examineExample(i1)首先检查，如果它超出tolerance而违背KKT条件，那么它就成为第一个乘子；
    
89. //然后，寻找第二个乘子(位置为i2),通过调用takeStep(i1,i2)来优化这两个乘子。
    
90. int examineExample(int i1)
    
91. {
    
92.     float y1,alph1,E1,r1;
    
93.     y1=target[i1];
    
94.     alph1=alph[i1];
    
95.     if(alph1>0&&alph1<C)
    
96.         E1=error_cache[i1];
    
97.     else
    
98.         E1=learned_func(i1)-y1;//learned_func为计算输出函数
    
99.     r1=y1*E1;
    
100.     //////违反KKT条件的判断
     
101.     if((r1>tolerance&&alph1>0)||(r1<-tolerance&&alph1<C))
     
102.     {
     
103.         /////////////使用三种方法选择第二个乘子
     
104.         //1：在non-bound乘子中寻找maximum fabs(E1-E2)的样本
     
105.         //2：如果上面没取得进展,那么从随机位置查找non-boundary 样本
     
106.         //3：如果上面也失败，则从随机位置查找整个样本,改为bound样本
     
107.         if (examineFirstChoice(i1,E1))//1
     
108.              {
     
109.                      return 1;
     
110.               }
     
111. 
     
112.               if (examineNonBound(i1))//2
     
113.               {
     
114.                  return 1;
     
115.               }
     
116. 
     
117.               if (examineBound(i1))//3
     
118.              {
     
119.                  return 1;
     
120.               }
     
121.     }
     
122.     ///没有进展
     
123.           return 0;
     
124. }
     
125. 
     
126. ////////////1：在non-bound乘子中寻找maximum fabs(E1-E2)的样本
     
127. int examineFirstChoice(int i1,float E1)
     
128. {
     
129.    int k,i2;
     
130.    float tmax;
     
131.    for(i2=-1,tmax=0.0,k=0;k<end_support_i;k++)//*******************************end_support_i
     
132.        if(alph[k]>0&&alph[k]<C)
     
133.        {
     
134.            float E2,temp;
     
135.            E2=error_cache[k];
     
136.            temp=fabs(E1-E2);
     
137.            if(temp>tmax)
     
138.            {
     
139.                tmax=temp;
     
140.                i2=k;
     
141.            }
     
142.        }
     
143.    if(i2>=0)
     
144.    {
     
145.        if(takeStep(i1,i2))
     
146.            return 1;
     
147.    }
     
148.    return 0;
     
149. }
     
150. ///////////////    2：如果上面没取得进展,那么从随机位置查找non-boundary 样本
     
151. int examineNonBound(int i1)
     
152. {
     
153.        int k,k0 = rand()%end_support_i;
     
154.         int i2;
     
155.        for (k = 0; k < end_support_i; k++)
     
156.        {
     
157.               i2 = (k + k0) % end_support_i;//从随机位开始
     
158. 
     
159.               if (alph[i2] > 0.0 && alph[i2] < C)
     
160.               {
     
161.                  if (takeStep(i1, i2))
     
162.                  {
     
163.                     return 1;
     
164.               }
     
165.               }
     
166.        }
     
167.        return 0;
     
168. }
     
169. ////////////3：如果上面也失败，则从随机位置查找整个样本,(改为bound样本)
     
170. int examineBound(int i1)
     
171. {
     
172.        int k,k0 = rand()%end_support_i;
     
173.         int i2;
     
174.        for (k = 0; k < end_support_i; k++)
     
175.        {
     
176.               i2 = (k + k0) % end_support_i;//从随机位开始
     
177. 
     
178.               //if (alph[i2]= 0.0 || alph[i2]=C)//修改******************
     
179.               {
     
180.                  if (takeStep(i1, i2))
     
181.                  {
     
182.                     return 1;
     
183.               }
     
184.               }
     
185.        }
     
186.        return 0;   
     
187. }
     
188. ///////////takeStep()
     
189. //用于优化两个乘子，成功，返回1，否则，返回0
     
190. int takeStep(int i1,int i2)
     
191. {
     
192.     int y1,y2,s;
     
193.     float alph1,alph2;//两个乘子的旧值
     
194.     float a1,a2;//两个乘子的新值
     
195.     float E1,E2,L,H,k11,k22,k12,eta,Lobj,Hobj,delta_b;
     
196.    
     
197.     if(i1==i2) return 0;//不会优化两个同一样本
     
198.     //给变量赋值
     
199.     alph1=alph[i1];
     
200.     alph2=alph[i2];
     
201.     y1=target[i1];
     
202.     y2=target[i2];
     
203.     if(alph1>0&&alph1<C)
     
204.         E1=error_cache[i1];
     
205.     else
     
206.         E1=learned_func(i1)-y1;//learned_func(int)为非线性的评价函数，即输出函数
     
207.     if(alph2>0&&alph2<C)
     
208.         E2=error_cache[i2];
     
209.     else
     
210.         E2=learned_func(i2)-y2;
     
211.     s=y1*y2;
     
212.     //计算乘子的上下限
     
213.     if(y1==y2)
     
214.     {
     
215.         float gamma=alph1+alph2;
     
216.         if(gamma>C)
     
217.         {
     
218.             L=gamma-C;
     
219.             H=C;
     
220.         }
     
221.         else
     
222.         {
     
223.             L=0;
     
224.             H=gamma;
     
225.         }
     
226.     }
     
227.     else
     
228.     {
     
229.         float gamma=alph1-alph2;
     
230.         if(gamma>0)
     
231.         {
     
232.             L=0;
     
233.             H=C-gamma;
     
234.         }
     
235.         else
     
236.         {
     
237.             L=-gamma;
     
238.             H=C;
     
239.         }
     
240.     }//计算乘子的上下限
     
241.     if(L==H) return 0;
     
242.     //计算eta
     
243.     k11=kernel_func(i1,i1);//kernel_func(int,int)为核函数
     
244.     k22=kernel_func(i2,i2);
     
245.     k12=kernel_func(i1,i2);
     
246.     eta=2*k12-k11-k22;
     
247.     if(eta<0)
     
248.     {
     
249.         a2=alph2+y2*(E2-E1)/eta;//计算新的alph2
     
250.         //调整a2，使其处于可行域
     
251.         if(a2<L)     a2=L;
     
252.         if(a2>H)    a2=H;
     
253.     }
     
254.     else//此时，得分别从端点H,L求目标函数值Lobj,Hobj，然后设a2为求得最大目标函数值的端点值
     
255.     {
     
256.         float c1=eta/2;
     
257.         float c2=y2*(E2-E1)-eta*alph2;
     
258.         Lobj=c1*L*L+c2*L;
     
259.         Hobj=c1*H*H+c2*H;
     
260.         if(Lobj>Hobj+eps)//eps****************
     
261.             a2=L;
     
262.         else if(Lobj<Hobj-eps)
     
263.             a2=H;
     
264.         else
     
265.             a2=alph2;//加eps的目的在于，使得Lobj与Hobj尽量分开，如果，很接近，就认为没有改进(make progress)
     
266.     }
     
267.     if(fabs(a2-alph2)<eps*(a2+alph2+eps))
     
268.         return 0;
     
269.     a1=alph1-s*(a2-alph2);//计算新的a1
     
270.     if(a1<0)//调整a1,使其符合条件*****??????????????????????????????????????????
     
271.     {
     
272.         a2+=s*a1;
     
273.         a1=0;
     
274.     }
     
275.     else if(a1>C)
     
276.     {
     
277.         float t=a1-C;
     
278.         a2+=s*t;
     
279.         a1=C;
     
280.     }
     
281.     //更新阀值b
     
282.     {
     
283.         float b1,b2,bnew;
     
284.         if(a1>0&&a1<C)
     
285.             bnew=b+E1+y1*(a1-alph1)*k11+y2*(a2-alph2)*k12;
     
286.         else{
     
287.             if(a2>0&&a2<C)
     
288.                 bnew=b+E2+y1*(a1-alph1)*k12+y2*(a2-alph2)*k22;
     
289.             else{
     
290.                 b1=b+E1+y1*(a1-alph1)*k11+y2*(a2-alph2)*k12;
     
291.                 b2=b+E2+y1*(a1-alph1)*k12+y2*(a2-alph2)*k22;
     
292.                 bnew=(b1+b2)/2;
     
293.             }
     
294.         }
     
295.         delta_b=bnew-b;
     
296.         b=bnew;
     
297.     }
     
298.    
     
299.     //对于线性情况，要更新权向量，这里不用了
     
300.     //更新error_cache，对取得进展的a1,a2,所对应的i1,i2的error_cache[i1]=error_cache[i2]=0
     
301.     {
     
302.         float t1=y1*(a1-alph1);
     
303.         float t2=y2*(a2-alph2);
     
304.         for(int i=0;i<end_support_i;i++)
     
305.         if(0<alph[i]&&alph[i]<C)
     
306.             error_cache[i]+=t1*kernel_func(i1,i)+t2*(kernel_func(i2,i))-delta_b;
     
307.             error_cache[i1]=0.0;   
     
308.             error_cache[i2]=0.0;
     
309.     }
     
310.     alph[i1]=a1;//store a1,a2 in the alpha array
     
311.     alph[i2]=a2;
     
312.     return 1;//说明已经取得进展
     
313. }
     
314. 
     
315. //////////learned_func(int)
     
316. //评价分类学习函数
     
317. float learned_func(int k)
     
318. {
     
319.     float s=0.0;
     
320.     for(int i=0;i<end_support_i;i++)
     
321.         if(alph[i]>0)
     
322.             s+=alph[i]*target[i]*kernel_func(i,k);
     
323.         s-=b;
     
324.     return s;
     
325. }
     
326. //计算点积函数dot_product_func(int,int)
     
327. float dot_product_func(int i1,int i2)
     
328. {
     
329.     float dot=0;
     
330.     for(int i=0;i<d;i++)
     
331.         dot+=dense_points[i1][i]*dense_points[i2][i];   
     
332.     return dot;
     
333. }
     
334. //径向机核函数RBF:kernel_func(int,int)
     
335. float kernel_func(int i1,int i2)
     
336. {
     
337.     float s=dot_product_func(i1,i2);
     
338.     s*=-2;
     
339.     s+=precomputed_self_dot_product[i1]+precomputed_self_dot_product[i2];  
     
340.     return exp(-s/two_sigma_squared);
     
341. }
     
342. //初始化initialize()
     
343. void initialize()
     
344. {
     
345.     //初始化阀值b为0
     
346.     b=0.0;
     
347.     //初始化alph[]为0
     
348.     for(int i=0;i<end_support_i;i++)
     
349.     alph[i]=0.0;
     
350.    
     
351.     //设置样本值矩阵
     
352.     setX();
     
353.     //设置目标值向量
     
354.     setT();
     
355.     //设置预计算点积
     
356.     {for(int i=0;i<N;i++)
     
357.         precomputed_self_dot_product[i]=dot_product_func(i,i);//注意，这是对训练样本的设置，对于测试样本也应考虑?????????????????
     
358.     }
     
359. }
     
360. //计算误差率error_rate()
     
361. float error_rate()
     
362. {    ofstream to("d:\\smo_test.txt");
     
363.     int tp=0,tn=0,fp=0,fn=0;
     
364.     float ming=0,te=0,total_q=0,temp=0;
     
365.     for(int i=first_test_i;i<N;i++)
     
366.     {    temp=learned_func(i);
     
367.         if(temp>0&&target[i]>0)
     
368.             tp++;
     
369.         else if(temp>0&&target[i]<0)
     
370.             fp++;
     
371.         else if(temp<0&&target[i]>0)
     
372.             fn++;
     
373.         else if(temp<0&&target[i]<0)
     
374.             tn++;
     
375.         to<<i<<"  实际输出"<<temp<<endl;
     
376.     }
     
377.     total_q=(float)(tp+tn)/(float)(tp+tn+fp+fn);//总精度
     
378.     ming=(float)tp/(float)(tp+fn);
     
379.     te=(float)tp/(float)(tp+fp);
     
380.     to<<"---------------测试结果-----------------"<<endl;
     
381.     to<<"tp="<<tp<<"   tn="<<tn<<"  fp="<<fp<<"  fn="<<fn<<endl;
     
382.     to<<"ming="<<ming<<"  te="<<te<<"  total_q="<<total_q<<endl;
     
383.     return (1-total_q);
     
384. }
     
385. //计算样本X[]
     
386. void setX(){
     
387.     ifstream ff("D:\\17_smo.txt");////////////////////////
     
388.     if(!ff)
     
389.     {
     
390.         cerr<<"error!!"<<endl;
     
391.         exit(1);
     
392.     }
     
393.     int i=0,j=0;
     
394.     char ch;
     
395.     while(ff.get(ch))
     
396.     {
     
397.         if(isspace(ch)) continue;
     
398.         if(isdigit(ch))
     
399.         {
     
400.             dense_points[i][j]=(int)ch-48;
     
401.             j++;
     
402.             if(j==d)
     
403.             {
     
404.                 j=0;
     
405.                 i++;
     
406.             }
     
407.         }
     
408.     }
     
409. }
     
410. //set targetT[]
     
411. void setT()
     
412. {
     
413.     for(int i=0;i<6045;i++)
     
414.         target[i]=1;
     
415.     for(i=6045;i<12090;i++)
     
416.         target[i]=-1;
     
417.     for(i=12090;i<13097;i++)
     
418.         target[i]=1;
     
419.     for(i=13097;i<14104;i++)
     
420.         target[i]=-1;
     
421. }

复制代码

风花雪月

另外有人将其改为IDL语言，http://hi.baidu.com/hilbertspace ... a738c550da4b70.html

代码如下：
实现了硬最大线性间隔分界(最简单的情况)

1. ;;;;;SVM的SMO优化实现
   
2. ;;;;借用研学论坛上luke的代码
   
3. ;;Introduction: x: 二维样本数据.
   
4. ;;;;;;;;;;;;;;;;n: 数据规模
   
5. ;;;;;;;;;;;;;;;;y; 数据标志
   
6. ;;;;;;;;;;;;;;;;w; 边界线斜率
   
7. ;;;;;;;;;;;;;;;;b: 边界线截距
   
8. ;;;;;;;;;;;;;;;;所以 y=wx+b 是边界线
   
9. 
   
10. 
    
11. PRO svmlineartest
    
12. 
    
13. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
14. ;generate the training set;;;;;;;;
    
15. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
16. x=fltarr(2,40)
    
17. y=intarr(40)
    
18. x[0:1,0:19]=randomu(1.5,40)
    
19. x[0:1,20:39]=randomu(2.5,40)+0.5
    
20. y[0:19]=-1
    
21. y[20:39]=1
    
22. n=40
    
23. plot,x[0,n/2:n-1],x[1,n/2:n-1],psym=1,xrange=[0,2],yrange=[0,2]
    
24. oplot,x[0,0:n/2-1],x[1,0:n/2-1],psym=1,color=255
    
25. 
    
26. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
27. ;;;finished generate the training set;;;;
    
28. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
29. 
    
30. 
    
31. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
32. ;;;;;;minimise w^2, subject to: yi*(wxi+b) >= 1, then it can be converted to
    
33. ;;;;;;sum wi*alphi=0,   & ,   w- sum yi*alphi*xi =0 , then get
    
34. ;;;;;;maxmize L= sum alphi-0.5* sum (yi*yj*alphi*alphj*<xi.*xj> subject to: sum yi*alphi =0, ai>=0
    
35. ;;;;;;sumj alphi*yi*yj*xi,*xj=2
    
36. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
37. 
    
38. w=fltarr(2)
    
39. 
    
40. w=svmw(x,y,n)
    
41. END
    
42. 
    
43. Function svmw,x,y,n
    
44. 
    
45. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
46. ;;;;;;;intialization;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
47. ;;;;;;;使用指针方便全局操作;;;;;;;;;;;;;;;;;;;;;;
    
48. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
49. 
    
50. alph=fltarr(n)
    
51. C=100
    
52. tol=0.001
    
53. eps=1e-3
    
54. error_cache=fltarr(n)
    
55. alphptr=ptr_new(alph)
    
56. b=0.0
    
57. bptr=ptr_new(b)
    
58. errorptr=ptr_new(error_cache)
    
59. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
60. ;;;;;;;Finish intialization;;;;;;;;;;;;;;;;;
    
61. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
62. 
    
63. numChange=0
    
64. examineAll=1
    
65. 
    
66. 
    
67. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
68. ;;;;;;;以下两层循环，开始时检查所有样本，选择不符合KKT条件的两个乘子进行优化
    
69. ;;;;;;;选择成功，返回1，否则，返回
    
70. ;;;;;;;成功了，numChanged必然>0，从第二遍循环时，就不从整个样本中去寻找不符合KKT条件的两个乘子进行优化，
    
71. ;;;;;;;而是从边界的乘子中去寻找，因为非边界样本需要调整的可能性更大，边界样本往往不被调整而始终停留在边界上。
    
72. ;;;;;;;如果，没找到，再从整个样本中去找，直到整个样本中再也找不到需要改变的乘子为止，此时，算法结束。
    
73. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
74. 
    
75. while ((numChange gt 0) or examineAll) do begin
    
76. 
    
77.    numChanged =0
    
78.    if (examineAll) then begin
    
79.      for i=0,n-1 do begin
    
80.           numChanged+=examineExample(i,x,y,alphptr,n,bptr,errorptr)
    
81.      endfor
    
82.    endif
    
83. 
    
84.    if (examineAll eq 1) then examineAll=0   else if (numChanged eq 0) then examineAll=1
    
85. 
    
86. endwhile
    
87. 
    
88. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
89. ;;;;计算w,b并画出图;;;;;;;;;;;;;;
    
90. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
91. w=fltarr(2)
    
92. w[0]= total((*alphptr)*x[0,*]*y)
    
93. w[1]= total((*alphptr)*x[1,*]*y)
    
94. b=-(max(w[0]*x[0,0:n/2-1]+w[1]*x[1,0:n/2-1])+min(w[0]*x[0,n/2:n-1]+w[1]*x[1,n/2:n-1]))/2
    
95. print,w
    
96. print,b
    
97. 
    
98. x1new=findgen(100)*0.01*2
    
99. ynew=-(w[0]*x1new+b)/w[1]
    
100. oplot, x1new,ynew
     
101. End
     
102. 
     
103. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
104. ;;;;;;;;;;takestep用于优化两个乘子，成功，返回1，否则，返回0;;;;
     
105. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
106. 
     
107. Function takestep,i1,i2,x,y,alphptr,n,bptr,errorptr
     
108. 
     
109. 
     
110. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
111. ;;;;;;;;;;取出乘子和x,y,定义C值;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
112. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
113. alph=*alphptr
     
114. b=*bptr
     
115. C=100
     
116. tol=0.001
     
117. eps=1e-3
     
118. error_cache=*errorptr
     
119. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
120. ;;;;;;;;;;完成取出乘子和x,y;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
121. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
122. 
     
123. 
     
124.     if (i1 eq i2) then return, 0
     
125. 
     
126.     alph1 = alph[i1]
     
127.     y1 = y[i1]
     
128.     alph2 = alph[i2]
     
129.     y2 = y[i2]
     
130. s=y1*y2
     
131. if ((alph1 gt 0) and (alph1 lt C)) then E1=error_cache[i1] else E1=learned_func(i1,x,y,alph,n)-y1;//learned_func(int)为非线性的评价函数，即输出函数
     
132. if ((alph2 gt 0) and (alph2 lt C)) then E2=error_cache[i2] else E2=learned_func(i2,x,y,alph,n)-y2;
     
133. 
     
134. if y1 ne y2 then begin
     
135.     L=   0 > (alph[i2]-alph[i1])
     
136.     H=   C < (C+alph[i2]-alph[i1])
     
137. endif
     
138. if y1 eq y2 then begin
     
139.     L=   0 > (alph[i2]+alph[i1]-C)
     
140.     H=   C < (alph[i2]+alph[i1])
     
141. endif
     
142. 
     
143.     if (L eq H) then return, 0
     
144. 
     
145.      k11=kernel_func(i1,i1,x,y,n);//kernel_func(int,int)为核函数
     
146. k22=kernel_func(i2,i2,x,y,n);
     
147. k12=kernel_func(i1,i2,x,y,n)
     
148.     eta = 2*k12-k11-k22
     
149. 
     
150.     if (eta lt 0) then begin
     
151. 
     
152.         a2=alph2 - y2*(E1-E2)/eta
     
153.          if (a2 lt L) then a2=L   else if (a2 gt H) then a2=H
     
154. 
     
155.     endif else begin
     
156.       c1=eta/2;
     
157.    c2=y2*(E2-E1)-eta*alph2;
     
158.         Lobj=c1*L*L+c2*L
     
159.         Hobj=c1*H*H+c2*H
     
160. 
     
161.          if (Lobj gt (Hobj+eps))    then begin
     
162.            a2=L
     
163.          endif else if (Lobj lt Hobj-eps) then begin
     
164.              a2=H
     
165.          endif else begin
     
166.              a2=alph2
     
167.          endelse
     
168.    endelse
     
169. 
     
170.    if (abs(a2-alph2) lt eps*(a2+alph2+eps))   then return, 0
     
171.    a1=alph1+s*(alph2-a2)
     
172. 
     
173.    if (a1 lt 0) then begin
     
174.    a2+=s*a1
     
175.    a1=0
     
176.    endif else if (a1 gt C) then begin
     
177.     t=a1-C
     
178.    a2+=s*t
     
179.    a1=C
     
180. endif
     
181. 
     
182. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
183. ;;;;;;;;;;;;更新阈值b;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
184. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
185. 
     
186. 
     
187. if (a1 gt 0) and (a1 lt C) then begin
     
188.     bnew=b+E1+y1*(a1-alph1)*k11+y2*(a2-alph2)*k12;
     
189. endif else if (a2 gt 0) and (a2 lt C) then begin
     
190.     bnew=b+E2+y1*(a1-alph1)*k12+y2*(a2-alph2)*k22
     
191. endif else begin
     
192.      b1=b+E1+y1*(a1-alph1)*k11+y2*(a2-alph2)*k12
     
193.      b2=b+E2+y1*(a1-alph1)*k12+y2*(a2-alph2)*k22
     
194.      bnew=(b1+b2)/2
     
195. endelse
     
196. 
     
197. delta_b=bnew-b
     
198. b=bnew
     
199. *bptr=b
     
200. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
201. ;;;;;;;;;;;;完成更新阈值b;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
202. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
203. 
     
204. 
     
205. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
206. ;Update weight vector to reflect change in a1&a2,if linear SVM;;;;;;;
     
207. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
208. 
     
209. 
     
210. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
211. ;;;;;;;;;;;;更新error cache;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
212. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
213.       t1=y1*(a1-alph1);
     
214.    t2=y2*(a2-alph2);
     
215.    for i=0,n-1 do begin
     
216.    if (0 lt alph[i] and alph[i] lt C) then begin
     
217.     error_cache[i]+=t1*kernel_func(i1,i,x,y,n)+t2*(kernel_func(i2,i,x,y,n))-delta_b;
     
218.     error_cache[i1]=0.0;
     
219.     error_cache[i2]=0.0;
     
220.    endif
     
221.    endfor
     
222.    *errorptr=error_cache
     
223. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
224. ;;;;;;;;;;;;完成更新error cache;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
225. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
226. 
     
227. 
     
228. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
229. ;;;;;;;;;;;;更新lagrange乘子;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
230. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
231.         alph[i1]=a1;Store a1 in the alpha array
     
232.         alph[i2]=a2;Store a2 in the alpha array
     
233.         *alphptr=alph
     
234. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
235. ;;;;;;;;;;;;完成更新lagrange乘子;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
236. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
237. 
     
238.         return, 1
     
239. 
     
240. End
     
241. 
     
242. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
243. ;;;;;;;;;examineExample程序;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
244. ;;;;;;;;;假定第一个乘子ai(位置为i1)，examineExample(i1)首先检查;;;;;;;;;;;;;;;;;;;;;;;;;
     
245. ;;;;;;;;;如果它超出tolerance而违背KKT条件，那么它就成为第一个乘子;;;;;;;;;;;;;;;;;;;;;;;
     
246. ;;;;;;;;;然后，寻找第二个乘子(位置为i2),通过调用takeStep(i1,i2)来优化这两个乘子;;;;;;;;;
     
247. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
248. 
     
249. Function examineExample,i1,x,y,alphptr,n,bptr,errorptr
     
250. 
     
251. alph=(*alphptr)
     
252.    tol=0.001
     
253.      C=100
     
254.      eps=1e-3
     
255.      error_cache=*errorptr
     
256. 
     
257. 
     
258.    y1=y[i1]
     
259.    alph1=alph[i1]
     
260. 
     
261. if ((alph1 gt 0) and (alph1 lt C)) then E1=error_cache[i1] else E1=learned_func(i1,x,y,alph,n)-y1;//learned_func为计算输出函数
     
262. r1=y1*E1
     
263. 
     
264. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
265. ;;;;;;;;;;;;;;;;;违反KKT条件的判断;;;;;;;;;;;;;;;;;;
     
266. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
267. 
     
268. if (((r1 lt -tol) and (alph1 lt C)) or (( r1 gt tol) and (alph1 gt 0))) then begin
     
269.     ;if(size(where (alph ne 0)) gt 1 ) then begin
     
270. 
     
271. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
272. ;;;;;;;;;;;;;;;;;;;;;;;使用三种方法选择第二个乘子;;;;;;;;;;;;;;;;;;;;;;;;
     
273. ;;;;;;;;;;1：在non-bound乘子中寻找maximum fabs(E1-E2)的样本;;;;;;;;;;;;;;
     
274. ;;;;;;;;;;2：如果上面没取得进展,那么从随机位置查找non-boundary 样本;;;;;;
     
275. ;;;;;;;;;;3：如果上面也失败，则从随机位置查找整个样本,改为bound样本;;;;;;
     
276. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
277.       if (examineFirstChoice(i1,E1,x,y,alphptr,n,bptr,errorptr)) then return, 1
     
278.       if (examineNonBound(i1,x,y,alphptr,n,bptr,errorptr)) then   return, 1
     
279.       if (examineBound(i1,x,y,alphptr,n,bptr,errorptr)) then return, 1
     
280. 
     
281.     endif
     
282. ;endif
     
283.     return, 0
     
284. End
     
285. 
     
286. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
287. ;;;;;;1：在non-bound乘子中寻找maximum fabs(E1-E2)的样本;;;;;;;;;;;;;;;;;;
     
288. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
289. Function examineFirstChoice,i1,E1,x,y,alphptr,n,bptr,errorptr
     
290. C=100
     
291.      alph=(*alphptr)
     
292.      error_cache=*errorptr
     
293.     for k=0,n-1 do begin
     
294. 
     
295.          i2=-1
     
296.          tmax=0.0
     
297. 
     
298.      if(alph[k] gt 0) and (alph[k] lt C) then begin
     
299. 
     
300.       E2=error_cache[k];
     
301.       temp=abs(E1-E2);
     
302.       if (temp gt tmax) then begin
     
303. 
     
304.        tmax=temp;
     
305.        i2=k;
     
306.       endif
     
307.       endif
     
308.      endfor
     
309. 
     
310.       if (i2 ge 0) then begin
     
311. 
     
312.      if (takeStep(i1,i2,x,y,alphptr,n,bptr,errorptr)) then return, 1
     
313. 
     
314.     endif
     
315. 
     
316.     return, 0
     
317. 
     
318. end
     
319. 
     
320. 
     
321. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
322. ;;;;;;2：如果上面没取得进展,那么从随机位置查找non-boundary 样本;;;;;;;;;;
     
323. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
324. function examineNonBound,i1,x,y,alphptr,n,bptr,errorptr
     
325. C=100
     
326. alph=(*alphptr)
     
327. error_cache=*errorptr
     
328.      k0 = randomu(seed,1) mod n-1
     
329. 
     
330.      for k = 0, n-1 do begin
     
331. 
     
332.          i2 = (k + k0) mod n-1;//从随机位开始
     
333. 
     
334.          if (alph[i2] gt 0.0) and (alph[i2] lt C) then begin
     
335. 
     
336.             if (takeStep(i1, i2,x,y,alphptr,n,bptr,errorptr)) then begin
     
337. 
     
338.                return, 1;
     
339. 
     
340.          endif
     
341.          endif
     
342.       endfor
     
343.      return,0
     
344. 
     
345. end
     
346. 
     
347. 
     
348. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
349. ;;;;;;如果上面也失败，则从随机位置查找整个样本,(改为bound样本);;;;;;;;;;;
     
350. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
351. function examineBound,i1,x,y,alphptr,n,bptr,errorptr
     
352. C=100
     
353.    alph=(*alphptr)
     
354.    error_cache=*errorptr
     
355.         k0 = randomu(seed,1) mod (n-1);
     
356. 
     
357.      for k = 0, n-1 do begin
     
358. 
     
359.          i2 = (k + k0) mod (n-1);//从随机位开始
     
360. 
     
361.         if (alph[i2] eq 0.0 )or (alph[i2] eq C) then begin
     
362. 
     
363.             if (takeStep(i1, i2,x,y,alphptr,n,bptr,errorptr)) then   return, 1
     
364.             endif
     
365.       endfor
     
366.      return,0
     
367. end
     
368. 
     
369. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
370. ;;;;;;compute the learn function(y=wx+b) //test;;;;;;;;;;;;;;;;;;
     
371. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
372. 
     
373. Function learned_func ,i,x,y,alphlinshi,n
     
374. 
     
375. s=fltarr(n)
     
376. 
     
377. for j=0,n-1 do begin
     
378.    s[j]=alphlinshi[i]*alphlinshi[j]*y[i]*y[j]*(x[*,i] ## transpose(x[*,j]))
     
379. endfor
     
380. 
     
381. learned_func=total(s)
     
382. return, learned_func
     
383. 
     
384. End
     
385. 
     
386. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
387. ;;;;;;;;;;;;核函数kernel_func(int,int);;;;;;;;;;;;;;;;;;;;;;;;;;;
     
388. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
389. 
     
390. function kernel_func,i1,i2,x,y,n
     
391. func=(x[*,i1] ## transpose(x[*,i2]))
     
392. return,func
     
393. ;float s=dot_product_func(i1,i2);
     
394. ;s*=-2;
     
395. ;s+=precomputed_self_dot_product[i1]+precomputed_self_dot_product[i2];
     
396. ;return exp(-s/two_sigma_squared);
     
397. end

复制代码

jytan

太感谢了。

Marina_sun

  非常感谢 很有借鉴性

傻B小蜗牛

什么都不说了，感谢的一塌糊涂

novelbean

谢谢楼主无私分享，很好很强大！！