在上一篇博客中主要整理了支持向量机的数学原理，本文主要致力于用SMO算法来求解支持向量机最后的凸二次规划问题。
SMO算法
SMO算法的全称是Sequential minimal optimization，又名序列最小最优化算法，是用于快速求解支持向量机凸二次规划问题的算法。其基本思路是：如果所有变量的解都满足此最优化问题的KKT条件，那么这个最优化问题的解就得到了，因为KKT条件是该最优化问题的充分必要条件。否则，选择两个变量，固定其他变量，针对这两个变量构建一个二次规划问题，进行求解。SMO不断将原问题分解为子问题进行求解。整个SMO算法包括两个部分，求解两个变量二次规划的解析方法和选择变量的启发式方法。
两个变量二次规划问题求解
假设选择两个变量$\alpha_1$和$\alpha_2$，其他变量$\alpha_i$是固定的，于是SMO的最优化问题的子问题可以写成：
其中$\varsigma$是常数，$K_{ji} = K(x_i, x_j)$，目标函数中省略了不含$\alpha_1$和$\alpha_2$的常数项。因为$\alpha_1 y_1 + \alpha_2 y_2 = \varsigma$是常数，所以可以化为$\alpha_1 + \alpha_1\alpha_2y_2 = y_1\varsigma = C$，所以得到$\alpha_1 = C - \alpha_1\alpha_2y_2$，将此式子代入上面的目标函数，就将目标函数转为了关于单变量$\alpha_2$的优化问题。用转化后的目标函数对$\alpha_2$求偏导，可以得到
其中
还有一个需要注意的地方是，$\alpha_2^{new}$必须落在区间$[L, H]$之内，即
如果$y_1 \not = y_2$，则
如果$y_1 = y_2$，则
至于为什么$L$和$H$是这样得到的，其实画个图就可以出来了。然后根据$\alpha_2^{new}$的值可以求得$\alpha_1^{new}$的值：
变量的选择
选择第一个变量的过程称为外层循环，外层循环在训练样本中选取违反KKT条件最严重的样本点，具体的，检查训练样本点是否满足KKT条件，即
第二个变量选择使得$|E_1 - E_2|$最大的。
在每次完成两个变量的优化后，都要重新计算阈值$b$。当$0 < \alpha_i^{new} < C$时，由KKT条件可知：
于是可以得到
又因为
所以
同理
如果$\alpha_1^{new}$和$\alpha_2^{new}$同时满足$0 < \alpha_i^{new} < C$，那么$b_1^{new} = b_2^{new}$。如果$\alpha_1^{new}$和$\alpha_2^{new}$是0或者C，那么$\alpha_1^{new}$，$\alpha_2^{new}$和它们之间的数都符合，这时选择它们的中点。在每次完成两个变量的优化之后，还需要更新$E_i$，并且要用$b^{new}$来更新，
其中$S$是所有支持向量$x_j$的集合。
SMO算法实现
下面用python来实现SMO算法求解支持向量机。本文实现的是一个带高斯核的支持向量机，代码架构大致如下：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38class SVM:
def __init__(self):
pass
def _rbf_kernel(self, x1, x2, gamma):
"""径向基核函数"""
distance = np.linalg.norm(x1 - x2) ** 2
return np.exp(-gamma * distance)
def _cal_E_i(self, i):
"""计算Ei，注意只用到支持向量"""
pass
def _select_j(self, i):
"""SMO算法的第二层循环，根据选择的第一个点来选择第二个点"""
pass
def _update(self, i, j):
"""两个变量的二次规划算法"""
pass
def fit(self, X, y, C=1.0, gamma='auto', kernel='rbf', n_iteration=1000, tol=1e-3):
"""
对数据进行拟合
Args:
X: 特征矩阵
y: 标签矩阵
C: 对误差的惩罚项
gamma: rbf核函数的参数
n_iteration: 迭代次数
tol: 停止训练的误差精度
"""
pass
def predict(self, X):
"""预测"""
pass
对于支持向量机的求解，也就是fit函数，其代码主要分为以下几个步骤：
- 外层循环，选取第一个点
- 内层循环，根据第一个点选取第二个点
- 根据选取的两个点进行二变量问题的求解，更新参数
外层循环本文采取简单的遍历，内层循环采取上文所说的选择使得$|E_1 - E_2|$最大的样本，具体实现：1
2
3
4
5
6
7
8
9
10
11
12def _select_j(self, i):
"""SMO算法的第二层循环，根据选择的第一个点来选择第二个点"""
max_delta_E, j = 0, -1
n, _ = self.X.shape
for k in range(n):
if k == i:
continue
delta_E = abs(self.E[i, 0] - self.E[k, 0])
if delta_E > max_delta_E:
max_delta_E = delta_E
j = k
return j
二变量规划问题的求解，则按照上文所说的步骤一步步来，具体实现：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45def _update(self, i, j):
"""两个变量的二次规划算法"""
alpha_i_old, alpha_j_old = self.alpha[i, 0], self.alpha[j, 0]
if self.y[i, 0] != self.y[j, 0]:
L = max(0, self.alpha[j, 0] - self.alpha[i, 0])
H = min(self.C, self.C + self.alpha[j, 0] - self.alpha[i, 0])
else:
L = max(0, self.alpha[i, 0] + self.alpha[j, 0] - self.C)
H = min(self.C, self.alpha[i, 0] + self.alpha[j, 0])
if L == H:
return 0
eta = self.K[i, i] + self.K[j, j] - 2 * self.K[i, j]
if eta < 0:
return 0
alpha_j_new = alpha_j_old + self.y[j, 0] * (self.E[i, 0] - self.E[j, 0]) / eta
# 更新第一个变量
if alpha_j_new > H:
self.alpha[j, 0] = H
elif alpha_j_new < L:
self.alpha[j, 0] = L
else:
self.alpha[j, 0] = alpha_j_new
if abs(self.alpha[j, 0] - alpha_j_old) < 1e-5: # 更新量太小
return 0
# 更新第二个变量
self.alpha[i, 0] = alpha_i_old + self.y[i, 0] * self.y[j, 0] * (alpha_j_old - self.alpha[j, 0])
# 更新b
b1 = self.b - self.E[i, 0] - self.y[i, 0] * self.K[i, i] * (self.alpha[i, 0] - alpha_i_old) \
- self.y[j, 0] * self.K[j, i] * (self.alpha[j, 0] - alpha_j_old)
b2 = self.b - self.E[j, 0] - self.y[i, 0] * self.K[i, j] * (self.alpha[i, 0] - alpha_i_old) \
- self.y[j, 0] * self.K[j, j] * (self.alpha[j, 0] - alpha_j_old)
if self.alpha[i, 0] > 0 and self.alpha[i, 0] < self.C:
self.b = b1
elif self.alpha[j, 0] > 0 and self.alpha[j, 0] < self.C:
self.b = b2
else:
self.b = (b1 + b2) / 2
# 更新E
self.E[i, 0] = self._cal_E_i(i)[0]
self.E[j, 0] = self._cal_E_i(j)[0]
return 1
最后fit函数的具体实现为：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42def fit(self, X, y, C=1.0, gamma='auto', kernel='rbf', n_iteration=1000, tol=1e-3):
"""
对数据进行拟合
Args:
X: 特征矩阵
y: 标签矩阵
C: 对误差的惩罚项
gamma: rbf核函数的参数
n_iteration: 迭代次数
tol: 停止训练的误差精度
"""
n, m = X.shape
if gamma == 'auto': # 初始化gamma
self.gamma = 1 / m
else:
self.gamma = gamma
# 计算Gram相似性矩阵
self.K = np.zeros((n, n))
for i in range(n):
for j in range(i, n):
self.K[i, j] = self.K[j, i] = self._rbf_kernel(X[i], X[j], self.gamma)
self.X, self.y = X, y
self.C, self.tol = C, tol
# alpha和b都初始化为0
self.alpha, self.b = np.zeros((n, 1)), 0
# 初始化E
self.E = (np.dot((self.alpha * y).T, self.K) + self.b).reshape(-1, 1) - y
for iteration in range(n_iteration):
pairs_changed = 0
for i in range(n):
Ei = self.E[i, 0]
# 选取不满足KKT条件的第一个点
if ((y[i, 0] * Ei < -tol) and (self.alpha[i, 0] < C)) or ((y[i, 0] * Ei > tol) and (self.alpha[i, 0] > 0)):
j = self._select_j(i)
if j != -1:
pairs_changed += self._update(i, j)
if pairs_changed == 0:
print(f'no more pairs to change, stop at iteration {iteration + 1}')
break
最后用数据测试本文支持向量机的实现，在默认参数的情况下，训练集和测试集的Precision都能达到0.98以上，完整代码在这里。
总结
本文实现了SMO算法求解带高斯核的支持向量机，采取的两个变量的选择方法是《统计学习方法》一书中提到的简化版，所以还有待改进。另外测试算法性能的时候发现，参数C和参数gamma对性能的影响还是挺大的，它们的默认值参考了sklearn的支持向量机实现，分别是1.0和1/n_features。参数C是对支持向量机对误差的容忍度，比较好理解，而gamma直观来说会影响高斯核的分布情况，但通常情况下为什么要设置成这样一个值，还有待探究。