SVM算法之SMO算法

Halcom

SVM算法之SMO算法：百度网盘视频链接：
视频链接：http://pan.baidu.com/s/1kV5fPg7
录制的视频是算法底层原理讲解，底层代码实现，方便大家真正掌握算法实质，开发出更加出色的算法。录制视频的初衷是：避免读者朋友利用大把时间学习已有常见算法，本系列视频旨在让读者朋友快速深入了解这些常见算法原理，有更多的时间去研究更加高大上算法（价值）。

具体链接在halcom.cn论坛，联系人QQ：3283892722
该论坛是一个学习交流平台，我会逐一的和大家分享学习。
欢迎大家录制视频，并提交给我，我来设置视频，你可在论坛进行打赏分享。
视频专用播放器：http://halcom.cn/forum.php?mod=viewthread&tid=258&extra=page%3D1

运行环境：win7+32bit+matlab2014a
x：400x2
y：400x1

1. % SVM using Sequential Minimal Optimization (SMO)
   
2. clear all;clc;
   
3. data=csvread('LinearlySeprerableData.csv');
   
4. data(:,1:end-1)=zscore(data(:,1:end-1));
   
5. x = data(:,1:end-1);
   
6. y = data(:,end);
   
7. N = length(y);
   
8. C = 0.5;       % Concluded after Cross-Validation
   
9. tol = 10e-5;
   
10. alpha = zeros(N,1);
    
11. bias = 0;
    
12. 
    
13. %  SMO Algorithm
    
14. while (1)
    
15.     changed_alphas=0;
    
16.     N=size(y,1)
    
17.     for i=1:N
    
18.         Ei=sum(alpha.*y.*K(x,x(i,:),'l'))-y(i);
    
19.         if ((Ei*y(i)<-tol) && alpha(i)<C)||(Ei*y(i) > tol && (alpha(i) > 0))
    
20.             for j=[1:i-1,i+1:N]
    
21.                 Ej=sum(alpha.*y.*(x*x(j,:)'))-y(j);
    
22.                 alpha_iold=alpha(i);
    
23.                 alpha_jold=alpha(j);
    
24.                
    
25.                 if y(i)~=y(j)
    
26.                     L=max(0,alpha(j)-alpha(i));
    
27.                     H=min(C,C+alpha(j)-alpha(i));
    
28.                 else
    
29.                     L=max(0,alpha(i)+alpha(j)-C);
    
30.                     H=min(C,alpha(i)+alpha(j));
    
31.                 end
    
32.                
    
33.                 if (L==H)
    
34.                     continue
    
35.                 end
    
36.                
    
37.                 eta = 2*x(j,:)*x(i,:)'-x(i,:)*x(i,:)'-x(j,:)*x(j,:)';
    
38.                
    
39.                 if eta>=0
    
40.                     continue
    
41.                 end
    
42.                
    
43.                 alpha(j)=alpha(j)-( y(j)*(Ei-Ej) )/eta;
    
44.                 if alpha(j) > H
    
45.                     alpha(j) = H;
    
46.                 end
    
47.                 if alpha(j) < L
    
48.                     alpha(j) = L;
    
49.                 end
    
50.                
    
51.                 if norm(alpha(j)-alpha_jold,2) < tol
    
52.                     continue
    
53.                 end
    
54.                
    
55.                 alpha(i)=alpha(i)+y(i)*y(j)*(alpha_jold-alpha(j));
    
56.                 b1 = bias - Ei - y(i)*(alpha(i)-alpha_iold)*x(i,:)*x(i,:)'...
    
57.                     -y(j)*(alpha(j)-alpha_jold)*x(i,:)*x(j,:)';
    
58.                 b2 = bias - Ej - y(i)*(alpha(i)-alpha_iold)*x(i,:)*x(j,:)'...
    
59.                     -y(j)*(alpha(j)-alpha_jold)*x(j,:)*x(j,:)';
    
60.                
    
61.                
    
62.                 if 0<alpha(i)<C
    
63.                     bias=b1;
    
64.                 elseif 0<alpha(j)<C
    
65.                     bias=b2;
    
66.                 else
    
67.                     bias=(b1+b2)/2;
    
68.                 end
    
69.                 changed_alphas=changed_alphas+1;
    
70.             end
    
71.         end
    
72.     end
    
73.     if changed_alphas==0
    
74.         break
    
75.     end
    
76.     x=x((find(alpha~=0)),:);
    
77.     y=y((find(alpha~=0)),:);
    
78.     alpha=alpha((find(alpha~=0)),:);
    
79. end
    
80. % Weights
    
81. % W=sum(alpha.*y.*x)
    
82. W(1)=sum(alpha.*y.*x(:,1));
    
83. W(2)=sum(alpha.*y.*x(:,2));
    
84. % Bias
    
85. bias =mean( y - x*W')
    
86. % Support Vectors
    
87. disp('Number of support Vectors : ')
    
88. disp(N)
    
89. Xsupport=x;
    
90. Ysupport=y;
    
91. % Accuracy and F-measure
    
92. x=data(:,1:end-1);
    
93. y=data(:,end);
    
94. fx=sign(W*x'+bias)';
    
95. [~, Accuracy, F_measure ] = confusionMatrix( y,fx )
    
96. % Plotting the Decision Boundry
    
97. hold on
    
98. scatter(x(y==1,1),x(y==1,2),'b')
    
99. scatter(x(y==-1,1),x(y==-1,2),'r')
    
100. scatter(Xsupport(Ysupport==1,1),Xsupport(Ysupport==1,2),'.b')
     
101. scatter(Xsupport(Ysupport==-1,1),Xsupport(Ysupport==-1,2),'.r')
     
102. 
     
103. x1 = -2:0.01:2;
     
104. fn=((-bias-W(1)*x1)/W(2));
     
105. fn1=((-1-bias-W(1)*x1)/W(2));
     
106. fn2=((1-bias-W(1)*x1)/W(2));
     
107. plot(x1,fn,'r','Linewidth',2);
     
108. plot(x1,fn1,'Linewidth',1);
     
109. plot(x1,fn2,'Linewidth',1);
     
110. hold off

复制代码

参考链接：

【1】SVR回归（线性，多项式、RBF）预测模型
【2】Python底层SVR代码分享：http://pan.baidu.com/s/1kU773lt
【3】MATLAB底层SVR代码分享：http://pan.baidu.com/s/1pLG3dBd
【4】MATLAB底层SVC、SVR代码分享：http://pan.baidu.com/s/1dFEL6vB
【5】http://blog.csdn.net/yclzh0522/article/details/6900707
【6】http://blog.csdn.net/luoshixian099/article/details/51227754
【7】http://www.cnblogs.com/biyeymyhj ... /07/17/2591592.html
【8】支持向量机序列讲解http://blog.pluskid.org/?page_id=683
【9】python实现，理论与代码相结合的解说http://www.cnblogs.com/vivounicorn/archive/2011/06/01/2067496.html
【10】http://www.cnblogs.com/vivounicorn/archive/2011/01/13/1934296.html

Halcom

核函数如下：

1. function ker=K(Xtrain,x_i,k)
   
2. 
   
3. if k=='g'
   
4.     for i=1:size(Xtrain,1)
   
5.         ker(i,1)=exp(-norm(Xtrain(i,:)-x_i)); %gaussian Kernel
   
6.     end
   
7. elseif k=='l'
   
8.     for i=1:size(Xtrain,1)
   
9.         ker(i,1)=Xtrain(i,:)*x_i'; %linear Kernel
   
10.     end
    
11. elseif k=='p'
    
12.     for i=1:size(Xtrain,1)
    
13.         ker(i,1)=(Xtrain(i,:)*x_i').^3; %poly3 Kernel
    
14.     end
    
15. end
    
16. 
    
17. end

复制代码