1 前备知识
在这里简略讲一下使用方法,具体原理和推导公式不展开讲了。
1.1 拉格朗日乘子法
拉格朗日乘子法就是求函数
f(x1,x2,...)在约束条件
g(x1,x2,...)=0下的极值的方法。其主要思想是将约束条件函数与原函数联立,从而求出使原函数取得极值的各个变量的解。
首先看下面的例题:
min f=2x12+3x22+x32s.t. 2x1+x2−1=0 2x2+x3−2=0
第一步将每个约束条件都分配一个乘子
αi,在将目标函数和所有的约束函数相加,得到函数:
L=f+i=1∑mgiαi
其中每个约束条件
gi的右边都是0,所以
∑i=1mgi=0.
L=(2x12+3x22+x32)+α1(2x1+x2−1)+α2(2x2+x3−2)
第二步对
xi求偏导:
⎩⎨⎧∂x1∂L=4x1+2α1∂x2∂L=6x2+α1+2α2∂x3∂L=2x3+α2
令偏导数等于0,用
αi表示
x:
⎩⎨⎧x1=−2α1x2=−6α1+2α2x3=−2α2
将所得
x代入约束条件
g中,求得
α:
{α1=−2/5α2=−72/45
得到
α的值,代入上式得到
x的最优解。
1.2 KKT条件
我们可以发现,1.1讲的拉格朗日乘子法中,它的约束条件都是等式,那么对于约束条件是不等式的应该怎么办呢?
对于一个新的极值问题:
min f=x12−2x1+x22+5s.t. x1+10x2>10 10x1−x2<10
为了统一,首先将约束条件都转化为小于号:
min f=x12−2x1+x22+5s.t. 10−x1−10x2<0 10x1−x2−10<0
依旧是分配乘子并求和:
L=f+i=1∑mgiαi+i=1∑mhiβi
其中
gi是不等式约束条件,
hi是等式约束条件。(此例中没有等式)
L=(x12−2x1+x22+5)+α1(10−x1−10x2)+α2(10x1−x2−10)
KKT条件就是最优值,KKT条件为:
-
L对每个
x求偏导等于
0;
-
h(x)=0;
-
gi(x)<=0
-
αi>=0
-
∑αigi(x)=0
可以发现,将3、4、5合并就是:
αigi(x)=0
对于上例题,接下来的操作就是:
一、
L对每个
x求偏导等于
0求出
x的表达式。
二、将
x的表达式代入
αigi(x)=0,求出
α。
三、将
α代回,求出
x。
2 SVM
2.1 简介
支持向量机(support vector machines, SVM)是一种二分类问题模型。
它的目标是找到一个尽可能正确分类,且“确信度”尽可能高的超平面。
其中“确信度”指的是:正确分类的样本点,距离超平面越远,该样本点的确信度就越高。(我对这个样本点分类正确的信任程度)
换而言之,就是该超平面的鲁棒性要好,泛化能力要强。
对于线性可分支持向量机,分类超平面为:
w∗⋅x+b∗=0
相应的分类决策函数
f(x)=sign(w∗⋅x+b∗)
称为线性可分支持向量机。
2.2 函数间隔与几何间隔
函数间隔和几何间隔是用来描述计算“确信度”的。
2.1.1 函数间隔
在超平面
w⋅x+b=0确定的情况下,
∣w⋅xi+b∣的值可以作为衡量样本点
xi确信度的一个指标,
∣w⋅xi+b∣就是该样本点的函数间隔。
定义(函数间隔):对于给定的训练数据集
T和超平面
(w,b),定义超平面
(w,b)关于样本点
(xi,yi)的函数间隔为:
γi^=∣w⋅xi+b∣
超平面
(w,b)关于样本点
(xi,yi)的函数间隔为所有样本点函数间隔的最小值:
γ^=min γi^, i=1,2,...,N
函数间隔虽然可以表示预测的确信度,但是当
w和
b成比例增加时,超平面没有改变,但函数间隔却成倍增加。
例如超平面
λw∗⋅x+λb∗=0与
w∗⋅x+b∗=0等价,但是
γi^=λ∣w⋅xi+b∣。
为了解决这个问题,引入了几何间隔。
2.2.2 几何间隔
在二维坐标系中,点是样本,线是分离超平面,那么点到线的距离就是几何间隔。
点到线的距离公式:
d=A2+B2
∣Ax0+By0+C∣
扩展到多维坐标系:
d=∣∣w∣∣2∣w⋅x+b∣
其中
∣∣w∣∣2=∑i=1mw
为L2范数。
我们记
γi=∣∣w∣∣2∣w⋅xi+b∣
同样:
γ=min γi, i=1,2,...,N
2.2.3 间隔最大化
回想一下,我们SVM的目标是什么来着?
是寻找一个“确信度”尽可能高的超平面,也就是一个几何间隔尽可能大的超平面。
那么我们可以得到一个约束最优化问题:
使几何间隔最大化的
w和
b,并且满足约束条件所有样本的几何间隔大于
γ.
w,bmax γ s.t. ∣∣w∣∣2yi(w⋅xi+b)≥γ, i=1,..,N
又因为函数间隔和几何间隔的关系:
γ=∣∣w∣∣2γ^
上述问题可以化为:
w,bmax ∣∣w∣∣2γ^ s.t. yi(w⋅xi+b)≥γ^, i=1,..,N
函数间隔
γ^的取值并不影响最优化问题的解,则可以令
γ^=1。
这个其实也好理解,我们看上式,
w和
b就是两个参数,令
γ^=1,就相当于将
w′=γ^w和
b′=γ^b,对于该约束最优化问题的解没有任何影响。就好比解方程时等式两边同时除以一个数。
另外,最大化
∣∣w∣∣21和最小化
21∣∣w∣∣2是等价的。
所以上面的问题就变成了:
w,bmin 21∣∣w∣∣2 s.t. yi(w⋅xi+b)−1≥0, i=1,..,N
这就是SVM的基本型(对于线性可分问题),后面主要就是这个约束问题的求解。
2.3 对偶问题
我们可以发现上面的约束最优化问题本身就是一个凸二次规划问题,所以我们可以使用更高效的方法去求解,也就是使用拉格朗日乘子法得到其“对偶问题”。
首先稍微做一下调整:
w,bmin 21∣∣w∣∣2 s.t. 1−yi(w⋅xi+b)≤0, i=1,..,N
定义拉格朗日乘子
αi,根据
L=f+i=1∑mgiαi+i=1∑mhiβi
得到:
L=21∣∣w∣∣2+i=1∑mαi[1−yi(w⋅xi+b)]
令
L对
w和
b求偏导等于零,得:
⎩⎨⎧w=∑i=1mαiyixi0=∑i=1mαiyi
虽然把
b给消去了,但是并没有影响,后续过程中原式的
b也会消去。
将
w=∑i=1mαiyixi代入
L中,即可用
αi将
w和
b替换掉,得到对偶问题:
αmax i=1∑mαi−21i=1∑mj=1∑mαiαjyiyjxiTxjs.t. i=1∑mαiyi=0 , αi≥0, i=1,...,m
另外,需要注意,这里的约束问题需要用到KKT条件:
-
L对
w和
b求偏导等于
0;
-
yi(w⋅xi+b)−1≥0
-
αi>=0
-
αi[yi(w⋅xi+b)−1]=0
-
h(x)=0; (这里没有等式约束,可以忽略)
整理一下,可简化为两种情况:
- 当
αi=0时,
yi(w⋅xi+b)−1≥0
- 当
0<αi<ξi时,
yi(w⋅xi+b)−1=0
(
ξi为松弛变量,下面会讲。)
到这里,只要解出
αi,就可以得到
w和
b,从而得到分离超平面。
那么如何去求解
αi呢?一种比较高效的方法就是SMO算法。
在此之前需要先讲一下松弛变量和核函数。
2.4 松弛变量和核函数
2.4.1 线性SVM与松弛变量
在此之前所讲的都是数据集时线性可分的,但是还有的数据集是线性不可分的,这两种情况合起来就是线性SVM。
如下图,线性不可分的数据集:
线性不可分意味着某些样本点
(xi,yi)不能满足函数间隔大于等于1的约束条件。为了解决这个问题,给每一个样本点
(xi,yi)引进一个松弛变量
ξi≥0,使函数间隔加上松弛变量大于等于1。这样,约束条件变为:
yi(w⋅xi+b)−1+ξi≥0
同时,对每个松弛变量
ξi,需要支付一个代价
ξi,则目标函数就由原来的
21∣∣w∣∣2变为:
21∣∣w∣∣2+Ci=1∑Nξi
其中
C>0称为惩罚参数,用来控制松弛变量的代价高低。
所以对于线性不可分的SVM的基本模型就是:
w,b,ξmin 21∣∣w∣∣2+Ci=1∑Nξi s.t. yi(w⋅xi+b)ξi≥1, ξi≥0, i=1,..,N
2.4.2 非线性SVM与核函数
2.4.2.1 核技巧和核函数
简单说就是,对于非线性问题,可以将样本通过一个
φ(x)从原始空间映射到一个更高维的特征空间,使得样本在这个特征空间内线性可分。
如下图,对于异或问题,就是一个非线性问题,原始问题是在一个二维空间中,当我们将样本特征空间做一个映射,提升到三维空间中,就能容易找到一个分离超平面。
非线性SVM的基本模型为:
w,b,ξmin 21∣∣w∣∣2 s.t. yi(w⋅φ(xi)+b)≥1, i=1,..,N
其对偶问题为:
αmax i=1∑mαi−21i=1∑mj=1∑mαiαjyiyjφ(xi)Tφ(xj)s.t. i=1∑mαiyi=0 , αi≥0, i=1,...,m
特征空间的维数可能非常高。如果支持向量机的求解只用到内积运算,而在低维输入空间又存在某个函数
K(xi,xj) ,它恰好等于在高维空间中这个内积,即
K(xi,xj)=<φ(xi)⋅φ(xj)>=φ(xi)Tφ(xj)那么支持向量机就不用计算复杂的非线性变换,而由这个函数
K(xi,xj) 直接得到非线性变换的内积,使大大简化了计算。这样的函数
K(xi,xj) 称为核函数。
则该对偶问题可以改写为:
αmax i=1∑mαi−21i=1∑mj=1∑mαiαjyiyjK(xi,xj)s.t. i=1∑mαiyi=0 , αi≥0, i=1,...,m
2.4.2.2 常用的核函数
3 SMO算法
序列最小最优化(sequential minimal optimization,SMO)算法,可以高效地实现支持向量机问题。SMO算法在这里用来更新优化
αi的值。
算法的基本思想就是,每次挑选出两个变量(假设为
α1和
α2),固定其它的变量(
αi),每次只更新
α1和
α2,循环迭代多次,尽可能接近最优解。
SMO算法其实就做了两件事:
- 变量的挑选方法:挑选出
α1和
α2。
- 两个变量二次规划的求解方法:更新
α1和
α2。
3.1 两个变量二次规划的求解方法
假设选择的两个变量为
α1和
α2,其它变量
αi,(i=3,...,N)是固定的,则SMO的最优化问题的子问题可以写成:
α1,α2min W(α1,α2)=21K11α12+21K22α22+y1y2K12α1α2−(α1+α2)+y1α1i=3∑NyiαiKi1+y2α2i=3∑NyiαiKi2s.t. α1y1+α2y2=−i=3∑Nyiαi=δ0≤αi≤C, i=1,2
其中
δ是常数,目标函数中省略了不含
α1,α2的常数项。
接下来,我们从约束条件入手:
α1y1+α2y2=δ0≤αi≤C, i=1,2
我们假设考虑为变量
α2的最优化问题。
假设问题的初始可行解为
α1old和
α2old,更新后的解为
α1new和
α2new,并记未经剪辑时
α2的最优解为
α2new_unc。(未经剪辑就是不一定满足
0≤α2new_unc≤C)
我们先求
α2new_unc,再对其约束得到
α2new。
我们假设最优值
α2new必须满足:
L≤α2new≤H
如上图所示,分两种情况讨论,L、H的值就是线段
lj和边界相交的点,可以求出:
-
y1=y2:
L=max(0,α2old+α1old−C),H=min(C,α2old+α1old)
-
y1=y2
L=max(0,α2old−α1old),H=min(C,C+α2old−α1old)
得到
L,H后,我们先放一放,先去求
α2new_unc的值:
记
g(x)=i=1∑NαiyiK(xi,x)+b Ei=g(xi)−yi
g(x)为预测值,
Ei为预测值与真实值之差。
则:
α2new_unc=α2old+y2ηE1−E2
其中,
η=K11+K22−2K12=∣∣ϕ(x1)−ϕ(x2)∣∣2
再求
α2new:
α2new=⎩⎨⎧H ,α2new_unc>Hα2new_unc ,L≤α2new_unc≤HL ,α2new_unc<L
根据
α1y1+α2y2=δ
得到:
α1new=α1old+y1y2(α2old−α2new)
于是得到新的
α1和
α2。
3.2 变量的挑选方法
SMO算法要挑选的两个变量,一个(
α1)是违反KKT条件的,另一个(
α2)的选择标准是希望能使
α2有足够大的变化。
3.2.1 第一个变量的选择
SMO称选择第一个变量的过程称为外层循环。外层循环在训练样本中选取违反KKT条件的样本点,并将其对应的
αi作为第一个变量。
KKT条件:
- 当
αi=0时,
yi(w⋅xi+b)−1≥0
- 当
0<αi<ξi时,
yi(w⋅xi+b)−1=0
一般把松弛变量统一为一个量,记为
C。则:
当
0<αi<C时,
yi(w⋅xi+b)−1=0
在检验选取过程中,外层循环首先遍历符合
0<αi<C条件的,再遍历符合
αi=0条件的,检验他们是否满足KKT条件,将第一个不满足KKT条件的的作为
α1。
3.2.2 第二个变量的选择
SMO称选择第二个变量的过程称为内层循环。我们的选择标准是希望能使
α2有足够大的变化。
根据公式:
α2new_unc=α2old+ηyi(E1−E2)
可知
α2new是依赖于
∣E1−E2∣的,所以:
- 当
E1≥0时,选择所有样本点中最小的
Ei作为
E2;
- 当
E1<0时,选择所有样本点中最大的
Ei作为
E2;
同时将挑选出的
Ei相对应的
αi作为第二个变量(
α2)
3.2.3 计算并更新阈值
b和差值
Ei
4 Python代码实现SVM
import numpy as np
class SVM:
def init_args(self, max_iter, features, labels):
self.max_iter = max_iter
self.m, self.n = features.shape
self.X = features
self.Y = labels
self.b = 0.0
self.alpha = np.ones(self.m)
self.E = [self.calc_E(i) for i in range(self.m)]
self.C = 1.0
def kernel(self, x1, x2):
sum = 0
for i in range(self.n):
sum += x1[i]*x2[i]
return sum
def calc_g(self, i):
g = self.b
for j in range(self.m):
g += self.alpha[j]*self.Y[j]*self.kernel(self.X[i], self.X[j])
return g
def calc_E(self, i):
return self.calc_g(i) - self.Y[i]
def judge_KKT(self, i):
if self.alpha[i]==0 and self.Y[i]*self.calc_g(i)>=1:
return True
elif 0<self.alpha[i]<self.C and self.Y[i]*self.calc_g(i)==1:
return True
return False
def get_alpha(self):
for i in range(self.m):
if self.judge_KKT(i) == False:
E1 = self.E[i]
if E1 >= 0:
j = min(range(self.m), key=lambda index : self.E[index])
else:
j = max(range(self.m), key=lambda index : self.E[index])
return i, j
def train(self, max_iter, features, labels):
self.init_args(max_iter, features, labels)
for i in range(self.max_iter):
i1, i2 = self.get_alpha()
if self.Y[i1] == self.Y[i2]:
L = max(0, self.alpha[i2]+self.alpha[i1]-self.C)
H = min(self.C, self.alpha[i2]+self.alpha[i1])
else:
L = max(0, self.alpha[i2]-self.alpha[i1])
H = min(self.C, self.alpha[i2]+self.alpha[i1]+self.C)
eta = self.kernel(self.X[i1], self.X[i1]) + self.kernel(self.X[i2], self.X[i2]) - 2*self.kernel(self.X[i1], self.X[i2])
alpha2_new_unc = self.alpha[i2] + self.Y[i2] * (self.E[i1] - self.E[i2]) / eta
if alpha2_new_unc > H:
alpha2_new = H
elif L <= alpha2_new_unc <= H:
alpha2_new = alpha2_new_unc
elif alpha2_new_unc < L:
alpha2_new = L
alpha1_new = self.alpha[i1] + self.Y[i1] * self.Y[i2] * (self.alpha[i2] - alpha2_new)
b1_new = -self.E[i1] - self.Y[i1] * self.kernel(self.X[i1], self.X[i1]) * (alpha1_new-self.alpha[i1]) - self.Y[i2] * self.kernel(self.X[i2], self.X[i1]) * (alpha2_new-self.alpha[i2])+ self.b
b2_new = -self.E[i2] - self.Y[i1] * self.kernel(self.X[i1], self.X[i2]) * (alpha1_new-self.alpha[i1]) - self.Y[i2] * self.kernel(self.X[i2], self.X[i2]) * (alpha2_new-self.alpha[i2])+ self.b
if 0 < alpha1_new < self.C:
b_new = b1_new
elif 0 < alpha2_new < self.C:
b_new = b2_new
else:
b_new = (b1_new + b2_new) / 2
self.alpha[i1] = alpha1_new
self.alpha[i2] = alpha2_new
self.b = b_new
self.E[i1] = self.calc_E(i1)
self.E[i2] = self.calc_E(i2)
print("Train: {0} iterations have been done.".format(self.max_iter))
from sklearn.svm import SVC
svc = SVC()
svc.fit(X, Y)
svc.score(Xt, Yt)
参考: