求一完整的SVM分类器的程序，matlab编写的。感激不尽。。。。

这个是非线性svm的：

1.命令函数部分：

clear%清屏

clc

X =load('data.txt')

n = length(X)%总样本数量

y = X(:,4)%类别标志

X = X(:,1:3)

TOL = 0.0001%精度要求

C = 1%参数，对损失函数的权重

b = 0%初始设置截距b

Wold = 0%未更新a时的W(a)

Wnew = 0%更新a后的W(a)

for i = 1 : 50%设置类别标志为1或者-1

y(i) = -1

end

a = zeros(n,1)%参数a

for i = 1 : n%随机初始化a,a属于[0,C]

a(i) = 0.2

end

%为简化计算，减少重复计算进行的计算

K = ones(n,n)

for i = 1 :n%求出K矩阵，便于之后的计算

for j = 1 : n

K(i,j) = k(X(i,:),X(j,:))

end

end

sum = zeros(n,1)%中间变量，便于之后的计算，sum(k)=sigma a(i)*y(i)*K(k,i)

for k = 1 : n

for i = 1 : n

sum(k) = sum(k) + a(i) * y(i) * K(i,k)

end

end

while 1%迭代过程

%启发式选点

n1 = 1%初始化，n1,n2代表选择的2个点

n2 = 2

%n1按照第一个违反KKT条件的点选择

while n1 <= n

if y(n1) * (sum(n1) + b) == 1 &&a(n1) >= C &&a(n1) <= 0

break

end

if y(n1) * (sum(n1) + b) >1 &&a(n1) ~= 0

break

end

if y(n1) * (sum(n1) + b) <1 &&a(n1) ~=C

break

end

n1 = n1 + 1

end

%n2按照最大化|E1-E2|的原则选取

E1 = 0

E2 = 0

maxDiff = 0%假设的最大误差

E1 = sum(n1) + b - y(n1)%n1的误差

for i = 1 : n

tempSum = sum(i) + b - y(i)

if abs(E1 - tempSum)>maxDiff

maxDiff = abs(E1 - tempSum)

n2 = i

E2 = tempSum

end

end

%以下进行更新

a1old = a(n1)

a2old = a(n2)

KK = K(n1,n1) + K(n2,n2) - 2*K(n1,n2)

a2new = a2old + y(n2) *(E1 - E2) / KK%计算新的a2

%a2必须满足约束条件

S = y(n1) * y(n2)

if S == -1

U = max(0,a2old - a1old)

V = min(C,C - a1old + a2old)

else

U = max(0,a1old + a2old - C)

V = min(C,a1old + a2old)

end

if a2new >V

a2new = V

end

if a2new <U

a2new = U

end

a1new = a1old + S * (a2old - a2new)%计算新的a1

a(n1) = a1new%更新a

a(n2) = a2new

%更新部分值

sum = zeros(n,1)

for k = 1 : n

for i = 1 : n

sum(k) = sum(k) + a(i) * y(i) * K(i,k)

end

end

Wold = Wnew

Wnew = 0%更新a后的W(a)

tempSum = 0%临时变量

for i = 1 : n

for j = 1 : n

tempSum= tempSum + y(i )*y(j)*a(i)*a(j)*K(i,j)

end

Wnew= Wnew+ a(i)

end

Wnew= Wnew - 0.5 * tempSum

%以下更新b：通过找到某一个支持向量来计算

support = 1%支持向量坐标初始化

while abs(a(support))<1e-4 &&support <= n

support = support + 1

end

b = 1 / y(support) - sum(support)

%判断停止条件

if abs(Wnew/ Wold - 1 ) <= TOL

break

end

end

%输出结果：包括原分类，辨别函数计算结果，svm分类结果

for i = 1 : n

fprintf('第%d点:原标号 ',i)

if i <= 50

fprintf('-1')

else

fprintf(' 1')

end

fprintf('判别函数值%f 分类结果',sum(i) + b)

if abs(sum(i) + b - 1) <0.5

fprintf('1\n')

else if abs(sum(i) + b + 1) <0.5

fprintf('-1\n')

else

fprintf('归类错误\n')

end

end

end

2.名为f的功能函数部分:

function y = k(x1,x2)

y = exp(-0.5*norm(x1 - x2).^2)

end

3.数据：

0.8871 -0.34918.3376 0

1.25191.20836.5041 0

-1.19251.93381.8790 0

-0.12772.43712.6971 0

1.96973.09066.0391 0

0.76030.82411.5323 0

1.63823.55164.4694 0

1.3438 -0.45395.9366 0

-1.3361 -2.02011.6393 0

-0.38863.30418.0450 0

-0.67806.0196 -0.4084 0

0.3552 -0.10511.2458 0

1.65604.07860.8521 0

0.81173.54516.8925 0

1.4773 -1.93403.9256 0

-0.0732 -0.95260.4609 0

0.15214.37112.2600 0

1.48200.74930.3475 0

0.61404.52618.3776 0

0.57213.34603.7853 0

0.52694.14524.3900 0

1.7879 -0.53902.5516 0

0.98855.76250.1832 0

-0.33182.4373 -0.6884 0

1.35785.47093.4302 0

2.7210 -1.12684.7719 0

0.5039 -0.10252.3650 0

1.11071.68853.7650 0

0.78621.35877.3203 0

1.0444 -1.58413.6349 0

1.77951.72764.9847 0

0.67101.4724 -0.5504 0

0.23030.2720 -1.6028 0

1.7089 -1.73994.8882 0

1.00590.55575.1188 0

2.30500.85452.8294 0

1.95550.98980.3501 0

1.71411.54133.8739 0

2.27495.32804.9604 0

1.61710.52703.3826 0

3.6681 -1.84094.8934 0

1.19641.87811.4146 0

0.77882.10480.0380 0

0.79165.09063.8513 0

1.08071.88495.9766 0

0.63402.60303.6940 0

1.9069 -0.06097.4208 0

1.65994.94098.1108 0

1.37630.88993.9069 0

0.84851.46886.7393 0

3.67926.10924.9051 1

4.38127.21486.1211 1

4.39713.41397.7974 1

5.07167.7253 10.5373 1

5.30788.81386.1682 1

4.14485.51562.8731 1

5.36096.04584.0815 1

4.74526.63521.3689 1

6.02746.5397 -1.9120 1

5.31743.01346.7935 1

7.24593.69703.1246 1

6.10078.10875.5568 1

5.99246.92385.7938 1

6.02635.33337.5185 1

3.64708.09156.4713 1

3.65437.22647.5783 1

5.01146.53353.5229 1

4.43487.4379 -0.0292 1

3.60873.73513.0172 1

3.53745.53547.6578 1

6.00482.0691 10.4513 1

3.14234.00035.4994 1

3.40127.15368.3510 1

5.54715.1372 -1.5090 1

6.50895.49118.0468 1

5.45836.76745.9353 1

4.17272.97983.6027 1

5.16728.41364.8621 1

4.88083.55141.9953 1

5.49384.19983.2440 1

5.45425.88034.4269 1

4.87433.96418.1417 1

5.97626.77112.3816 1

6.69457.28581.8942 1

4.73015.76521.6608 1

4.70845.36233.2596 1

6.04083.31387.7876 1

4.60248.35170.2193 1

4.70546.6633 -0.3492 1

4.71395.63626.2330 1

4.0850 10.71183.3541 1

6.10886.16354.2292 1

4.98365.40426.7422 1

6.13876.19492.5614 1

6.07007.03733.3256 1

5.68815.13639.9254 1

7.20582.35704.7361 1

4.29727.32454.7928 1

4.77948.12353.1827 1

3.92826.4092 -0.6339 1

我先直观地阐述我对SVM的理解，这其中不会涉及数学公式，然后给出Python代码。

SVM是一种二分类模型，处理的数据可以分为三类：

线性可分，通过硬间隔最大化，学习线性分类器

近似线性可分，通过软间隔最大化，学习线性分类器

线性不可分，通过核函数以及软间隔最大化，学习非线性分类器

线性分类器，在平面上对应直线；非线性分类器，在平面上对应曲线。

硬间隔对应于线性可分数据集，可以将所有样本正确分类，也正因为如此，受噪声样本影响很大，不推荐。

软间隔对应于通常情况下的数据集（近似线性可分或线性不可分），允许一些超平面附近的样本被错误分类，从而提升了泛化性能。

如下图：

实线是由硬间隔最大化得到的，预测能力显然不及由软间隔最大化得到的虚线。

对于线性不可分的数据集，如下图：

我们直观上觉得这时线性分类器，也就是直线，不能很好的分开红点和蓝点。

但是可以用一个介于红点与蓝点之间的类似圆的曲线将二者分开，如下图：

我们假设这个黄色的曲线就是圆，不妨设其方程为x^2+y^2=1，那么核函数是干什么的呢？

我们将x^2映射为X，y^2映射为Y，那么超平面变成了X+Y=1。

那么原空间的线性不可分问题，就变成了新空间的（近似）线性可分问题。

此时就可以运用处理（近似）线性可分问题的方法去解决线性不可分数据集的分类问题。

---------------------------------------------------------------------------------------------------------------------------

以上我用最简单的语言粗略地解释了SVM，没有用到任何数学知识。但是没有数学，就体会不到SVM的精髓。因此接下来我会用尽量简洁的语言叙述SVM的数学思想，如果没有看过SVM推导过程的朋友完全可以跳过下面这段。

对于求解（近似）线性可分问题：

由最大间隔法，得到凸二次规划问题，这类问题是有最优解的（理论上可以直接调用二次规划计算包，得出最优解）

我们得到以上凸优化问题的对偶问题，一是因为对偶问题更容易求解，二是引入核函数，推广到非线性问题。

求解对偶问题得到原始问题的解，进而确定分离超平面和分类决策函数。由于对偶问题里目标函数和分类决策函数只涉及实例与实例之间的内积，即<xi,xj>。我们引入核函数的概念。

拓展到求解线性不可分问题：

如之前的例子，对于线性不可分的数据集的任意两个实例：xi,xj。当我们取某个特定映射f之后，f(xi)与f(xj)在高维空间中线性可分，运用上述的求解（近似）线性可分问题的方法，我们看到目标函数和分类决策函数只涉及内积<f(xi),f(xj)>。由于高维空间中的内积计算非常复杂，我们可以引入核函数K(xi,xj)=<f(xi),f(xj)>，因此内积问题变成了求函数值问题。最有趣的是，我们根本不需要知道映射f。精彩！

我不准备在这里放推导过程，因为已经有很多非常好的学习资料，如果有兴趣，可以看：CS229 Lecture notes

最后就是SMO算法求解SVM问题，有兴趣的话直接看作者论文：Sequential Minimal Optimization:A Fast Algorithm for Training Support Vector Machines

我直接给出代码：SMO+SVM

在线性可分数据集上运行结果：

图中标出了支持向量这个非常完美，支持向量都在超平面附近。

在线性不可分数据集上运行结果（200个样本）：

核函数用了高斯核，取了不同的sigma

sigma=1，有189个支持向量，相当于用整个数据集进行分类。

sigma=10，有20个支持向量，边界曲线能较好的拟合数据集特点。

我们可以看到，当支持向量太少，可能会得到很差的决策边界。如果支持向量太多，就相当于每次都利用整个数据集进行分类，类似KNN。

拉格朗日

function y=lagrange(x0,y0,x)

n=length(x0)m=length(x)

for i=1:m

z=x(i)

s=0.0

for k=1:n

p=1.0

for j=1:n

if j~=k

p=p*(z-x0(j))/(x0(k)-x0(j))

end

end

s=p*y0(k)+s

end

y(i)=s

end

SOR迭代法的Matlab程序

function [x]=SOR_iterative(A,b)

% 用SOR迭代求解线性方程组,矩阵A是方阵

x0=zeros(1,length(b))% 赋初值

tol=10^(-2)% 给定误差界

N=1000% 给定最大迭代次数

[n,n]=size(A)% 确定矩阵A的阶

w=1% 给定松弛因子

k=1

% 迭代过程

while k=N

x(1)=(b(1)-A(1,2:n)*x0(2:n)')/A(1,1)

for i=2:n

x(i)=(1-w)*x0(i)+w*(b(i)-A(i,1:i-1)*x(1:i-1)'-A(i,i+1:n)*x0(i+1:n)')/A(i,i)

end

if max(abs(x-x0))=tol

fid = fopen('SOR_iter_result.txt', 'wt')

fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n')

fprintf(fid,'迭代次数: %d次\n\n',k)

fprintf(fid,'x的值\n\n')

fprintf(fid, '%12.8f \n', x)

break

end

k=k+1

x0=x

end

if k==N+1

fid = fopen('SOR_iter_result.txt', 'wt')

fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n')

fprintf(fid,'迭代次数: %d次\n\n',k)

fprintf(fid,'超过最大迭代次数，求解失败！')

fclose(fid)

end

Matlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y'=f(x,y),使用差分概念。

(Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'）

Yn+1=Yn+h*f(Xn,Yn)

另外根据微分中值定理，存在0t1,使得

Yn+1=Yn+h*f(Xn+th,Y(Xn+th))

这里K＝f(Xn+th,Y(Xn+th))称为平均斜率，龙格库塔方法就是求得K的一种算法。

利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为O(h^5)的四阶龙格库塔公式：

K1＝f(Xn,Yn)

K2=f(Xn+h/2,Yn+(h/2)*K1)

K3=f(Xn+h/2,Yn+(h/2)*K2)

K4=f(Xn+h,Yn+h*K3)

Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6)

所以，为了更好更准确地把握时间关系，应自己在理解龙格库塔原理的基础上，编写定步长的龙格库塔函数，经过学习其原理，已经完成了一维的龙格库塔函数。

仔细思考之后，发现其实如果是需要解多个微分方程组，可以想象成多个微分方程并行进行求解，时间，步长都是共同的，首先把预定的初始值给每个微分方程的第一步，然后每走一步，对多个微分方程共同求解。想通之后发现，整个过程其实很直观，只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数：

function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）

n=floor((b-a)/h)%求步数

x(1)=a%时间起点

y(:,1)=y0%赋初值，可以是向量，但是要注意维数

for ii=1:n

x(ii+1)=x(ii)+h

k1=ufunc(x(ii),y(:,ii))

k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2)

k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2)

k4=ufunc(x(ii)+h,y(:,ii)+h*k3)

y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6

%按照龙格库塔方法进行数值求解

end

调用的子函数以及其调用语句：

function dy=test_fun(x,y)

dy = zeros(3,1)%初始化列向量

dy(1) = y(2) * y(3)

dy(2) = -y(1) + y(3)

dy(3) = -0.51 * y(1) * y(2)

对该微分方程组用ode45和自编的龙格库塔函数进行比较，调用如下：

[T,F] = ode45(@test_fun,[0 15],[1 1 3])

subplot(121)

plot(T,F)%Matlab自带的ode45函数效果

title('ode45函数效果')

[T1,F1]=runge_kutta1(@test_fun,[1 1 3],0.25,0,15)%测试时改变test_fun的函数维数，别忘记改变初始值的维数

subplot(122)

plot(T1,F1)%自编的龙格库塔函数效果

title('自编的 龙格库塔函数')

---