QP快速优化器以及其他solver

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

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

此外，优化算法是个无底洞，对数学要求很高，除非专门做优化算法本身的学者，不建议过多时间投入到此处。建议在学习完一些经典的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 \]

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

SMO 的核心理念是：

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

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

\[\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) \]

固定除 \(\alpha_i, \alpha_j\) 以外的变量后，这变成了一个关于 \(\alpha_i, \alpha_j\) 的二元二次函数。

由于

\[\alpha_i y_i + \alpha_j y_j = -\sum_{k \ne i,j} \alpha_k y_k = \text{常数} \]

故\(\alpha_i,\alpha_j\)更新前后，应当满足

\[\alpha_i^{old}y_i+\alpha_j^{old}y_j=\alpha_i y_i + \alpha_j y_j \]

令\(s = y_i y_j\)，可以将 \(\alpha_i\) 表示为：

\[\alpha_i = \alpha_i^{\text{old}} + s (\alpha_j^{\text{old}} - \alpha_j) \]

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

此时的目标函数事实上变成了关于\(\alpha_i\),\(\alpha_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{其他项} \]

代入：

\[\alpha_i = \alpha_i^{\text{old}} + s (\alpha_j^{\text{old}} - \alpha_j) \]

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

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

其中：

• \(\eta = 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）

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

• 若 \(y_i \ne y_j\)：

  \[L = \max(0, \alpha_j^{\text{old}} - \alpha_i^{\text{old}}),\quad H = \min(C, C + \alpha_j^{\text{old}} - \alpha_i^{\text{old}}) \]

• 若 \(y_i = y_j\)：

  \[L = \max(0, \alpha_i^{\text{old}} + \alpha_j^{\text{old}} - C),\quad H = \min(C, \alpha_i^{\text{old}} + \alpha_j^{\text{old}}) \]

然后进行裁剪：

\[\alpha_j^{\text{new}} = \text{clip}(\alpha_j^{\text{new, unclipped}}, L, H) \]

第四步：更新另一个变量 \(\alpha_i\)

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

\[\alpha_i^{\text{new}} = \alpha_i^{\text{old}} + s(\alpha_j^{\text{old}} - \alpha_j^{\text{new}}) \]

第五步：更新阈值 \(b\)

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

\[b_1 = b^{\text{old}} - E_i - y_i ( \alpha_i^{\text{new}} - \alpha_i^{\text{old}} ) K_{ii} - y_j ( \alpha_j^{\text{new}} - \alpha_j^{\text{old}} ) K_{ij} \]

\[b_2 = b^{\text{old}} - E_j - y_i ( \alpha_i^{\text{new}} - \alpha_i^{\text{old}} ) K_{ij} - y_j ( \alpha_j^{\text{new}} - \alpha_j^{\text{old}} ) K_{jj} \]

最终的更新方式：

• 如果 \(0 < \alpha_i < C\)：设 \(b = b_1\)
• 否则如果 \(0 < \alpha_j < C\)：设 \(b = b_2\)
• 否则取均值 \(b = (b_1 + b_2)/2\)

SOR（逐次超松弛）

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

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

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

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

\[A \mathbf{x} = \mathbf{b} \]

其中：

• \(A \in \mathbb{R}^{n \times n}\)：系数矩阵
• \(\mathbf{x} \in \mathbb{R}^n\)：待求解向量
• \(\mathbf{b} \in \mathbb{R}^n\)：常数向量

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

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

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

\[A\mathbf{x} = \mathbf{b} \]

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

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

\[A = D + L + U \]

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

将 \(A\mathbf{x} = \mathbf{b}\) 改写为：

\[D\mathbf{x} = \mathbf{b} - (L + U)\mathbf{x} \]

两边左乘 \(D^{-1}\)，得到不动点形式：

\[\mathbf{x} = D^{-1} \left( \mathbf{b} - (L + U)\mathbf{x} \right) \]

自然导出 Jacobi 迭代公式：

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

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

对于方程组第 \(i\) 行：

\[a_{i1}x_1 + \cdots + a_{ii}x_i + \cdots + a_{in}x_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_1^{(k+1)}, \dots, x_{i-1}^{(k+1)}\)，导致信息利用不充分。

其收敛性取决于迭代矩阵 \(M = D^{-1}(L + U)\) 的谱半径 \(\rho(M)\)：

• 若 \(\rho(M) < 1\)，则收敛
• 若 \(\rho(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 加入一个松弛因子 \(\omega \in (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) \]

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

通过松弛因子\(\omega\)引入惯性项，动态调整步长：

\[x_i^{(k+1)}=(1-\omega)x_i^{(k)}+\omega(Gauss-Seidel更新量) \]

\(\omega>1\):放大修正步长(超松弛),加速逼近解。\(\omega<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) \]

其中损失函数\(\xi\)为L1或L2损失：

• L1-SVM：\(\xi(w; x_i, y_i) = \max(1 - y_i w^T x_i, 0)\)

• L2-SVM：\(\xi(w; x_i, y_i) = \max(1 - y_i w^T x_i, 0)^2\)

对偶问题推导：引入拉格朗日乘子\(\alpha_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 \]

其中：

• \(\bar{Q} = Q + D\)，\(Q_{ij} = y_i y_j x_i^T x_j\)

• L1-SVM：\(D_{ii}=0\)，\(U=C\)

• L2-SVM：\(D_{ii}=1/(2C)\)，\(U=\infty\)

为了方便叙述，先定义这样一组向量

\[\boldsymbol{\alpha}^{k,i}=[\alpha_1^{k+1},\ldots,\alpha_{i-1}^{k+1},\alpha_i^k,\ldots,\alpha_l^k]^T \]

此时的\(\alpha\)表示正在进行第K+1轮迭代，按顺序已经迭代完了前i-1个分量，正准备迭代其中第i个分量。

DCDM的子问题定义：
固定除\(\alpha_i\)外的所有变量，优化单变量\(\alpha_i\)：

\[\min_{d} f(\alpha + d e_i) \quad \text{s.t.} \quad 0 \leq \alpha_i + d \leq 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{常数项} \]

其中梯度分量：
\(\nabla_i f(\alpha) = (\bar{Q} \alpha)_{i} - 1 = y_i w^T x_i - 1 + D_{ii} \alpha_i\) \(w = \sum_j y_j \alpha_j x_j\)。

闭式解推导：

1. 二次函数极小值：最优\(d\)为无约束解：

\[d^* = -\frac{\nabla_i f(\alpha)}{\bar{Q}_{ii}} \]

1. 投影到可行域：

\[\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）区间上)

此时问题变成了如何高效求解\(\nabla_i f(\alpha)\)，

自然地想法是直接求解\(\nabla_i f(\boldsymbol{\alpha})=(\bar{Q}\boldsymbol{\alpha})_i-1\)=

\(\sum_{j=1}^l\bar{Q}_{ij}\alpha_j-1\),

其中\(\bar{Q}_{ij}=x_i^\top x_j\),

时间复杂度是\(O(l\bar{n})\),

\(\bar{n}\)是平均非零特征数，这样的计算开销是巨大的。

（可能会疑惑为何不预先计算出整个矩阵Q？事实上是因为存储成本太高！若样本数\(l=10^5\), 那么Q存储\(10^{10}\)个浮点数，粗略记每个单精度浮点数4字节，此时内存开销已经来到将近40GB，完全不可行）

因此我们计算分量式子

\[\nabla_i f(\boldsymbol{\alpha})=y_i\boldsymbol{w}^T\boldsymbol{x}_i-1+D_{ii}\alpha_i \]

当然，如果按照定义式子\(\boldsymbol{w}=\sum_{j=1}^ly_j\alpha_j\boldsymbol{x}_j\)求解\(w\),复杂度就回到了\(O(l\bar{n})\),为此作者也对\(w\)在每一轮进行迭代

\[\boldsymbol{w}\leftarrow\boldsymbol{w}+(\alpha_i-\bar{\alpha}_i)y_i\boldsymbol{x}_i \]

根据当前 \(\alpha_i\) 的位置，PG的计算分为以下三种情况：

\(\alpha_i\) 的位置	PG的计算公式	物理意义
\(0 < \alpha_i < U\)（内部）	\(PG = G\)	变量未触及边界，直接使用原始梯度 \(G\) 的方向更新。
\(\alpha_i = 0\)（下边界）	\(PG = \min(G, 0)\)	仅允许梯度方向为负（\(G < 0\)）时更新，避免 \(\alpha_i < 0\)。
\(\alpha_i = U\)（上边界）	\(PG = \max(G, 0)\)	仅允许梯度方向为正（\(G > 0\)）时更新，避免 \(\alpha_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 \succeq 0\)，\(A \in \mathbb{R}^{m \times 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 \in \mathbb{R}^m\) 为对偶变量，\(\rho > 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 + \rho A^T A)^{-1} \left( q + A^T (\rho z^k - y^k) \right) \]

(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 + \rho (A x^{k+1} - z^{k+1}) \]

完整算法流程

• 输入：初始 \(x^0, z^0, y^0\)，参数 \(\rho > 0\)，容差 \(\epsilon\)。
  迭代直到收敛：
• 1. 更新 \(x^{k+1} = -(P + \rho A^T A)^{-1} (q + A^T (\rho 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 + \rho (A x^{k+1} - z^{k+1})\)
• 4. 检查终止条件 \(\|r^k\| \leq \epsilon\), \(\|s^k\| \leq \epsilon\)

\(\nabla_i f(\alpha) = (\bar{Q} \alpha)i - 1 = y_i w^T x_i - 1 + D{ii} \alpha_iw = \sum_j y_j \alpha_j x_j\)