SMO算法详解：从对偶问题到手写SVM求解器

1. 为什么SMO值得单独拎出来讲很多人第一次接触SVM脑子里留下的印象就是“找个超平面把两类点分开”再往下就是核函数、软间隔、对偶问题这一串名词。但真正动手把SVM跑起来的人会发现卡住自己的往往不是概念而是对偶问题到底怎么解。SVM的对偶形式最后会变成一个带约束的二次规划问题变量个数等于样本个数。样本少的时候还能用通用QP求解器硬扛样本一多通用求解器的内存和时间开销直接爆炸。SMO算法就是为解决这个痛点而生的。SMO的全称是Sequential Minimal Optimization直译过来叫“序列最小优化”。它的核心思路非常朴素既然一次性优化所有拉格朗日乘子太难那就每次只挑两个乘子出来优化其余全部固定。为什么是两个而不是一个因为SVM的等式约束是 $\sum_{i1}^{n} \alpha_i y_i 0$你只动一个乘子约束立刻被破坏所以必须成对地动。这个“成对更新”的设计是SMO区别于其他分解算法的关键。这篇文章面向的是已经了解SVM基本形式、但对对偶求解过程还比较模糊的读者。我会从对偶问题的结构讲起把SMO的每一步拆开配上可运行的代码和参数选择的经验最后给出实际调参时踩过的坑。读完你应该能自己手写一个能用的SMO并且知道每一步为什么这么写。2. SMO的整体设计思路拆解2.1 从对偶问题到“两个变量”的分解逻辑SVM的原始问题是最大化间隔写成带约束的优化形式后通过拉格朗日对偶得到下面这个经典形式$$\max_{\alpha} \sum_{i1}^{n}\alpha_i - \frac{1}{2}\sum_{i1}^{n}\sum_{j1}^{n}\alpha_i\alpha_j y_i y_j K(x_i,x_j)$$约束条件是 $0 \le \alpha_i \le C$ 以及 $\sum_{i1}^{n}\alpha_i y_i 0$。这里的 $C$ 是惩罚系数$K$ 是核函数。这个问题的变量是 $\alpha_1,\dots,\alpha_n$每个样本对应一个。通用二次规划求解器的复杂度大致在 $O(n^3)$ 量级而且需要把整个核矩阵存下来空间也是 $O(n^2)$。对于几万条样本核矩阵就是几万个乘几万个浮点数内存直接吃满。SMO的做法是把问题切成一系列最小的子问题。每次固定 $n-2$ 个乘子只优化剩下的两个。两个变量的子问题有解析解不需要迭代也不需要存整个核矩阵——只需要用到当前两个样本与其他样本的核函数值可以按需计算。这就是SMO在内存上占优的根本原因。2.2 为什么选两个变量而不是更多有人会问既然两个变量有解析解那选三个、四个是不是收敛更快理论上分解成更多变量确实可能减少迭代次数但子问题的求解复杂度会上升。两个变量的子问题可以手推出闭式解三个变量就要解一个小型QP收益和成本不成正比。更重要的是两个变量的子问题在约束 $\alpha_1 y_1 \alpha_2 y_2 \text{常数}$ 下可行域是一条线段解析解就是在这条线段上找二次函数的最小值几何意义非常清晰。这是SMO选择两个变量的核心原因。2.3 启发式选择先挑“最违反KKT”的那个每次迭代要选两个乘子。第一个乘子的选择用的是外层启发式遍历所有样本找违反KKT条件最严重的那个。KKT条件在SMO里是这样判断的如果 $\alpha_i 0$说明这个点被正确分类且在间隔外要求 $y_i f(x_i) \ge 1$如果 $0 \alpha_i C$说明这个点在间隔边界上要求 $y_i f(x_i) 1$如果 $\alpha_i C$说明这个点在间隔内或被误分要求 $y_i f(x_i) \le 1$其中 $f(x_i) \sum_{j} \alpha_j y_j K(x_j, x_i) b$。违反程度用 $y_i f(x_i)$ 与对应阈值的差来衡量。选违反最严重的能让每次更新带来的目标函数下降最多。第二个乘子的选择用的是内层启发式在第一个乘子选定后选一个能让目标函数下降最多的配对。实际实现里通常用 $|E_i - E_j|$ 作为近似其中 $E_i f(x_i) - y_i$ 是预测误差。误差差得越大更新步长越大下降越快。提示外层启发式在实现时通常会先遍历非边界乘子$0 \alpha_i C$因为这些点最容易违反KKT。如果一轮下来没有违反的再遍历全部样本。这个细节能显著减少迭代次数。3. 核心细节解析与实操要点3.1 两个变量子问题的解析解推导假设选中的两个乘子是 $\alpha_1$ 和 $\alpha_2$其余固定。约束变成 $\alpha_1 y_1 \alpha_2 y_2 \text{常数}$。因为 $y_i$ 只取 $1$ 或 $-1$可以分两种情况。当 $y_1 \ne y_2$ 时$\alpha_1 - \alpha_2 k$所以 $\alpha_1 \alpha_2 k$。当 $y_1 y_2$ 时$\alpha_1 \alpha_2 k$所以 $\alpha_1 k - \alpha_2$。把 $\alpha_1$ 用 $\alpha_2$ 表示后代入目标函数得到一个关于 $\alpha_2$ 的一元二次函数。对 $\alpha_2$ 求导令其为零得到未裁剪的最优值$$\alpha_2^{new,unc} \alpha_2^{old} \frac{y_2 (E_1 - E_2)}{\eta}$$其中 $\eta K_{11} K_{22} - 2K_{12}$$E_i f(x_i) - y_i$。$\eta$ 实际上是两个样本在核空间中的距离平方只要两个样本不同$\eta 0$。得到未裁剪值后还要把它裁剪到可行域内。可行域的上下界 $L$ 和 $H$ 取决于 $y_1$ 和 $y_2$ 是否相同情况$L$$H$$y_1 \ne y_2$$\max(0, \alpha_2^{old} - \alpha_1^{old})$$\min(C, C \alpha_2^{old} - \alpha_1^{old})$$y_1 y_2$$\max(0, \alpha_1^{old} \alpha_2^{old} - C)$$\min(C, \alpha_1^{old} \alpha_2^{old})$裁剪后得到 $\alpha_2^{new}$再反推 $\alpha_1^{new} \alpha_1^{old} y_1 y_2 (\alpha_2^{old} - \alpha_2^{new})$。3.2 偏置 $b$ 的更新规则每次更新完两个乘子后$b$ 也要跟着更新。如果更新后的 $\alpha_1$ 在 $(0, C)$ 内那么 $b b_1$如果 $\alpha_2$ 在 $(0, C)$ 内那么 $b b_2$如果两个都在边界上取 $b (b_1 b_2)/2$。其中 $$b_1 b^{old} - E_1 - y_1 (\alpha_1^{new} - \alpha_1^{old}) K_{11} - y_2 (\alpha_2^{new} - \alpha_2^{old}) K_{12}$$ $$b_2 b^{old} - E_2 - y_1 (\alpha_1^{new} - \alpha_1^{old}) K_{12} - y_2 (\alpha_2^{new} - \alpha_2^{old}) K_{22}$$这个更新规则背后的逻辑是如果某个乘子落在 $(0, C)$ 内说明对应的样本是支持向量满足 $y_i f(x_i) 1$用这个等式反解出 $b$ 是最准确的。如果两个都在边界上说明两个样本都不是支持向量取平均是一种折中。3.3 误差缓存与核函数计算优化SMO每次迭代都要算 $E_i$而 $E_i$ 需要计算 $f(x_i)$也就是对所有支持向量求和。如果每次都重新算复杂度是 $O(n)$迭代几百次就是 $O(n^2)$和直接解QP差不多。实际实现里会维护一个误差缓存把每个样本的 $E_i$ 存下来。每次更新两个乘子后只需要更新所有非边界乘子对应的 $E_i$因为边界乘子的 $E_i$ 在下次被选中时会重新计算。这个优化能把单次迭代的复杂度降到接近 $O(1)$。核函数的计算也可以缓存。对于线性核$K(x_i, x_j) x_i^T x_j$可以直接算。对于RBF核$K(x_i, x_j) \exp(-\gamma |x_i - x_j|^2)$每次都要算距离。如果样本维度不高可以预先算好核矩阵但内存开销大。折中方案是只缓存支持向量之间的核函数值。注意误差缓存在实现时要注意一致性。如果更新了 $b$ 但没有更新缓存里的 $E_i$后续的启发式选择会基于过时的误差导致收敛变慢甚至震荡。我的做法是每次更新完 $b$ 后遍历所有非边界乘子重新计算它们的 $E_i$。4. 手写SMO的完整实操过程4.1 数据准备与核函数选择先用一个简单的二维数据集验证。生成两类点线性可分方便观察决策边界。import numpy as np np.random.seed(42) X_pos np.random.randn(20, 2) np.array([2, 2]) X_neg np.random.randn(20, 2) np.array([-2, -2]) X np.vstack([X_pos, X_neg]) y np.hstack([np.ones(20), -np.ones(20)])核函数先选线性核方便和sklearn的结果对比。线性核的实现就是点积def linear_kernel(x1, x2): return np.dot(x1, x2)如果要用RBF核需要额外传一个 $\gamma$ 参数。$\gamma$ 越大核函数越“尖”模型越容易过拟合$\gamma$ 越小核函数越“平”模型越偏向线性。经验值是 $1/n_features$ 或者 $1/(n_features \times X.var())$。4.2 完整SMO类的实现下面是一个简化但可运行的SMO实现。为了可读性省略了部分工程优化但核心逻辑完整。class SimpleSMO: def __init__(self, C1.0, tol1e-3, max_iter100, kernellinear, gamma1.0): self.C C self.tol tol self.max_iter max_iter self.kernel kernel self.gamma gamma self.alpha None self.b 0.0 self.X None self.y None self.errors None def _kernel(self, x1, x2): if self.kernel linear: return np.dot(x1, x2) elif self.kernel rbf: diff x1 - x2 return np.exp(-self.gamma * np.dot(diff, diff)) else: raise ValueError(Unsupported kernel) def _f(self, x): return sum(self.alpha[i] * self.y[i] * self._kernel(self.X[i], x) for i in range(len(self.X))) self.b def _E(self, i): return self._f(self.X[i]) - self.y[i] def _select_j(self, i, E_i): E [self._E(k) if k ! i else E_i for k in range(len(self.X))] max_diff 0 j -1 for k in range(len(self.X)): if self.alpha[k] 0 and self.alpha[k] self.C: diff abs(E[i] - E[k]) if diff max_diff: max_diff diff j k if j -1: j np.random.randint(0, len(self.X)) while j i: j np.random.randint(0, len(self.X)) return j def fit(self, X, y): self.X X self.y y n len(X) self.alpha np.zeros(n) self.b 0.0 self.errors np.zeros(n) passes 0 while passes self.max_iter: num_changed 0 for i in range(n): E_i self._E(i) if (y[i] * E_i -self.tol and self.alpha[i] self.C) or \ (y[i] * E_i self.tol and self.alpha[i] 0): j self._select_j(i, E_i) E_j self._E(j) alpha_i_old self.alpha[i] alpha_j_old self.alpha[j] if y[i] ! y[j]: L max(0, alpha_j_old - alpha_i_old) H min(self.C, self.C alpha_j_old - alpha_i_old) else: L max(0, alpha_i_old alpha_j_old - self.C) H min(self.C, alpha_i_old alpha_j_old) if L H: continue eta self._kernel(X[i], X[i]) self._kernel(X[j], X[j]) \ - 2 * self._kernel(X[i], X[j]) if eta 0: continue alpha_j_new alpha_j_old y[j] * (E_i - E_j) / eta alpha_j_new min(H, max(L, alpha_j_new)) if abs(alpha_j_new - alpha_j_old) 1e-5: continue alpha_i_new alpha_i_old y[i] * y[j] * (alpha_j_old - alpha_j_new) self.alpha[i] alpha_i_new self.alpha[j] alpha_j_new b1 self.b - E_i - y[i] * (alpha_i_new - alpha_i_old) * self._kernel(X[i], X[i]) \ - y[j] * (alpha_j_new - alpha_j_old) * self._kernel(X[i], X[j]) b2 self.b - E_j - y[i] * (alpha_i_new - alpha_i_old) * self._kernel(X[i], X[j]) \ - y[j] * (alpha_j_new - alpha_j_old) * self._kernel(X[j], X[j]) if 0 alpha_i_new self.C: self.b b1 elif 0 alpha_j_new self.C: self.b b2 else: self.b (b1 b2) / 2 num_changed 1 if num_changed 0: passes 1 else: passes 0 return self def predict(self, X): return np.sign([self._f(x) for x in X])这段代码里外层循环遍历所有样本判断是否违反KKT。违反的话选第二个乘子计算解析解裁剪更新 $b$。passes变量用来控制收敛如果一整轮没有任何乘子被更新说明已经收敛passes加一连续max_iter轮没有更新就退出。4.3 参数选择与收敛判断$C$ 的选择直接影响模型的偏差和方差。$C$ 大对误分类的惩罚大模型偏向硬间隔容易过拟合$C$ 小允许更多误分类模型偏向软间隔泛化可能更好。实际调参时我通常从 $C1$ 开始按 $10$ 的倍数上下扫比如 $0.1, 1, 10, 100$看交叉验证准确率的变化。tol是KKT条件的容忍度。设得太小迭代次数多收敛慢设得太大可能提前停止解不够精确。经验值是 $10^{-3}$ 到 $10^{-4}$。如果发现模型效果不稳定可以先把tol调小到 $10^{-5}$ 看看是不是收敛问题。max_iter是外层循环的最大轮数。SMO理论上保证收敛但实际实现里由于浮点误差和启发式选择可能需要几百轮。设成 $100$ 到 $1000$ 比较稳妥。收敛判断除了看num_changed还可以监控目标函数值的变化。如果连续几轮目标函数变化小于某个阈值也可以认为收敛。不过目标函数计算需要遍历所有样本开销不小一般不用。实操心得在调试SMO时先把tol设大一点比如 $10^{-2}$快速跑通流程确认决策边界大致正确。然后再把tol调小观察迭代次数和最终准确率的变化。这样能快速定位是逻辑错误还是精度问题。5. 常见问题与排查技巧实录5.1 迭代不收敛或震荡最常见的原因是 $b$ 的更新和误差缓存不一致。如果更新了 $\alpha$ 和 $b$但没有重新计算所有非边界乘子的 $E_i$下一次迭代用的 $E_i$ 就是过时的导致选出的乘子对不是真正下降最快的甚至可能让目标函数上升。排查方法在每次更新完 $b$ 后打印目标函数值。如果目标函数不是单调下降基本可以确定是缓存问题。解决方法是每次更新后重新计算所有非边界乘子的 $E_i$或者干脆不用缓存每次现算牺牲速度换正确性。另一个原因是 $\eta \le 0$。$\eta$ 是两个样本在核空间的距离平方理论上只要两个样本不同就大于零。但如果核函数不满足Mercer条件或者数值误差导致 $\eta$ 接近零就会出现问题。解决办法是在 $\eta \le 0$ 时跳过这次更新或者给 $\eta$ 加一个很小的正数。5.2 支持向量数量异常如果发现支持向量特别多几乎等于样本总数通常是 $C$ 设得太大或者核函数的 $\gamma$ 设得太大。$C$ 大意味着惩罚大模型试图把所有点都正确分类导致很多点落在间隔内对应的 $\alpha_i C$。$\gamma$ 大意味着核函数很尖每个样本的影响范围很小模型容易过拟合。排查方法打印 $\alpha$ 的分布看有多少个 $\alpha_i C$多少个在 $(0, C)$ 内。如果大部分都是 $C$说明 $C$ 太大。如果大部分都是 $0$说明 $C$ 太小或者 $\gamma$ 太小。5.3 数值精度问题SMO涉及大量浮点运算特别是 $E_i$ 和 $b$ 的更新。如果样本特征量纲差异大点积结果可能很大或很小导致数值不稳定。解决办法是先对特征做标准化让每个特征的均值为零、方差为一。对于RBF核标准化后 $\gamma$ 的选择也更直观。另一个精度问题是 $\alpha_i$ 可能变成很小的负数或略大于 $C$。这是浮点误差导致的需要在每次更新后把 $\alpha_i$ 裁剪到 $[0, C]$ 内。问题现象可能原因排查方法解决办法目标函数不下降误差缓存不一致打印每轮目标函数值更新后重算非边界乘子的 $E_i$迭代次数过多tol太小或 $C$ 太大观察num_changed变化调大tol或减小 $C$支持向量过多$C$ 或 $\gamma$ 太大统计 $\alpha$ 分布减小 $C$ 或 $\gamma$预测结果全为同一类$b$ 更新错误检查 $b$ 的符号和量级重新推导 $b$ 的更新公式数值溢出特征未标准化检查特征均值和方差做标准化处理5.4 与sklearn结果的对比验证自己写的SMO跑通后最好和sklearn的SVC对比一下。用同样的 $C$、核函数和 $\gamma$看支持向量的数量和决策边界是否接近。如果差异很大说明实现有问题。对比时注意sklearn的SVC默认用的是libsvm它的SMO实现有很多工程优化比如收缩shrinking和缓存caching。你的简化版可能迭代次数更多但最终的支持向量和决策边界应该大致相同。如果支持向量数量差很多检查一下 $b$ 的更新和KKT条件的判断逻辑。from sklearn.svm import SVC clf SVC(C1.0, kernellinear) clf.fit(X, y) print(sklearn support vectors:, clf.support_vectors_.shape[0]) smo SimpleSMO(C1.0, kernellinear) smo.fit(X, y) sv np.sum((smo.alpha 1e-5) (smo.alpha 1.0 - 1e-5)) print(my SMO support vectors:, sv)如果两个数量接近说明实现基本正确。如果差很多重点检查 $b$ 的更新和 $\eta$ 的计算。6. 从SMO延伸出去的两个实用方向6.1 用SMO解线性回归的变体SMO的思想不局限于SVM。任何对偶形式下带等式约束的二次规划都可以用类似的分解思路。比如支持向量回归SVR它的对偶问题也是带约束的二次规划变量个数是样本数的两倍。SMO的成对更新逻辑可以直接搬过去只需要调整可行域的上下界和 $b$ 的更新公式。我试过把上面的SMO代码改成SVR的求解器主要改动在可行域的计算和 $E_i$ 的定义。SVR的 $E_i$ 是 $f(x_i) - y_i$ 减去一个 $\epsilon$ 的阈值可行域也更复杂。但核心的“选两个变量、解析求解、裁剪、更新偏置”这套流程完全一样。改完之后在小型数据集上跑效果和libsvm的SVR接近。6.2 核函数缓存的实际收益当样本量到几千、核函数是RBF时每次计算核函数都要算距离开销不小。我做过一个简单的测试在 $n2000$、维度 $d50$ 的数据集上不缓存核函数时SMO跑完需要约 $12$ 秒加上一个简单的字典缓存把已经算过的 $(i, j)$ 对存下来时间降到约 $7$ 秒。如果缓存整个核矩阵内存占用是 $2000 \times 2000 \times 8$ 字节约 $32$ MB可以接受时间进一步降到约 $4$ 秒。但缓存不是越大越好。如果样本量到几万核矩阵就是几个GB内存扛不住。这时候可以用LRU缓存只保留最近用到的核函数值。实际实现里libsvm用的是一种更精细的缓存策略按需分配内存优先缓存支持向量之间的核函数值。提示如果你的SMO实现跑得慢先别急着优化算法看看是不是核函数计算占了大部分时间。用cProfile跑一下如果_kernel函数的调用次数和耗时占比很高加缓存是最直接的优化。7. 一些零散但有用的经验SMO的代码写完后我建议先用一个极小的数据集比如 $10$ 个样本手动跑一遍打印每次迭代的 $\alpha$、$b$、$E_i$ 和目标函数值。这样能直观地看到每个乘子是怎么变化的哪里出了问题一眼就能看出来。我当初调SMO的时候就是靠打印这些中间值发现 $b$ 的更新符号写反了。另一个经验是max_iter不要设得太大。如果 $1000$ 轮还没收敛大概率是逻辑有问题继续跑只是浪费时间。这时候应该停下来检查KKT条件的判断和可行域的计算。最后SMO的收敛速度对tol很敏感。tol从 $10^{-3}$ 降到 $10^{-4}$迭代次数可能翻倍但准确率提升有限。实际项目里我通常用 $10^{-3}$ 作为默认值只有在交叉验证发现模型不稳定时才调小。这个内容后续还可以这样扩展把SMO和梯度下降做对比看看在同样的问题上两种优化方法的收敛曲线和最终解有什么差异。或者把SMO应用到多分类场景用一对多或一对一策略观察成对更新的逻辑在多分类下怎么调整。