支撑向量机 SVM 学习笔记（Matlab代码）

 

 

% 支持向量机SVM分类算法

clear all
% ------------------------------------------------------------%
% 构造两类训练数据集     x2=aa*x1+bb+(-)b1

aa=3;
bb=6;
b1=0.2;       

x1(:,1) = -1:0.1:1;

n = length(x1(:,1));

x1(:,2) = aa.*x1(:,1) + bb + b1 + abs(randn(n,1));
y1 = ones(n,1);

x2(:,1) =  -1:0.1:1;
x2(:,2) = aa.*x2(:,1) + bb - b1 - abs(randn(n,1));

y2 = -ones(n,1);


figure;
plot(x1(:,1),x1(:,2),'bx',x2(:,1),x2(:,2),'k.');
hold on;

X = [x1;x2]; % 训练样本
Y = [y1;y2]; % 训练目标,n×1的矩阵,n为样本个数,值为+1或-1
% ------------------------------------------------------------%

tic
% 解二次优化方城
n = length(Y);
H = (Y*Y').*(X*X'); % liner kernel 
f = -ones(n,1);
A = [];
b = [];
Aeq = Y';
beq = 0;
lb = zeros(n,1);
ub = 100*ones(n,1);
a0 = zeros(n,1);
options = optimset;
options.LargeScale = 'off';
options.Display = 'off';
[a,fval,eXitflag,output,lambda] = quadprog(H,f,A,b,Aeq,beq,lb,ub,a0,options);
eXitflag

time=toc


   
%先复制上面的代码 运行后  再复制下面的代码 运行 呵呵 （对Matlab不熟，不知道为什么一次执行和分开执行看到的结果不同）

% 以下是分类平面：
Y2=a.*Y;
 W(1)=sum(Y2.*(X(:,1)));
 W(2)=sum(Y2.*(X(:,2)));

aLarge=find(a>0.1);
j=aLarge(1);

S(:,1)=Y.*a.*X(:,1);
S(:,2)=Y.*a.*X(:,2);

S2=S*(X(j,:)');

b=Y(j)-sum(S2);


 xx1=x1(:,1);
 xx2=-(W(1)*xx1+b)/W(2);
 plot(xx1,xx2);              

线性核  分界线

 

 

径向基核  分界线             可以看出非线性核可以拟合非线性曲线（对应 核空间中的线性曲线）

 

 

 

代码如下  还有部分小bug  例如有时候随机生成的数据 对应奇异矩阵  解二次规划 解不出来    此时可以重新运行

 

代码

 1function Result=CalcKernel(kernel,X,Y)
 2
 3
 4dimensionX=size(X);
 5m=dimensionX(1);
 6
 7
 8dimensionY=size(Y);
 9n=dimensionY(1);
10
11
12
13Result=zeros(m,n);
14
15if strcmp(kernel.type,'linear')==1
16    
17    for row=1:m
18        for col=1:n
19            Result(row,col)=X(row,:)*Y(col,:)';
20        end
21    end
22end
23
24
25if strcmp(kernel.type,'q-order polynomial kernel')==1
26    
27    q=kernel.value;
28    
29    for row=1:m
30        for col=1:n         
31            Result(row,col)=(X(row,:)*Y(col,:)' +1)^q;
32        end
33    end
34end
35
36
37if strcmp(kernel.type,'Radial basis kernel')==1
38    
39    q=2*kernel.value*kernel.value;
40    
41    for row=1:m
42        for col=1:n         
43            t=X(row,:)-Y(col,:);
44            Result(row,col)=exp(-t*t'/q);
45        end
46    end
47end
48
49

 

 

 

 1function result=CalcValueBySVM(svm, x)
 2
 3
 4result=svm.b;
 5
 6n=size(svm.y);
 7
 8for j=1:n
 9    result=result+ svm.a(j)*svm.y(j)*CalcKernel(svm.ker, x, svm.x(j,:));
10end
11
12
13

 

 

 1
 2function SVM_Test
 3
 4
 5clear all
 6% ------------------------------------------------------------%
 7% 构造两类训练数据集    以 y=ax^2+ b*x+c 为分界线   shift   
 8%  if a==0  then it is a line
 9
10a=0;
11b=3;
12c=0.5;
13
14shift=0.5;
15
16b1=0.5;       
17
18xx=-1:0.1:1;
19
20n = length(xx);
21
22x1=zeros(n,2);
23x2=zeros(n,2);
24
25x1(:,1) = xx';
26x2(:,1)=xx';
27
28
29
30x1(:,2)=a*x1(:,1).*x1(:,1)+b*x1(:,1)+c + shift+   abs(randn(n,1)) ;
31x2(:,2)=a*x2(:,1).*x2(:,1)+b*x2(:,1)+c - shift-  abs(randn(n,1)) ;
32
33y1 = ones(n,1);
34y2 = -ones(n,1);
35
36figure;
37plot(x1(:,1),x1(:,2),'bx',x2(:,1),x2(:,2),'k.');
38hold on;
39
40X = [x1;x2]; % 训练样本
41Y = [y1;y2]; % 训练目标,n×1的矩阵,n为样本个数,值为+1或-1
42
43
44
45% ------------------------------------------------------------%
46
47%ker=struct('type','linear');
48ker=struct('type','Radial basis kernel','value',1);
49
50svm=SVM_DivideTwoClass(ker,X,Y)
51
52
53
54
55
56%plot the  decision function curve
57
58for i=1:n
59
60    tx1=x1(i,1);    
61    tx2=a*tx1*tx1+b*tx1+c;    
62      
63    min=1000;
64    
65    for step=-2:0.01:2  
66           
67        tx=[tx1,tx2+step];      
68        
69        f=CalcValueBySVM(svm,tx);
70        
71        if(abs(f)<min)
72            min=abs(f);            
73            decisionCurve(i,:)=tx;
74        end
75    end
76end
77
78
79plot(decisionCurve(:,1),decisionCurve(:,2));    
80
81
82
83 
84 
85 for i=1:2*n
86     f=CalcValueBySVM(svm,X(i,:));
87    
88     YY(i)=sign(f);        
89 end
90 YY=YY'
91
92

 

 补充以前漏的，呵呵，刚刚调试出来的

%---------------- author: feathersky----------------

function svm=SVM_DivideTwoClass(ker, X,Y)

% 解二次优化方城
n = length(Y);
 %H = (Y*Y').*(X*X'); % liner kernel
H = (Y*Y').*CalcKernel(ker,X,X); %  kernel
f = -ones(n,1);
A = [];
b = [];
Aeq = Y';
beq = 0;
lb = zeros(n,1);
ub = 1000*ones(n,1);
a0 = zeros(n,1);
options = optimset;
options.LargeScale = 'off';
options.Display = 'off';
[a,fval,eXitflag,output,lambda] = quadprog(H,f,A,b,Aeq,beq,lb,ub,a0,options);
eXitflag

svm.ker=ker;
svm.x=X;
svm.y=Y;
svm.a=a;



% aLarge=find(a>0.01);    %找到第一个不等于0的a ,   
% j=aLarge(1);
% svm.b=Y(j)-CalcKernel(ker,X(j,:),w);


% ------------------------------------------------------------%%
% 求 b
epsilon = 1e-8; % 如果小于此值则认为是0
i_sv = find(a>epsilon); % 支持向量下标
tmp = (Y.*a)'*CalcKernel(ker,X,X(i_sv,:)); % 行向量
b = 1./Y(i_sv)-tmp';
b = mean(b);
svm.b=b
fprintf('Construct function Y = sign(tmp+b):')

% ------------------------------------------------------------%