求解大规模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\)。
闭式解推导：
二次函数极小值：最优\(d\)为无约束解：
\[d^* = -\frac{\nabla_i f(\alpha)}{\bar{Q}_{ii}}
\]
投影到可行域：
\[\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\)。
迭代直到收敛：
更新 \(x^{k+1} = -(P + \rho A^T A)^{-1} (q + A^T (\rho z^k - y^k))\)
更新 \(z^{k+1} = \Pi_{[l, u]} \left( A x^{k+1} + \frac{y^k}{\rho} \right)\)
更新 \(y^{k+1} = y^k + \rho (A x^{k+1} - z^{k+1})\)
检查终止条件 \(\|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\)