针对 “支持向量机” 将分为如下几个部分进行介绍:html
- 支持向量机 01 - 线性可分支持向量机和线性支持向量机
- 支持向量机 02 - SMO(序列最小化)(本篇)
- 支持向量机 03 - 非线性支持向量机
- 支持向量机 04 - SVR(支持向量回归)
本内容将介绍 SMO(序列最小化)算法,包含详细公式推导以及 Python 代码实现。
python
在学习本内容前,须要对支持向量机基本理论知识有必定了解,若是您还不了解,能够参阅 支持向量机 01 - 线性可分支持向量机和线性支持向量机。
web
咱们知道,支持向量机的学习问题能够形式化为求解凸二次规划问题。这样的凸二次规划问题具备全局最优解,而且有许多最优化算法能够用于这一问题的求解。可是当训练样本容量很大时,这些算法每每变得很是低效,以至没法使用。算法
SMO 算法是一种启发式算法,其基本思路是:若是全部变量的解都知足此最优化问题的 KKT 条件,那么这个最优化问题的解就获得了。由于 KKT 条件是该最优化问题的充要条件。不然,选择两个变量,固定其余变量,针对这两个变量构建一个二次规划问题。这个二次规划问题关于这两个变量的解应该更接近原始二次规划问题的解,由于这会使得原始二次规划问题的目标函数值变得更小。重要的是,这时子问题能够经过解析方法求解,这样就能够大大提升整个算法的计算速度。子问题有两个变量,一个是违反 KKT 条件最严重的那一个,另外一个由约束条件自动肯定。如此,SMO 算法将原问题不断分解为子问题并对子问题求解,进而达到求解原问题的目的。
缓存
1、SMO 算法描述
整个 SMO 算法包括两个部分:求解两个变量二次规划的解析方法和选择变量的启发式方法。数据结构
1.1 两个变量二次规划的求解方法
1.1.1 问题转化
在 支持向量机 01 - 线性可分支持向量机和线性支持向量机 中讲解了 SVM 的基本原理,了解到 SMO 算法要解以下凸二次规划的对偶问题app
αmax(i=1∑mαi−21i=1∑mj=1∑mαiαjy(i)y(j)x(i)⋅x(j))s.t.i=1∑mαiy(i)=0αi≥0,i=1,2,⋯,m(1)dom
添加一个负号,将最大值问题转换成最小值问题机器学习
αmin(21i=1∑mj=1∑mαiαjy(i)y(j)x(i)⋅x(j)−i=1∑mαi)s.t.i=1∑mαiy(i)=0αi≥0,i=1,2,⋯,m(2)ide
1.1.2 转换为二元函数
假设选择的两个变量为
α1 和
α2,其余变量
αi(i=3,4,⋯,m) 是固定的,式(2)可写为
W(α1,α2)=21i=1∑mj=1∑mαiαjy(i)y(j)x(i)⋅x(j)−i=1∑mαi=21i=1∑m(j=1∑2αiαjy(i)y(j)x(i)⋅x(j)+j=3∑mαiαjy(i)y(j)x(i)⋅x(j))−αi−α2−i=3∑mαi=21i=1∑2(j=1∑2αiαjy(i)y(j)x(i)⋅x(j)+j=3∑mαiαjy(i)y(j)x(i)⋅x(j))+21i=3∑m(j=1∑2αiαjy(i)y(j)x(i)⋅x(j)+j=3∑mαiαjy(i)y(j)x(i)⋅x(j))−αi−α2−i=3∑mαi=21i=1∑2j=1∑2αiαjy(i)y(j)x(i)⋅x(j)+i=1∑2j=3∑mαiαjy(i)y(j)x(i)⋅x(j)+21i=3∑mj=3∑mαiαjy(i)y(j)x(i)⋅x(j)−α1−α2−i=3∑mαi=21α12x(1)⋅x(1)+21α22x(2)x(2)+y(1)y(2)α1α2x(1)⋅x(2)+y(1)α1j=3∑mαjy(j)x(1)⋅x(j)+y(2)α2j=3∑mαjy(j)x(2)⋅x(j)−α1−α2+21i=3∑mj=3∑mαiαjy(i)y(j)x(i)⋅x(j)−i=3∑mαi(3)
为了描述方便,咱们定义以下符号:
f(x(i))=j=1∑mαjy(j)x(j)⋅x(i)+b=j=1∑mαjy(j)x(i)⋅x(j)+b(4)
vi=j=3∑mαjy(j)x(i)⋅x(j)=f(x(i))−j=1∑2αjy(j)x(i)⋅x(j)−b(5)
根据式(4)和(5),目标函数式(3)可写为
W(α1,α2)=21α12x(1)⋅x(1)+21α22x(2)⋅x(2)+y(1)y(2)α1α2x(1)⋅x(2)+y(1)α1v1+y(2)α2v2−α1−α2+constant(6)
其中
constant=21i=3∑mj=3∑mαiαjy(i)y(j)x(i)⋅x(j)−i=3∑mαi(7)
由于
constant 部分对
α1 和
α2 来讲,属于常数项,在求导的时候,直接变为 0,因此咱们不须要关心这部分的内容。
1.1.3 转换为一元函数
由于有
∑i=1mαiy(i)=0,可得
α1y(1)+α2y(2)=−i=3∑mαiy(i)(8)
对变量
α1 和
α2 来讲,
αi 和
y(i)(
i=3,4,⋯,m )可看做常数项,所以有
α1y(1)+α2y(2)=B(9)
其中,
B 为一个常数。将等式(9)两边同时乘以
y(1),得
α1=By(1)−α2y(1)y(2)=(B−α2y(2))y(1)(10)
将其带入式(6),获得只含变量
α2 的目标函数
W(α2)=21(B−α2y(2))2x(1)⋅x(1)+21α22x(2)⋅x(2)+y(2)(B−α2y(2))α2x(1)⋅x(2)+(B−α2y(2))v1+y(2)α2v2−(B−α2y(2))y(1)−α2+constant(11)
1.1.4 求解
α2new,unclipped
将式(11)对
α2 求导,得
∂α2∂W(α2)=x(1)⋅x(1)α2+x(2)⋅x(2)α2−2x(1)⋅x(2)α2−Bx(1)⋅x(1)y(2)+Bx(1)⋅x(2)y(2)+y(1)y(2)−v1y(2)+v2y(2)−1(12)
令其为
0,得
α2new=x(1)⋅x(1)+x(2)⋅x(2)−2x(1)⋅x(2)y(2)(y(2)−y(1)+B(x(1)⋅x(1)−x(1)⋅x(2))+v1−v2)(13)
定义以下符号
Ei=f(x(i))−y(i)(14)
η=x(1)⋅x(1)+x(2)⋅x(2)−2x(1)⋅x(2)(15)
其中
Ei 为偏差项,
η 为学习速率。
因为已知
B=α1oldy(1)+α2oldy(2)(16)
vi=f(x(i))−j=1∑2αjy(j)x(i)⋅x(j)−b(17)
将式(16)和(17)带入式(13), 将
α2new,unclipped 化简后得
α2new,unclipped=α2old+ηy(2)(E1−E2)(18)
1.1.5 求解
α2new(对原始
α2new,unclipped 解进行修剪)
上面求出的
α2new,unclipped 解是没有通过修剪的。咱们知道
αi 存在如下约束条件
⎩⎪⎨⎪⎧0<αi<Cα1y(1)+α2y(2)=B(19)
可知
α2new 的取值须要知足必定条件,假设为
L≤α2new≤H(20)
图-1 在二维平面上直观地表达了式(19)中的约束条件,从而可知
- 当
y(1)̸=y(2) 时,存在
α2old−α1old=k,得
L=max(0,α2old−α1old),H=min(C,C+α2old−α1old)(21)
- 当
y(1)=y(2) 时,存在
α2old+α1old=k,得
L=max(0,α2old+α1old−C),H=min(C,α2old+α1old)(22)
图-1
则通过修剪后,
α2new 的解为
α2new=⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧H,α2new,unclipped,L,α2new,unclipped>HL≤α2new,unclipped≤Hα2new,unclipped<L(23)
1.1.6 求解
α1new
由于存在如下等式条件
α1oldy(1)+α2oldy(2)=α1newy(1)+α2newy(2)(24)
从而得出
α1new=α1old+y(1)y(2)(α2old−α2new)(25)
1.1.7 求解阈值
bnew
在求得
α1new 和
α2new 后,须要根据它们更新
b。
-
若是
0<α1new<C
由 KKT 条件可知(在 支持向量机 01 - 线性可分支持向量机和线性支持向量机 中 2.2 节有进行介绍),当
0<αi<C 时,这个样本点是支持向量。所以,知足
y(1)(wTx(1)+b)=1(26)
上面公式两边同时乘以
y(1) 得
wTx(1)+b=y(1)(27)
将
w=∑i=1mαiy(i)x(i) 代入,得
i=1∑mαiy(i)x(i)⋅x(1)+b=y(1)(28)
由于咱们是根据
α1new 和
α2new 的值更新
b,整理可得:
b1new=y(1)−i=3∑mαiy(i)x(i)⋅x(1)−α1newy(1)x(1)⋅x(1)−α2newy(2)x(2)⋅x(1)(29)
其中有
y(1)−i=3∑mαiy(i)x(i)⋅x(1)=y(1)−i=1∑mαiy(i)x(i)⋅x(1)−bold+α1oldy(1)x(1)⋅x(1)+α2oldy(2)x(2)⋅x(1)+bold=y(1)−f(x(1))+α1oldy(1)x(1)⋅x(1)+α2oldy(2)x(2)⋅x(1)+bold=−E1+α1oldy(1)x(1)⋅x(1)+α2oldy(2)x(2)⋅x(1)+bold(30)
将式(30)代入式(29)得
b1new=bold−E1−y(1)(α1new−α1old)x(1)⋅x(1)−y(2)(α2new−α2old)x(2)x(1)(31)
-
若是
0<α2new<C
按照上面的步骤,一样可求得
b2new
b2new=bold−E2−y(1)(α1new−α1old)x(1)⋅x(2)−y(2)(α2new−α2old)x(2)x(2)(32)
若是
α1new 和
α2new 同时知足条件
0<αinew<C,则
bnew=b1new=b2new;若是
α1new 和
α2new 是
0 或者
C,那么
b1new 和
b2new 以及它们之间的数都是符合 KKT 条件的阈值,这是选择它们的重点做为
bnew。则
bnew 为
bnew=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧b1new,b2new,2(b1new+b2new),0<α1new<C0<α2new<Cotherwise(33)
1.1.8 变量更新算法的步骤
根据上面讲解的内容,咱们来整理一下变量更新算法的步骤:
E1=f(x(1))−y(1)=i=1∑mαiy(i)x(i)⋅x(1)+b−y(1)
E2=f(x(2))−y(2)=i=1∑mαiy(i)x(i)⋅x(2)+b−y(2)
- 步骤 2:计算
α2new 取值范围
⎩⎪⎪⎨⎪⎪⎧if(y(1)̸=y(2))L=max(0,α2old−α1old),H=min(C,C+α2old−α1old)if(y(1)=y(2))L=max(0,α2old+α1old−C),H=min(C,α2old+α1old)
η=x(1)⋅x(1)+x(2)⋅x(2)−2x(1)x(2)
- 步骤 4:求解
α2new,unclipped
α2new,unclipped=α2old+ηy(2)(E1−E2)
- 步骤 5:求解
α2new
α2new=⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧H,α2new,unclipped,L,α2new,unclipped>HL≤α2new,unclipped≤Hα2new,unclipped<L
- 步骤 6:求解
α1new
α1new=α1old+y(1)y(2)(α2old−α2new)
- 步骤 7:求解
b1new 和
b2new
b1new=bold−E1−y(1)(α1new−α1old)x(1)⋅x(1)−y(2)(α2new−α2old)x(2)⋅x(1)
b2new=bold−E2−y(1)(α1new−α1old)x(1)⋅x(2)−y(2)(α2new−α2old)x(2)⋅x(2)
- 步骤 8:求解
bnew
bnew=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧b1new,b2new,2(b1new+b2new),0<α1new<C0<α2new<Cotherwise
前面进行介绍时,咱们有假设选择两个变量
α1 和
α2,可是在实际实现算法时,咱们应该如何选择这两个变量呢?下面就来介绍一下。
1.2 变量
α1 和
α2 的选择方法
1.2.1 第 1 个变量的选择
SMO 称选择第 1 个变量的过程为外层循环。外层循环在训练样本中选取违反 KKT 条件最严重的样本点,并将其对应的变量做为第 1 个变量。具体地,检验训练样本点是否知足 KKT 条件,即(在 支持向量机 01 - 线性可分支持向量机和线性支持向量机 中 2.2 节有进行介绍)
αi=0⇔y(i)f(x(i))≥10<αi<C⇔y(i)f(x(i))=1αi=C⇔y(i)f(x(i))≤1(34)
从而可得知,如下几种状况将会不知足 KKT 条件
y(i)f(x(i))≥1andαi>0y(i)f(x(i))=1and(αi=0orαi=C)y(i)f(x(i))≤1andαi<C(35)
在检验过程当中,外层循环首先遍历全部知足条件
0<αi<C 的样本点,即在间隔边界上的支持向量点,检验它们是否知足 KKT 条件。若是这些样本点都知足 KKT 条件,那么遍历整个训练集,检验它们是否知足 KKT 条件。
1.2.2 第 2 个变量的选择
SMO 称选择第 2 个变量的过程为内循环。假设在外层循环中已经找到第 1 个变量
α1,如今要在内层循环中找第 2 个变量
α2。第 2 个变量选择的标准是但愿能使
α2 有足够大的变化。
由式(18)可知,
α2new 是依赖于
∣E1−E2∣ 的,为了加快计算速度,一种简单的作法是选择
α2,使其对应的
∣E1−E2∣ 最大。为了节省计算时间,能够将全部
Ei 值保存在一个列表中。
在特殊状况下,若是内层循环经过以上方法选择的
α2 不能使目标函数有足够的降低(等价于
∣α2new−α2∣ 很小),那么能够采用如下启发式规则继续选择
α2。遍历在间隔边界上的支持向量点,依次将其对应的变量做为
α2 试用,知道目标函数有足够的降低。若仍然找不到合适的
α2,那么遍历训练数据集;若仍是找不到合适的
α2,则放弃第 1 个
α2,再经过外层循环寻求另外的
α1。
由 Osuna 定理可知,只需选取的
α1 和
α2 中有一个不知足 KKT 条件,目标函数就会在迭代后减少。
1.3 SMO 算法的步骤
咱们总结一下 SMO 算法的步骤:
- 步骤 1:初始化
α 和
b,好比初始化
α 为零向量,
b 为 0。(注意这里的
α 是一个列向量)
- 步骤 2:选取优化变量
α1 和
α2,而后更新
α1、
α2 和
b。
- 步骤 3:若是知足中止条件(即前面提到的“全部变量的解都知足此最优化问题的 KKT 条件”)则结束;不然,跳到步骤 2。
2、Python 代码实现
如下是 Python 3 的代码实现:
import numpy as np
import random
import matplotlib.pyplot as plt
class OptStruct:
""" 数据结构,存储须要操做的数据 """
def __init__(self, input_mat, label_mat, c, toler):
""" :param input_mat: 样本特征值 :param label_mat: 样本标签值 :param c: 惩罚因子 :param toler: 容错率 """
self.x = input_mat
self.label_mat = label_mat
self.c = c
self.toler = toler
self.m = np.shape(input_mat)[0]
self.alphas = np.mat(np.zeros((self.m, 1)))
self.b = 0
self.e_cache = np.mat(np.zeros((self.m, 2)))
class SVM:
def __init__(self):
self.alphas = None
self.b = None
self.w = None
pass
@staticmethod
def calc_ek(os, k):
""" 计算偏差 E 值 :param os: :param k: :return: """
fxk = float(np.multiply(os.alphas, os.label_mat).T *\
(os.x * os.x[k, :].T)) + os.b
ek = fxk - float(os.label_mat[k])
return ek
@staticmethod
def select_j_rand(i, m):
""" 随机选择第二个 alpha :param i: 第一个 alpha 对应的索引值 :param m: alpha 参数的个数(即全部训练样本的总数) :return: 第二个 alpha 对应的索引值 """
j = i
while j == i:
j = int(random.uniform(0, m))
return j
@staticmethod
def select_j(i, os, ei):
""" 选择第二个 alpha :param i: 第一个 alpha 的索引值 :param os: 数据结构 :param ei: 第一个 alpha 对应样本点的偏差 :return: 返回第二个 alpha 的索引值 和 对应的偏差值 """
max_k = -1
max_delta_e = 0
ej = 0
os.e_cache[i] = [1, ei]
valid_e_cache_list = np.nonzero(os.e_cache[:, 0].A)[0]
if (len(valid_e_cache_list)) > 1:
for k in valid_e_cache_list:
if k == i:
continue
ek = SVM.calc_ek(os, k)
delta_e = abs(ei - ek)
if delta_e > max_delta_e:
max_k = k
max_delta_e = delta_e
ej = ek
return max_k, ej
else:
j = SVM.select_j_rand(i, os.m)
ej = SVM.calc_ek(os, j)
return j, ej
@staticmethod
def update_ek(os, k):
""" 更新偏差缓存 :param os: 数据结构 :param k: 须要更新偏差项对应的索引值 """
ek = SVM.calc_ek(os, k)
os.e_cache[k] = [1, ek]
@staticmethod
def clip_alpha(alpha_j, high, low):
""" 修剪 alpha_j :param alpha_j: alpha_j 的未修剪的值 :param high: alpha_j 的上限 :param low: alpha_j 的下限 :return: alpha_j 修剪后的值 """
if alpha_j > high:
alpha_j = high
elif alpha_j < low:
alpha_j = low
return alpha_j
def inner_l(self, i, os):
""" 选择第二个 alpha,并更新 alpha 和 b :param i: 第一个 alpha 对应样本点的索引值 :param os: 数据结构 :return: 是否有更新 alpha 对 """
ei = SVM.calc_ek(os, i)
if ((os.label_mat[i] * ei < -os.toler) and (os.alphas[i] < os.c)) or \
((os.label_mat[i] * ei > os.toler) and (os.alphas[i] > 0)):
j, ej = SVM.select_j(i, os, ei)
alpha_iold = os.alphas[i].copy()
alpha_jold = os.alphas[j].copy()
if os.label_mat[i] != os.label_mat[j]:
low = max(0, os.alphas[j] - os.alphas[i])
high = min(os.c, os.c + os.alphas[j] - os.alphas[i])
else:
low = max(0, os.alphas[j] + os.alphas[i] - os.c)
high = min(os.c, os.alphas[j] + os.alphas[i])
if low == high:
print('low == high')
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
if eta >= 0:
print('eta >= 0')
return 0
os.alphas[j] -= os.label_mat[j] * (ei - ej) / eta
os.alphas[j] = SVM.clip_alpha(os.alphas[j], high, low)
SVM.update_ek(os, j)
if abs(os.alphas[j] - alpha_jold) < 0.00001:
print('j not moving enough')
return 0
os.alphas[i] += os.label_mat[j] * os.label_mat[i] * \
(alpha_jold - os.alphas[j])
SVM.update_ek(os, i)
b1 = os.b - ei - os.label_mat[i] * (os.alphas[i] - alpha_iold) * \
os.x[i, :] * os.x[i, :].T