重剑无锋：超大规模QP快速优化器以及其他solver

求解大规模QPP问题，几乎是是每一种基于核方法的机器学习算法最终要面对的问题，求解器的发展也有及其长的历史，本文记录笔者学习过程中遇到的几种代表性，重要的方法。

接着，本文将继续推导其他不限于QPP的经典solver的实现细节。包括内点法等经典方法，以及现代神经网络训练使用到的adam，Nadam等主流优化器。书山有路勤为径，学海无涯苦作舟，机器学习理论与方法正在不断快速迭代，知识也应当是不断更新的，笔者也会持续更新此文。

此外，优化算法是个无底洞，对数学要求很高，除非专门做优化算法本身的学者，不建议过多时间投入到此处。建议在学习完一些经典的solver后，直接使用现成的库即可。比如scikit-learn，libsvm，torch，cvxopt，osqp等。

SMO

1996年，John Platt 发布一个称为 SMO（Sequential Minimal Optimization，序列最小优化）的强大算法。（SVM）中，我们最终要解决的优化问题是一个二次规划问题（QP）：

$$ \max_{\boldsymbol\alpha} \quad W(\boldsymbol\alpha) = \sum_{i=1}^n \alpha_i - \frac{1}{2} \sum_{i,j} \alpha_i \alpha_j y_i y_j K(\mathbf{x}_i, \mathbf{x}_j) $$

$$ \text{s.t. } \sum_{i=1}^n \alpha_i y_i = 0,\quad 0 \le \alpha_i \le C $$

这是一个约束优化问题，变量是 α₁, …, α_n，它们的个数和样本数量一样。求解上述问题需要专用的二次规划器（如内点法、拉格朗日对偶算法）。当样本数 n 较大（比如上万），存储和运算开销会非常高。所以我们需要一种更高效、结构化的优化方法。

SMO 的核心理念是：

“每次只优化两个变量（两个 α），将高维约束问题分解成一系列二元子问题，这些子问题可以解析求解。”

为什么是两个？因为我们有一个等式约束：

$$ \sum_{i=1}^n \alpha_i y_i = 0 $$

优化两个变量时，可以直接通过代数法消去其中一个，从而将约束转化为一维问题。

第一步：固定其他变量，优化两个变量

目标函数（对偶）：

$$ W(\alpha) = \sum_{i=1}^n \alpha_i - \frac{1}{2} \sum_{i,j=1}^n \alpha_i \alpha_j y_i y_j K(x_i, x_j) $$

固定除 α_i, α_j 以外的变量后，这变成了一个关于 α_i, α_j 的二元二次函数。

由于 α_iy_i + α_jy_j = −∑_{k ≠ i, j}α_ky_k = 常数 故α_i, α_j更新前后，应当满足 α_i^oldy_i + α_j^oldy_j = α_iy_i + α_jy_j 令s = y_iy_j，可以将 α_i 表示为：

α_i = α_i^old + s(α_j^old − α_j)

第二步：构造二次目标函数

此时的目标函数事实上变成了关于α_i,α_j的二元二次函数。 $$ W(\alpha) = \text{常数} + \alpha_i + \alpha_j - \frac{1}{2} \left[ \alpha_i^2 K_{ii} + \alpha_j^2 K_{jj} + 2 \alpha_i \alpha_j y_i y_j K_{ij} \right] - \text{其他项} $$

代入：

α_i = α_i^old + s(α_j^old − α_j)

将上述整理，写为标准一维二次函数：

$$ W(a_j) = -\frac{1}{2} \eta a_j^2 + A a_j + \text{常数} $$

其中：

• η = K_ii + K_jj − 2K_ij
• A = y_j(E_i − E_j)
• *E_i = f(x_i) − y_i，E_j = f(x_j) − y_j：预测误差

$$ \frac{dW}{d\alpha_j} = y_j (E_i - E_j) - \eta \cdot (\alpha_j - \alpha_j^{\text{old}}) \Rightarrow $$

令导数为 0，得到解析更新公式：

$$ \alpha_j^{\text{new, unclip}} = \alpha_j^{\text{old}} + \frac{y_j (E_i - E_j)}{\eta} \tag{2} $$

第三步：边界裁剪（clip）

20250505163359

为满足盒约束，根据 y_i, y_j 的关系，推导出合法区间 [L, H]。：

• 若 y_i ≠ y_j：

  L = max (0, α_j^old − α_i^old),  H = min (C, C + α_j^old − α_i^old)

• 若 y_i = y_j：

  L = max (0, α_i^old + α_j^old − C),  H = min (C, α_i^old + α_j^old)

然后进行裁剪：

α_j^new = clip(α_j^{new, unclipped}, L, H)

第四步：更新另一个变量 α_i

用约束表达式反求另一个变量：

α_i^new = α_i^old + s(α_j^old − α_j^new)

第五步：更新阈值 b

b的更新是一种启发式方法，不做过多赘述，有两种情况下可选的更新策略：

b₁ = b^old − E_i − y_i(α_i^new − α_i^old)K_ii − y_j(α_j^new − α_j^old)K_ij

b₂ = b^old − E_j − y_i(α_i^new − α_i^old)K_ij − y_j(α_j^new − α_j^old)K_jj

最终的更新方式：

• 如果 0 < α_i < C：设 b = b₁
• 否则如果 0 < α_j < C：设 b = b₂
• 否则取均值 b = (b₁ + b₂)/2

SOR（逐次超松弛）

因为解大规模QPP等同于解大规模线性方程组，因此解方程组的方法也可以用于解QPP。

我们从最基本的开始讲解，循序渐进地理解从雅可比（Jacobi）迭代法发展到 SOR（Successive Over-Relaxation，逐次超松弛法）的来龙去脉。

一、目标：求解线性方程组

我们的问题是： 给定一个线性方程组：

Ax = b

其中：

• A ∈ ℝ^n × n：系数矩阵
• x ∈ ℝⁿ：待求解向量
• b ∈ ℝⁿ：常数向量

当 A 较大或稀疏时，直接解法（如高斯消元、LU 分解）代价太高，所以我们使用迭代法。

二、雅可比迭代法（Jacobi Method）

在数值线性代数中，我们希望解线性方程组：

Ax = b

当矩阵 A 维度较大或稀疏时，直接求逆代价过高，因此考虑构造一个迭代过程，令解 x 逐步逼近。

这里我们将矩阵 A 分解为三部分：

A = D + L + U

• D：对角线上的元素（对角矩阵）
• L：下三角部分（不含对角）
• U：上三角部分（不含对角）

将 Ax = b 改写为：

Dx = b − (L + U)x

两边左乘 D⁻¹，得到不动点形式：

x = D⁻¹(b − (L + U)x)

自然导出 Jacobi 迭代公式：

$$ \boxed{ \mathbf{x}^{(k+1)} = D^{-1} \left( \mathbf{b} - (L + U)\mathbf{x}^{(k)} \right) } $$

同理可以使用分量形式表达

对于方程组第 i 行：

a_i1x₁ + ⋯ + a_iix_i + ⋯ + a_inx_n = b_i

解出 x_i：

$$ x_i = \frac{1}{a_{ii}} \left( b_i - \sum_{j \ne i} a_{ij}x_j \right) $$

引入迭代思想，即将未知分量 x_j 全部替换为上一步的近似值 x_j^(k)，得到：

$$ x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j \ne i} a_{ij}x_j^{(k)} \right) $$

Jacobi 的主要问题是收敛慢。在一次迭代中，更新 x_i^(k + 1) 时，不利用已经更新的 x₁^(k + 1), …, x_i − 1^(k + 1)，导致信息利用不充分。

其收敛性取决于迭代矩阵 M = D⁻¹(L + U) 的谱半径 ρ(M)：

• 若 ρ(M) < 1，则收敛
• 若 ρ(M) > 1，则发散

例如：

$$ A = \begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix} $$

Jacobi 的迭代矩阵谱半径约为 1.08 > 1，因此不收敛。

✅ 小结

Jacobi 方法从线性方程组中构造出一个显式的不动点形式，每一步使用上次的全部值进行更新，适合并行计算。但由于更新中未充分利用最新值，收敛性较差，通常作为后续 Gauss-Seidel 或 SOR 方法的引子。

高斯-赛德尔（Gauss-Seidel）迭代法

Jacobi 有个问题：每次迭代都等到所有 x_j^(k) 用完才算完，不充分利用已更新的结果。

Gauss-Seidel 改进：一边算一边更新。

更新公式：

$$ x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(k)} \right) $$

也就是：前面的 x_j 用新值，后面的用旧值。

优点：通常比 Jacobi 收敛快 缺点：仍可能收敛慢，特别是矩阵条件数不好时

逐次超松弛法（SOR）

Gauss-Seidel 其实是 Jacobi 的“小步走”，SOR 就是加速这个过程的“大步走”。

在 Gauss-Seidel 的基础上，SOR 加入一个松弛因子 ω ∈ (0, 2)，做“加权平均”：

$$ x_i^{(k+1)} = (1 - \omega) x_i^{(k)} + \frac{\omega}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(k)} \right) $$

• 若 ω = 1，就退化为 Gauss-Seidel；
• 若 ω ∈ (1, 2)，为超松弛（加速）；
• 若 ω ∈ (0, 1)，为欠松弛（更稳定）。

通过松弛因子ω引入惯性项，动态调整步长：

x_i^(k + 1) = (1 − ω)x_i^(k) + ω(Gauss − Seidel更新量)

ω > 1:放大修正步长(超松弛),加速逼近解。ω < 1:缩减不稳定更新(低松弛),抑制震荡。

DCDM对偶坐标下降方法

源自ICML一篇论文(https://www.csie.ntu.edu.tw/~cjlin/papers/cddual.pdf)[https://www.csie.ntu.edu.tw/~cjlin/papers/cddual.pdf]

对偶坐标下降方法（DCDM）的数学推导全流程

原始SVM问题： $$ \min_w \frac{1}{2} w^T w + C \sum_{i=1}^l \xi(w; x_i, y_i) $$ 其中损失函数ξ为L1或L2损失： - L1-SVM：ξ(w; x_i, y_i) = max (1 − y_iw^Tx_i, 0)

• L2-SVM：ξ(w; x_i, y_i) = max (1 − y_iw^Tx_i, 0)²

对偶问题推导：引入拉格朗日乘子α_i，构造对偶函数。通过KKT条件，原始问题的对偶形式为

$$ \min_\alpha \frac{1}{2} \alpha^T \bar{Q} \alpha - e^T \alpha \quad \text{s.t.} \quad 0 \leq \alpha_i \leq U $$ 其中： - Q̄ = Q + D，Q_ij = y_iy_jx_i^Tx_j

• L1-SVM：D_ii = 0，U = C

• L2-SVM：D_ii = 1/(2C)，U = ∞

为了方便叙述，先定义这样一组向量 α^k, i = [α₁^k + 1, …, α_i − 1^k + 1, α_i^k, …, α_l^k]^T 此时的α表示正在进行第K+1轮迭代，按顺序已经迭代完了前i-1个分量，正准备迭代其中第i个分量。

DCDM的子问题定义： 固定除α_i外的所有变量，优化单变量α_i： min_df(α + de_i)  s.t.  0 ≤ α_i + d ≤ U 其中$f(\alpha) = \frac{1}{2} \alpha^T \bar{Q} \alpha - e^T \alpha$。

目标函数泰勒展开后(二阶近似)是一个关于d的二次函数： $$ f(\alpha + d e_i) = \frac{1}{2} \bar{Q}_{ii} d^2 + \nabla_i f(\alpha) d + \text{常数项} $$ 其中梯度分量：∇_if(α) = (Q̄α)_i − 1 = y_iw^Tx_i − 1 + D_iiα_i，w = ∑_jy_jα_jx_j。

闭式解推导： 1. 二次函数极小值：最优d为无约束解： $$ d^* = -\frac{\nabla_i f(\alpha)}{\bar{Q}_{ii}} $$ 2. 投影到可行域：

$$\alpha_i^{k,i+1}=\min\left(\max\left(\alpha_i^{k,i}-\frac{\nabla_if(\boldsymbol{\alpha}^{k,i})}{\bar{Q}_{ii}},0\right),U\right)$$

(即使用clip函数限制在（0，U）区间上)

此时问题变成了如何高效求解∇_if(α)，自然地想法是直接求解$\nabla_if(\boldsymbol{\alpha})=(\bar{Q}\boldsymbol{\alpha})_i-1=\sum_{j=1}^l\bar{Q}_{ij}\alpha_j-1$,其中Q̄_ij = x_i^⊤x_j,时间复杂度是O(ln̄),n̄是平均非零特征数，这样的计算开销是巨大的。（可能会疑惑为何不预先计算出整个矩阵Q？事实上是因为存储成本太高！若样本数l = 10⁵, 那么Q存储10¹⁰个浮点数，粗略记每个单精度浮点数4字节，此时内存开销已经来到将近40GB，完全不可行）

因此我们计算分量式子 ∇_if(α) = y_iw^Tx_i − 1 + D_iiα_i 当然，如果按照定义式子$\boldsymbol{w}=\sum_{j=1}^ly_j\alpha_j\boldsymbol{x}_j$求解w,复杂度就回到了O(ln̄),为此作者也对w在每一轮进行迭代 w ← w + (α_i − ᾱ_i)y_ix_i

20250505215415

根据当前 α_i 的位置，PG的计算分为以下三种情况：

αi 的位置	PG的计算公式	物理意义
0 < αi < U（内部）	PG = G	变量未触及边界，直接使用原始梯度 G 的方向更新。
αi = 0（下边界）	PG = min (G, 0)	仅允许梯度方向为负（G < 0）时更新，避免 αi < 0。
αi = U（上边界）	PG = max (G, 0)	仅允许梯度方向为正（G > 0）时更新，避免 αi > U。

OSQP（https://osqp.org/）

原始论文地址(https://web.stanford.edu/~boyd/papers/pdf/osqp.pdf)

介绍OSQP之前，需要先介绍ADMM（Alternating Direction Method of Multipliers）方法。ADMM是一个非常经典的优化算法，OSQP就是基于ADMM的一个变种。

ADMM

Boyd等人的推导：（https://stanford.edu/~boyd/papers/pdf/admm_talk.pdf）

ADMM 将优化问题分解为两个交替更新的子问题： $$ \begin{aligned} \min_{x,z} \quad & f(x) + g(z) \\ \text{s.t.} \quad & A x + B z = c \end{aligned} $$ 写出增广拉格朗日函数 $$\begin{aligned}L_\rho(\mathbf{x},\mathbf{z},\boldsymbol{\lambda})&=Q_\rho(\mathbf{x},\mathbf{z})+\boldsymbol{\lambda}^\top(\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{z}-\mathbf{c})\\&=f(\mathbf{x})+g(\mathbf{z})+\frac{\rho}{2}\|\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{z}-\mathbf{c}\|_2^2+\boldsymbol{\lambda}^\top(\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{z}-\mathbf{c})\end{aligned}$$ 迭代公式： $$ \begin{aligned} x^{k+1} &= \arg\min_x \left( f(x) + \frac{\rho}{2} \| A x + B z^k - c + \lambda^k \|^2 \right) \\ z^{k+1} &= \arg\min_z \left( g(z) + \frac{\rho}{2} \| A x^{k+1} + B z - c + \lambda^k \|^2 \right) \\ \lambda^{k+1} &= \lambda^k + \rho (A x^{k+1} + B z^{k+1} - c) \end{aligned} $$

OSQP演变

从ADMM到OSQP的数学推导

考虑凸二次规划问题： $$ \begin{aligned} \min_{x} \quad & \frac{1}{2}x^T P x + q^T x \\ \text{s.t.} \quad & l \leq A x \leq u \end{aligned} $$ 其中 P ≽ 0，A ∈ ℝ^m × n。引入辅助变量 z，将原问题改写为： $$ \begin{aligned} \min_{x,z} \quad & \frac{1}{2}x^T P x + q^T x + I_{[l, u]}(z) \\ \text{s.t.} \quad & A x = z \end{aligned} $$ 其中 I_[l, u](z) 是指示函数： $$ I_{[l, u]}(z) = \begin{cases} 0 & \text{if } l \leq z \leq u \\ +\infty & \text{otherwise} \end{cases} $$

构造增广拉格朗日函数： $$ \mathcal{L}_\rho(x, z, y) = \frac{1}{2}x^T P x + q^T x + I_{[l, u]}(z) + y^T (A x - z) + \frac{\rho}{2} \|A x - z\|^2 $$ 其中 y ∈ ℝ^m 为对偶变量，ρ > 0 为惩罚参数。采用ADMM迭代步骤交替更新 x、z、y：

(1) x-子问题 固定 z^k, y^k，求解： $$ x^{k+1} = \arg\min_x \left( \frac{1}{2}x^T P x + q^T x + y^{kT} (A x - z^k) + \frac{\rho}{2} \|A x - z^k\|^2 \right) $$ 展开并整理目标函数： $$ x^{k+1} = \arg\min_x \frac{1}{2}x^T (P + \rho A^T A) x + \left( q + A^T ( \rho z^k - y^k ) \right)^T x $$ 其解析解为： x^k + 1 = −(P + ρA^TA)⁻¹(q + A^T(ρz^k − y^k))

(2) z-子问题 固定 x^k + 1, y^k，求解： $$ z^{k+1} = \arg\min_z \left( I_{[l, u]}(z) + \frac{\rho}{2} \| z - (A x^{k+1} + \frac{y^k}{\rho}) \|^2 \right) $$ 等价于投影到区间 ([l, u])： $$ z^{k+1} = \Pi_{[l, u]} \left( A x^{k+1} + \frac{y^k}{\rho} \right) $$

(3) 对偶变量更新 y^k + 1 = y^k + ρ(Ax^k + 1 − z^k + 1)

完整算法流程 - 输入：初始 x⁰, z⁰, y⁰，参数 ρ > 0，容差 ϵ。
迭代直到收敛： - 1. 更新 x^k + 1 = −(P + ρA^TA)⁻¹(q + A^T(ρz^k − y^k)) - 2. 更新 $z^{k+1} = \Pi_{[l, u]} \left( A x^{k+1} + \frac{y^k}{\rho} \right)$ - 3. 更新 y^k + 1 = y^k + ρ(Ax^k + 1 − z^k + 1) - 4. 检查终止条件 ∥r^k∥ ≤ ϵ, ∥s^k∥ ≤ ϵ