第9.9节 从零实现SVM分类算法

各位朋友大家好，欢迎来到月来客栈，我是掌柜空字符。

本期推送内容目录如下，如果本期内容对你有所帮助，欢迎点赞、转发支持掌柜！

• 9.9 从零实现SVM分类算法
• 9.9.1 常见核函数实现
• 9.9.2 SMO求解过程实现
• 9.9.3 SVM二分类代码实现
• 9.9.4 SVM多分类代码实现
• 9.9.5 小结

9.9 从零实现SVM分类算法

经过前面几节内容的介绍我们现在已经清楚了 SVM 的基本原理，并且根据第9.8.5节内容的讲解，对于SVM的求解过程也有了一定的认识，对于整个SVM内容的介绍就只差最后一步编码实现了。下面，笔者将根据前面介绍的各个求解公式来一步一步介绍如何实现一个简单SVM分类器。

9.9.1 常见核函数实现

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

1. 线性核实现

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

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

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

2. 高斯核实现

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

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行用来计算高斯核函数里的参数（参考的是sklearn中的做法）；第6~7行是返回经高斯核函数变换后的结果。

9.9.2 SMO求解过程实现

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

根据式(9.24)可知，预测函数的编码实现为：

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()为核函数。

进一步，根据式(9.147)和(9-149)可知，和以及的编码实现为：

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)

同时，根据式(9.151)和(9-152)可知，计算和的编码实现如下：

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

此时根据式(9.153)和式(9.154)便可以编码实现的计算，如下：

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)

最后，根据式(9.160)到式(9.162)可以编码实现的计算，如下：

 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

在实现上述相关计算函数后，便可以第9.8.5节中的伪代码来实现SVM参数的求解过程。由于这部分代码较长，所以这里笔者分成块进行介绍。

 1 def smo(C, tol, max_passes, data_x, data_y, kernel):
 2     m, n = data_x.shape
 3     b, passes = 0., 0
 4     alphas = np.zeros(shape=(m))
 5     alphas_old = np.zeros(shape=(m))
 6     while passes < max_passes:
 7         num_changed_alphas = 0
 8         for i in range(m):
 9             x_i, y_i, alpha_i = data_x[i], data_y[i], alphas[i]
10             f_x_i = f_x(data_x, data_y, alphas, x_i, b, kernel)
11             E_i = compute_E_i(f_x_i, y_i)
12             if ((y_i * E_i < -tol and alpha_i < C) or (y_i * E_i > tol and alpha_i > 0.)):
13                 j = select_j(i, m)
14                 x_j, y_j, alpha_j = data_x[j], data_y[j], alphas[j]
15                 f_x_j = f_x(data_x, data_y, alphas, x_j, b, kernel)
16                 E_j = compute_E_i(f_x_j, y_j)
17                 alphas_old[i], alphas_old[j] = alpha_i, alpha_j
18                 L, H = compute_L_H(C, alpha_i, alpha_j, y_i, y_j)

在上述代码中，第2~5行用来初始化需要进行求解的相关参数；第8行开始计算训练集中每个样本点对应的参数；第9~11行是依次去每个样本点，然后计算其预测值及预测值与真实值之间的误差；第13~16行是取另一个待优化的参数以及其对应的样本点，并计算真实值与预测值之间的误差；第18行是计算得到当前对应的约束范围。

 1                 if L == H:
 2                     continue
 3                 eta = compute_eta(x_i, x_j, kernel)
 4                 if eta <= 0:
 5                     continue
 6                 alpha_j = compute_alpha_2(alpha_j, E_i, E_j, y_j, eta)
 7                 alpha_j = clip_alpha_2(alpha_j, H, L)
 8                 alphas[j] = alpha_j
 9                 if np.abs(alpha_j - alphas_old[j]) < 10e-5:
10                     continue
11                 alpha_i = compute_alpha_1(alpha_i, y_i, y_j, alpha_j, alphas_old[j])
12                 alphas[i] = alpha_i
13                 b1 = compute_b1(b, E_i, y_i, alpha_i, alphas_old[i],
14                                 x_i, y_i, alpha_j, alphas_old[j], x_j, kernel)
15                 b2 = compute_b2(b, E_j, y_i, alpha_i, alphas_old[i],
16                                 x_i, x_j, y_j, alpha_j, alphas_old[j], kernel)
17                 b = clip_b(alpha_i, alpha_j, b1, b2, C)
18                 num_changed_alphas += 1
19         if num_changed_alphas == 0:
20             passes += 1
21         else:
22             passes = 0
23     return alphas, b

在上述代码中，第3行是根据当前这一组样本点来计算得到；第6~12行是计算得到当前这一组样本点对应的参数和，并更新到参数列表中；第13~17是根据计算得到偏置；第19~22行是用来统计不再发生变化时的次数；第23行是返回最终计算得到的参数和偏置。

到此，对于利用SMO算法来求解SVM中参数的代码就介绍完了。

9.9.3 SVM二分类代码实现

在完成SMO算法的求解编码过程后，便可以来编码实现一个基础的SVM二分类器。首先，定义一个类并完成对应的初始化方法，代码如下：

 1 class SVM(object):
 2     def __init__(self, C, tol, kernel='rbf', max_passes=20):
 3         self.C = C 
 4         self.tol = tol  
 5         self.max_passes = max_passes 
 6         self.alphas = [] 
 7         self.bias = [] 
 8         if kernel == 'rbf':
 9             self.kernel = kernel_rbf
10         elif kernel == 'linear':
11             self.kernel = kernel_linear
12         else:
13             raise ValueError(f"核函数{kernel}未实现")

在上述代码中，第3~5行是上面smo函数中对应的参数，这里就不再赘述；第6~7行用来保存每个二分类器计算得到的参数，因为在多分类问题中采用的是ovr策略，同时bias用来保存每个分类器对应的偏置；第8~13行则是取对应的核函数。

进一步，实现二分类器的拟合过程，代码如下：

1     def _fit_binary(self, X, y):
2         alphas, bias = smo(C=self.C, tol=self.tol,
3                            max_passes=self.max_passes,
4                            data_x=X, data_y=y, kernel=self.kernel)
5         self.alphas.append(alphas)
6         self.bias.append(bias)

在上述代码中，第1行X是原始的训练特征，y是每个样本对应的标签值（只含有 -1 和 +1）；第2~4行是根据SMO算法来求解得到当前二分类器对应的参数；第5~6行是分别将当前二分类器对应的参数放入到列表中。