10.10 从零实现SVM分类算法#

经过前面几节内容的介绍我们现在已经清楚了 SVM 的基本原理，并且根据第10.9.5节内容的讲解，对于SVM的求解过程也有了一定的认识，对于整个SVM内容的介绍就只差最后一步编码实现了。下面，我们将根据前面介绍的各个求解公式来一步一步介绍如何实现一个简单SVM分类器，完整示例代码可参见 AllBooKCode/Chapter10/C16_svm_impl.py 文件。

10.10.1 常见核函数实现#

根据10.8.4节内容介绍可知，常见的核函数有线性核函数、多项式核函数、高斯核函数和Sigmoid核函数等，这里我们以典型的线性核与高斯核函数为例进行实现，剩余的几种各位读者可以自行实现并验证。

1. 线性核实现

线性核的实现比较简单，根据式(10-127)可知实现如下：

1 def kernel_linear(X, x):
2     return np.dot(X, x)

在上述代码中，第1行中X为支持向量形状为[m,n]或者[n,]，x为待预测样本形状为[n,]。第2行是返回线性组合后的结果形状为[n,]。

2. 高斯核实现

高斯核的作用是将低维特征空间映射到无穷维的特征空间，根据式(10-125)可知实现如下：

1 def kernel_rbf(X, x):
2     if X.ndim > 1:
3         sigma = 1 / X.shape[1]
4     else:
5         sigma = 1 / X.shape[0]
6     k = -0.5 * np.sum((X - x) ** 2, axis=-1) / sigma ** 2
7     return np.exp(k)

在上述代码中，第2~5行用来计算高斯核函数里的参数$\sigma$（这里参考的是sklearn中的做法）。第6~7行是返回经高斯核函数变换后的结果。

10.10.2 SMO求解过程实现#

在介绍完核函数的实现部分后再来SMO算法的求解实现过程。首先，需要根据第10.9.3和第10.9.4小节中的内容来实现相关辅助函数。

根据式(10-120)可知，预测函数$f(x)$的编码实现为：

1 def f_x(X, y, alphas, x, b, kernel):
2     k = kernel(X, x)
3     r = alphas * y * k
4     return np.sum(r) + b

在上述代码中，第1行X和y为训练集，alphas为求解得到的参数，x为预测样本，b为偏置，kernel()为核函数。

进一步，根据式(10-147)和式(10-149)可知，$\eta$和$E_i$以及$\alpha_2^{new}$的编码实现为：

1 def compute_eta(x_1, x_2, kernel):
2     return kernel(x_1, x_1) - 2 * kernel(x_1, x_2) + kernel(x_2, x_2)
3 
4 def compute_E_i(f_x_i, y_i):
5     return f_x_i - y_i
6 
7 def compute_alpha_2(alpha_2, E_1, E_2, y_2, eta):
8     return alpha_2 + (y_2 * (E_1 - E_2) / eta)

同时，根据式(10-151)和式(10-152)可知，计算$L$和$H$的编码实现如下：

1 def compute_L_H(C, alpha_1, alpha_2, y_1, y_2):
2     L = np.max((0., alpha_2 - alpha_1))
3     H = np.min((C, C + alpha_2 - alpha_1))
4     if y_1 == y_2:
5         L = np.max((0., alpha_1 + alpha_2 - C))
6         H = np.min((C, alpha_1 + alpha_2))
7     return L, H

此时根据式(10-153)和式(10-154)便可以编码实现$\alpha$的计算，如下：

1 def clip_alpha_2(alpha_2, H, L):
2     if alpha_2 > H:
3         return H
4     if alpha_2 < L:
5         return L
6     return alpha_2
7 
8 def compute_alpha_1(alpha_1, y_1, y_2, alpha_2, alpha_2_old):
9     return alpha_1 + y_1 * y_2 * (alpha_2_old - alpha_2)

最后，根据式(10-160)~式(10-162)可以编码实现$b$的计算，如下：

 1 def compute_b1(b, E_1, y_1, alpha_1, alpha_1_old,
 2                x_1, y_2, alpha_2, alpha_2_old, x_2, kernel):
 3     p1 = b - E_1 - y_1 * (alpha_1 - alpha_1_old) * kernel(x_1, x_1)
 4     p2 = y_2 * (alpha_2 - alpha_2_old) * kernel(x_1, x_2)
 5     return p1 - p2
 6 
 7 def compute_b2(b, E_2, y_1, alpha_1, alpha_1_old,
 8                x_1, x_2, y_2, alpha_2, alpha_2_old, kernel):
 9     p1 = b - E_2 - y_1 * (alpha_1 - alpha_1_old) * kernel(x_1, x_2)
10     p2 = y_2 * (alpha_2 - alpha_2_old) * kernel(x_2, x_2)
11     return p1 - p2
12 
13 def clip_b(alpha_1, alpha_2, b1, b2, C):
14     if alpha_1 > 0 and alpha_1 < C:
15         return b1
16     if alpha_2 > 0 and alpha_2 < C:
17         return b2
18     return (b1 + b2) / 2