Matlab实现DCDM训练L2-SVM：对偶坐标下降法原理与实践

前阵子要在 Matlab 里实现一个大间距线性分类器的训练模块，数据是高维稀疏的文本特征，样本量大概几十万，特征维度有几百万。我一开始图省事直接写全批量梯度下降，结果被两个问题折腾得够呛：一是步长真是玄学，给大了直接发散，给小了每一轮都在原地磨蹭；二是算完一次全量梯度只更新一次权重向量，收敛慢到让人怀疑人生。后来把目标函数写成了对偶问题，改用 DCDM（Dual Coordinate Descent Method，对偶坐标下降法）来训练，整个求解过程突然就清爽了：不需要调学习率，不需要做线搜索，每一轮只碰一个对偶变量，但权重向量会跟着同步更新，代码写起来也特别直白。

这篇文章就把我这套实现从头到尾拆开讲一遍，包括损失函数怎么选、拉格朗日对偶怎么推、单变量更新公式怎么来的、Matlab 代码怎么写，以及我实际调参和排坑时踩过的几个坑。适合正在做线性模型但不想只调第三方库 API 的开发者，也适合想彻底搞懂“对偶坐标下降为什么快”的初学者。

1. DCDM是什么：一个适合大规模线性模型的优化器

1.1 从坐标下降到对偶坐标下降

坐标下降法的思路很朴素：每次只优化一个变量，其他变量保持不动，一个变量一个变量地扫过去。原始问题里的权重向量w在特征维度很高时，直接对w做坐标下降是很痛苦的，因为样本的损失函数是全局耦合的，你动一个w_j，所有样本的预测值都会变，目标函数不是一个简单的单变量二次函数。

但如果先做对偶变换，把优化对象换成拉格朗日乘子α，情况就完全不同。在不少线性模型里，对偶问题的每个单变量子问题都长得很规矩：一个一元二次函数加上一个非负约束。这个子问题可以直接用闭式解算出来，不需要任何步长参数。这就是 DCDM 的核心思想：把“对w做坐标下降”变成“对α做坐标下降”。

这个转换的好处还在于，对偶变量的个数等于样本数n，而不是特征维度d。对文本分类这种d上千万但n只有几万的场景，对偶空间的计算代价会舒服很多。

1.2 它能解决什么问题

DCDM 最适合的场景是“线性模型 + 中等规模以上的样本量 + 稀疏特征”。文本分类、点击率预估、基因表达数据建模，都属于这一类。它的单轮复杂度跟算一次全量梯度差不多，但一轮迭代内部会更新很多次对偶变量，所以实际收敛速度通常比朴素梯度法快很多。

另外一个很实际的价值是：DCDM 是确定性优化算法，不涉及随机采样，也不需要维护优化器的动量或者自适应学习率状态。这对工程实现非常友好。如果你只是在 Matlab 里做一个原型验证，或者需要把某个定制损失函数改成自己的变体，手写这样一套求解器往往比继续跟黑盒库搏斗更可控。

提示：这里的 DCDM 是一族方法的总称。不同论文里可能会用来处理 SVM、逻辑回归或者各种正则化线性模型。本文的样例实现以 L2-SVM（二次合页损失）为例，因为它对偶形式最简洁，闭式更新公式也最容易理解。

2. 原理拆解：从原始问题到对偶坐标下降

2.1 先定一个具体的原始问题：L2-SVM

假设我们有n个样本x_i，标签y_i ∈ {+1, -1}，特征维度是d。我用的是 L2-SVM 原始问题：

min_w 0.5 * ||w||^2 + (C / 2) * Σ_i max(0, 1 - y_i * x_i^T * w)^2

这里的C是正则化系数。第一项限制模型复杂度，第二项是二次合页损失。和标准合页损失max(0, 1 - y_i * x_i^T * w)相比，平方项让损失函数光滑，对偶问题的约束也会简单一些。

需要注意一个细节：我这里写的是(C / 2) * loss^2，而不是很多人习惯写的C * loss^2。这个系数最后会影响对偶问题里对角项的系数。如果你用的是别的版本，推导结果里会出现1 / (2C)或者1 / C的差别，别对着公式硬套，最好自己从头推一遍。

2.2 拉格朗日对偶：把主变量消掉

引入松弛变量ξ_i，原始问题可以改写成约束形式：

min_{w,ξ} 0.5 * ||w||^2 + (C / 2) * Σ_i ξ_i^2

s.t. y_i * x_i^T * w ≥ 1 - ξ_i

拉格朗日函数写出来：

L = 0.5 * ||w||^2 + (C / 2) * Σ_i ξ_i^2 - Σ_i α_i * (y_i * x_i^T * w - 1 + ξ_i)

其中α_i ≥ 0。对w和ξ_i求导并让导数为零：

w = Σ_i α_i * y_i * x_i

ξ_i = α_i / C

把这两个关系代回拉格朗日函数，得到对偶问题的最大化形式。为了符合“最小化”的习惯，取个负号，就得到：

min_α F(α) = 0.5 * α^T * Q * α + 0.5 / C * ||α||^2 - Σ_i α_i

Q_ij = y_i * y_j * x_i^T * x_j

约束只有α_i ≥ 0，没有上界。这一点和标准合页 SVM 不同，二次合页损失的对偶变量没有0 ≤ α_i ≤ C这个框，公式写起来会更省事。

2.3 单变量子问题的闭式解

固定其他α_j，只考虑α_i，目标函数在α_i上就是一个一元二次函数。因为这个坐标子问题里含有二次项，极值点是解析可求的。

定义当前样本对偶预测值：

z_i = y_i * x_i^T * w

因为w = Σ_j α_j * y_j * x_j，所以z_i实际上就等于(Q * α)_i。

固定其他变量时，关于α_i的目标可以写成：

0.5 * (Q_ii + 1 / C) * α_i^2 + (Σ_{j≠i} Q_ij * α_j - 1) * α_i + const

令A_i = Q_ii + 1 / C = ||x_i||^2 + 1 / C，那么新的α_i就是：

α_i_new = max(0, α_i_old - g_i / A_i)

其中g_i = z_i - 1 + α_i_old / C

这个形式特别像“梯度下降”，但分母A_i是二阶信息，而不是人为设定的学习率。max(0, ·)只是把结果投影到非负约束上。由于A_i > 0，这个更新保证目标函数不上升。

求出α_i_new之后，对照最原始的关系式w = Σ_i α_i * y_i * x_i，权重向量只需要做增量更新：

w ← w + (α_i_new - α_i_old) * y_i * x_i

这个增量更新非常关键。它意味着我们不需要每一轮都重新累加所有支持向量，只要碰过哪一个样本，就顺手把它的贡献加进w。整个算法的代价因此维持在O(nnz)量级，而不是O(n^2)。

3. Matlab实现与关键代码

3.1 数据准备与预计算

写代码前先把数据规范好。标签必须转成+1 / -1，如果原始标签是0 / 1，记得做一下y = 2 * y - 1。

特征矩阵X建议用稀疏矩阵存储。高维文本特征里绝大多数位置是零，稀疏存储能节省大量内存，迭代速度也快很多。

function [w, alpha, info] = dcdm_l2svm(X, y, C, opts) % DCDM_L2SVM -- 对偶坐标下降法训练 L2-SVM % X: n x d 矩阵，支持稀疏矩阵 % y: n x 1 标签向量，只能取 +1 或 -1 % C: 正则化系数，必须大于 0 % opts.maxiter: 最大迭代轮数 % opts.tol: 对偶间隙停止阈值 % opts.verbose: 是否打印每轮日志 if nargin < 4 opts = struct(); end if ~isfield(opts, 'maxiter'), opts.maxiter = 200; end if ~isfield(opts, 'tol'), opts.tol = 1e-6; end if ~isfield(opts, 'verbose'), opts.verbose = false; end [n, d] = size(X); alpha = zeros(n, 1); w = zeros(d, 1); xnorm = full(sum(X.^2, 2)); Acoef = xnorm + 1 / C; history = [];

xnorm就是每个样本的||x_i||^2，它是坐标子问题分母里的主要部分。一次算好能省掉循环里重复计算的开销。

3.2 核心循环代码

进入主循环之后，我每一轮都会先把样本顺序随机打乱，再逐个更新对偶变量。随机轮转比固定顺序更不容易出现周期性的停滞，这是我在高相关特征上对比过的结果。

for iter = 1:opts.maxiter max_delta = 0; perm = randperm(n); for r = 1:n i = perm(r); yi = y(i); [idx, ~, vals] = find(X(i, :)); if isempty(idx) z = 0; else z = yi * (w(idx)' * vals); end g = z - 1 + alpha(i) / C; alpha_new = max(0, alpha(i) - g / Acoef(i)); delta = alpha_new - alpha(i); if abs(delta) > 0 alpha(i) = alpha_new; w(idx) = w(idx) + delta * yi * vals; max_delta = max(max_delta, abs(delta)); end end

这里我没有直接用X(i, :) * w，而是先find出行里的非零元素下标，再对w(idx)做点乘。原因是 Matlab 对稀疏行的处理有时会隐式生成临时数组，在n很大时这个临时开销会被放大很多倍。用find + vals的方式可以精确控制访问范围。

w(idx) = w(idx) + delta * yi * vals这一步相当于把w = w + delta * yi * x_i只作用在非零维度上，不会碰那些本来就不需要变化的零位置。

3.3 停机准则：用对偶间隙判断收敛

训练循环的停止条件不能只看max_delta，因为坐标下降有时候会出现一种情况：大部分α_i已经稳定，但一小部分还在缓慢变化。更可靠的指标是对偶间隙。

原始目标J(w)和对偶目标F(α)之间有一个天然的间距，最优解处两者相等。因为我的每一次坐标更新都保证对偶目标单调下降，原始目标在任何一组可行(w, α)下都大于等于对偶目标，所以gap = primal - dual是非负的，越接近零说明解越接近最优。

Xw = X * w; primal = 0.5 * (w' * w) + 0.5 * C * sum(max(0, 1 - y .* Xw).^2); dual = 0.5 * sum(alpha .* (y .* Xw)) + 0.5 / C * sum(alpha.^2) - sum(alpha); gap = primal - dual; if gap < 0 && opts.verbose fprintf('warning: negative gap %.3e\n', gap); end gap = max(gap, 0); info.gap(iter) = gap; info.primal(iter) = primal; if opts.verbose fprintf('iter %3d primal %.6f dual %.6f gap %.3e max_delta %.3e\n', ... iter, primal, dual, gap, max_delta); end if iter > 2 && gap < opts.tol break; end end end

计算dual时用到的y .* Xw就是Q * α的等价形式。这是对偶问题里最贵的部分，但它本身也只要一次稀疏矩阵乘法，和求一次原始损失差不多。相比起“看梯度模长是不是小于阈值”，对偶间隙有一个绝对尺度的解释：小于1e-6说明当前原始损失离最优损失的距离已经被压到这个量级以下。

注意：由于浮点误差，偶尔会出现gap为极小负值的情况。这只是数值扰动，不代表理论上出了问题，我在代码里直接max(gap, 0)处理掉了，但日志里仍然会打印 warning，方便排查。

3.4 完整调用示例

调用方式很简单，核心参数就一个C。

rng(42); n = 2000; d = 50; X = randn(n, d); w_true = zeros(d, 1); w_true(randperm(d, 10)) = randn(10, 1); y = sign(X * w_true + 0.1 * randn(n, 1)); y(y == 0) = 1; opts = struct('maxiter', 100, 'tol', 1e-8, 'verbose', true); [w, alpha, info] = dcdm_l2svm(X, y, 1, opts);

训练结束后，w可以直接用来做预测。alpha里非零元素对应的样本就是支持向量，只不过 L2-SVM 的支持向量概念不像标准软间隔 SVM 那么直观，但用来分析哪些样本在边界附近仍然有价值。

4. 实验效果与调参心得

4.1 在合成数据上的收敛表现

上面那段代码生成的合成数据里，真实权重只有 10 个非零维度，其余都是噪声。用C = 1跑下来，我的日志里对偶间隙变化大致是这样：

这个收敛速度在稠密数据上已经很快了。换成稀疏高维文本特征时，每轮时间主要花在稀疏矩阵乘法和行扫描上，但轮数基本还在同一个数量级。

这里想强调一个经验：不要一开始就把tol设成1e-10。对大多数应用，1e-5就已经能保证模型性能差别可以忽略。对偶间隙压得太低，后面的迭代基本是在给浮点噪声买单，收益很小。

4.2 初始化与C值的影响

代码里默认从α = 0、w = 0开始。这是最省心的初始化，因为对偶问题在这个位置一定是可行的，不会出现需要一个不可行起始点的情况。

如果你已经有了一组旧权重，比如在上一批数据上训练过，或者想扫描多个C值，可以用旧解做热启动：直接把旧α和旧w填进去继续跑。DCDM 对热启动天然友好，因为它本身就是在“接着上一步”微调每个变量。我在做C值网格搜索时，通常按从大到小或者从小到大的顺序逐个训练，每个新解都比从零开始快一半以上。

C值影响的是正则化和损失的权衡。C越大，模型越看重训练样本的拟合，α的更新步长分母1 / C越小，每次变化更剧烈，收敛激进但容易过拟合。C越小，正则化越强，坐标子问题越“钝”，更新更保守，训练更稳。实际使用时我会先取C = 1跑通流程，再根据验证集调几个数量级的网格，比如0.01, 0.1, 1, 10。

4.3 和梯度类方法的对比

我用同一个合成数据集对比过全批量梯度下降和 DCDM。梯度下降的关键问题是步长。步长太大，损失容易震荡；步长太小，前几十轮几乎看不出进展。即便我手工给梯度下降配了一个还算合适的步长，DCDM 仍然更快到达同等训练损失。

这个对比并不严谨，毕竟梯度下降的步长和初始化都影响结果。但方向性结论是明确的：在二次合页损失这种光滑凸问题上，DCDM 的单变量闭式更新等价于自带二阶缩放信息的精准坐标步，而梯度下降每一步只用了局部一阶信息，步长全靠拍脑袋，差距在早期迭代尤其明显。

5. 实战中遇到的坑与排查

5.1 特征尺度不一致的坑

有一次我把某个实验数据直接扔进训练函数，发现训练损失下降得特别慢，后面一查，某个特征的范围是零到一千万，其他特征是零到一。这个特征让xnorm里出现很大的离群值，导致对应样本的坐标步分母特别大，一次更新只挪动一点点。

问题的本质是||x_i||^2在不同样本之间差异太悬殊。解决办法也很常规：训练前先对特征做标准化，让每个特征有差不多的尺度。对稀疏矩阵不要做全局z-score，这会破坏稀疏性，更好的做法是只对非零位置做缩放，或者直接把每个特征缩放到统一区间。特征尺度一致之后，DCDM 的收敛速度会稳定很多。

5.2 稀疏矩阵行访问的坑

Matlab 里循环访问稀疏矩阵的每一行，写法不同性能差异非常大。我最早用的是xi = X(i, :); z = yi * (xi * w);，在二三十万样本的数据上，一轮循环要跑十几秒。改成find(X(i, :))拿到非零下标和值，再对w(idx)做点乘之后，同一条循环只有原来的大约四分之一时间。

还有一个隐藏问题：在更新w时，如果直接写w = w + delta * yi * full(X(i, :))'，每循环一个样本就会生成一个d维的临时稠密向量。d是百万级时，这个操作会引发大量内存分配。用w(idx) = w(idx) + delta * yi * vals就不会有这个问题。

5.3 数值误差导致的负间隙

我在长期训练末尾观察到过对偶间隙出现-1e-12这种值。一开始我还以为是推导出错，后来发现单纯是浮点累加误差。primal和dual在收敛时都是大数相减，精度稍微波动，间隙就可能变成负的。

处理办法有两个：一个是设置gap = max(gap, 0)，不把它当作严重问题；另一个是在打印日志时保留更多有效数字，观察负值的绝对量级。如果负到1e-8以上，那就不是浮点误差，而是实现里有 bug 了。

5.4 常见问题速查表

6. 后续扩展方向

6.1 把L2-SVM替换成逻辑回归

如果你的业务更习惯用逻辑回归而不是 SVM，DCDM 的大框架完全不变，需要替换的只是单变量子问题。逻辑回归的对偶问题里会出现log(1 + exp(·))项，坐标子问题不再是一元二次函数，求闭式解需要解带对数项的标量方程。常见做法是在该坐标方向上做几次一维牛顿迭代，把步长限制在可行域内，剩下的结构（更新w、计算对偶间隙、随机轮转）全部照搬。

这种替换最大的好处是工程代码复用度很高。我后来扩展逻辑回归版本时，核心循环没改，只改了单变量求解函数和损失计算函数，大约两三个小时就完成了原型验证。

6.2 偏置项和分块更新

标准 SVM 的偏置项b不应该被正则化。如果直接给X加一列常数特征，然后用这套代码训练，等于让b也参与了 L2 正则化，模型严格来说变成了带偏置的岭分类器，而不是教科书式的 SVM。如果你对偏置项有严格要求，需要在原始问题里单独处理b，这样对偶问题会多一个Σ_i α_i * y_i = 0的约束，单变量更新就要成对地做，这和 SMO 算法的思路就比较接近了。

工程上另一个简单的折中方案是把常数特征那一列单独留出来，训练后再根据验证集手动调一个b。虽然不够优雅，但在多数任务里完全够用。

6.3 面向更大规模数据的收敛技巧

当n继续增大，比如超过一百万行，逐样本扫描一轮都会变得很慢。此时可以用“收缩”策略：训练过程中有一部分α_i会长期停留在零附近，这些样本对应的坐标子问题实际上可以跳过。做法是每跑几轮就用 KKT 条件筛选一遍候选集，只对可能大于零的变量继续做坐标更新，每隔一定轮数再重新做一次全量检查，防止误伤。

我实际使用过这个技巧，在百万级稀疏数据上能把单轮时间压缩到原来的三分之一左右，而且没有引入明显的精度损失。需要留意的是收缩不能做得太激进，否则容易把一些缓慢变化的变量提前排除，导致对偶间隙停在阈值以上。

另外，DCDM 的随机轮转顺序在单线程上已经足够好，如果想进一步压榨性能，可以考虑按块并行：把样本分块，每个工作线程负责一块的坐标更新，块与块之间用共享w同步。不过这会引入并发安全的问题，复杂度会上升很多，我建议先跑通单线程版本再考虑。

这一版 DCDM 实现后来被我沿用到好几个内部数据集的实验里，最大的收获其实不是性能，而是调试体验：对偶间隙一个数字就能告诉我当前解离最优还差多远，比盯着梯度模长靠谱得多。如果你也在写自己的线性模型训练器，我强烈建议把对偶间隙这个指标留在代码里，它会在你判断“模型到底收没收敛”的时候帮你省下大量时间。