版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/eeeee123456/article/details/80108096
引言:
支持向量机(support vector machines,SVM)是一种二值分类模型。分类学习最基本的想法就是基于训练集在样本空间中找到一个分隔超平面(separating hyperplane),将不同类别的样本分开。但能将训练样本分开的分隔面可能有很多,我们应该努力去找哪一个呢?直观上看,应该去找位于两类训练样本‘正中间’的分隔面,即图1中的红线。
图1 线性可分支持向量机
如何寻找这条红线?我们希望找到离分隔面最近的点,确保它们离分隔面的距离尽可能远。这里点到分隔面的距离称为间隔(margin)。我们要做的就是为红线设一线性方程,用点到线的距离公式计算间隔,求出能使间隔最大化的线性方程即可求出分隔超平面。
由于寻找分隔超平面就是寻找能使间隔最大化的函数的值,介绍理论知识部分会用到求函数极值的许多数学知识,均在‘数学预备知识’中介绍,此部分只给出结论,不予证明,有此数学基础可跳过。另外,现实中很多数据并不像图1那样清晰可分,更多时候是如图2或图3所示,它们分别对应支持向量机中由简至繁的模型:线性可分支持向量机、线性支持向量机、非线性支持向量机。
图2 线性支持向量机
图3 非线性支持向量机
一、数学预备知识
无约束条件下求函数极值
1.函数极值必要条件
设函数
z=f(x,y)
在点
(x0,y0)
具有偏导数,且在点
(x0,y0)
处有极值,则有
fx(x0,y0)=0
fy(x0,y0)=0
以上关于二元函数的概念,可推广到
n
元函数。
例1 某厂要用钢板做成一个体积为
2m3
的长方体有盖水箱,问当长、宽、高各取怎样的尺寸时,才能使用料最省?
解:
设水箱的长、宽各为
x
、
y
,则其高应为
2xy
,此水箱所用材料的面积
A=2(xy+x⋅2xy+y⋅2xy)=2(xy+2y+2x) x>0,y>0
下面求使上述函数取得最小值的点(x,y),求其对x、y的一阶偏导数得
Ax=2(y−2x2)=0
Ay=2(x−2y2)=0
解得
x=2‾√3
,
y=2‾√3
。
根据题意可知,水箱所用材料面积的最小值一定存在,并在开区域
D={(x,y)|x>0,y>0}
内取得,又知函数在
D
内只有唯一的驻点
(2‾√3,2‾√3)
,因此可判定,当长、宽、高各取
2‾√3
、
2‾√3
、
2xy=2‾√3
时,使用料最省。
求等式约束条件下函数极值,可用拉格朗日乘子法求解
2.拉格朗日乘子法
要找函数
z=f(x,y)
在约束条件
s.t. φ(x,y)=0 (1)
下的可能极值点,可以先作拉格朗日函数
L(x,y,λ)=f(x,y)+λφ(x,y)
其中
λ
为参数,求其对
x
与
y
的一阶偏导数,并使之为零,然后与方程(1)联立起来:
fx(x,y)+λφx(x,y)=0
fy(x,y)+λφy(x,y)=0
φ(x,y)=0
由这方程组解出
x
、
y
及
λ
,这样得到的
(x,y)
就是函数
f(x,y)
在附加条件
φ(x,y)=0
下的可能极值点。
λ
称为拉格朗日乘子。
这方法还可以推广到自变量多于两个而条件多于一个的情形。
例如要求函数
u=f(x,y,z,t)
在约束条件
s.t. φ1(x,y,z,t)=0,φ2(x,y,z,t)=0 (2)
下的可能极值点,可以先作拉格朗日函数
L(x,y,z,t,λ1,λ2)=f(x,y,z,t)+λ1φ1(x,y,z,t)+λ2φ2(x,y,z,t)
其中
λ1
、
λ2
均为参数,求其对
x
与
y
的一阶偏导数,并使之为零,然后与(2)中2个方程联立起来求解,这样得出的
(x,y,z,t)
就是函数
f(x,y,z,t)
在附加条件(2)下的可能极值点。
至于如何确定所求得的点是否为极值点,在实际问题中往往可根据问题本身的性质来判定。
例2 求表面积为
a2
体积为最大的长方体的体积。
解:
设长方体的三棱长为
x,y,z
,则问题就是在条件
φ(x,y,z)=2xy+2xz+2yz−a2=0下
,求函数
V=xyz (x>0,y>0,z>0)
的最大值。作拉格朗日函数
L(x,y,z,λ)=xyz+λ(2xy+2xz+2yz−a2)
求其对
x,y,z
的偏导数,并使之为零,得到
yz+2λ(y+z)=0
xz+2λ(x+z)=0
xy+2λ(x+y)=0
再与约束条件
φ(x,y,z)
联立求解,得
x=y=z=6√6a
,这是唯一可能的极值点,因为由问题本身可知最大值一定存在,所以最大值就在这个可能的极值点处取得。
固表面积为
a2
体积为最大的长方体的体积为
xyz=6√36a3
。
求不等式约束条件下函数极值,可用K-T条件求解
3.库恩-塔克条件(K-T条件)
定义: 线性规划的一般形式
设
X∗
是非线性规划(4)式的极小点,且
X∗
点的所有起作用约束的梯度线性无关,则存在向量
T=(λ1,λ2,⋯,λn)T
,使下述条件成立:
条件(5)式常简称为K-T条件,满足这个条件的点
X∗
称为库恩-塔克点(或K-T点)。
为了得出非线性规划(3)式的库恩-塔克条件,我们用
代替约束条件
hi(X)=0
,这样即可得出(5)求解。
上式中,
λ1,λ2,...,λn
称为广义拉格朗日(Lagrange)乘子。
例3 用K-T条件解下述非线性规划问题
min f(x)=(x−3)2
s.t. 0≤x≤5
解:
先将该非线性规划问题写成以下形式
min f(x)=(x−3)2
s.t. g1(x)=x≥0
g2(x)=5−x≥0
写出其目标函数和约束函数的梯度:
∇f(x)=2(x−3)
∇g1(x)=1,∇g2(x)=−1
对第一个和第二个约束条件分别引入广义拉格朗日乘子
λ1
和
λ2
,则可写出该问题的K-T条件如下:
2(x−3)−λ1+λ2=0
λ1x=0
λ2(5−x)=0
λ1,λ2≥0
解得
λ1=0
,
λ2=0
,
x=3
。
固函数在约束条件下的最小值为
f(3)=0
。
例4 写出下述非线性规划问题的K-T条件
min f(x1,x2)=(x1−3)2+(x2−2)2
s.t. 4−x1−x2≥0
x1,x2≥0
解:
先将该非线性规划问题写成以下形式
min f(x1,x2)=(x1−3)2+(x2−2)2
s.t. g1(x1,x2)=4−x1−x2≥0
g2(x1,x2)=x1≥0
g3(x1,x2)=x2≥0
写出其目标函数和约束函数的梯度:
∇f(x1,x2)=(2x1−6,2x2−4)T
∇g1(x1,x2)=(−1,−1)T,∇g2(x1,x2)=(1,0)T,∇g3(x1,x2)=(0,1)T
对三个约束条件分别引入广义拉格朗日乘子
λ1
、
λ2
和
λ3
,则可写出该问题的K-T条件如下:
(2x1−6,2x2−4)T−λ1(−1,−1)T−λ2(1,0)T−λ3(0,1)T=0
分解为:
2x1−6+λ1−λ2=0
2x2−4+λ1−λ3=0
λ1(4−x1−x2)=0
λ2x1=0
λ3x2=0
λ1,λ2,λ3≥0
求不等式约束条件下函数极值,若原始问题过于复杂,可转换成对偶问题
4.拉格朗日对偶问题
1)原始问题
假设
f(x)
、
ci(x)
、
hj(x)
是定义在
Rn
上的连续可微函数,考虑约束条件
称此约束最优化问题为原始最优化问题或原始问题。
作广义拉格朗日函数
L(x,α,β)=f(x)+∑i=1kαici(x)+∑j=1lβjhj(x)
其中
x=(x1,x2,...,xn)T∈Rn
,
αi
和
βj
是拉格朗日乘子且
αi≥0
。考虑
x
的函数:
θp(x)=maxα,β,αi≥0L(x,α,β)
这里下标P表示原始问题,则
θp(x)=maxα,β,αi≥0f(x)+∑i=1kαici(x)+∑j=1lβjhj(x)
因为
αi≥0
,
ci(x)≤0
,
hj(x)=0
,解得
θP(x)=f(x)
。
考虑极小化问题
minxθp(x)=minxmaxα,β,αi≥0L(x,α,β)
它是与原始最优化问题(6)等价的,即它们有相同的解。问题
minxmaxα,β,αi≥0L(x,α,β)
称为广义拉格朗日的极小极大问题。
为了方便,定义原始问题的最优值
p∗=minxθp(x)
称为原始问题的值。
2)对偶问题
θD(α,β)=minxL(x,α,β)
再考虑极大化
θD(α,β)=minxL(x,α,β)
,即
maxα,β,αi≥0θD(α,β)=maxα,β,αi≥0minxL(x,α,β)
问题
maxα,β,αi≥0minxL(x,α,β)
称为广义拉格朗日的极大极小问题。
将拉格朗日函数的极大极小问题表示为约束最优化问题:
称为原始问题的对偶问题。
为了方便,定义对偶问题的最优值
d∗=maxα,β,αi≥0θD(α,β)
称为对偶问题的值。
3)原始问题与对偶问题的关系
对原始问题(6)和对偶问题(7),假设函数
f(x)
和
ci(x)
是凸函数,
hj(x)
是仿射函数,并且不等式约束
ci(x)
是严格可行的,即存在
x
,对所有
i
有
ci(x)<0
,则存在
x∗
、
α∗
、
β∗
,使得
x∗
是原始问题的解,
α∗
、
β∗
是对偶问题的解,并且
p∗=d∗=L(x∗,α∗,β∗)
x∗
和
α∗
、
β∗
分别是原始问题和对偶问题的解的充分必要条件是
x∗
、
α∗
、
β∗
满足K-T条件。
例5 求下列非线性规划问题的对偶问题
min x12+x22
s.t. 4−x1−x2≤0
解:
该问题的拉格朗日函数为
L(x,α)=x12+x22+α(4−x1−x2), α≥0
设
f(x)=x12+x22
,其中
x=(x1,x2)T
。则原始问题为:
minx f(x)=minx θP(x)=minxmaxα,α≥0L(x,α)
对偶问题为:
maxα,α≥0θD(α)=maxα,α≥0minxL(x,α)
解
minxL(x,α)=minxx12+x22+α(4−x1−x2)
,求其对
x1
、
x2
的一阶偏导数得
Lx1=2x1−α=0
Lx2=2x2−α=0
解得
x1=α2
,
x2=α2
,代入得
minxL(x,α)=−α22+4α
求得对偶问题为
maxα,α≥0θD(α)=maxα,α≥0−α22+4α
二、线性可分支持向量机
1.基本型
给定训练集
D={(x1,y1),(x2,y2),⋯,(xN,yN)}
,
yi∈{−1,1}
,分类学习最基本的想法就是基于训练集
D
在样本空间中找到一个分隔超平面,将不同类别的样本分开。
在样本空间中,分隔超平面可通过如下线性方程来描述:
wTx+b=0 (8)
其中
w=(w1,w2,⋯,wn)T
为法向量,决定了分隔面的方向,
b
为位移项,决定了分隔面与原点之间的距离。样本空间中任意点x到分隔面的距离可写为
r=|wTx+b|||w||
证明:
考虑二维空间,则(8)可表示为
w1x1+w2x2+b=0
,点
(x1,x2)
到线的距离为
r=|w1x1+w2x2+b|w12+w22‾‾‾‾‾‾‾‾‾‾√
考虑三维空间,则(8)可表示为
w1x1+w2x2+w3x3+b=0
,点
(x1,x2,x3)
到面的距离为
r=|w1x1+w2x2+w3x3+b|w12+w22+w32‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√
以此类推,可得
n
维空间,点
(x1,x2,…,xn
)到分隔面的距离为
r=|w1x1+w2x2+⋯+wnxn+b|w12+w22+⋯+wn2‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√=|wTx+b|||w||
假设分隔面能将训练样本正确分类,即对于
(xi,yi)∈D
,若
yi=+1
,则有
wTx+b>0
;若
yi=−1
,则有
wTx+b<0
。令
如图4所示,距离分隔面最近的这几个训练样本点(位于黑线上的点)使式(9)的等号成立,它们被称为“支持向量”,两个异类支持向量到超平面的距离(两条黑线的间隔)之和为
γ=2||w||
它被称为间隔。
证明:
假设
ξ>0
且
ξ2=1
。
若
yi=+1
,则有
wTx+b>0
等价于
wTx+b≥ξ
,等号两边同时乘以
ξ
得
ξwTx+ξb≥1
;
若
yi=−1
,则有
wTx+b<0
等价于
wTx+b≤−ξ
,等号两边同时乘以
ξ
得
ξwTx+ξb≤−1
;
即若分隔面能将训练样本正确分类,则总存在缩放变换
ξw→w′
和
ξb→b′
使式(9)成立。
图4
欲找到具有“最大间隔”的分隔面,也就是要找到能满足式(9)中约束的参数
w
和
b
,使
得
γ
最大,即
显然为了最大化间隔,仅需要最大化
||w||−1
,这等价于最小化
||w||2
,于是式(10)可重写为
这就是支持向量机的基本型。
2.对偶问题
基本型的对偶问题为:
证明:
基本型的拉格朗日函数为
L(w,b,α)=||w||22+∑i=1nαi(1−yi(wTxi+b)), i=1,2,...,n
其中
α=(α1,α2,...,αn)
,设
f(w,b)=||w||22
,则原始问题为
minw,b f(w,b)=minw,b θP(w,b)=minw,bmaxα,α≥0L(w,b,α)
对偶问题为:
maxα,α≥0θD(α)=maxα,α≥0minw,bL(w,b,α)
要求对偶问题,先求
minw,bL(w,b,α)
,即求
L(w,b,α)
对
w
和
b
的一阶偏导数并使之为零
Lw=w−∑i=1nαiyixi=0
Lb=−∑i=1nαiyi=0
解得
w=∑i=1nαiyixi (13)
0=∑i=1nαiyi (14)
将(13)代入拉格朗日函数
L(w,b,α)
,并利用(14)得
已知
minw,bL(w,b,α)
则可得对偶问题(12)
3.利用对偶问题求解
w
和
b
设
α∗=(α1∗,α2∗,...,αn∗)T
是对偶问题(12)的解,则存在下标
j
,使得
αj∗>0
,并可按下式求得原始问题(11)的解
w∗
、
b∗
得
w∗=∑i=1nαiyixi (15)
b∗=yi−∑i=1nαiyixiTxj (16)
证明:
w∗
、
b∗
和
α∗
分别是原始问题和对偶问题的解的充分必要条件是
w∗
、
b∗
、
α∗
满足K-T条件。
先将(11)式写成以下形式:
minw,b f(w,b)=||w||22
s.t. gi(w,b)=yi(wTxi+b)−1≥0, i=1,2,...,n
写出其目标函数和约束函数的梯度:
∇f(w,b)=(||w||,0)T,∇gi(w,b)=(xiyi,yi)T
对
n
个约束条件分别引入广义拉格朗日乘子
α1,α2,...,αn
,则可写出该问题的K-T条件如下:
(||w||,0)T−∑i=1nαi(xiyi,yi)T=0
将上式分解为:
||w||=∑i=1nαixiyi=0
−∑i=1nαiyi=0
αi(yi(wTxi+b)−1)=0
αi≥0, i=1,2,...,n
由此解得式(15)
其中至少有一个
αj∗>0
(用反证法,假设
α∗=0
,由式(15)可知
w∗=0
,而
w∗=0
不是原始问题(11)的解,产生矛盾),对此
j
有
yj(w∗xj+b∗)−1=0(17)
注意到
yj2=1
,将(17)两边同时乘以
yj
,将式(15)代入即得式(16)
综上所述,对于给定的线性可分训练数据集,可以先求对偶问题(12)的解
α∗
,再利用式(15) (16)求得原始问题的解
w∗
和
b∗
,从而得到分隔超平面的线性方程。这种算法称为线性可分支持向量机的对偶算法,是线性可分支持向量机学习的基本算法。
三、线性支持向量机
1.基本型
在前面的讨论中,我们一直假定训练样本是线性可分的,即存在一个分隔面能将不同类的样本完全划分开。然而在现实任务中,往往存在一些异常点。线性不可分意味着某些样本点
(xi,yi)
不能满足函数间隔大于等于1的约束条件,如图5绿色框中的两个样本点。为了解决这个问题,可以对每个样本点
(xi,yi)
引进一个松弛变量
ξi≥0
,使函数间隔加上松弛变量大于等于1,这样,约束条件变为
yi(wTxi+b)+ξi≥1
目标函数由原来的
minw,b||w||22
变成
minw,b,ξ ||w||22+C∑i=1nξi (18)
这里
C>0
称为惩罚参数,一般由应用问题决定,
C
值大时对误分类的惩罚增大,
C
值小时对误分类的惩罚减小。最小化目标函数(18)包含两层含义:使
||w||22
尽量小即间隔尽量大,同时使误分类点的个数尽量小。至此,线性支持向量机基本型如下:
图5
2.对偶问题
基本型的对偶问题为:
证明:
基本型的拉格朗日函数为
L(w,b,ξ,α,μ)=||w||22+C∑i=1nξi+∑i=1nαi(1−yi(wTxi+b)−ξ)−∑i=1nμiξi
其中
αi≥0
,
μi≥0
,设
f(w,b,ξ)=||w||22+C∑ni=1ξi
,则原始问题为
minw,b,ξ f(w,b,ξ)=minw,b,ξ θP(w,b,ξ)=minw,b,ξ maxα,μL(w,b,ξ,α,μ)
对偶问题为:
maxα,μθD(α,μ)=maxα,μminw,b,ξL(w,b,ξ,α,μ)
解
minw,b,ξL(w,b,ξ,α,μ)
,求
L(w,b,ξ,α,μ)
对
w
、
b
和
ξ
的一阶偏导数得
Lw=w−∑i=1nαixiyi=0
Lb=−∑i=1nαiyi=0
Lξ=C−αi−μi=0
解得
w=∑i=1nαixiyi (21)
∑i=1nαiyi=0 (22)
C−αi−μi=0 (23)
将(21)代入
minw,b,ξL(w,b,ξ,α,μ)
,并利用(22)(23)求得对偶问题
maxα,μ∑i=1nαi−12∑i=1n∑j=1nαiαjyiyjxiTxj
s.t. ∑i=1nαiyi=0
C−αi−μi=0 (24)
αi≥0, i=1,2,…,n
μi≥0, i=1,2,…,n
利用等式(24)得
μi=C−αi
,即可消去
μi
,并将(24)的约束写成
0≤αi≤C
从而得到最终的对偶问题(20)
3.利用对偶问题求解
w
和
b
设
α∗=(α1∗,α2∗,...,αn∗)T
是对偶问题(20)的解,若存在
α∗的
一个分量
αj∗
,满足
0<αj∗<C
,则原始问题(19)的解
w∗
、
b∗
可按下式求得
w∗=∑i=1nαixiyi (25)
b∗=yi=∑i=1nαiyixiTxj (26)
证明:
w∗
、
b∗
和
α∗
分别是原始问题和对偶问题的解的充分必要条件是
w∗
、
b∗
、
α∗
满足K-T条件。
先将式(19)写成以下形式:
minw,b,ξ f(w,b,ξ)=||w||22+C∑i=1nξi
s.t. gi(w,b,ξi)=yi(wTxi+b)+ξi−1≥0
hi(w,b,ξi)=ξi≥0
写出其目标函数和约束函数的梯度:
∇f(w,b,ξ)=(||w||,0,C)T
∇gi(w,b,ξi)=(xiyi,yi,1)T,∇hi(w,b,ξi)=(0,0,1)T
对
2n
个约束条件分别引入广义拉格朗日乘子
α1,α2,...,αn
和
μ1,μ2,…,μn
,则可写出该问题的K-T条件如下:
(||w||,0,C)T−∑i=1nαi(xiyi,yi,1)T=∑i=1nμi(0,0,1)T=0
将上式分解为:
||w||−∑i=1nαixiyi=0
−∑i=1nαiyi=0
C−αi−μi=0
αi(yi(wTxi+b)+ξi−1)=0, i=1,2,…,n
μiξi=0, i=1,2,…,n
αi≥0, i=1,2,...,n
μi≥0, i=1,2,...,n
由此解得式(25)
选择一个
αj∗
满足
0<αj∗<C
,对此有
yj(w∗xj+b∗)−1=0 (27)
注意到
yj2=1
,将(27)两边同时乘以
yj
,将式(25)代入即得式(26)
综上所述,对于给定的训练数据集,可以先求对偶问题(20)的解
α∗
,再利用式(25)(26)求得原始问题的解
w∗
和
b∗
,从而得到分隔超平面的线性方程。
四、非线性支持向量机
1、非线性分类问题
一般来说,对给定的一个训练数据集,如果能用一个超曲面将数据集正确分开,则称这个问题为非线性分类问题。
在图3中,数据点处于一个圆中,无法用直线(线性模型)将数据正确分开,但可以用一条椭圆曲线(非线性模型)将它们正确分开。
非线性问题往往不好求解,所以希望能用线性分类问题的方法解决,所采用的方法是进行一个非线性变换,将非线性问题变换为线性问题,通过解变换后的线性问题的方法求解原来的非线性问题。
2、核函数的定义
设X是输入空间,又称H为特征空间(希尔伯特空间),如果存在一个从X到H的映射
φ(x):X→H
使得对所有z,x∈X,函数K(x,z)满足条件
K(x,z)=φ(x)·φ(z)
则称K(x,z)为核函数,φ(x)为映射函数,式中φ(x)·φ(z)为φ(x)和φ(z)的内积。
核技巧的想法是,在学习与预测中只定义核函数K(x,z),而不显式地定义映射函数φ。通常,直接计算K(x,z)比较容易,而通过φ(x)和φ(z)计算K(x,z)并不容易。
例6
如下图,左边椭圆方程为
w1(x(1))2+w2(x(2))2+b=0
定义映射
y=φ(x)=((x(1))2,(x(2))2)T
则可将椭圆变换为下图中右边的直线
w1y(1)+w2y(2)+b=0
根据定义可得对所有z,x∈X,核函数为
K(x,z)=φ(x)·φ(z)=(x(1))2(z(1))2+(x(2))2(z(2))2
3、核技巧在支持向量机中的应用
我们注意到在线性支持向量机的对偶问题中,无论是目标函数还是决策函数都只涉及输入实例与实例之间的内积。在对偶问题(20)的目标函数中的内积
xTixj
可以用核函数
K(xi,xj)=xTixj
来代替。
4、常用核函数
五、序列最小最优化(Sequential Minimal Optimization,SMO)算法
如何求解式(20)的解α*呢?人们通过利用问题本身的特性,提出了很多高效算法,SMO就是其中一个著名的代表,这种算法1998年由Platt提出。
将对偶问题(20)的目标函数求极大转换为求极小,即可得到SMO算法要解的如下问题:
其中
K(xi,xj)=xiTxj
,在这个问题中,变量是拉格朗日乘子,一个变量
αi
对应一个样本点
(xi,yi)
,变量的总数等于训练样本容量
N
。
1.SMO算法思路
SMO是一种启发式算法,该算法并不直接求解对偶问题,而是从K-T条件入手,如果所有变量的解
α∗
都满足此最优化问题的K-T条件,那么这个最优化问题的解就得到了,因为K-T条件是该最优化问题的充分必要条件。具体做法如下:
首先,初始化一个
α
,此
α
满足对偶问题的约束条件(29)(30),若
α
也刚好满足K-T条件,得解。然而,很多时候我们需要不断优化
α
,直到
α
即满足(29)(30)也满足K-T条件,以此来求得最优解。
如何优化
α
呢?
先选择两个待优化的变量
αiold
和
αjold
(其中一个是违反K-T条件最严重的,另一个由约束条件自动确定),其他变量固定,从而构造一个新的二次规划问题,问题就变成简单的求二次函数极值。我们不太可能求解一次就让优化后的
αnew
满足K-T条件,但是关于优化后的
αinew
和
αjnew
应该更接近原始问题的解,因为这会使得原始二次规划问题的目标函数值变得更小。
如此,SMO算法将原问题不断分解为子问题并对子问题求解,进而达到求解原问题的目的。
2.构造两个变量的二次规划问题
不失一般性,假设选择的两个变量是
α1
和
α2
,其他变量
αi
(
i=3,4,…,n
)是固定的。于是(28)-(30)的子问题可以写成:
其中
Kij=K(xi,xj)
,
i,j=1,2,…,n
,
ξ
是常数。
3.两个变量二次规划问题的求解方法
假设问题的初始解为
αiold
、
αjold
,最优解为
αinew
、
αjnew
,并且假设未经剪辑(忽略约束条件(33))时
α2
的最优解为
α2new,unc
。
为了叙述简单,记
g(x)=∑j=1nαjyjK(x,xj)+b
令
Ei=g(xi)−yi=(∑j=1nαiyiK(xi,xj)+b)−yi, i=1,2 (34)
当
i=1,2
时,
Ei
为函数
g(x)
对输入
xi
的预测值与真实输出
yi
之差。
则
α2new,unc=α2old+y2(E1−E2)K11+K22−2K12 (35)
求
α2new
得
其中当
y1
、
y2
异号
L=max(0,α2old−α1old),H=min(C,C+α2old−α1old)
当
y1
、
y2
同号
L=max(0,α2old+α1old−C),H=min(C,α2old+α1old)
由
α2new
求得
α1new
是
α1new=α1old+y1y2(α2old−α2new) (36)
证明:
①求未经剪辑(忽略约束条件(33))时
α2
的最优解
α2new,unc
引进记号
vi=∑j=3nαjyjK(xi,xj)=g(xi)−∑j=12αjyjK(xi,xj)+b−b
目标函数可写成
W(α1,α2)=12K11α12+12K22α22+y1y2K12α1α2−(α1+α+2)+y1v1α1+y2v2α2
由
α1y1=ξ−α2y2
,等式两边同时乘以
y1
,由
yi2=1
可将
α1
表示为
α1=(ξ−α2y2)y1
代入式(31),得到只有α2的函数:
对α2求导数得
令其为0,得到
将
ξ=α1oldy1+α2oldy2
代入,得到
解得
α2new,unc=α2old+y2(E1−E2K11+K22−2K12
②求经剪辑(考虑约束条件(33))
α2new
的取值范围
假设取值范围为
L≤α2new≤H
考虑约束条件
0≤αi≤C,i=1,2
时
α2
的最优解,首先分析约束条件
α1y1+α2y2=−∑i=3nyiαi=ξ
将
α1
当做自变量,
α2
为因变量,则上式变成
α2=−y1y2α1+y2ξ (37)
当
y1
、
y2
异号,令
k=y2ξ
,则式(37)为直线方程
α2=α1+k
,图像如图6。
A线为
k>0
时,
α2new
的最大值为
C
,最小值为
k=α2old−α1old
;
B线为
k<0
时,
α2new
的最大值为
C+k=C+α2old−α1old
,最小值为0;
综上,
α2new
需满足(
α2
的取值即要在直线上,又要在方框内):
max(0,α2old−α1old)≤α2new≤min(C,C+α2old−α1old)
得
L=max(0,α2old−α1old)
,
H=min(C,C+α2old−α1old)
。
图6
当
y1
、
y2
同号,令
k=y2ξ
,则式(37)为直线方程
α2=−α1+k
,图像如图7。
A线为
k>C
时,
α2new
的最大值为
C
,最小值为
−Ck=α2old+α1old−C
;
B线为
k<C
时,
α2new
的最大值为
k=α2old+α1old
,最小值为0;
综上,
α2new
需满足(α2的取值即要在直线上,又要在方框内):
max(0,α2old+α1old−C)≤α2new≤min(C,α2old+α1old)
得
L=max(0,α2old+α1old−C)
,
H=min(C,α2old+α1old)
。
图7
③
α2new
最终解
因
α2new
需满足
α2new=α2new,unc
,又需满足
L≤α2new≤H
,综上可得:
④解
α1new
由式(32)得
α1oldy1+α2oldy2=ξ=α1newy1+α2newy2
上式解得
α1new=α1old+y1y2(α2old−α2new)
4.两个变量(
α1old
和
α2old
)的选择方法
①第1个变量的选择
SMO称选择第1个变量的过程为外层循环,外层循环在训练样本中选取违反K-T条件最严重的样本点,并将其对应的变量作为第1个变量,具体地,检验训练样本点(xi,yi)是否满足K-T条件,即
证明:
在“线性支持向量机”利用对偶问题求w和b那一部分,我们求得如下K-T条件为:
将
g(xi)
代入式(42),则式(42)可写成
αi(yig(xi)+ξi−1)=0
。
②第2个变量的选择
SMO称选择第2个变量的过程为内层循环。假设在外层循环中已经找到第1个变量
α1old
,现在要在内层循环中找第2个变量
α2old
,第2个变量的选择标准是希望能使
α2
有足够大的变化。
由式(35)可知,
α2new
是依赖于
|E1−E2|
的,为了加快计算速度,一种简单的做法是选择
α2
,使其对应的
|E1−E2|
最大。因为
α1old
已定,
E1
也确定了。如果
E1
为正,选择最小的
Ei
作为
E2
;如果
E1
为负,选择最大的
Ei
作为
E2
。为了节省计算时间,将所有
Ei
值保存在一个列表中。
在特殊情况下,如果内层循环通过以上方法选择的
α2new
不能使目标函数有足够的下降,那么采用以下启发式规则继续选择
α2new
。遍历在间隔边界上的支持向量点,依次将其对应的变量作为
α2new
试用,直到目标函数有足够的下降。若找不到合适的
α2new
,那么遍历训练数据集;若仍找不到合适的
α2new
,则放弃第1个
α1old
,再通过外层循环寻求另外的
α1old
。
5.计算b和差值
Ei
在每次完成两个变量的优化后,都要重新计算
b
值,因为
b
值关系到下一次优化时
Ei
的计算。
bnew
值的选择是使
α1new
、
α2new
满足K-T条件。
当
0<α1new<C
b1new=−E1−y1K11(α1new−α1old)−y2K21(α2new−α2old)+bold
当
0<α2new<C
b2new=−E2−y1K12(α1new−α1old)−y2K22(α2new−α2old)+bold
如果
α1new
、
α2new
同时满足条件
0<αinew<C
,
i=1,2
,那么
bnew=b1new=b2new
如果
α1new
、
α2new
其中一个满足
0<αinew<C
,另一个不满足
0<αinew<C
,这里假设
α2new
满足
0<α2new<C
,另一个即
α1new
是0或者C,那么
bnew=b2new
如果
α1new
、
α2new
是0或者C,那么
b1new
和
b2new
以及它们之间的数都符合K-T条件,这时选择它们的中点作为
bnew
,那么
bnew=b1new+b2new2
在每次完成两个变量的优化之后,还必须更新对应的
Ei
值,并将它们保存在列表中。
Ei
值的更新要用到
bnew
值,以及所有支持向量对应的
αj
,得
Einew=∑SαjyjK(xi,xj)+bnew−yi
其中,
S
是所有支持向量
xj
的集合。
6.SMO算法
输入:训练数据集
D={(x1,y1),(x2,y2),⋯,(xN,yN)}
,
yi∈−1,1
,
i=1,2,…,N
,精度
ε
,惩罚参数
C
输出:近似解
α∗
①取初始值
α(0)=0
,令
k=0
②选取优化变量
α1(k)
、
α2(k)
,解析求解两个变量的最优化问题,求得最优解
α1(k+1)
、
α2(k+1)
,更新
α
为
α(k+1)
。
③若在精度ε范围内满足停机条件
则转④;否则令
k=k+1
,转②
④取
α∗=α(k+1)
五、SMO代码实现
以下代码来自Peter Harrington《Machine Learing in Action》。
通过一个外循环来选择第一个alpha值,并且其选择过程会在两种方式之间进行交替:一种方式是在所有数据集上进行单遍扫描,另一种方式是在alpha值大于0和小于C的数据集上进行单遍扫描。第二个alpha值的选择与上述第二个变量基本一致。停机条件为无alpha可优化或到达最大循环次数。
主函数smoP需要五个参数:数据集、类别标签、常数C、精度ε和取消前最大的循环次数。
代码如下(保存为smo.py):
# -- coding: utf-8 --
from numpy import *
def loadDataSet(fileName):
dataMat = []; labelMat = []
fr = open(fileName)
for line in fr.readlines():
lineArr = line.strip().split('\t')
dataMat.append([float(lineArr[0]), float(lineArr[1])])
labelMat.append(float(lineArr[2]))
return dataMat,labelMat
def selectJrand(i,m):
j=i
while (j==i):
j = int(random.uniform(0,m))
return j
def clipAlpha(aj,H,L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
class optStruct:
def __init__(self,dataMatIn, classLabels, C, toler):
self.X = dataMatIn
self.labelMat = classLabels
self.C = C
self.tol = toler
self.m = shape(dataMatIn)[0]
self.alphas = mat(zeros((self.m,1))) #初始化alphas为0
self.b = 0 #初始化b为0
self.eCache = mat(zeros((self.m,2))) #初始化Ek为0
def calcEk(oS, k):
fXk = float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T)) + oS.b
Ek = fXk - float(oS.labelMat[k]) #对应式(34)Ei=g(xi)+b-yi
return Ek
def selectJ(i, oS, Ei): #返回使|E1-E2|值最大的E2及E2的位置j
maxK = -1
maxDeltaE = 0
Ej = 0
oS.eCache[i] = [1,Ei] #更新位置i的的E1,用1标示该值为更新过的Ei
validEcacheList = nonzero(oS.eCache[:,0].A)[0] #返回eCache数组第一列不为0的值,即更新过的Ei的索引
if (len(validEcacheList)) > 1: #如果eCache的Ek个数大于1,返回使|E1-E2|值最大的E2及E2的位置j
for k in validEcacheList:
if k == i: continue
Ek = calcEk(oS, k)
deltaE = abs(Ei - Ek)
if (deltaE > maxDeltaE):
maxK = k; maxDeltaE = deltaE; Ej = Ek
return maxK, Ej
else: #如果eCache的Ek都是没有更新过的,则随机选择一个与E1不同位置的E2并返回E2的位置j
j = selectJrand(i, oS.m)
Ej = calcEk(oS, j)
return j, Ej
def updateEk(oS, k):
Ek = calcEk(oS, k)
oS.eCache[k] = [1,Ek]
def innerL(i, oS):
Ei = calcEk(oS, i) #计算E1
if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)):
# (yiEi<0 => yig(xi)<1 and alpha<C) or (yiEi>0 => yig(xi)>1 and alpha>0) 选取违背K-T条件的alpha1
j,Ej = selectJ(i, oS, Ei) #选择smo的第二个变量
alphaIold = oS.alphas[i].copy() #alpha1old
alphaJold = oS.alphas[j].copy() #alpha2old
if (oS.labelMat[i] != oS.labelMat[j]): #y1、y2异号
L = max(0, oS.alphas[j] - oS.alphas[i]) #L=max(0,alpha2old-alpha1old)
H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i]) #H=min(0,C+alpha2old-alpha1old)
else: #y1、y2同号
L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C) #L=max(0,alpha2old+alpha1old-C)
H = min(oS.C, oS.alphas[j] + oS.alphas[i]) #H=min(0,alpha2old+alpha1old)
if L==H: print "L==H"; return 0 #异常情况,返回
eta = 2.0 * oS.X[i,:]*oS.X[j,:].T - oS.X[i,:]*oS.X[i,:].T - oS.X[j,:]*oS.X[j,:].T #2K12-K11-K22
if eta >= 0: print "eta>=0"; return 0 #异常情况,返回
oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta #式(35)计算alpha2new,unc
oS.alphas[j] = clipAlpha(oS.alphas[j],H,L) #返回alpha2new
updateEk(oS, j) #更新E2的值
if (abs(oS.alphas[j] - alphaJold) < 0.00001): print "j not moving enough"; return 0 #alpha2基本不变,返回
oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j]) #式(36)求alpha1new
updateEk(oS, i) #更新E1的值
#b1new
b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T
#b2new
b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[j,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].T
#根据b值的计算结论更新b值
if (0 < oS.alphas[i]) and (oS.C > oS.alphas[i]): oS.b = b1
elif (0 < oS.alphas[j]) and (oS.C > oS.alphas[j]): oS.b = b2
else: oS.b = (b1 + b2)/2.0
return 1
else: return 0
def smoP(dataMatIn, classLabels, C, toler, maxIter):
#样本X,样本Y,惩罚参数C,精度e,最大循环次数maxIter
oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler)
iter = 0
entireSet = True
alphaPairsChanged = 0
while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
#循环次数大于maxIter或alphaPairsChanged连续2次等于0,则退出循环
alphaPairsChanged = 0
if entireSet:
for i in range(oS.m): #循环所有alphas,作为smo的第一个变量
alphaPairsChanged += innerL(i,oS) #inner函数:如果有任意一对alphas值发生改变,返回1
print "fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
iter += 1
else:
nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
for i in nonBoundIs: #循环大于0小于C的alphas,作为smo的第一个变量
alphaPairsChanged += innerL(i,oS)
print "non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)
iter += 1
if entireSet: entireSet = False
elif (alphaPairsChanged == 0): entireSet = True
print "iteration number: %d" % iter
return oS.b,oS.alphas
def calcWs(alphas,dataArr,classLabels): #利用式(25)alphas[i]*y[i]*xi求w的值
X = mat(dataArr)
labelMat = mat(classLabels).transpose()
m,n = shape(X)
w = zeros((n,1))
for i in range(m):
w += multiply(alphas[i]*labelMat[i],X[i,:].T)
return w
数据集如下(保存为textSet.txt):
3.542485 1.977398 -1
3.018896 2.556416 -1
7.551510 -1.580030 1
2.114999 -0.004466 -1
8.127113 1.274372 1
7.108772 -0.986906 1
8.610639 2.046708 1
2.326297 0.265213 -1
3.634009 1.730537 -1
0.341367 -0.894998 -1
3.125951 0.293251 -1
2.123252 -0.783563 -1
0.887835 -2.797792 -1
7.139979 -2.329896 1
1.696414 -1.212496 -1
8.117032 0.623493 1
8.497162 -0.266649 1
4.658191 3.507396 -1
8.197181 1.545132 1
1.208047 0.213100 -1
1.928486 -0.321870 -1
2.175808 -0.014527 -1
7.886608 0.461755 1
3.223038 -0.552392 -1
3.628502 2.190585 -1
7.407860 -0.121961 1
7.286357 0.251077 1
2.301095 -0.533988 -1
-0.232542 -0.547690 -1
3.457096 -0.082216 -1
3.023938 -0.057392 -1
8.015003 0.885325 1
8.991748 0.923154 1
7.916831 -1.781735 1
7.616862 -0.217958 1
2.450939 0.744967 -1
7.270337 -2.507834 1
1.749721 -0.961902 -1
1.803111 -0.176349 -1
8.804461 3.044301 1
1.231257 -0.568573 -1
2.074915 1.410550 -1
-0.743036 -1.736103 -1
3.536555 3.964960 -1
8.410143 0.025606 1
7.382988 -0.478764 1
6.960661 -0.245353 1
8.234460 0.701868 1
8.168618 -0.903835 1
1.534187 -0.622492 -1
9.229518 2.066088 1
7.886242 0.191813 1
2.893743 -1.643468 -1
1.870457 -1.040420 -1
5.286862 -2.358286 1
6.080573 0.418886 1
2.544314 1.714165 -1
6.016004 -3.753712 1
0.926310 -0.564359 -1
0.870296 -0.109952 -1
2.369345 1.375695 -1
1.363782 -0.254082 -1
7.279460 -0.189572 1
1.896005 0.515080 -1
8.102154 -0.603875 1
2.529893 0.662657 -1
1.963874 -0.365233 -1
8.132048 0.785914 1
8.245938 0.372366 1
6.543888 0.433164 1
-0.236713 -5.766721 -1
8.112593 0.295839 1
9.803425 1.495167 1
1.497407 -0.552916 -1
1.336267 -1.632889 -1
9.205805 -0.586480 1
1.966279 -1.840439 -1
8.398012 1.584918 1
7.239953 -1.764292 1
7.556201 0.241185 1
9.015509 0.345019 1
8.266085 -0.230977 1
8.545620 2.788799 1
9.295969 1.346332 1
2.404234 0.570278 -1
2.037772 0.021919 -1
1.727631 -0.453143 -1
1.979395 -0.050773 -1
8.092288 -1.372433 1
1.667645 0.239204 -1
9.854303 1.365116 1
7.921057 -1.327587 1
8.500757 1.492372 1
1.339746 -0.291183 -1
3.107511 0.758367 -1
2.609525 0.902979 -1
3.263585 1.367898 -1
2.912122 -0.202359 -1
1.731786 0.589096 -1
2.387003 1.573131 -1
运行命令如下:
通过loadDataSet函数转换数据集,运行smoP函数获取b及alphas值。
最后通过calcWs函数计算w的值。
可利用数据验证,若该值大于0属于1类,若该值小于0则属于-1类。
以上全部内容参考书籍如下:
李航《统计学习方法》
周志华《机器学习》
Peter Harrington《Machine Learing in Action》
《高等数学第六版下册(同济大学)》
《运筹学 第三版》清华大学出版社