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在实现上述相关计算函数后，便可以第10.9.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行开始计算训练集中每个样本点对应的参数$\alpha_i$。第9~11行是依次取每个样本点，然后计算其预测值及预测值与真实值之间的误差。第13~16行是取另一个待优化的参数$\alpha_j$以及其对应的样本点，并计算真实值与预测值之间的误差。第18行是计算得到当前$\alpha_i,\alpha_j$对应的约束范围。
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行是根据当前这一组样本点来计算得到$\eta$。第6~12行是计算得到当前这一组样本点对应的参数$\alpha_i$和$\alpha_j$，并更新到参数列表中。第13~17行是根据$\alpha_i,\alpha_j$计算得到偏置$b$。第19~22行是用来统计$\alpha_i$不再发生变化时的次数。第23行是返回最终计算得到的参数和偏置。
到此，对于利用SMO算法来求解SVM中参数$\alpha$的代码就介绍完了。
10.10.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行用来保存每个二分类器计算得到的$\alpha$参数，因为在多分类问题中采用的是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行是分别将当前二分类器对应的参数放入到列表中。
在拟合得到分类器对应的参数后，便可以用其对新样本进行预测。这里需要先实现一个方法来完成对单个样本的预测，示例代码如下：
1 def predict_one_sample(self, x, y, alphas, bias):
2 y_score = np.sum(y * self.kernel(self._X, x) * alphas) + bias
3 return y_score进一步，实现对多个样本点的预测，示例代码如下：
1 def _predict_binary(self, X, y, alphas, bias):
2 y_scores = []
3 for x in X:
4 y_scores.append(self.predict_one_sample(x, y, alphas, bias))
5 return y_scores10.10.4 SVM多分类代码实现#
在完成SVM二分类的代码实现后，下一步便可以通过ovr的策略（这部分内容可以参见3.2节内容）来实现多分类支持向量机。在二分类的基础上，定义一个类方法来完成多个二分类器的拟合过程，示例代码如下：
1 def fit(self, X, y):
2 self._X = X
3 labelbin = LabelBinarizer(neg_label=-1)
4 Y = labelbin.fit_transform(y)
5 self.classes_ = labelbin.classes_
6 if Y.shape[1] == 1: #
7 Y = np.concatenate((-1 * Y, Y), axis=1)
8 self.n_classes = Y.shape[1]
9 self.Y = Y
10 for c in range(self.n_classes):
11 self._fit_binary(self._X, Y[:, c])
12 self.alphas = np.vstack((self.alphas))
13 self.bias = np.array(self.bias) 在上述代码中，第3~4行用于将原始的类别标签转化为one-hot形式的标签值，且正样本用1表示负样本用-1表示，最终将会得到一个形状为[m,c]的标签矩阵（m表示样本数量，c表示类别数量）。第5行代码表示得到原始的分类标签。
例如对于标签y=[0,0,1,1,2]来说，其转换后的结果为：
1 y = [0, 0, 1, 1, 2]
2 labelbin = LabelBinarizer(neg_label=-1)
3 Y = labelbin.fit_transform(y)
4 print(Y)
5 [[ 1 -1 -1]
6 [ 1 -1 -1]
7 [-1 1 -1]
8 [-1 1 -1]
9 [-1 -1 1]]
10 print(labelbin.classes_)
11 [0 1 2]同时，第6~7行判断当数据集为二分类时fit_transform处理后的结果并不是one-hot形式，因此需要得到上述类似的标签矩阵。第10~11是分别为每个类别拟合得到一个SVM二分类器。第12~13行是保存得到每个分类器对应的参数。
在完成多分类SVM的拟合过程后，便可以借助上面实现的二分类预测方法来完成多分类样本的预测过程，示例代码如下：
1 def predict(self, X, return_prob=False):
2 all_y_scores = []
3 for c in range(self.n_classes):
4 y_scores = self._predict_binary(X, self.Y[:, c], self.alphas[c], self.bias[c])
5 all_y_scores.append(y_scores)
6 all_y_scores = np.vstack((all_y_scores)).transpose() # [n_samples,n_classes]
7 prob = np.exp(all_y_scores) / np.sum(np.exp(all_y_scores), 1, keepdims=True)
8 y_pred = np.argmax(prob, axis=-1)
9 if return_prob:
10 return y_pred, prob
11 return y_pred在上述代码中，第3~5行是使用每个二分类器来对新输入的所有样本进行预测。第6~8行则是根据返回得到的预测值来计算得到每个样本最终对应的类别标签。这里需要注意的是，由于ovr策略是为每一个类别都构建一个二分类器，因此在根据式(10-120)来判断样本类别时不再是通过计算结果是否大于0来决定，而是通过每个分类器对应计算得到的最大值来决定。
在完成整个SVM分类器的实现后，便可以进行建模使用，示例代码如下：
1 def load_data():
2 x, y = load_iris(return_X_y=True)
3 x_train, x_test, y_train, y_test = train_test_split(x, y,
4 test_size=0.3, random_state=2022)
5 ss = StandardScaler()
6 x_train = ss.fit_transform(x_train)
7 x_test = ss.transform(x_test)
8 return x_train, x_test, y_train, y_test
9
10 def test_iris_classification():
11 x_train, x_test, y_train, y_test = load_data()
12 model = SVM(C=1., tol=0.001, max_passes=20, kernel='rbf')
13 model.fit(x_train, y_train)
14 y_pred = model.predict(x_test)
15 print("手动实现准确率：", accuracy_score(y_pred, y_test))
16
17 model = SVC(C=1, kernel='rbf')
18 model.fit(x_train, y_train)
19 y_pred = model.predict(x_test)
20 print("sklean上准确率：", accuracy_score(y_pred, y_test))
21
22 if __name__ == '__main__':
23 test_iris_classification()在上述代码运行结束后，便会得到如下所示的输出结果：
1 手动实现准确率： 0.978
2 sklean上准确率： 0.978到此，对于整个支持向量机从零实现的部分就介绍完了。
10.10.5 小结#
在本节中，我们首先介绍了两种常用核函数的实现方法；然后介绍了SMO算法中各个求解公式的实现过程以及其自身的编码实现过程；接着介绍了如何从零实现原始的二分类SVM分类器以及对新样本的预测过程；最后，基于ovr的思想实现了多分类的SVM分类器，并对整个模型的使用方法进行了示例。
总结一下，在本章中我们首先介绍了SVM的基本思想，即最大化分类间隔，接着详细介绍了SVM的原理及推导过程，进一步又介绍了SVM中的线性不可分问题及核函数的应用，然后介绍了SVM中软间隔的由来、拉格朗日乘数法和对偶优化问题，并分别介绍了SVM中硬间隔和软间隔的优化求解问题，和一种高效的凸二次规划求解方法SMO算法，最后还一步一步从零介绍了如何实现一个简单的SVM分类器。
引用#
[1] 李航.统计学习方法［M］.2版.北京： 清华大学出版社,2019.
[2] Andrew Ng,Machine Learning,Stanford University,CS229,Spring 2019.
[3] PEDREGOSA.scikitlearn: Machine Learning in Python［J］.JMLR 12,2011: 28252830.
[4] https://en.wikipedia.org/wiki/Lagrange_multiplier
[5] 徐小湛.高等数学学习手册［M］.北京： 科学出版社,2005.
[6] John C. Platt. Sequential Minimal Optimization: A Fast Algorithm for Training Support Vector Machines,Microsoft Research Technical Report MSRTR9814.
[7] https://en.wikipedia.org/wiki/Sequential_minimal_optimization
[8] https://en.wikipedia.org/wiki/Coordinate_descent
[9] 邓乃扬,田英杰. 数据挖掘中的新方法：支持向量机［M］. 北京： 科学出版社,2004.