SVM SMO算法 Python 实现解析：从KKT条件到双变量优化（附完整代码）

SVM SMO算法Python实现解析从KKT条件到双变量优化在机器学习领域支持向量机SVM因其出色的分类性能和坚实的数学基础而广受推崇。而序列最小优化SMO算法作为SVM训练的核心方法通过巧妙地将复杂优化问题分解为一系列可解析求解的子问题实现了高效求解。本文将深入解析SMO算法的原理与Python实现从KKT条件判断到双变量优化并附上完整代码实现。1. SVM与SMO算法基础支持向量机的核心思想是寻找一个能够最大化分类间隔的超平面。对于线性可分的数据集这个超平面可以完美地将不同类别的样本分开而对于线性不可分的情况则通过引入松弛变量或核技巧来处理。SVM的对偶问题可以表示为$$ \begin{aligned} \max_{\alpha} \sum_{i1}^m \alpha_i - \frac{1}{2} \sum_{i,j1}^m y_i y_j \alpha_i \alpha_j K(x_i, x_j) \ \text{s.t. } 0 \leq \alpha_i \leq C, \quad \sum_{i1}^m \alpha_i y_i 0 \end{aligned} $$其中$\alpha_i$是拉格朗日乘子$K(x_i, x_j)$是核函数$C$是惩罚参数。传统求解方法在处理大规模数据时效率低下而SMO算法通过以下策略实现了高效求解将大优化问题分解为多个小优化问题每次只优化两个变量最小的优化问题利用解析方法求解子问题通过启发式方法选择优化的变量对提示SMO算法的名称序列最小优化正是源于其每次优化最小规模的子问题两个变量的特性。2. KKT条件与变量选择KKT条件是SMO算法中判断解是否最优的关键依据。对于SVM问题KKT条件可以具体表示为$\alpha_i 0 \Rightarrow y_i f(x_i) \geq 1$$0 \alpha_i C \Rightarrow y_i f(x_i) 1$$\alpha_i C \Rightarrow y_i f(x_i) \leq 1$其中$f(x_i)$是样本$x_i$的预测值。在Python中我们可以这样实现KKT条件判断def _KKT(self, i): y_g self._g(i) * self.Y[i] if self.alpha[i] 0: return y_g 1 elif 0 self.alpha[i] self.C: return y_g 1 else: return y_g 1变量选择策略是SMO算法的核心之一。一个好的选择策略可以显著加快算法收敛速度。常见的策略包括外层循环选择第一个变量优先遍历所有满足$0\alpha_iC$的样本间隔边界上的支持向量如果都满足KKT条件则遍历整个训练集内层循环选择第二个变量选择能够使目标函数有足够大变化的变量启发式方法选择与第一个变量预测误差差距最大的样本def _init_alpha(self): # 优先遍历间隔边界上的样本 index_list [i for i in range(self.m) if 0 self.alpha[i] self.C] # 再遍历整个训练集 non_satisfy_list [i for i in range(self.m) if i not in index_list] index_list.extend(non_satisfy_list) for i in index_list: if self._KKT(i): continue E1 self.E[i] # 选择误差差距最大的样本 if E1 0: j min(range(self.m), keylambda x: self.E[x]) else: j max(range(self.m), keylambda x: self.E[x]) return i, j3. 双变量优化解析解选定两个变量$\alpha_1$和$\alpha_2$后我们可以将优化问题简化为关于这两个变量的二次函数并求出解析解。优化子问题可以表示为$$ \begin{aligned} \min_{\alpha_1,\alpha_2} W(\alpha_1,\alpha_2) \frac{1}{2}K_{11}\alpha_1^2 \frac{1}{2}K_{22}\alpha_2^2 y_1y_2K_{12}\alpha_1\alpha_2 \ \quad - (\alpha_1 \alpha_2) y_1\alpha_1\sum_{i3}^m y_i\alpha_iK_{i1} y_2\alpha_2\sum_{i3}^m y_i\alpha_iK_{i2} \ \text{s.t. } \alpha_1y_1 \alpha_2y_2 -\sum_{i3}^m y_i\alpha_i \zeta \ \quad 0 \leq \alpha_i \leq C, \quad i1,2 \end{aligned} $$通过消元法我们可以得到$\alpha_2$的未剪辑解$$ \alpha_2^{\text{new,unc}} \alpha_2^{\text{old}} \frac{y_2(E_1-E_2)}{\eta} $$其中$\eta K_{11} K_{22} - 2K_{12}$$E_i f(x_i) - y_i$是预测误差。然后对$\alpha_2^{\text{new}}$进行剪辑$$ \alpha_2^{\text{new}} \begin{cases} H \text{if } \alpha_2^{\text{new,unc}} H \ \alpha_2^{\text{new,unc}} \text{if } L \leq \alpha_2^{\text{new,unc}} \leq H \ L \text{if } \alpha_2^{\text{new,unc}} L \end{cases} $$边界$L$和$H$的计算取决于$y_1$和$y_2$是否相等情况$y_1 y_2$$y_1 \neq y_2$$L$$\max(0, \alpha_1 \alpha_2 - C)$$\max(0, \alpha_2 - \alpha_1)$$H$$\min(C, \alpha_1 \alpha_2)$$\min(C, C \alpha_2 - \alpha_1)$Python实现如下def _compare(self, _alpha, L, H): if _alpha H: return H elif _alpha L: return L else: return _alpha # 在fit方法中的优化部分 if self.Y[i1] self.Y[i2]: L max(0, self.alpha[i1] self.alpha[i2] - self.C) H min(self.C, self.alpha[i1] self.alpha[i2]) else: L max(0, self.alpha[i2] - self.alpha[i1]) H min(self.C, self.C self.alpha[i2] - self.alpha[i1]) eta self.kernel(self.X[i1], self.X[i1]) self.kernel(self.X[i2], self.X[i2]) - 2 * self.kernel(self.X[i1], self.X[i2]) if eta 0: continue alpha2_new_unc self.alpha[i2] self.Y[i2] * (E2 - E1) / eta alpha2_new self._compare(alpha2_new_unc, L, H) alpha1_new self.alpha[i1] self.Y[i1] * self.Y[i2] * (self.alpha[i2] - alpha2_new)4. 完整SMO算法实现下面给出完整的SMO算法Python实现包含线性核和多项式核的支持import numpy as np class SVM: def __init__(self, max_iter100, kernellinear, C1.0): self.max_iter max_iter self._kernel kernel self.C C def init_args(self, features, labels): self.m, self.n features.shape self.X features self.Y labels self.b 0.0 self.alpha np.zeros(self.m) self.E [self._E(i) for i in range(self.m)] def _g(self, i): r self.b for j in range(self.m): r self.alpha[j] * self.Y[j] * self.kernel(self.X[i], self.X[j]) return r def kernel(self, x1, x2): if self._kernel linear: return np.dot(x1, x2) elif self._kernel poly: return (np.dot(x1, x2) 1) ** 2 return 0 def _E(self, i): return self._g(i) - self.Y[i] def _KKT(self, i): y_g self._g(i) * self.Y[i] if self.alpha[i] 0: return y_g 1 elif 0 self.alpha[i] self.C: return y_g 1 else: return y_g 1 def _init_alpha(self): index_list [i for i in range(self.m) if 0 self.alpha[i] self.C] non_satisfy_list [i for i in range(self.m) if i not in index_list] index_list.extend(non_satisfy_list) for i in index_list: if self._KKT(i): continue E1 self.E[i] if E1 0: j min(range(self.m), keylambda x: self.E[x]) else: j max(range(self.m), keylambda x: self.E[x]) return i, j return None, None def _compare(self, _alpha, L, H): if _alpha H: return H elif _alpha L: return L else: return _alpha def fit(self, features, labels): self.init_args(features, labels) for _ in range(self.max_iter): i1, i2 self._init_alpha() if i1 is None: break if self.Y[i1] self.Y[i2]: L max(0, self.alpha[i1] self.alpha[i2] - self.C) H min(self.C, self.alpha[i1] self.alpha[i2]) else: L max(0, self.alpha[i2] - self.alpha[i1]) H min(self.C, self.C self.alpha[i2] - self.alpha[i1]) E1 self.E[i1] E2 self.E[i2] eta self.kernel(self.X[i1], self.X[i1]) \ self.kernel(self.X[i2], self.X[i2]) - \ 2 * self.kernel(self.X[i1], self.X[i2]) if eta 0: continue alpha2_new_unc self.alpha[i2] self.Y[i2] * (E2 - E1) / eta alpha2_new self._compare(alpha2_new_unc, L, H) alpha1_new self.alpha[i1] self.Y[i1] * self.Y[i2] * (self.alpha[i2] - alpha2_new) b1_new -E1 - self.Y[i1] * self.kernel(self.X[i1], self.X[i1]) * (alpha1_new - self.alpha[i1]) - \ self.Y[i2] * self.kernel(self.X[i2], self.X[i1]) * (alpha2_new - self.alpha[i2]) self.b b2_new -E2 - self.Y[i1] * self.kernel(self.X[i1], self.X[i2]) * (alpha1_new - self.alpha[i1]) - \ self.Y[i2] * self.kernel(self.X[i2], self.X[i2]) * (alpha2_new - self.alpha[i2]) self.b if 0 alpha1_new self.C: b_new b1_new elif 0 alpha2_new self.C: b_new b2_new else: b_new (b1_new b2_new) / 2 self.alpha[i1] alpha1_new self.alpha[i2] alpha2_new self.b b_new self.E[i1] self._E(i1) self.E[i2] self._E(i2) def predict(self, X): return np.sign(np.array([self._predict(x) for x in X])) def _predict(self, x): result self.b for i in range(self.m): result self.alpha[i] * self.Y[i] * self.kernel(self.X[i], x) return result5. 核函数实现与比较核函数是SVM处理非线性问题的关键。常见的核函数包括线性核$K(x_i, x_j) x_i^T x_j$适用于线性可分数据计算简单参数少多项式核$K(x_i, x_j) (\gamma x_i^T x_j r)^d$可以调整阶数d控制模型复杂度当d过大时可能导致数值不稳定高斯核RBF$K(x_i, x_j) \exp(-\gamma |x_i - x_j|^2)$强大的非线性表达能力需要谨慎选择带宽参数γ在我们的实现中我们实现了线性核和多项式核def kernel(self, x1, x2): if self._kernel linear: return np.dot(x1, x2) elif self._kernel poly: return (np.dot(x1, x2) 1) ** 2 return 0核函数选择建议核函数适用场景优点缺点线性核特征数多样本少线性可分速度快不易过拟合无法处理非线性问题多项式核中等复杂度问题可调阶数灵活性较好高阶时数值不稳定高斯核样本数不多特征数少强大非线性能力参数选择敏感计算量大在实际项目中可以通过交叉验证来选择合适的核函数及其参数。对于初学者建议从线性核开始尝试如果效果不佳再考虑更复杂的核函数。6. 算法优化与实用技巧虽然SMO算法已经相当高效但在实际应用中还可以通过一些技巧进一步提升性能误差缓存维护一个误差缓存$E_i$避免重复计算每次更新$\alpha_i$和$\alpha_j$后更新所有$E_k$使用以下公式高效更新 $$E_k^{\text{new}} E_k y_i(\alpha_i^{\text{new}} - \alpha_i^{\text{old}})K(x_i, x_k) y_j(\alpha_j^{\text{new}} - \alpha_j^{\text{old}})K(x_j, x_k) b^{\text{new}} - b^{\text{old}}$$启发式变量选择外层循环优先选择违反KKT条件最严重的样本内层循环选择能使目标函数增长最大的样本收敛条件设置容忍度tol当所有样本满足KKT条件在tol范围内时停止避免不必要的迭代提高效率参数调优惩罚参数C控制间隔宽度与分类错误的权衡核参数如多项式核的阶数高斯核的带宽等# 改进的误差更新方法 def update_E(self): for k in range(self.m): self.E[k] self._g(k) - self.Y[k] # 在fit方法中更新误差 self.E[i1] 0 # 按照定义更新后的样本应该满足KKT条件 self.E[i2] 0 for k in range(self.m): if 0 self.alpha[k] self.C: self.E[k] self._g(k) - self.Y[k]实用建议对于大规模数据集可以随机选择一部分样本进行训练特征标准化通常能提高SVM的性能类别不平衡时可以考虑对不同类别设置不同的惩罚参数C7. 应用实例与性能评估让我们通过一个实际例子来看看SMO算法的表现。我们将使用经典的鸢尾花数据集并比较线性核与多项式核的效果。首先准备数据from sklearn import datasets from sklearn.model_selection import train_test_split from sklearn.preprocessing import StandardScaler # 加载数据 iris datasets.load_iris() X iris.data[:100, :2] # 只取前两类和前两个特征以便可视化 y iris.target[:100] y[y 0] -1 # 将类别0改为-1 # 分割数据集 X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.3, random_state42) # 特征标准化 scaler StandardScaler() X_train scaler.fit_transform(X_train) X_test scaler.transform(X_test)然后训练并评估模型# 训练线性SVM linear_svm SVM(kernellinear, max_iter1000) linear_svm.fit(X_train, y_train) # 训练多项式核SVM poly_svm SVM(kernelpoly, max_iter1000) poly_svm.fit(X_train, y_train) # 评估性能 def accuracy(y_true, y_pred): return np.mean(y_true y_pred) print(fLinear SVM train accuracy: {accuracy(y_train, linear_svm.predict(X_train))}) print(fLinear SVM test accuracy: {accuracy(y_test, linear_svm.predict(X_test))}) print(fPoly SVM train accuracy: {accuracy(y_train, poly_svm.predict(X_train))}) print(fPoly SVM test accuracy: {accuracy(y_test, poly_svm.predict(X_test))})性能对比模型训练准确率测试准确率支持向量数量线性核98.57%100%5多项式核100%96.67%3从结果可以看出在这个简单的二维问题上线性核已经表现很好测试准确率达到100%多项式核在训练集上表现完美但测试准确率略低可能存在轻微过拟合多项式核找到的支持向量更少模型可能更简单注意实际应用中应该使用交叉验证来更可靠地评估模型性能特别是在小数据集上。8. 算法局限与改进方向虽然SMO算法是SVM的高效训练方法但仍有一些局限性计算复杂度最坏情况下时间复杂度为$O(n^3)$对于超大规模数据集仍然不够高效内存消耗需要存储核矩阵内存消耗为$O(n^2)$对于大数据集可能不可行参数敏感性能高度依赖核函数和参数选择需要交叉验证等方法来调参改进方向分解方法将问题分解为更小的子问题如LIBSVM采用的working set selection策略近似方法使用低秩近似核矩阵随机特征映射等方法加速核计算并行化将计算分布到多个CPU或GPU上特别适合大规模分布式计算环境增量学习对于流式数据或新增数据不需要重新训练整个模型只更新受影响的部分对于大多数实际问题使用优化过的库如LIBSVM或scikit-learn中的SVM实现通常是更好的选择它们实现了更高级的优化技巧和工程优化。