凸优化算法路线图——从梯度下降到内点法

前置自测

📋 前置自测（答不出 ≥ 2 题 → 先回 10_凸分析基础、20_共轭函数与proximal算子 与 30_凸优化问题与对偶理论 复习）

1. 什么是 \(L\)-光滑性？写出其等价的梯度 Lipschitz 条件和二次上界不等式。
2. \(\mu\)-强凸函数的定义是什么？它如何保证极小点唯一？
3. 写出 proximal 算子 \(\text{prox}_{\alpha g}(v)\) 的定义，并给出 \(g = \|\cdot\|_1\) 时的闭式解。
4. 什么是 Slater 条件？它在对偶理论中起什么作用？
5. 条件数 \(\kappa = L/\mu\) 的物理含义是什么？它如何影响优化算法的收敛速度？

---

本章目标

学完本章后，你将能够：(1) 对任意凸优化问题，根据其结构（光滑性、强凸性、约束类型、规模）选择最优算法并给出收敛率的理论依据；(2) 独立推导梯度下降、加速梯度法、FISTA、ADMM 的收敛性证明；(3) 理解内点法的多项式时间保证与 self-concordance 理论；(4) 掌握 OSQP、HPIPM、acados 等机器人实战求解器的底层算法原理与选型依据。

---

1. 全景地图与算法演进脉络 ⭐

1.1 动机：为什么需要一张"算法路线图"？

假设你正在为一台四足机器人设计模型预测控制器（MPC）。控制器每 20 毫秒需要求解一个包含 200 个决策变量和 500 个约束的二次规划（QP）。你面前摆着至少五种可选算法：梯度下降、加速梯度法、ADMM、内点法、active-set 方法。选哪一个？为什么？ 如果选错了，你的机器人要么因为求解超时而摔倒，要么因为精度不够而抖动。

这不是一个可以靠"试试看"回答的问题。每种算法都有严格的收敛率理论，这些理论告诉我们：在给定的问题结构（光滑性、强凸性、约束数量、稀疏模式）下，每种算法到达 \(\epsilon\)-最优解需要多少次迭代，每次迭代的计算代价是多少。算法选择 = 识别问题结构 + 查收敛率表 + 匹配计算预算——这正是本章要教会你的能力。

1.2 如果不了解算法理论会怎样

本质洞察：凸优化算法的核心叙事是"结构换速度"——你对问题结构知道得越多（光滑、强凸、可分、稀疏），就能选择越快的算法。不了解理论的工程师只能在黑箱求解器之间盲目切换。

工程实践中最常见的三类失败：

失败 1：选了"万能"算法却太慢。 用通用内点法求解一个 \(10^5\) 变量的 LASSO 问题。IPM 每步需要 \(O(n^3)\) 的 Newton 系统求解，\(n = 10^5\) 意味着单步就要几分钟——而 FISTA（利用 \(\ell_1\) 的 prox 结构）可以在几秒内到达工程精度。

失败 2：追求理论最优阶却忽略常数。 FISTA 理论上比 ADMM 快一阶（\(O(1/k^2)\) vs \(O(1/k)\)），但在 MPC 的时序求解场景中，ADMM 的 warm-start 效率让它实际更快。理论阶次是渐近行为，实际性能还取决于常数、warm-start、因子缓存。

失败 3：不识别问题结构，暴力通用化。 把具有时域耦合（block-tridiagonal）结构的 MPC QP 交给通用 QP 求解器，复杂度是 \(O(N^3 T^3)\)（\(N\) 状态维度，\(T\) 时域步数）。而 HPIPM 利用 Riccati 递归结构，复杂度降为 \(O(N^3 T)\)——快了 \(T^2\) 倍。

1.3 历史演进脉络

凸优化算法的发展可以用"三次加速革命"来概括：

时代	核心突破	代表方法	收敛率（光滑凸）
1847	Cauchy 最速下降	梯度下降	\(O(1/k)\)
1964	Polyak 动量	Heavy-ball	二次上 \(O((\sqrt{\kappa}-1)/(\sqrt{\kappa}+1))^k\)
1983	Nesterov 加速	AGD	\(O(1/k^2)\)——理论最优
1984	Karmarkar 内点	IPM	\(O(\sqrt{m}\log(1/\epsilon))\) Newton 步
2004	近端算子 + 加速	FISTA	非光滑复合 \(O(1/k^2)\)
2011	算子分裂 + 分布式	ADMM (Boyd)	\(O(1/k)\) 但常数极小 + 可分布
2020	结构化 IPM + Riccati	HPIPM/acados	机器人 kHz 级 MPC

本质洞察：Nesterov 1983 年的加速不是"又一个技巧"，而是一个理论极限——他同时证明了 \(O(1/k^2)\) 是**不可再改进的**（黑箱下界定理）。此后所有"加速"都是在**结构假设**（强凸、可 prox、有限和、稀疏）上换来的，黑箱世界已经被封顶。

1.4 知识体系结构

本章的知识树如下：

凸优化算法
├── 一阶无约束方法
│   ├── 梯度下降（§2）⭐⭐ ——— 一切的起点
│   ├── Nesterov 加速（§3）⭐⭐⭐ ——— 理论最优
│   ├── 近端梯度法 ISTA/FISTA（§4）⭐⭐ ——— 非光滑复合
│   └── 其他一阶法（§7）⭐⭐⭐
│       ├── Frank-Wolfe ——— 投影贵时的替代
│       ├── Mirror Descent ——— 几何自适应
│       └── 次梯度法 ——— 非光滑兜底
├── 分裂方法
│   └── ADMM（§5）⭐⭐ ——— 可分结构 + 分布式
├── 二阶 / 内点法
│   └── IPM（§6）⭐⭐⭐ ——— 高精度 + 约束
├── 随机方法
│   └── SGD / SVRG / SAGA（§7.4）⭐⭐⭐
└── 收敛率下界与 PEP（§8）⭐⭐⭐⭐ ——— 理论天花板

阅读建议：§2-§5 是主干（⭐⭐），必须逐行推导；§6 是高精度与约束场景的核心工具；§7-§8 是研究级拓展。

⚠️ 常见陷阱

💡 概念误区：认为"越新的算法越好" - 新手想法："ADMM 是 2011 年的，FISTA 是 2009 年的，所以 ADMM 比 FISTA 好。" - 实际上：算法没有"更好"之说，只有"更适合"。FISTA 对 \(f+g\)（\(g\) 有廉价 prox）结构最优；ADMM 对两块都非光滑的可分结构最优；IPM 对中小规模高精度约束问题最优。 - 正确思维：先识别问题结构，再查收敛率表。

🧠 思维陷阱：混淆"迭代复杂度"与"总计算时间" - 新手想法："FISTA 只需 \(O(\sqrt{L/\epsilon})\) 次迭代，IPM 需要 \(O(\sqrt{m}\log(1/\epsilon))\) 次 Newton 步，所以 FISTA 更快。" - 实际上：FISTA 每步 \(O(n)\)（一次梯度 + 一次 prox），IPM 每步 \(O(n^3)\)（一次 Newton 系统）。总时间 = 迭代数 × 每步代价。对 \(n=100\) 的 QP，IPM 可能 20 步（每步 \(O(10^6)\)）就结束，FISTA 可能需要 \(10^4\) 步（每步 \(O(100)\)）。 - 分水岭：中小规模高精度 → IPM；大规模中精度 → 一阶法。

练习

1. 画出上述知识树的详细版本，为每个节点标注"适用的问题结构"和"收敛率"。
2. 对比梯度下降与内点法在以下三个场景中的总计算时间：(a) \(n=50\) 的 QP，\(\epsilon = 10^{-8}\)；(b) \(n = 10^5\) 的 LASSO，\(\epsilon = 10^{-3}\)；(c) \(n = 500\) 的 MPC QP（带 warm-start），\(\epsilon = 10^{-4}\)。

---

2. 梯度下降完整理论 ⭐⭐

2.1 动机：最简单的优化算法

考虑无约束凸优化问题 \(\min_{x \in \mathbb{R}^n} f(x)\)，其中 \(f\) 是凸且可微的。最自然的想法是：沿着目标函数下降最快的方向走一步。负梯度 \(-\nabla f(x)\) 正是这个方向（在欧氏范数意义下的最速下降方向）。

\[x_{k+1} = x_k - \alpha \nabla f(x_k)\]

这就是梯度下降法（Gradient Descent, GD）。Cauchy 在 1847 年首次提出这个想法。但一个关键问题立刻出现：步长 \(\alpha\) 该选多大？

• 太大：跨过最优点，来回振荡，甚至发散
• 太小：每步前进极少，需要天文数字般的迭代次数
• 刚好：存在一个"甜点"步长，使收敛最快

本质洞察：梯度下降的收敛率完全由**条件数** \(\kappa = L/\mu\) 决定（\(L\) 是 Lipschitz 常数，\(\mu\) 是强凸参数）。条件数衡量的是目标函数的"椭球性"——\(\kappa\) 越大，等高线越像扁长的椭球，梯度方向与最优方向的夹角越大，收敛越慢。

2.2 下降引理——一阶方法的基石

定理 1（下降引理 / Descent Lemma）。 若 \(f\) 是 \(L\)-光滑的（即 \(\nabla f\) 满足 \(\|\nabla f(x) - \nabla f(y)\| \le L\|x - y\|\)），则对任意 \(x, y \in \mathbb{R}^n\)：

\[f(y) \le f(x) + \langle \nabla f(x), y - x \rangle + \frac{L}{2}\|y - x\|^2\]

为什么这个引理如此重要？ 它提供了一个可计算的二次上界——给定当前点 \(x\) 和梯度 \(\nabla f(x)\)，我们可以构造一个抛物面，保证 \(f\) 在任何地方都不超过这个抛物面。梯度下降法本质上就是在每步最小化这个二次上界。

完整证明（从 \(L\)-光滑性出发）：

Step 1：定义辅助函数 \(g(t) = f(x + t(y-x))\)，则 \(g(0) = f(x)\)，\(g(1) = f(y)\)。

Step 2：由微积分基本定理： $\(f(y) - f(x) = g(1) - g(0) = \int_0^1 g'(t)\,dt = \int_0^1 \langle \nabla f(x + t(y-x)),\, y-x \rangle\,dt\)$

Step 3：加减 \(\langle \nabla f(x), y-x \rangle\)： $\(f(y) - f(x) = \langle \nabla f(x), y-x \rangle + \int_0^1 \langle \nabla f(x + t(y-x)) - \nabla f(x),\, y-x \rangle\,dt\)$

Step 4：对积分项用 Cauchy-Schwarz 和 \(L\)-Lipschitz： $\(|\langle \nabla f(x + t(y-x)) - \nabla f(x), y-x \rangle| \le \|\nabla f(x + t(y-x)) - \nabla f(x)\| \cdot \|y-x\| \le Lt\|y-x\|^2\)$

Step 5：积分： $\(f(y) - f(x) \le \langle \nabla f(x), y-x \rangle + L\|y-x\|^2 \int_0^1 t\,dt = \langle \nabla f(x), y-x \rangle + \frac{L}{2}\|y-x\|^2 \quad \blacksquare\)$

直觉理解：\(L\)-光滑性意味着 \(f\) 的曲率不超过 \(L\)——函数不会"弯得太快"。因此，一阶 Taylor 展开加上 \(\frac{L}{2}\|y-x\|^2\) 的修正就足以把 \(f\) 盖住。

2.3 凸情形 \(O(1/k)\) 收敛率

定理 2。 若 \(f\) 凸且 \(L\)-光滑，步长 \(\alpha = 1/L\)，则梯度下降满足：

\[f(x_k) - f^* \le \frac{2L\|x_0 - x^*\|^2}{k}\]

完整证明（5 步，无跳步）：

Step 1：单步下降。令 \(y = x_{k+1} = x_k - \frac{1}{L}\nabla f(x_k)\)，代入下降引理： $\(f(x_{k+1}) \le f(x_k) + \langle \nabla f(x_k), -\frac{1}{L}\nabla f(x_k) \rangle + \frac{L}{2}\|\frac{1}{L}\nabla f(x_k)\|^2\)$ $\(= f(x_k) - \frac{1}{L}\|\nabla f(x_k)\|^2 + \frac{1}{2L}\|\nabla f(x_k)\|^2 = f(x_k) - \frac{1}{2L}\|\nabla f(x_k)\|^2\)$

这告诉我们每步至少下降 \(\frac{1}{2L}\|\nabla f(x_k)\|^2\)。

Step 2：利用凸性连接到最优值。凸性给出 \(f^* \ge f(x_k) + \langle \nabla f(x_k), x^* - x_k \rangle\)，即： $\(f(x_k) - f^* \le \langle \nabla f(x_k), x_k - x^* \rangle \le \|\nabla f(x_k)\| \cdot \|x_k - x^*\|\)$

由 Cauchy-Schwarz，\(\|\nabla f(x_k)\|^2 \ge \frac{(f(x_k) - f^*)^2}{\|x_k - x^*\|^2}\)。

Step 3：合并 Step 1 和 Step 2。令 \(\delta_k = f(x_k) - f^*\)，\(R = \|x_0 - x^*\|\)： $\(\delta_{k+1} \le \delta_k - \frac{1}{2L}\|\nabla f(x_k)\|^2 \le \delta_k - \frac{\delta_k^2}{2L R^2}\)$

（这里用了 \(\|x_k - x^*\| \le \|x_0 - x^*\| = R\)，由步长 \(1/L\) 的非扩张性保证。）

Step 4：取倒数。令 \(a_k = 1/\delta_k\)： $\(a_{k+1} = \frac{1}{\delta_{k+1}} \ge \frac{1}{\delta_k - \frac{\delta_k^2}{2LR^2}} = \frac{1}{\delta_k} \cdot \frac{1}{1 - \frac{\delta_k}{2LR^2}} \ge a_k + \frac{1}{2LR^2}\)$

（最后一步用了 \(\frac{1}{1-t} \ge 1+t\) 对 \(t \in [0,1)\)。）

Step 5：递归得 \(a_k \ge a_0 + \frac{k}{2LR^2}\)，即 \(\frac{1}{\delta_k} \ge \frac{k}{2LR^2}\)（忽略 \(a_0\) 项），因此： $\(f(x_k) - f^* = \delta_k \le \frac{2LR^2}{k} = \frac{2L\|x_0 - x^*\|^2}{k} \quad \blacksquare\)$

工程解读：要达到 \(\epsilon\)-最优，需要 \(k \ge \frac{2LR^2}{\epsilon} = O(L/\epsilon)\) 次迭代。对光滑凸函数，GD 的收敛率是**次线性**的——\(1/k\) 意味着精度每提高一倍需要加倍迭代次数。

2.4 强凸情形线性收敛

定理 3。 若 \(f\) 还是 \(\mu\)-强凸的，步长 \(\alpha = 1/L\)，则：

\[\|x_k - x^*\|^2 \le \left(1 - \frac{\mu}{L}\right)^k \|x_0 - x^*\|^2\]

完整证明：

Step 1：强凸性给出更强的下界：\(f(y) \ge f(x) + \langle \nabla f(x), y-x \rangle + \frac{\mu}{2}\|y-x\|^2\)。

Step 2：取 \(y = x^*\)，利用 \(\nabla f(x^*) = 0\)（最优性条件）： $\(f(x) - f^* \le \langle \nabla f(x), x - x^* \rangle - \frac{\mu}{2}\|x - x^*\|^2\)$

Step 3：利用强凸-光滑互补不等式（co-coercivity）：对 \(\mu\)-强凸且 \(L\)-光滑的 \(f\)： $\(\langle \nabla f(x) - \nabla f(y), x - y \rangle \ge \frac{\mu L}{\mu + L}\|x-y\|^2 + \frac{1}{\mu + L}\|\nabla f(x) - \nabla f(y)\|^2\)$

取 \(y = x^*\)：\(\langle \nabla f(x), x - x^* \rangle \ge \frac{\mu L}{\mu + L}\|x - x^*\|^2 + \frac{1}{\mu + L}\|\nabla f(x)\|^2\)

Step 4：计算 \(\|x_{k+1} - x^*\|^2\)： $\(\|x_{k+1} - x^*\|^2 = \|x_k - \frac{1}{L}\nabla f(x_k) - x^*\|^2 = \|x_k - x^*\|^2 - \frac{2}{L}\langle \nabla f(x_k), x_k - x^* \rangle + \frac{1}{L^2}\|\nabla f(x_k)\|^2\)$

Step 5：代入 co-coercivity（Step 3 的结果）： $\(\le \|x_k - x^*\|^2 - \frac{2}{L}\left[\frac{\mu L}{\mu+L}\|x_k-x^*\|^2 + \frac{1}{\mu+L}\|\nabla f(x_k)\|^2\right] + \frac{1}{L^2}\|\nabla f(x_k)\|^2\)$

\[= \left(1 - \frac{2\mu}{\mu+L}\right)\|x_k-x^*\|^2 + \left(\frac{1}{L^2} - \frac{2}{L(\mu+L)}\right)\|\nabla f(x_k)\|^2\]

第二项的系数 \(= \frac{\mu+L - 2L}{L^2(\mu+L)} = \frac{\mu-L}{L^2(\mu+L)} \le 0\)（因为 \(\mu \le L\)），所以：

\[\|x_{k+1} - x^*\|^2 \le \left(1 - \frac{2\mu}{\mu+L}\right)\|x_k - x^*\|^2 \le \left(1 - \frac{\mu}{L}\right)\|x_k - x^*\|^2 \quad \blacksquare\]

（最后一步用了 \(1 - \frac{2\mu}{\mu+L} = \frac{L-\mu}{L+\mu} \le 1 - \frac{\mu}{L}\)。更紧的率是 \(\left(\frac{\kappa-1}{\kappa+1}\right)^2\)，通过精确步长 \(\alpha = \frac{2}{\mu+L}\) 可达。）

工程解读：线性收敛意味着每步误差乘以固定常数 \(1 - 1/\kappa\)。达到 \(\epsilon\) 精度需要 \(k = O(\kappa \log(1/\epsilon))\) 步。当 \(\kappa = 10^4\)（机器人逆动力学 Hessian 的典型值）时，\(\log(10^{-6})/\log(1-10^{-4}) \approx 1.4 \times 10^5\) 次迭代——这就是加速法的出场动机。

2.5 步长选择：回溯线搜索

实践中 \(L\) 常常未知或难以精确估计。Armijo 回溯线搜索提供了一种自适应选择步长的策略：

Armijo 条件：\(f(x - \alpha \nabla f(x)) \le f(x) - c_1 \alpha \|\nabla f(x)\|^2\)，其中 \(c_1 \in (0,1)\)（通常取 \(c_1 = 10^{-4}\)）。

回溯策略：从 \(\alpha_0 = 1\) 开始，若不满足 Armijo 则 \(\alpha \leftarrow \rho \alpha\)（\(\rho = 0.5\) 或 \(0.9\)），直到满足。

# Armijo 回溯线搜索的 Python 实现
def armijo_backtrack(f, grad_f, x, d, c1=1e-4, rho=0.5, alpha0=1.0):
    """
    参数：
        f: 目标函数
        grad_f: 梯度函数
        x: 当前点
        d: 下降方向（通常为 -grad_f(x)）
        c1: Armijo 参数（充分下降条件）
        rho: 步长缩减系数
        alpha0: 初始步长
    返回：满足 Armijo 条件的步长 alpha
    """
    alpha = alpha0
    fx = f(x)
    gx = grad_f(x)
    slope = gx @ d  # 方向导数，必须 < 0（下降方向）

    while f(x + alpha * d) > fx + c1 * alpha * slope:
        alpha *= rho  # 步长缩小为原来的 rho 倍

    return alpha

为什么 Armijo 总是终止？ 由 Taylor 展开，当 \(\alpha \to 0\) 时 \(f(x + \alpha d) \approx f(x) + \alpha \nabla f(x)^T d\)。因为 \(d\) 是下降方向（\(\nabla f(x)^T d < 0\)），对足够小的 \(\alpha\) 必然满足 \(f(x+\alpha d) \le f(x) + c_1 \alpha \nabla f(x)^T d\)。

2.6 与机器人学的联系

回顾凸分析：梯度下降虽然简单，但它是所有复杂优化器的"原子操作"。在机器人学中：

• SLAM 后端的 Gauss-Newton：本质上是对 \(J^TJ\) 近似 Hessian 后的 Newton 步，当残差小时每步等价于在二次近似上做"一步梯度下降"
• RL 的策略梯度：REINFORCE 算法是目标函数 \(J(\theta) = \mathbb{E}[R]\) 上的随机梯度下降
• MPC 的预条件迭代：OSQP 的 ADMM 子步骤中涉及梯度下降型更新

⚠️ 常见陷阱

⚠️ 编程陷阱：步长 \(1/L\) 中 \(L\) 的估计错误 - 错误做法：凭经验设 \(\alpha = 0.01\) - 现象：收敛极慢或振荡发散 - 根本原因：\(L\) 是目标函数梯度的 Lipschitz 常数，取决于 Hessian 的最大特征值。对二次函数 \(f(x) = \frac{1}{2}x^TQx\)，\(L = \lambda_{\max}(Q)\) - 正确做法：使用 Armijo 回溯自适应选择步长；或用 \(L\) 的上界估计

💡 概念误区：认为梯度方向是"到最优点"的方向 - 新手想法："负梯度指向最优点的方向" - 实际上：负梯度是**局部最速下降方向**（在欧氏范数意义下），但只有对球形等高线（\(\kappa = 1\)）才指向最优点。对椭球形等高线，梯度方向可能与最优方向几乎垂直——这正是条件数大时收敛慢的几何原因 - 正确理解：梯度 = 等高面的法向量，最优方向 = 当前点到最优点的连线。两者夹角由条件数决定

🧠 思维陷阱：认为"收敛率 \(O(1/k)\) 已经够好了" - 对 \(\epsilon = 10^{-6}\)，\(O(1/k)\) 需要 \(10^6\) 次迭代 - 加速到 \(O(1/k^2)\) 只需要 \(10^3\) 次迭代——快了 1000 倍 - 再到线性收敛 \(O((1-c)^k)\)，只需要 \(\frac{6}{c}\log 10 \approx 14/c\) 次

练习

1. （手推） 对二次函数 \(f(x) = \frac{1}{2}(x_1^2 + 100x_2^2)\)（\(\kappa = 100\)），从 \(x_0 = (1, 1)\) 出发，手工计算 GD 前 5 步的轨迹，画出在等高线上的路径。
2. （编程） 实现带 Armijo 回溯的梯度下降，在 Rosenbrock 函数 \(f(x,y) = (1-x)^2 + 100(y-x^2)^2\) 上测试，记录迭代次数与步长的变化。
3. （证明） 证明：对 \(L\)-光滑凸函数，\(\frac{1}{2L}\|\nabla f(x)\|^2 \le f(x) - f^*\)。（提示：从下降引理出发，令 \(y = x - \frac{1}{L}\nabla f(x)\)。）

---

3. Nesterov 加速梯度法 ⭐⭐⭐

3.1 动机：能否打破 \(O(1/k)\) 的壁垒？

梯度下降的 \(O(1/k)\) 收敛率看起来是"合理"的——每步利用一次梯度信息，似乎无法更快。但 Nesterov 在 1983 年（年仅 27 岁！）证明了一个惊人的结果：只需稍加修改迭代格式，就能达到 \(O(1/k^2)\)——平方加速。更惊人的是，他同时证明了这是**不可再改进的理论极限**。

本质洞察：Nesterov 加速的本质不是"走得更远"，而是"在动量方向上向前看一步再计算梯度"。这个小小的改变——先外推再求梯度（lookahead）——把收敛率从 \(O(1/k)\) 提升到 \(O(1/k^2)\)。直觉上，这相当于从"沿最陡方向爬山"升级为"带惯性的智能滑行"。

3.2 Nesterov 加速梯度法（AGD）的三种等价形式

形式 1（动量形式，最常用）： $\(y_k = x_k + \frac{k-1}{k+2}(x_k - x_{k-1})\)$ $\(x_{k+1} = y_k - \frac{1}{L}\nabla f(y_k)\)$

形式 2（三序列形式，证明最方便）： $\(y_k = (1-\tau_k)x_k + \tau_k v_k\)$ $\(x_{k+1} = y_k - \frac{1}{L}\nabla f(y_k)\)$ $\(v_{k+1} = v_k - \frac{1}{\tau_k L}\nabla f(y_k)\)$

其中 \(\tau_k = 2/(k+2)\)。

形式 3（ODE 极限，物理直觉最强）： $\(\ddot{X}(t) + \frac{3}{t}\dot{X}(t) + \nabla f(X(t)) = 0\)$

这是一个带时变阻尼的二阶 ODE（Su-Boyd-Candes 2016）。\(3/t\) 阻尼系数随时间衰减——起初阻尼大（稳定启动），后来阻尼小（允许加速）。

3.3 与 Heavy-Ball 的关键区别

Polyak Heavy-ball 方法（1964）： $\(x_{k+1} = x_k - \alpha \nabla f(x_k) + \beta(x_k - x_{k-1})\)$

看起来和 Nesterov AGD 几乎一样——都用了"动量"\(\beta(x_k - x_{k-1})\)。但两者有本质区别：

特性	Heavy-ball	Nesterov AGD
梯度计算点	当前点 \(x_k\)	外推点 \(y_k = x_k + \beta(x_k - x_{k-1})\)
参数 \(\beta\)	固定常数	时变 \(\frac{k-1}{k+2}\)（关键！）
一般凸保证	无 \(O(1/k^2)\) 保证	\(O(1/k^2)\)
强凸保证	\((1-2\sqrt{\mu/L})^k\)（仅二次）	\((1-\sqrt{\mu/L})^k\)（一般强凸）
反例	Lessard-Recht-Packard 2016	无（已证最优）

⚠️ 这是初学者最大的陷阱之一：Heavy-ball 在二次函数上表现完美，但在一般凸函数上**没有加速保证**。Lessard 等人 2016 年用 IQC（积分二次约束，控制论工具）框架构造了 heavy-ball 不收敛的反例。永远不要把 heavy-ball 当成 Nesterov 的等价替代。

3.4 收敛率证明（Lyapunov 方法）

定理 4（Nesterov 加速 \(O(1/k^2)\)）。 对 \(L\)-光滑凸 \(f\)，AGD 满足： $\(f(x_k) - f^* \le \frac{2L\|x_0 - x^*\|^2}{(k+1)^2}\)$

证明策略：构造 Lyapunov 函数（能量函数）\(E_k\)，证明 \(E_k\) 单调不增。

Step 1：定义 Lyapunov 函数： $\(E_k = t_k^2 (f(x_k) - f^*) + \frac{L}{2}\|v_k - x^*\|^2\)$

其中 \(t_0 = 1\)，\(t_{k+1}\) 满足 \(t_{k+1}^2 - t_{k+1} = t_k^2\)（即 \(t_k \approx k/2\)）。

Step 2：证明 \(E_{k+1} \le E_k\)。利用下降引理： $\(f(x_{k+1}) \le f(y_k) - \frac{1}{2L}\|\nabla f(y_k)\|^2\)$

以及凸性：\(f(y_k) \le (1-\tau_k)f(x_k) + \tau_k f(x^*) + \tau_k \langle \nabla f(y_k), v_k - x^* \rangle\)

Step 3：将上述两式合并，利用 \(v_{k+1} = v_k - \frac{1}{\tau_k L}\nabla f(y_k)\)： $\(t_{k+1}^2 (f(x_{k+1}) - f^*) + \frac{L}{2}\|v_{k+1} - x^*\|^2 \le t_k^2(f(x_k) - f^*) + \frac{L}{2}\|v_k - x^*\|^2\)$

即 \(E_{k+1} \le E_k\)。

Step 4：由 \(E_k \le E_0 = f(x_0) - f^* + \frac{L}{2}\|x_0 - x^*\|^2 \le \frac{L}{2}\|x_0-x^*\|^2 + \frac{L}{2}\|x_0-x^*\|^2 = L\|x_0-x^*\|^2\)（利用了 \(f(x_0)-f^* \le \frac{L}{2}\|x_0-x^*\|^2\)）。

又 \(E_k \ge t_k^2(f(x_k) - f^*)\)，\(t_k \ge (k+1)/2\)，所以： $\(f(x_k) - f^* \le \frac{E_k}{t_k^2} \le \frac{L\|x_0-x^*\|^2}{(k+1)^2/4} = \frac{4L\|x_0-x^*\|^2}{(k+1)^2} \quad \blacksquare\)$

（精确常数取决于 \(t_k\) 的初始化方式，标准结果为 \(\frac{2L\|x_0-x^*\|^2}{(k+2)^2}\)。）

3.5 强凸情形：\(\sqrt{\kappa}\) 替代 \(\kappa\)

对 \(\mu\)-强凸且 \(L\)-光滑的 \(f\)，AGD 以固定动量 \(\beta = \frac{\sqrt{L}-\sqrt{\mu}}{\sqrt{L}+\sqrt{\mu}} = \frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\) 给出：

\[f(x_k) - f^* \le L\|x_0 - x^*\|^2 \left(1 - \sqrt{\frac{\mu}{L}}\right)^k = L\|x_0-x^*\|^2 \left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^k\]

达到 \(\epsilon\) 精度需要 \(O(\sqrt{\kappa}\log(1/\epsilon))\) 步——比 GD 的 \(O(\kappa\log(1/\epsilon))\) 快了 \(\sqrt{\kappa}\) 倍。

数值对比：\(\kappa = 10^4\)，\(\epsilon = 10^{-6}\)： - GD：\(\kappa \cdot 14 = 1.4 \times 10^5\) 步 - AGD：\(\sqrt{\kappa} \cdot 14 = 1400\) 步——快了 100 倍！

3.6 加速的三种理解视角

视角 1：Estimate Sequence（Nesterov 原始）。构造一族递减的二次上界 \(\phi_k(x)\) 使得 \(\min \phi_k \ge f(x_k) - \epsilon_k\)。这是最严谨但最晦涩的路线。

视角 2：ODE 极限（Su-Boyd-Candes 2016）。加速梯度法是 ODE \(\ddot{X} + \frac{3}{t}\dot{X} + \nabla f = 0\) 的离散化。物理类比：在势能场中滑动的粒子，阻尼随时间减小。

ODE 参数	优化对应	物理类比
\(X(t)\)	迭代点 \(x_k\)	粒子位置
\(\dot{X}(t)\)	动量 \(x_k - x_{k-1}\)	粒子速度
\(\nabla f(X)\)	梯度力	势能场的力
\(\frac{3}{t}\dot{X}\)	时变阻尼	空气阻力（随时间减弱）

视角 3：PEP（Performance Estimation Problem, Drori-Teboulle 2014）。把"找一阶方法最坏收敛率"本身写成一个 SDP（半定规划），数值求解得到紧的上下界。PEP 的意义在于它把算法分析从"天才的灵光一现"变成了"SDP 数值计算"——你可以自动验证任意新算法的最坏性能。

3.7 自适应重启策略

AGD 在实践中有一个问题：动量积累可能导致振荡（像一个停不下来的摆锤）。解决方案是**重启**（restart）：当检测到"目标值增大"或"梯度与动量方向夹角 > 90°"时，重置动量为零。

重启准则（O'Donoghue-Candes 2015）： - 目标重启：\(f(x_{k+1}) > f(x_k)\) 时重启 - 梯度重启：\(\langle \nabla f(y_k), x_k - x_{k-1} \rangle > 0\) 时重启

重启后的 AGD 在强凸场景下自动适应未知的 \(\mu\)，不需要知道强凸参数。

⚠️ 常见陷阱

⚠️ 编程陷阱：动量参数 \(\beta_k = (k-1)/(k+2)\) 中 \(k\) 从 0 还是 1 开始 - 不同论文/库的惯例不同。PyTorch 的 nesterov=True 用的是固定 \(\beta\)（需要用户提供），不是时变的 \(\frac{k-1}{k+2}\) - 错误的 \(\beta\) 初始化可能导致首步出错（\(k=0\) 时 \(\beta = -1/2\)，不合理） - 正确做法：\(k\) 从 1 开始，\(\beta_1 = 0\)（第一步无动量，等价于纯 GD）

💡 概念误区：认为"加速 = 动量 = Polyak momentum" - 深度学习社区常说的"momentum SGD"是 Polyak heavy-ball 的随机版，不是 Nesterov AGD - 两者在一般非凸/凸下有本质区别 - Nesterov AGD 需要在**外推点**计算梯度，而 momentum SGD 在当前点计算

🧠 思维陷阱：认为"\(O(1/k^2)\) 不可改进"意味着没有更快的方法 - Nesterov 下界是在**黑箱一阶 oracle** 模型下证明的 - 如果你知道更多结构（强凸、有限和、可 prox），可以利用这些结构做得更好 - 例如：CG 在二次函数上 \(n\) 步精确收敛，不违反下界（因为它利用了二次结构）

练习

1. （手推） 对 \(f(x) = \frac{1}{2}(x_1^2 + 100x_2^2)\)，从 \(x_0 = (1,1)\) 出发，手工执行 AGD 5 步，并与 GD 的轨迹对比。
2. （编程） 实现带自适应重启的 FISTA，在 \(\ell_1\) 正则 logistic 回归上测试，绘制 \(f(x_k) - f^*\) 的对数图，验证 \(O(1/k^2)\) 的斜率。
3. （概念） 解释为什么 Conjugate Gradient 在 \(n\) 步内精确求解 \(n\) 维二次规划，而这不违反 Nesterov 下界。（提示：CG 利用了函数的二次结构，不是"黑箱 oracle"。）

---

4. 近端梯度法与 FISTA ⭐⭐

4.1 动机：当目标函数"部分不可微"时

许多机器人与机器学习问题有如下结构：

\[\min_x \underbrace{f(x)}_{\text{光滑，如最小二乘}} + \underbrace{g(x)}_{\text{非光滑，如 } \|\cdot\|_1 \text{ 正则化}}\]

例如： - LASSO（\(\ell_1\) 稀疏恢复）：\(f = \frac{1}{2}\|Ax - b\|^2\)，\(g = \lambda\|x\|_1\) - 鲁棒估计：\(f\) 是数据项，\(g\) 是指示函数 \(I_{\mathcal{C}}\)（约束投影） - 全变分去噪：\(f = \frac{1}{2}\|x - y\|^2\)，\(g = \lambda\|\nabla x\|_1\)

梯度下降对 \(g\) 无能为力（\(g\) 不可微，没有传统梯度）。但如果 \(g\) 的 proximal 算子**可以高效计算，我们就能用**近端梯度法（Proximal Gradient Method, PGM）。

4.2 Proximal 算子回顾

回顾凸分析章节中定义的 proximal 算子：

\[\text{prox}_{\alpha g}(v) = \arg\min_x \left\{ g(x) + \frac{1}{2\alpha}\|x - v\|^2 \right\}\]

直觉：在"最小化 \(g\)"和"不要离 \(v\) 太远"之间找平衡。

\(g(x)\)	\(\text{prox}_{\alpha g}(v)\)	计算代价
\(\lambda\|x\|_1\)	\(\text{sign}(v) \odot \max(\|v\| - \alpha\lambda, 0)\)（软阈值）	\(O(n)\)
\(I_{\mathcal{C}}(x)\)	\(\text{proj}_{\mathcal{C}}(v)\)（投影）	取决于 \(\mathcal{C}\)
\(\lambda\|x\|_2\)	\(\max(1 - \frac{\alpha\lambda}{\|v\|}, 0) \cdot v\)（块阈值）	\(O(n)\)
\(\lambda\|X\|_*\)（核范数）	\(U \text{diag}(\text{soft}(\sigma, \alpha\lambda)) V^T\)	\(O(n^3)\)（SVD）

4.3 ISTA（Iterative Shrinkage-Thresholding Algorithm）

ISTA = 梯度下降 + proximal 步：

\[x_{k+1} = \text{prox}_{\alpha g}\left(x_k - \alpha \nabla f(x_k)\right)\]

直觉：先沿 \(f\) 的梯度方向走一步（处理光滑部分），然后用 prox 算子"修正"到满足 \(g\) 结构的位置（处理非光滑部分）。

收敛率：\(F(x_k) - F^* \le \frac{L\|x_0 - x^*\|^2}{2k}\)，即 \(O(1/k)\)——和纯梯度下降一样。

4.4 FISTA（Fast ISTA）——Beck & Teboulle 2009

FISTA 是在 ISTA 外层加上 Nesterov 动量：

\[y_k = x_k + \frac{t_k - 1}{t_{k+1}}(x_k - x_{k-1})$$ $$x_{k+1} = \text{prox}_{\frac{1}{L} g}\left(y_k - \frac{1}{L}\nabla f(y_k)\right)\]

其中 \(t_1 = 1\)，\(t_{k+1} = \frac{1 + \sqrt{1 + 4t_k^2}}{2}\)。

定理 6（FISTA \(O(1/k^2)\)）。 对 \(F = f + g\)，\(f\) \(L\)-光滑凸，\(g\) 闭凸 proper： $\(F(x_k) - F^* \le \frac{2L\|x_0 - x^*\|^2}{(k+1)^2}\)$

证明核心思想：构造 Lyapunov 函数 \(E_k = t_k^2(F(x_k) - F^*) + \frac{L}{2}\|v_k - x^*\|^2\)，证明其单调不增。

Step 1：由 proximal 梯度步的下降性质（proximal descent lemma）： $\(F(x_{k+1}) \le F(y_k) - \frac{1}{2L}\|\nabla f(y_k)\|^2 + \langle \nabla f(y_k), y_k - x_{k+1} \rangle - \frac{L}{2}\|x_{k+1} - y_k\|^2\)$

更精确地（利用 prox 的最优性条件）： $\(F(x_{k+1}) \le F(u) + \langle \nabla f(y_k), y_k - u \rangle - \frac{L}{2}\|x_{k+1} - y_k\|^2 + \frac{L}{2}\|u - y_k\|^2 - \frac{L}{2}\|u - x_{k+1}\|^2\)$

对任意 \(u\)。

Step 2：取 \(u\) 的凸组合 \(u = (1-\tau_k)x_k + \tau_k x^*\)，利用 \(F\) 的凸性展开后合并，经过代数化简可得 \(E_{k+1} \le E_k\)。

Step 3：结论：\(F(x_k) - F^* \le \frac{E_k}{t_k^2} \le \frac{E_0}{t_k^2} \le \frac{2L\|x_0-x^*\|^2}{(k+1)^2}\)。

4.5 FISTA 在 LASSO 上的应用

LASSO 问题：\(\min_x \frac{1}{2}\|Ax - b\|^2 + \lambda\|x\|_1\)

• \(f(x) = \frac{1}{2}\|Ax - b\|^2\)，\(\nabla f(x) = A^T(Ax - b)\)，\(L = \|A^TA\| = \sigma_{\max}^2(A)\)
• \(g(x) = \lambda\|x\|_1\)，\(\text{prox}_{\alpha g}(v) = \text{soft}_{\alpha\lambda}(v)\)

import numpy as np

def fista_lasso(A, b, lam, L=None, max_iter=1000, tol=1e-6):
    """
    FISTA 求解 LASSO: min 0.5*||Ax-b||^2 + lam*||x||_1

    参数：
        A: m x n 观测矩阵
        b: m 维观测向量
        lam: 正则化参数
        L: Lipschitz 常数（若 None 则自动计算）
    """
    m, n = A.shape
    if L is None:
        # L = 最大奇异值的平方 = A^T A 的最大特征值
        L = np.linalg.norm(A, ord=2) ** 2

    # 软阈值函数（ell_1 的 prox 算子）
    def soft_threshold(v, threshold):
        return np.sign(v) * np.maximum(np.abs(v) - threshold, 0)

    x = np.zeros(n)       # 当前迭代点
    x_prev = x.copy()     # 上一步迭代点
    t = 1.0               # FISTA 的 t 序列

    for k in range(max_iter):
        # Nesterov 动量：先外推到 y_k
        t_new = (1 + np.sqrt(1 + 4 * t**2)) / 2
        y = x + (t - 1) / t_new * (x - x_prev)

        # 梯度步 + proximal 步
        grad = A.T @ (A @ y - b)           # f 的梯度
        x_new = soft_threshold(y - grad / L, lam / L)  # prox 步

        # 收敛检查
        if np.linalg.norm(x_new - x) < tol:
            break

        # 更新
        x_prev = x
        x = x_new
        t = t_new

    return x

4.6 prox 算子的计算代价决定算法选择

反事实推理：如果 \(g\) 的 prox 也很贵怎么办？

\(g\)	prox 代价	建议替代方案
\(\ell_1\)（软阈值）	\(O(n)\)	FISTA 完美适用
指示函数 \(I_{\mathcal{C}}\)（投影）	取决于 \(\mathcal{C}\)	简单集（盒、球）→ FISTA；复杂集 → ADMM 或 Frank-Wolfe
核范数	\(O(\min(m,n) \cdot mn)\)（SVD）	大规模 → Frank-Wolfe（避免 SVD）
Total Variation	\(O(n)\)（fused LASSO 的线性时间算法）	ADMM 分裂

关键洞察：FISTA 的每步代价 = 一次梯度计算 + 一次 prox 评估。如果 prox 评估本身就是 \(O(n^3)\)（如核范数），那么 FISTA 的每步代价主导项变成了 prox——此时不如用 Frank-Wolfe（每步只需一个线性最小化 oracle，\(O(n^2)\) 或更低）。

⚠️ 常见陷阱

⚠️ 编程陷阱：FISTA 的 \(t\) 序列初始化 - \(t_1\) 必须是 1（不是 0）。如果 \(t_0 = 0\) 会导致除以零 - 某些实现用 \(t_{k+1}^2 = t_k^2 + t_{k+1}\)（等价形式），注意哪种形式被使用

💡 概念误区：认为 prox 越复杂的问题就不能用近端方法 - prox 的"代价"是可以被利用的。例如核范数的 prox 需要 SVD，但如果问题本身就需要反复做 SVD（如低秩矩阵补全），那么 prox 的代价是"免费的" - 关键是**prox 代价相对于梯度代价**的比值

练习

1. （推导） 推导 \(\ell_2\) 范数（非平方）\(g(x) = \lambda\|x\|_2\) 的 proximal 算子闭式解。
2. （编程） 在上述 FISTA-LASSO 代码中加入 Lyapunov 函数 \(E_k\) 的计算（需要知道 \(f^*\)，可用 CVXPY 预先求得），画出 \(E_k\) 随 \(k\) 的变化曲线，验证单调不增。
3. （跨章综合） 结合凸分析章节的 Moreau 分解定理 \(x = \text{prox}_f(x) + \text{prox}_{f^*}(x)\)，证明 \(\text{prox}_{\lambda\|\cdot\|_1}(v) = v - \lambda \cdot \text{prox}_{\frac{1}{\lambda}\|\cdot\|_\infty}(v/\lambda)\)。

---

5. ADMM（交替方向乘子法）⭐⭐

5.1 动机：当 prox 不可行、但问题可"拆分"时

FISTA 要求 \(g\) 有廉价的 prox。但如果问题有如下结构：

\[\min_{x,z} f(x) + g(z) \quad \text{s.t. } Ax + Bz = c\]

其中 \(f\) 和 \(g\) 各自有廉价的 prox（或闭式解），但合在一起没有——此时 ADMM 是最自然的选择。

核心思想：不直接求解耦合问题，而是**交替**处理 \(x\) 和 \(z\)——一个变量固定时解另一个，通过拉格朗日乘子传递耦合信息。

5.2 ADMM 算法推导

增广 Lagrangian： $\(L_\rho(x, z, y) = f(x) + g(z) + y^T(Ax + Bz - c) + \frac{\rho}{2}\|Ax + Bz - c\|^2\)$

ADMM 三步迭代：

\[x^{k+1} = \arg\min_x L_\rho(x, z^k, y^k) \quad \text{（x-update：固定 z, y，关于 x 最小化）}$$ $$z^{k+1} = \arg\min_z L_\rho(x^{k+1}, z, y^k) \quad \text{（z-update：固定 x, y，关于 z 最小化）}$$ $$y^{k+1} = y^k + \rho(Ax^{k+1} + Bz^{k+1} - c) \quad \text{（对偶更新：梯度上升步）}\]

为什么叫"交替方向乘子法"？ - "交替方向"：先优化 \(x\)，再优化 \(z\)（交替） - "乘子法"：对偶变量 \(y\)（Lagrange 乘子）通过梯度上升更新

缩放形式（scaled form，用 \(u = y/\rho\) 替代 \(y\)，更简洁）： $\(x^{k+1} = \arg\min_x \left\{ f(x) + \frac{\rho}{2}\|Ax + Bz^k - c + u^k\|^2 \right\}\)$ $\(z^{k+1} = \arg\min_z \left\{ g(z) + \frac{\rho}{2}\|Ax^{k+1} + Bz - c + u^k\|^2 \right\}\)$ $\(u^{k+1} = u^k + Ax^{k+1} + Bz^{k+1} - c\)$

5.3 收敛性证明

定理 7（ADMM 收敛）。 对闭凸 proper 的 \(f, g\)，若增广 Lagrangian \(L_0\) 有鞍点，则 ADMM 满足： - 原残差 \(\|Ax^k + Bz^k - c\| \to 0\) - 对偶残差 \(\rho\|B(z^k - z^{k-1})\| \to 0\) - 目标值遍历（平均）收敛率 \(O(1/k)\)

证明思路（He-Yuan 2012 的变分不等式方法）：

定义 \(w = (x, z, y)\)，最优性 KKT 系统写成变分不等式 \(\langle F(w^*), w - w^* \rangle \ge 0\)。

关键 Lyapunov 函数： $\(V_k = \frac{1}{\rho}\|y^k - y^*\|^2 + \rho\|Bz^k - Bz^*\|^2\)$

证明 \(V_{k+1} \le V_k - \rho\|r^{k+1}\|^2 - \frac{1}{\rho}\|s^{k+1}\|^2\)（其中 \(r\) 是原残差，\(s\) 是对偶残差）。

由 \(V_k\) 单调下降且下有界，求和得 \(\sum_{k=0}^K (\rho\|r^k\|^2 + \frac{1}{\rho}\|s^k\|^2) \le V_0\)，因此 \(\min_{k \le K} (\|r^k\|^2 + \|s^k\|^2) = O(1/K)\)。

强凸加速：若 \(f\) 或 \(g\) 是 \(\mu\)-强凸的，ADMM 线性收敛（Deng-Yin 2016），率取决于 \(\rho\) 的选择。最优 \(\rho^* = \sqrt{\mu L}\)（Giselsson-Boyd 2017）。

5.4 ADMM 与 Douglas-Rachford Splitting 的等价性

本质洞察：ADMM 作用在原问题 \(\min f(x) + g(z)\) s.t. \(x = z\) 上，等价于 Douglas-Rachford splitting 作用在**对偶问题**上。这个等价性（Eckstein-Bertsekas 1992）解释了为什么 ADMM 不需要原问题可微——它本质上是在对偶空间做近端迭代。

Douglas-Rachford splitting：对 \(\min h_1(u) + h_2(u)\)： $\(\hat{u}^{k+1} = \text{prox}_{\alpha h_2}(u^k)\)$ $\(u^{k+1} = u^k + \text{prox}_{\alpha h_1}(2\hat{u}^{k+1} - u^k) - \hat{u}^{k+1}\)$

当 \(h_1 = f^*(-\cdot)\)，\(h_2 = g^*(\cdot)\)（共轭函数）时，DR 恰好给出 ADMM 的对偶迭代。

5.5 OSQP：ADMM 求解 QP 的工业标准

OSQP（Stellato 等 2020）求解标准 QP：\(\min \frac{1}{2}x^T P x + q^T x\) s.t. \(l \le Ax \le u\)。

ADMM 分裂策略：引入辅助变量 \(z = Ax\)： $\(\min \frac{1}{2}x^T P x + q^T x + I_{[l,u]}(z) \quad \text{s.t. } Ax = z\)$

x-update 归结为解线性系统： $\(\begin{bmatrix} P + \sigma I & A^T \\ A & -I/\rho \end{bmatrix} \begin{bmatrix} x^{k+1} \\ \nu^{k+1} \end{bmatrix} = \begin{bmatrix} \sigma x^k - q \\ z^k - u^k/\rho \end{bmatrix}\)$

OSQP 的三大工程优势：

1. 因子分解缓存：KKT 矩阵左端在 MPC 闭环中**不变**（\(P, A\) 固定），只需一次 LDL 分解，后续每步只做前后代入 \(O(n^2)\)
2. Warm-start：MPC 时序求解中，上一时刻的解是当前时刻的极佳初始点，ADMM 迭代次数从几十次降到 5-10 次
3. 不可行检测：ADMM 原/对偶残差可以自然判断问题是否不可行或无界，无需额外逻辑

// OSQP C API 使用示例（MPC 场景）
#include <osqp/osqp.h>

// 设置问题数据（只做一次）
OSQPSettings settings;
osqp_set_default_settings(&settings);
settings.warm_starting = 1;  // 开启 warm-start
settings.eps_abs = 1e-4;     // MPC 通常不需要太高精度
settings.eps_rel = 1e-4;
settings.max_iter = 200;     // 限制最大迭代以保证实时性

OSQPSolver *solver;
osqp_setup(&solver, P, q, A, l, u, m, n, &settings);

// MPC 闭环中（每个控制周期执行）
osqp_update_data_vec(solver, q_new, l_new, u_new);  // 只更新 q, l, u
osqp_warm_start(solver, x_prev, y_prev);            // 用上次解做 warm-start
osqp_solve(solver);                                  // 求解（通常 5-15 步 ADMM）

5.6 \(\rho\) 参数选择——ADMM 的"阿喀琉斯之踵"

ADMM 的收敛率和 \(\rho\) 的选择**高度相关**： - \(\rho\) 过小：原残差收敛慢（对约束的惩罚弱） - \(\rho\) 过大：对偶残差收敛慢（乘子更新步太大） - 最优 \(\rho\)：需要知道问题的光滑/强凸参数，通常未知

自适应策略（Boyd 2011, §3.4.1）： $\(\rho^{k+1} = \begin{cases} \tau \rho^k & \text{if } \|r^k\| > \mu \|s^k\| \\ \rho^k / \tau & \text{if } \|s^k\| > \mu \|r^k\| \\ \rho^k & \text{otherwise} \end{cases}\)$

OSQP 使用更精细的"residual balancing"策略。但注意：\(\rho\) 变化时 KKT 矩阵改变，需要重新分解——OSQP 通过稀疏更新（rank-1 update）减少代价。

⚠️ 常见陷阱

⚠️ 编程陷阱：忘记 ADMM 的 \(\rho\) 初始化 - 错误做法：\(\rho = 1\)（对所有问题） - 后果：可能需要数千次迭代或完全不收敛 - 正确做法：\(\rho\) 初始值与问题的 \(\|P\|\), \(\|A\|\) 相关；OSQP 默认自动 scaling（Ruiz equilibration）

💡 概念误区：把 ADMM 的 \(O(1/k)\) 理论收敛率当成实际表现 - 理论 \(O(1/k)\) 是遍历（平均）收敛，实际 last-iterate 收敛可能更慢（波动） - 但 warm-start + 自适应 \(\rho\) 让 ADMM 在 MPC 场景下实际非常快（5-20 步） - 理论慢不代表实际慢，实际快的原因是结构利用（warm-start + 因子缓存）

练习

1. （推导） 写出 ADMM 求解 LASSO \(\min \frac{1}{2}\|Ax-b\|^2 + \lambda\|z\|_1\) s.t. \(x = z\) 的完整三步迭代（x-update, z-update, u-update），明确每步的闭式解。
2. （编程） 用 OSQP Python 接口求解一个 100 维盒约束 QP，分别测试 \(\rho = 0.01, 1, 100\) 的迭代次数。
3. （工程判断） 你需要以 50Hz 求解一个 MPC QP（200 变量，500 约束）。OSQP 和 qpOASES（active-set）哪个更合适？给出理由。

---

6. 内点法（IPM）⭐⭐⭐

6.1 动机：当需要高精度与约束保证时

一阶方法（GD、FISTA、ADMM）在中精度（\(\epsilon \sim 10^{-3}\) 到 \(10^{-4}\)）下非常高效，但要达到高精度（\(\epsilon < 10^{-8}\)）时迭代次数急剧增加。而内点法（Interior Point Method, IPM）具有**多项式时间复杂度**：无论精度要求多高，Newton 步数只增长 \(O(\log(1/\epsilon))\)。

本质洞察：内点法把约束优化转化为一系列**无约束**问题（通过 log-barrier），然后用 Newton 法快速求解每个无约束子问题。Newton 法的二次收敛保证了每步能把误差平方缩小——这是 IPM 快的根本原因。

6.2 Log-Barrier 方法

考虑约束凸问题：\(\min f(x)\) s.t. \(f_i(x) \le 0\)，\(i = 1,\ldots,m\)。

Barrier 函数： $\(\phi(x) = -\sum_{i=1}^m \log(-f_i(x))\)$

当 \(x\) 靠近约束边界时 \(\phi(x) \to +\infty\)——像一堵无形的"墙"阻止迭代点越界。

Central path（中心路径）：对参数 \(t > 0\)，定义： $\(x^*(t) = \arg\min_x \left\{ tf(x) + \phi(x) \right\}\)$

当 \(t \to \infty\) 时，\(x^*(t) \to x^*\)（原问题最优解）。中心路径是从内部逐渐逼近最优解的一条曲线。

Barrier 方法（外层循环）： 1. 初始化可行内点 \(x^{(0)}\)，\(t = t_0\) 2. 用 Newton 法求解 \(x^*(t)\)（centering step） 3. 更新 \(t \leftarrow \mu \cdot t\)（\(\mu > 1\)，通常 \(\mu = 10\) 或 \(\mu = 2\)） 4. 重复直到对偶间隙 \(m/t < \epsilon\)

6.3 IPM 复杂度分析

定理 9。 对带 \(m\) 个不等式约束的凸程序，若 barrier \(\phi\) 是 \(\nu\)-self-concordant barrier，则 IPM 在 \(O(\sqrt{\nu}\log(m/\epsilon))\) 次 Newton 步内达到 \(\epsilon\)-最优。

对 log-barrier，\(\nu = m\)，所以复杂度为 \(O(\sqrt{m}\log(1/\epsilon))\)。

与一阶法的对比：

方法	迭代数	每步代价	总代价	适用场景
GD	\(O(L/\epsilon)\)	\(O(n)\)	\(O(Ln/\epsilon)\)	大规模低精度
FISTA	\(O(\sqrt{L/\epsilon})\)	\(O(n)\)	\(O(n\sqrt{L/\epsilon})\)	大规模中精度
IPM	\(O(\sqrt{m}\log(1/\epsilon))\)	\(O(n^3)\)	\(O(n^3\sqrt{m}\log(1/\epsilon))\)	中小规模高精度

分水岭：当 \(n^3\sqrt{m}\log(1/\epsilon) < n \cdot \sqrt{L/\epsilon}\) 时 IPM 更快，即 \(n^2\sqrt{m/L} < 1/\sqrt{\epsilon}\)。对典型 MPC QP（\(n = 200\), \(m = 500\), \(\epsilon = 10^{-6}\)），IPM 约需 30 次 Newton 步。

6.4 Self-Concordance：IPM 多项式时间的数学基石

为什么 Newton 法在 barrier 上总是收敛？ 普通 Newton 法对非凸函数可能发散。但 barrier 函数具有 self-concordance 性质，保证 Newton 步在 Dikin 椭球内的仿射不变性。

定义：函数 \(\phi: \mathbb{R} \to \mathbb{R}\) 是 self-concordant 的，若 \(|\phi'''(x)| \le 2(\phi''(x))^{3/2}\)。

直觉：self-concordance 约束了三阶导数——函数的"弯曲速率"被其"曲率本身"所控制。这意味着 Newton 方向在局部是可信的——二阶近似不会因三阶项而被严重破坏。

验证 \(\phi(x) = -\log(x)\) 是 self-concordant 的： - \(\phi'(x) = -1/x\)，\(\phi''(x) = 1/x^2\)，\(\phi'''(x) = -2/x^3\) - \(|\phi'''(x)| = 2/x^3\)，\(2(\phi''(x))^{3/2} = 2(1/x^2)^{3/2} = 2/x^3\) ✓

6.5 HPIPM：结构化 IPM 与 Riccati 递归

对 MPC 的 QP（OCP 结构），Hessian 呈 block-tridiagonal 形式：

\[\begin{bmatrix} Q_0 + A_0^T \Lambda_1 A_0 & A_0^T \Lambda_1 B_0 & & \\ B_0^T \Lambda_1 A_0 & R_0 + B_0^T \Lambda_1 B_0 & & \\ & & Q_1 + A_1^T \Lambda_2 A_1 & \cdots \\ & & \vdots & \ddots \end{bmatrix}\]

直接 IPM 求解线性系统需要 \(O((nT)^3)\)。但利用 OCP 结构，可以用**Riccati 递归**在 \(O(n^3 T)\) 内求解——快了 \(T^2\) 倍。

这就是 HPIPM（Frison-Diehl 2020）的核心：在 IPM 的每次 Newton 步内用 Riccati 递归代替通用线性求解器。Boston Dynamics Atlas、MIT Mini Cheetah、ANYmal 都使用这条路线。

6.6 IPM 在机器人 MPC 中的应用

机器人系统	QP 规模	求解器	求解时间	频率
MIT Cheetah 3	120 变量, 200 约束	HPIPM	0.2 ms	1 kHz
ANYmal	200 变量, 500 约束	HPIPM	0.5 ms	400 Hz
Atlas (Boston Dynamics)	~500 变量	商业 IPM	< 1 ms	333 Hz
四旋翼 NMPC	600 变量 (SQP+HPIPM)	acados	1 ms	100 Hz

⚠️ 常见陷阱

⚠️ 编程陷阱：IPM 在病态问题上的数值失败 - 中心路径参数 \(t \to \infty\) 时，barrier Hessian 的条件数也 \(\to \infty\) - Newton 系统变得病态，float64 精度不够 - 正确做法：使用 Mehrotra predictor-corrector（标准做法）、或加正则化项

💡 概念误区：把 IPM 和 SQP 混为一谈 - IPM 对 NLP：沿**中心路径**（barrier 的 Newton 步），保持严格可行性 - SQP 对 NLP：在当前点线性化约束，求解 QP 子问题，可能违反约束 - Ipopt 是 IPM（不是 SQP！虽然常被误称）；SNOPT 是 SQP

练习

1. （证明） 验证 \(\phi(x) = -\log(x)\) 满足 self-concordance 条件，并证明 \(\phi(x) = -\sum_{i=1}^m \log(-f_i(x))\) 在线性约束 \(f_i(x) = a_i^T x - b_i\) 下也是 self-concordant 的。
2. （编程） 用 CVXPY 求解一个 50 维 SOCP，对比 ECOS（IPM）和 SCS（ADMM）的求解时间与精度。
3. （工程判断） 对一个 \(n = 10^4\) 变量的 SDP，IPM 的 Newton 系统维度是 \(n^2 \times n^2 = 10^8 \times 10^8\)。为什么这不可行？应该用什么替代方案？

---

7. 其他一阶方法鸟瞰 ⭐⭐⭐

7.1 Frank-Wolfe（条件梯度法）

适用场景：约束集 \(\mathcal{X}\) 上的投影很贵，但**线性最小化** \(\min_{s \in \mathcal{X}} \langle c, s \rangle\) 很便宜。

算法： $\(s_k = \arg\min_{s \in \mathcal{X}} \langle \nabla f(x_k), s \rangle\)$ $\(x_{k+1} = x_k + \gamma_k (s_k - x_k), \quad \gamma_k = \frac{2}{k+2}\)$

收敛率：\(f(x_k) - f^* \le \frac{2LD^2}{k+2}\)，其中 \(D\) 是 \(\mathcal{X}\) 的直径。

独特优势： 1. 迭代解是极点的凸组合 → 自动稀疏 2. 不需要投影 → 对核范数球、单纯形等复杂集合极有用 3. 无需 Lipschitz 常数 \(L\)（通过线搜索自适应）

与机器人学的联系：核范数约束的低秩矩阵补全（如相机标定中的结构-运动恢复）可用 Frank-Wolfe 避免昂贵的 SVD 投影。

7.2 Mirror Descent

动机：梯度下降用欧氏距离衡量"近与远"，但对某些几何结构（如概率单纯形 \(\Delta_n\)），欧氏距离不是最自然的度量。

核心思想：用 Bregman 散度 \(D_\phi(x, y) = \phi(x) - \phi(y) - \langle \nabla\phi(y), x-y \rangle\) 替代 \(\frac{1}{2}\|x-y\|^2\)。

算法： $\(x_{k+1} = \arg\min_{x \in \mathcal{X}} \left\{ \langle \nabla f(x_k), x \rangle + \frac{1}{\eta} D_\phi(x, x_k) \right\}\)$

关键结果：在概率单纯形上，用负熵 \(\phi(x) = \sum x_i \log x_i\) 作为镜像映射时，Mirror Descent 的 regret 从投影梯度的 \(O(\sqrt{n})\) 降到 \(O(\sqrt{\log n})\)——指数级改善。

这正是在线学习中"指数加权"（Multiplicative Weights）算法和强化学习中 softmax 策略梯度的数学来源。

7.3 次梯度法

对完全非光滑的 Lipschitz 凸函数（如 \(\ell_1\) 范数本身），次梯度法是唯一的通用方法：

\[x_{k+1} = x_k - \alpha_k g_k, \quad g_k \in \partial f(x_k)\]

收敛率：\(f(\bar{x}_K) - f^* \le \frac{G \cdot R}{\sqrt{K}}\)（其中 \(G\) 是次梯度的界，\(R = \|x_0 - x^*\|\)）。

重要区别： - 次梯度法**不是**下降方法——\(f(x_{k+1})\) 可能比 \(f(x_k)\) 大 - 步长必须递减（\(\alpha_k = 1/\sqrt{k}\)）；固定步长只能到 neighborhood - \(O(1/\sqrt{k})\) 在黑箱非光滑模型下是最优的（不可改进）

7.4 随机梯度法（SGD）与方差缩减

对有限和问题 \(f(x) = \frac{1}{n}\sum_{i=1}^n f_i(x)\)（机器学习典型结构）：

方法	收敛率（强凸）	每步代价	总代价
GD	\(O(\kappa\log(1/\epsilon))\)	\(O(n)\)	\(O(n\kappa\log(1/\epsilon))\)
SGD	\(O(1/(\mu\epsilon))\)	\(O(1)\)	\(O(1/(\mu\epsilon))\)
SVRG	\(O((n+\kappa)\log(1/\epsilon))\)	\(O(1)\)（+ 偶尔 \(O(n)\)）	\(O((n+\kappa)\log(1/\epsilon))\)
Katyusha	\(O((n+\sqrt{n\kappa})\log(1/\epsilon))\)	\(O(1)\)	理论最优

SVRG（Johnson-Zhang 2013）的核心思想：SGD 慢是因为梯度估计有方差。SVRG 定期计算全梯度 \(\mu = \nabla f(\tilde{x})\)，然后用方差缩减估计器 \(\nabla f_i(x) - \nabla f_i(\tilde{x}) + \mu\) 替代 \(\nabla f_i(x)\)。当 \(x \approx \tilde{x}\) 时，方差趋零。

与 RL 的联系：RL 中 GAE（Generalized Advantage Estimation）的 baseline 结构与 SVRG 同构：\(\hat{g} = g - b + \mathbb{E}[b]\)，都是"减去已知量的估计以降低方差"。

⚠️ 常见陷阱

🧠 思维陷阱：认为 SGD 的 \(O(1/\sqrt{k})\) 比 GD 的 \(O(1/k)\) 慢 - 比较的基准不同！SGD 的每步只用**一个**样本（\(O(1)\) 代价），GD 每步用**全部** \(n\) 个样本（\(O(n)\) 代价） - 以**总梯度评估次数**为基准：SGD 到 \(\epsilon\) 需要 \(O(1/\epsilon^2)\) 次单样本梯度；GD 需要 \(O(n/\epsilon)\) 次单样本梯度 - 当 \(n\) 很大时 SGD 明显更快（\(1/\epsilon^2 < n/\epsilon\) 当 \(n > 1/\epsilon\)）

练习

1. （概念辨析） Frank-Wolfe 的 \(O(1/k)\) 收敛率看起来和 GD 一样。为什么人们还用 FW？给出至少两个 FW 优于投影梯度法的具体场景。
2. （推导） 证明：在概率单纯形 \(\Delta_n\) 上，使用负熵镜像映射的 Mirror Descent 更新等价于 \(x_{k+1,i} \propto x_{k,i} \exp(-\eta g_{k,i})\)（指数加权更新）。

---

8. 收敛率下界与 PEP ⭐⭐⭐⭐

8.1 Nesterov 下界定理

定理 5（Nesterov 下界）。 存在 \(L\)-光滑凸函数 \(f\)（无穷维"最坏二次"），使得**任何**仅使用 \(\{\nabla f(x_0), \ldots, \nabla f(x_{k-1})\}\) 线性组合生成迭代的黑箱一阶方法满足：

\[f(x_k) - f^* \ge \frac{3L\|x_0 - x^*\|^2}{32(k+1)^2}\]

核心证明技巧：构造"最坏"目标函数为三对角二次形式 \(f(x) = \frac{L}{4}(x_1^2 + \sum_{i=1}^{n-1}(x_i - x_{i+1})^2)\)。这个函数有如下性质：前 \(k\) 步迭代只能探索 \(\text{span}(e_1, \ldots, e_k)\)（因为梯度的 Krylov 子空间性质），而最优解在 \((1,1,\ldots,1)^T\) 方向——需要 \(n\) 步才能到达。

为什么 CG 不违反下界？ 共轭梯度法在 \(n\) 维二次上精确收敛于 \(n\) 步。这似乎"违反"了 \(O(1/k^2)\) 下界。答案：CG 利用了**二次结构**（它知道问题是二次的），而 Nesterov 下界针对的是"黑箱一阶方法"（不利用任何特殊结构）。下界只对不知道目标函数解析形式的方法有效。

8.2 PEP 框架

Performance Estimation Problem（Drori-Teboulle 2014）：把"对一阶方法求最坏收敛率"写成一个 SDP：

\[\max_{f \in \mathcal{F}} (f(x_K) - f^*) \quad \text{s.t. } x_{k+1} = \text{algorithm}(x_0, \ldots, x_k, \nabla f(x_0), \ldots, \nabla f(x_k))\]

其中 \(\mathcal{F}\) 是函数类（如 \(L\)-光滑凸）。

PEPit 工具：Python 包 PEPit 让你用几行代码数值验证任意一阶方法的最坏性能：

from PEPit import PEP
from PEPit.functions import SmoothConvexFunction

# 验证 GD 在 L-光滑凸上 10 步后的最坏性能
problem = PEP()
f = problem.declare_function(SmoothConvexFunction, L=1.0)
x_0 = problem.set_initial_point()
x_star = f.stationary_point()

# 10 步 GD
x = x_0
for _ in range(10):
    x = x - 1.0 * f.gradient(x)  # 步长 1/L = 1

# 设置目标：最大化 f(x_10) - f(x*)
problem.set_objective(f(x) - f(x_star))
problem.set_initial_condition((x_0 - x_star)**2 <= 1)

# 求解 SDP
pepit_value = problem.solve()
# 对比理论上界 L*R^2/(2*(k+1)) = 1/(22) ≈ 0.0455
print(f"PEP 最坏值: {pepit_value:.6f}")
print(f"理论上界:   {1/(2*11):.6f}")

练习

1. （概念） 为什么"黑箱一阶方法"的定义排除了 CG 和 Newton 法？给出精确的"黑箱"定义。
2. （编程） 用 PEPit 验证 Nesterov AGD 在 \(k=10\) 步后的最坏值，并与理论上界 \(\frac{2LR^2}{(k+2)^2}\) 对比。

---

9. 算法选择决策树与机器人应用 ⭐⭐

9.1 算法选择的三条核心准则

1. 精度 × 规模 的乘积 决定 IPM 还是一阶法
2. 是否可 warm-start 决定 ADMM/FISTA 能否发挥优势
3. 结构红利（稀疏、可分、OCP、有限和）永远优于通用算法

9.2 完整决策流程图

输入：凸优化问题
├── Q1: 有约束吗？
│   ├── 是 → Q2（规模 × 精度）
│   │   ├── n,m <= 1e3 且 eps <= 1e-8 → IPM (ECOS/Mosek)
│   │   ├── MPC 结构（block-tridiagonal）→ HPIPM (Riccati IPM)
│   │   ├── n,m in [1e3, 1e5] QP 且 eps ~ 1e-4 → OSQP (ADMM)
│   │   ├── n,m >= 1e5 → SCS (ADMM on cones)
│   │   └── 投影贵（核范数球等）→ Frank-Wolfe
│   └── 否 → Q3（结构）
│       ├── 纯光滑 → Nesterov AGD (+ restart)
│       ├── f + g (g 有 prox) → FISTA
│       ├── 两块非光滑 → ADMM / Douglas-Rachford
│       ├── 非光滑 Lipschitz → 次梯度法
│       └── 有限和 (ML) → SVRG/SAGA/Adam
└── [机器人专用]
    ├── 线性 MPC → OSQP (鲁棒) 或 HPIPM (极速)
    ├── 非线性 MPC → acados (SQP + HPIPM)
    ├── SLAM BA → Ceres (Gauss-Newton + 稀疏 Cholesky)
    ├── RL 策略 → PPO/SAC (SGD with Adam)
    └── 可微优化层 → cvxpylayers / JAXopt

9.3 三大机器人场景深度剖析

场景 A：线性 MPC（四足机器人平衡控制）

问题结构：\(\min \sum_{t=0}^{T-1} (x_t^T Q x_t + u_t^T R u_t)\) s.t. \(x_{t+1} = Ax_t + Bu_t\)，\(u_{\min} \le u_t \le u_{\max}\)。

这是一个带盒约束的 QP，Hessian 呈 block-diagonal + 耦合项的带状结构。

求解策略	优势	劣势	典型求解时间（\(n_x=12, n_u=4, T=20\)）
OSQP (ADMM)	鲁棒、不可行检测、warm-start	精度有限(\(10^{-4}\))	0.5 ms
HPIPM (Riccati IPM)	极快、高精度	设置复杂	0.1 ms
qpOASES (active-set)	精确解、小规模快	大规模退化	0.3 ms

场景 B：SLAM 后端（Bundle Adjustment）

问题结构：\(\min \sum_{(i,j)} \rho(\|z_{ij} - h(x_i, l_j)\|^2_{\Sigma_{ij}})\)，其中 \(x_i \in SE(3)\), \(l_j \in \mathbb{R}^3\)。

这是一个**非凸**最小二乘问题（旋转流形上），严格意义上不属于凸优化。但其核心求解步骤——Gauss-Newton/LM 的线性系统——利用了与凸 QP 类似的稀疏结构。

关键区别：SLAM 不能直接用凸优化算法（问题非凸），但凸优化理论仍然重要： - SE-Sync 把位姿图优化**松弛**为 SDP（凸），可给出全局最优性证明 - 局部搜索的 Gauss-Newton 步等价于在每步的二次近似（凸）上做一次精确求解

场景 C：RL 策略优化

问题结构：\(\min_\theta -J(\theta) = -\mathbb{E}_{\tau \sim \pi_\theta}[R(\tau)]\)（非凸！）

虽然是非凸，但凸优化的思想大量出现： - TRPO：KL 约束下的自然梯度 ≈ Fisher 信息矩阵预条件 ≈ Mirror Descent - PPO clip：近似 trust-region ≈ proximal 方法 - LP-based RL：约束 MDP 在占用度量空间写成 LP，对偶给出 Lagrange 乘子

⚠️ 常见陷阱

🧠 思维陷阱：认为"凸优化理论对非凸问题没用" - 非凸优化的局部搜索方法（GN、LM、SQP）每步都在求解一个**凸**子问题 - 凸松弛（SDP relaxation）可以给非凸问题提供全局最优性证明 - 凸优化的收敛率分析框架可以推广到"收敛到稳定点"（Ghadimi-Lan 2013）

练习

1. （工程选型） 你需要为 10kg 足式机器人做 50Hz 线性 MPC，200 变量 500 约束，精度 \(10^{-4}\)。应选 OSQP、HPIPM、qpOASES、Mosek 中哪个？给出理由。再给一个场景让你换另一款。
2. （跨章综合） 结合对偶理论章节的 KKT 条件，解释为什么 OSQP 的"原残差+对偶残差 → 0"等价于 KKT 条件被满足。
3. （研究级） 阅读 SE-Sync 论文摘要，解释它如何把非凸 SLAM 松弛为 SDP，以及"zero duality gap"的条件是什么。

---

10. 收敛率速查表 ⭐⭐

10.1 完整速查表

方法 \ 函数类	Lipschitz 凸	光滑凸(\(L\))	光滑强凸(\(L,\mu\))	复合 \(f+g\)	带约束凸(\(m\) 约束)
次梯度法	\(O(1/\epsilon^2)\)	—	—	—	—
GD	—	\(O(L/\epsilon)\)	\(O(\kappa\log(1/\epsilon))\)	—	—
Nesterov AGD	—	\(O(\sqrt{L/\epsilon})\)	\(O(\sqrt{\kappa}\log(1/\epsilon))\)	—	—
ISTA	—	—	—	\(O(L/\epsilon)\)	—
FISTA	—	—	—	\(O(\sqrt{L/\epsilon})\)	—
ADMM	\(O(1/\epsilon)\)	\(O(1/\epsilon)\)	\(O(\kappa\log(1/\epsilon))\)	\(O(1/\epsilon)\)	\(O(1/\epsilon)\)
Frank-Wolfe	—	\(O(L D^2/\epsilon)\)	\(O(L D^2/\epsilon)\)	—	紧集
Mirror Descent	\(O(1/\epsilon^2)\)	\(O(L/\epsilon)\)	—	—	—
IPM	—	—	—	—	\(O(\sqrt{m}\log(1/\epsilon))\)

10.2 下界（不可再改进）

• Lipschitz 凸 + 一阶方法：\(\Omega(1/\epsilon^2)\)
• 光滑凸 + 一阶方法：\(\Omega(\sqrt{L/\epsilon})\)
• 光滑强凸 + 一阶方法：\(\Omega(\sqrt{\kappa}\log(1/\epsilon))\)

---

本章常见误解汇总

误解	正确理解
"越新的算法越好"	算法没有"更好"之说，只有"更适合"。FISTA 对 \(f+g\)（\(g\) 有廉价 prox）结构最优；ADMM 对两块可分结构最优；IPM 对中小规模高精度约束问题最优
"收敛率 \(O(1/k^2)\) 总比 \(O(1/k)\) 好"	收敛率是渐近行为，实际性能还取决于常数、warm-start、因子缓存。MPC 场景下 ADMM 的 \(O(1/k)\) 加 warm-start 可能比 FISTA 的 \(O(1/k^2)\) 实际更快
"迭代复杂度低 = 总计算时间短"	总时间 = 迭代数 \(\times\) 每步代价。FISTA 每步 \(O(n)\)，IPM 每步 \(O(n^3)\)。\(n=100\) 的 QP 上 IPM 可能更快
"大问题一定用一阶法"	稀疏大规模 QP（如 SLAM 的 Hessian）中，稀疏 Cholesky（二阶）可能比 ADMM（一阶）更快。关键是稀疏性改变了复杂度
"Nesterov 加速 = Polyak 动量"	Heavy-ball 在二次函数上完美，但在一般凸函数上没有 \(O(1/k^2)\) 保证。两者的关键区别在于梯度计算点（当前点 vs 外推点）和参数（固定 vs 时变）
"SGD 的 \(O(1/\sqrt{k})\) 比 GD 的 \(O(1/k)\) 慢"	比较基准不同：SGD 每步 \(O(1)\)，GD 每步 \(O(n)\)。以总梯度评估次数为基准，当 \(n\) 很大时 SGD 明显更快
"\(O(1/k^2)\) 是不可改进的理论极限"	这是黑箱一阶 oracle 下的极限。如果知道更多结构（强凸、有限和、二次），可以利用结构做得更好
"凸优化理论对非凸问题没用"	非凸优化的局部搜索（GN、SQP）每步解凸子问题；凸松弛（SDP）给非凸问题提供全局最优性证明

---

本章小结

符号表

符号	含义	首次出现
\(L\)	梯度 Lipschitz 常数（光滑性参数）	\(\S\)2
\(\mu\)	强凸参数	\(\S\)2
\(\kappa = L/\mu\)	条件数	\(\S\)2
\(\alpha\), \(\eta\)	步长（学习率）	\(\S\)2
\(\beta_k\)	Nesterov 动量参数	\(\S\)3
\(\text{prox}_f(v)\)	\(f\) 的 proximal 算子	\(\S\)4
\(\rho\)	ADMM 惩罚参数	\(\S\)5
\(\phi(x)\)	log-barrier 函数	\(\S\)6
\(\nu\)	self-concordant barrier 参数	\(\S\)6
\(D\)	约束集直径	\(\S\)7 (FW)
\(D_\phi(x,y)\)	Bregman 散度	\(\S\)7 (Mirror Descent)

方法	核心思想	适用结构	收敛率	机器人应用
GD	沿梯度反方向	光滑凸	\(O(1/k)\) / 线性	所有方法的基础
Nesterov AGD	动量 + lookahead	光滑凸	\(O(1/k^2)\)（最优）	RL optimizer
FISTA	AGD + prox	\(f\)(光滑)\(+g\)(prox)	\(O(1/k^2)\)	稀疏恢复、去噪
ADMM	分裂 + 对偶上升	可分结构	\(O(1/k)\)	MPC(OSQP)、分布式
IPM	log-barrier + Newton	约束凸，中小规模	\(O(\sqrt{m}\log(1/\epsilon))\)	MPC(HPIPM)、高精度
FW	线性最小化 oracle	投影贵的紧集	\(O(1/k)\)	核范数优化
Mirror Descent	Bregman 散度	特殊几何	对几何自适应	在线学习、softmax RL

---

累积项目：本章新增模块

项目：手写优化器库

本章新增： - gradient_descent.py：带 Armijo 回溯的 GD - nesterov_agd.py：带自适应重启的 AGD - fista.py：FISTA for LASSO - admm_qp.py：ADMM 求解盒约束 QP - benchmark.py：四种方法在同一 QP 上的收敛曲线对比

---

延伸阅读

资源	难度	主题
Boyd & Vandenberghe《Convex Optimization》Ch 9-11	⭐⭐	GD + Newton + IPM
Nesterov《Lectures on Convex Optimization》2nd ed.	⭐⭐⭐	加速法理论最权威
Beck《First-Order Methods in Optimization》	⭐⭐	PGM/FISTA 最详尽
Boyd et al. ADMM 综述 2011	⭐⭐	ADMM 圣经
Bubeck 综述 arXiv:1405.4980	⭐⭐⭐	120 页浓缩全景
OSQP 论文 (Stellato 2020)	⭐⭐	ADMM for QP 工程
HPIPM 论文 (Frison 2020)	⭐⭐⭐	Riccati IPM
Distill "Why Momentum Really Works"	⭐	动量可视化
PEPit 文档	⭐⭐⭐⭐	自动验证收敛率

---

🔧 故障排查手册

症状	可能原因	排查步骤	相关章节
GD 收敛极慢（\(>10^5\) 步）	条件数大/\(L\) 估计错误	1. 计算 Hessian 特征值范围 2. 检查步长是否 \(\le 1/L\) 3. 改用 AGD	§2
AGD 振荡不收敛	动量累积过大	1. 检查 \(\beta_k\) 公式 2. 加入自适应重启 3. 检查 \(L\) 是否准确	§3
FISTA 收敛后又上升	非单调行为（正常）/步长错	1. 检查 \(1/L\) 步长 2. 监控 Lyapunov \(E_k\) 3. 用重启	§4
ADMM 不收敛（残差不降）	\(\rho\) 选择不当	1. 打印原/对偶残差 2. 开启自适应 \(\rho\) 3. 检查问题是否真的可行	§5
IPM 输出 NaN	初始点不可行/Hessian 奇异	1. 确保初始点严格可行 2. 加正则化 3. 检查约束是否一致	§6
OSQP warm-start 无效	问题结构变化	1. 检查 P, A 是否改变 2. 只更新 q, l, u 3. 确认 warm_start=True	§5.5
FW 收敛极慢	步长或 LMO 有误	1. 检查 \(\gamma_k = 2/(k+2)\) 2. 验证线性最小化 oracle 3. 检查集合直径 \(D\)	§7.1
Mirror Descent 输出非概率	镜像映射不匹配	1. 检查 \(\phi\) 选择 2. 验证指数归一化 3. 确认定义域约束	§7.2
SGD 方差过大	未用 variance reduction	1. 切换 SVRG/SAGA 2. 减小学习率 3. 增大 mini-batch	§7.4

---

11. 梯度下降的深层分析与变体 ⭐⭐

11.1 精确步长选择与 Wolfe 条件

Armijo 回溯（§2.5）是实践中最常用的策略，但理论上还有更精细的步长选择方法。

精确线搜索：\(\alpha_k = \arg\min_{\alpha \ge 0} f(x_k - \alpha \nabla f(x_k))\)。对二次函数 \(f(x) = \frac{1}{2}x^\top Ax + b^\top x\)，精确线搜索有闭式解：

\[\alpha_k = \frac{\|\nabla f(x_k)\|^2}{\nabla f(x_k)^\top A \nabla f(x_k)}\]

推导：设 \(g_k = \nabla f(x_k) = Ax_k + b\)。沿方向 \(d_k = -g_k\) 走步长 \(\alpha\)：

\[\phi(\alpha) = f(x_k + \alpha d_k) = \frac{1}{2}(x_k + \alpha d_k)^\top A (x_k + \alpha d_k) + b^\top(x_k + \alpha d_k)\]

对 \(\alpha\) 求导令零：\(d_k^\top A(x_k + \alpha d_k) + b^\top d_k = 0\)，即 \(d_k^\top(g_k + \alpha A d_k) = 0\)。

由 \(d_k = -g_k\)：\(-g_k^\top g_k + \alpha g_k^\top A g_k = 0\)，解得 \(\alpha_k = \frac{g_k^\top g_k}{g_k^\top A g_k}\)。

Wolfe 条件（比 Armijo 更强，用于拟 Newton 方法）：

\[f(x_k + \alpha d_k) \le f(x_k) + c_1 \alpha \nabla f(x_k)^\top d_k \quad \text{（充分下降，同 Armijo）}$$ $$\nabla f(x_k + \alpha d_k)^\top d_k \ge c_2 \nabla f(x_k)^\top d_k \quad \text{（曲率条件，}0 < c_1 < c_2 < 1\text{）}\]

曲率条件排除了步长过小的情况——它要求新点处的方向导数不能比初始处小太多（即函数已经"足够平缓"了）。

为什么 L-BFGS 需要 Wolfe？ Armijo 只保证充分下降，但 L-BFGS 还需要曲率信息来维护 Hessian 近似的正定性。满足 strong Wolfe 条件的步长保证了 \(s_k^\top y_k > 0\)（\(s_k = x_{k+1} - x_k\)，\(y_k = g_{k+1} - g_k\)），这是 BFGS 更新公式正定性的前提。

11.2 坐标下降 ⭐⭐

动机：当变量维度很高（\(n \gg 10^4\)）且函数关于各坐标可分或有特殊结构时，坐标下降（Coordinate Descent, CD）比全梯度下降更快。

算法：每步只更新一个坐标 \(i\)： $\(x_{k+1}^{(i)} = \arg\min_{t} f(x_1, \ldots, x_{i-1}, t, x_{i+1}, \ldots, x_n)\)$

适用条件：\(f\) 可分（或近似可分），每个坐标子问题有闭式解。

**LASSO 的坐标下降**是最经典的应用。对 \(\min \frac{1}{2}\|Ax-b\|^2 + \lambda\|x\|_1\)，第 \(i\) 个坐标的更新：

\[x_i \leftarrow \text{soft}_{\lambda/\|a_i\|^2}\left(\frac{a_i^\top(b - A_{-i}x_{-i})}{\|a_i\|^2}\right)\]

其中 \(a_i\) 是 \(A\) 的第 \(i\) 列，\(A_{-i}\)、\(x_{-i}\) 是去掉第 \(i\) 个坐标后的矩阵和向量。

收敛率：对光滑凸函数，CD 的收敛率为 \(O(n/k)\)——比 GD 的 \(O(1/k)\) 多一个 \(n\) 因子。但每步代价只有 \(O(1)\)（vs GD 的 \(O(n)\)），所以**总计算量相当**。CD 的优势在于避免了存储和计算全梯度。

反事实推理：如果不用坐标下降而用全梯度下降求解 \(n = 10^6\) 的 LASSO，每步需要 \(O(mn)\) 的矩阵-向量乘，而 CD 每步只需 \(O(m)\)。在 \(n\) 很大时，CD 的内存和计算优势是决定性的——这就是 scikit-learn、glmnet 等库的默认 LASSO 求解器使用 CD 的原因。

11.3 投影梯度法的精确收敛分析 ⭐⭐

投影梯度法（Projected Gradient Descent, PGD）是约束凸优化中最基本的方法。回顾凸分析基础（§1.15）中的投影算子 \(\Pi_C\)，PGD 每步迭代为：

\[x_{k+1} = \Pi_C\left(x_k - \frac{1}{L}\nabla f(x_k)\right)\]

定理 11.1（PGD 收敛率）：设 \(f\) 在 \(C\) 上 \(L\)-光滑凸，则 PGD 满足：

\[f(x_k) - f^* \le \frac{L\|x_0 - x^*\|^2}{2k}\]

证明关键步骤：利用投影的非扩张性 \(\|\Pi_C(u) - x^*\| \le \|u - x^*\|\)（因为 \(x^* \in C\)），descent lemma 的证明可以逐步推广。

具体地，令 \(y_k = x_k - \frac{1}{L}\nabla f(x_k)\)，\(x_{k+1} = \Pi_C(y_k)\)。

Step 1：由 descent lemma：\(f(x_{k+1}) \le f(x_k) + \nabla f(x_k)^\top(x_{k+1} - x_k) + \frac{L}{2}\|x_{k+1} - x_k\|^2\)。

Step 2：由 \(x_{k+1} = \Pi_C(y_k)\) 的最优性：\(\langle y_k - x_{k+1}, x - x_{k+1} \rangle \le 0\) 对所有 \(x \in C\)。

即 \(\langle x_k - \frac{1}{L}\nabla f(x_k) - x_{k+1}, x - x_{k+1} \rangle \le 0\)，展开后可得到与无约束情形类似的递推关系。

Step 3：经过代数化简（与 §2.3 的 GD 证明相同的技巧），得到 \(f(x_k) - f^* \le \frac{L\|x_0-x^*\|^2}{2k}\)。

投影代价的讨论：PGD 的总代价 = \(k\) 次梯度计算 + \(k\) 次投影。对简单集合（Box、\(\ell_2\) 球），投影是 \(O(n)\)；对复杂集合（半正定锥 \(\mathbb{S}^n_+\)），投影是 \(O(n^3)\)（需要特征分解）。当投影代价远大于梯度代价时，Frank-Wolfe（§7.1）是更好的选择。

⚠️ 常见陷阱

⚠️ 编程陷阱：坐标下降的循环方式 - 错误做法：按固定顺序 \(i = 1, 2, \ldots, n\) 循环更新 - 后果：对某些问题收敛极慢（坐标之间强耦合时） - 正确做法：随机坐标选择（Randomized CD）有更好的理论保证。也可以用 Gauss-Southwell 规则（选梯度分量最大的坐标优先更新）

💡 概念误区：认为投影梯度法只是"梯度下降+投影"的简单拼凑 - 新手想法："先梯度步，再投影，两步独立" - 实际上：投影梯度法等价于 \(\text{prox}_{\delta_C}(x_k - \eta\nabla f(x_k))\)——它是 proximal gradient 的特殊情形（当 \(g = \delta_C\) 时）。理解了这个等价性，就自然理解了为什么 PGD 和 FISTA 有相同的收敛率框架

练习

1. 推导二次函数上精确线搜索的最速下降法的收敛率，证明 \(f(x_k) - f^* \le \left(\frac{\kappa-1}{\kappa+1}\right)^{2k}(f(x_0) - f^*)\)。
2. 对 elastic net \(\min \frac{1}{2}\|Ax-b\|^2 + \lambda_1\|x\|_1 + \frac{\lambda_2}{2}\|x\|^2\)，写出坐标下降的第 \(i\) 个坐标更新公式。
3. 对半正定锥约束 \(\min f(X)\) s.t. \(X \succeq 0\)，实现投影梯度法并分析每步复杂度。

---

12. 内点法的完整理论 ⭐⭐⭐

12.1 Self-Concordance 的完整理论

回顾 §6.4 中引入的 self-concordance 概念。这里给出完整的多维推广和关键定理。

定义 12.1（多维 self-concordant 函数）：\(f: \mathbb{R}^n \to \mathbb{R}\) 是 self-concordant 的，若对所有 \(x \in \text{dom}(f)\) 和所有 \(h \in \mathbb{R}^n\)：

\[|D^3 f(x)[h, h, h]| \le 2 (h^\top \nabla^2 f(x) h)^{3/2}\]

其中 \(D^3 f(x)[h, h, h] = \frac{d^3}{dt^3} f(x + th)\big|_{t=0}\)。

直觉：self-concordance 要求函数的三阶变化率被二阶曲率所"控制"。这意味着二阶 Taylor 近似在 Newton 步范围内是可靠的——Newton 方向不会被三阶项严重破坏。

self-concordant barrier 的定义：\(\phi\) 是 \(\nu\)-self-concordant barrier，若 \(\phi\) 是 self-concordant 的，且满足：

\[\|\nabla\phi(x)\|_{x}^* \le \sqrt{\nu}\]

其中 \(\|g\|_x^* = \sqrt{g^\top [\nabla^2\phi(x)]^{-1} g}\) 是 \(g\) 在 Hessian 度量下的对偶范数，\(\nu\) 称为 barrier 参数。

为什么 \(\nu\) 重要：barrier 参数 \(\nu\) 直接决定了 IPM 的迭代次数。标准 log-barrier \(\phi(x) = -\sum_{i=1}^m \log(-f_i(x))\) 的 barrier 参数为 \(\nu = m\)（约束数）。Nesterov-Nemirovsky (1994) 的核心定理是：

定理 12.2（IPM 复杂度）：对 \(\nu\)-self-concordant barrier，path-following IPM 在 \(O(\sqrt{\nu}\log(\nu/\epsilon))\) 次 Newton 步内达到 \(\epsilon\)-最优。

各锥的最优 barrier 参数：

锥	标准 barrier	\(\nu\)	最优 barrier	\(\nu_{\text{opt}}\)
\(\mathbb{R}^n_+\)	\(-\sum\log x_i\)	\(n\)	同上	\(n\)
\(\text{SOC}_n\)	\(-\log(t^2 - \|x\|^2)\)	2	同上	2
\(\mathbb{S}^n_+\)	\(-\log\det X\)	\(n\)	同上	\(n\)

SOC 的 \(\nu = 2\) 与维度无关——这就是为什么 SOCP 比 LP 的复杂度更好（对大维度约束）。\(\mathbb{S}^n_+\) 的 \(\nu = n\) 是最优的（Nesterov-Nemirovsky 1994），不能更低。

12.2 Newton 步的局部二次收敛

定理 12.3（Newton 步的 self-concordant 收敛）：设 \(\phi\) 是 self-concordant 的，\(x\) 满足 \(\lambda(x) = \|\nabla\phi(x)\|_x^* < 1/4\)，则 Newton 步 \(x^+ = x - [\nabla^2\phi(x)]^{-1}\nabla\phi(x)\) 满足：

\[\lambda(x^+) \le \frac{2\lambda(x)^2}{(1 - \lambda(x))^2}\]

二次收敛的含义：当 \(\lambda(x) < 0.25\) 时，每步 Newton 把 \(\lambda\) 大约平方缩小——从 0.1 到 0.01 到 0.0001。这就是 IPM 的"内层循环"如此快的原因：通常 4-8 步 Newton 就从粗糙解达到机器精度。

proof sketch：利用 self-concordance 条件，可以证明 Newton decrement \(\lambda(x)\) 在 Newton 步后满足上述递推。关键是三阶项被 self-concordance 约束，使得二阶 Taylor 近似的误差可控。

12.3 Barrier Method 的外层复杂度

完整 barrier method 算法：

输入：self-concordant barrier phi，初始可行内点 x0，参数 mu > 1
1. 初始化 t = t0 > 0
2. while m/t > epsilon:
    a. Centering：用 Newton 法求解 x*(t) = argmin[t*f0(x) + phi(x)]
    b. 更新：t <- mu*t
3. 返回 x*(t)

外层迭代次数：从 \(t_0\) 到 \(t_f = m/\epsilon\) 需要 \(\lceil\log(m/(t_0\epsilon))/\log\mu\rceil\) 次。

选择 \(\mu = 1 + 1/\sqrt{\nu}\)（理论最优），得外层 \(O(\sqrt{\nu}\log(m/\epsilon))\) 次，每次内层 \(O(1)\) 次 Newton 步（warm-start 后），总 Newton 步 \(O(\sqrt{\nu}\log(m/\epsilon))\)。

反事实推理：如果 barrier 不满足 self-concordance 怎么办？Newton 法的局部收敛域依赖于函数的三阶行为——如果三阶导数不受控，Newton 步可能过大导致发散，或者需要 line search 回退，丧失二次收敛。Self-concordance 保证了 Newton 步的**仿射不变**收敛分析——不依赖具体坐标系或 Lipschitz 常数。

12.4 Mehrotra Predictor-Corrector

实际 IPM 求解器（ECOS、Mosek、HPIPM）几乎都使用 Mehrotra predictor-corrector 策略而非简单的 path-following：

Step 1（Predictor）：在当前互补对 \((x, s)\)（\(s\) 是对偶松弛变量）处计算仿射 Newton 方向 \((\Delta x^{\text{aff}}, \Delta s^{\text{aff}})\)，相当于 \(\mu = 0\) 的 Newton 步。

Step 2（Centering parameter）：计算仿射步后的互补 \(\mu^{\text{aff}} = (x + \Delta x^{\text{aff}})^\top(s + \Delta s^{\text{aff}})/m\)，令 \(\sigma = (\mu^{\text{aff}}/\mu)^3\)（经验值）。

Step 3（Corrector）：求解修正后的 Newton 系统，目标是 \(XS\mathbf{1} = \sigma\mu\mathbf{1}\)（而非纯 centering 的 \(\mu\mathbf{1}\)）。

为什么 predictor-corrector 更快：predictor 步利用了"如果约束都是 inactive 的话最优方向是什么"的信息，corrector 步修正了互补条件。实践中 predictor-corrector 的外层迭代次数通常比标准 path-following 少 2-3 倍。

⚠️ 常见陷阱

💡 概念误区：认为 self-concordance 只是数学美感 - self-concordance 不是"锦上添花"，它是 IPM 多项式时间保证的**唯一已知手段** - 没有 self-concordance，Newton 法的收敛分析依赖 Lipschitz 常数——而 barrier 函数在约束边界附近 Lipschitz 常数趋于无穷 - self-concordance 提供了"仿射不变"的收敛保证——这就是为什么 IPM 对约束的 scaling 不敏感

🧠 思维陷阱：认为 IPM 迭代次数与问题维度 \(n\) 无关 - IPM 的 Newton 步数 \(O(\sqrt{m}\log(1/\epsilon))\) 确实与 \(n\) 无关 - 但每步 Newton 需要解一个 \(n \times n\) 线性系统，代价 \(O(n^3)\)（一般）或 \(O(n \cdot \text{nnz})\)（稀疏） - **总计算量**还是依赖 \(n\) 的——只是通过 Newton 步数和每步代价的乘积体现

练习

1. 验证 \(\phi(x) = -\sum_{i=1}^m \log(b_i - a_i^\top x)\) 对线性约束是 self-concordant 的，并计算其 barrier 参数 \(\nu\)。
2. 对一个 \(m = 100\) 约束的 LP，估计 barrier method 需要多少次 Newton 步达到 \(\epsilon = 10^{-8}\)（取 \(\mu = 1 + 1/\sqrt{m}\)）。
3. 解释为什么 SOCP 的 IPM 比 LP 的 IPM 对维度更友好（提示：SOC barrier 的 \(\nu\) 与约束维度无关）。

---

13. Frank-Wolfe 完整理论与应用 ⭐⭐⭐

13.1 从投影到线性最小化

回顾 §7.1 中引入的 Frank-Wolfe（条件梯度法）。这里深入分析其收敛性和独特优势。

为什么需要 Frank-Wolfe？ 考虑核范数约束优化 \(\min f(X)\) s.t. \(\|X\|_* \le \tau\)。投影到核范数球需要完整 SVD（\(O(\min(m,n) \cdot mn)\)），而 Frank-Wolfe 的线性最小化 oracle 只需要最大奇异值对（可用幂迭代 \(O(mn)\) 或更快的随机方法），代价低得多。

算法详述：

\[s_k = \arg\min_{s \in \mathcal{X}} \langle \nabla f(x_k), s \rangle \quad \text{（线性最小化 oracle, LMO）}$$ $$x_{k+1} = x_k + \gamma_k (s_k - x_k) \quad \text{（步长 } \gamma_k = \frac{2}{k+2} \text{ 或线搜索）}\]

关键洞察：\(s_k\) 是约束集的极值点（线性函数在紧凸集上的最优解在极值点取到）。因此 \(x_k\) 是 \(k\) 个极值点的凸组合——迭代解自动稀疏。

13.2 收敛率证明

定理 13.1（Frank-Wolfe \(O(1/k)\)）：设 \(f\) 凸且 \(L\)-光滑，\(\mathcal{X}\) 紧凸且直径 \(D = \max_{x,y \in \mathcal{X}}\|x-y\|\)。则：

\[f(x_k) - f^* \le \frac{2LD^2}{k+2}\]

证明：

Step 1：由 descent lemma： $\(f(x_{k+1}) \le f(x_k) + \gamma_k \langle \nabla f(x_k), s_k - x_k \rangle + \frac{L\gamma_k^2}{2}\|s_k - x_k\|^2\)$

Step 2：由 \(s_k\) 是 LMO 的解：\(\langle \nabla f(x_k), s_k \rangle \le \langle \nabla f(x_k), x^* \rangle\)（因为 \(x^* \in \mathcal{X}\)）。

因此 \(\langle \nabla f(x_k), s_k - x_k \rangle \le \langle \nabla f(x_k), x^* - x_k \rangle \le -(f(x_k) - f^*)\)。

Step 3：\(\|s_k - x_k\| \le D\)。代入 Step 1：

\[f(x_{k+1}) - f^* \le (1 - \gamma_k)(f(x_k) - f^*) + \frac{L\gamma_k^2 D^2}{2}\]

Step 4：取 \(\gamma_k = 2/(k+2)\)，用归纳法证明 \(f(x_k) - f^* \le \frac{2LD^2}{k+2}\)。

归纳基础：\(k=0\) 时 \(f(x_0) - f^* \le LD^2/2\)（由 descent lemma）。

归纳步骤：设 \(f(x_k) - f^* \le \frac{2LD^2}{k+2}\)。则：

\[f(x_{k+1}) - f^* \le \frac{k}{k+2} \cdot \frac{2LD^2}{k+2} + \frac{2LD^2}{(k+2)^2} = \frac{2LD^2(k+1)}{(k+2)^2} \le \frac{2LD^2}{k+3} \quad \blacksquare\]

13.3 Frank-Wolfe 的独特优势与局限

特性	Frank-Wolfe	投影梯度法	FISTA
每步核心操作	LMO（线性最小化）	投影 \(\Pi_C\)	梯度 + prox
解的稀疏性	自动（极点组合）	无保证	取决于 \(g\)
无需 \(L\)	可用线搜索	需要 \(1/L\) 步长	需要 \(1/L\) 步长
加速可能	一般不可（\(O(1/k)\) 是最优）	\(O(1/k^2)\) via AGD	\(O(1/k^2)\)
约束集要求	紧凸 + LMO 便宜	投影便宜	\(g\) 有 prox

FW 不可加速：Jaggi (2013) 和 Lan (2013) 证明了对一般紧凸集，FW 的 \(O(1/k)\) 收敛率是最优的——不存在 \(O(1/k^2)\) 的 FW 变体。但在强凸场景下，away-step FW（Lacoste-Julien & Jaggi 2015）可以达到线性收敛。

FW 的典型应用场景：

问题	LMO 操作	代价	投影代价（对比）
核范数球	最大奇异值对	\(O(mn)\)	\(O(\min(m,n) \cdot mn)\) SVD
单纯形 \(\Delta_n\)	\(\arg\min_i c_i\)	\(O(n)\)	\(O(n\log n)\) 排序
\(\ell_1\) 球	$\pm e_{\arg\max	c_i	}$
谱面体 \(\{X \succeq 0, \text{tr}(X) \le 1\}\)	最大特征向量	\(O(n^2)\)	\(O(n^3)\) 特征分解
流多面体	最短路径	$O(	V

13.4 Frank-Wolfe 在矩阵补全中的应用

矩阵补全问题：已知矩阵 \(M\) 在观测集 \(\Omega\) 上的值，恢复低秩矩阵 \(X\)：

\[\min_X \sum_{(i,j) \in \Omega} (X_{ij} - M_{ij})^2 \quad \text{s.t. } \|X\|_* \le \tau\]

Frank-Wolfe 的每步：

1. 计算梯度 \(G_k\)：\(G_{ij} = \begin{cases} 2(X_{ij}^{(k)} - M_{ij}) & (i,j) \in \Omega \\ 0 & \text{otherwise} \end{cases}\)

2. LMO：\(s_k = -\tau u_1 v_1^\top\)，其中 \((u_1, v_1)\) 是 \(G_k\) 的最大奇异值对

3. 更新：\(X_{k+1} = (1-\gamma_k)X_k + \gamma_k s_k\)

关键优势：\(X_k\) 保持为秩-\(k\) 矩阵（\(k\) 个秩-1 矩阵的凸组合），避免了存储和分解完整矩阵。

⚠️ 常见陷阱

⚠️ 编程陷阱：FW 的步长与 LMO 方向 - FW 不是 \(x_{k+1} = x_k - \gamma_k \nabla f(x_k)\)（那是梯度下降） - FW 是 \(x_{k+1} = (1-\gamma_k) x_k + \gamma_k s_k\)（凸组合） - LMO 方向 \(s_k - x_k\) 不一定是下降方向——但步长 \(\gamma_k\) 保证了充分下降

练习

1. 对单纯形上的优化 \(\min f(x)\) s.t. \(x \in \Delta_n\)，写出 Frank-Wolfe 的完整一步迭代（LMO + 更新），并与 Mirror Descent 对比。
2. 证明 Frank-Wolfe 的 \(O(1/k)\) 对一般紧凸集是最优的（提示：构造目标函数使得每步的 FW gap \(\delta_k = \max_{s \in \mathcal{X}} \langle \nabla f(x_k), x_k - s \rangle\) 恰好是 \(O(1/k)\)）。

---

14. 一阶 vs 二阶方法的完整对比 ⭐⭐

14.1 复杂度的三个维度

算法比较不能只看"迭代次数"，必须同时考虑三个维度：

维度	含义	一阶法	二阶法
迭代复杂度	到 \(\epsilon\) 最优需要多少步	\(O(\sqrt{L/\epsilon})\)	\(O(\sqrt{m}\log(1/\epsilon))\)
每步算术复杂度	每步需要多少次浮点运算	\(O(n)\) or \(O(\text{nnz})\)	\(O(n^3)\) or \(O(n \cdot \text{nnz})\)
存储复杂度	需要多少内存	\(O(n)\)	\(O(n^2)\) (Hessian)

总复杂度 = 迭代复杂度 \(\times\) 每步算术复杂度。

14.2 分水岭分析

设问题有 \(n\) 个变量、\(m\) 个约束、精度要求 \(\epsilon\)。

一阶法总代价（以 FISTA 为例）：\(T_1 = O(\sqrt{L/\epsilon}) \cdot O(n) = O(n\sqrt{L/\epsilon})\)

二阶法总代价（以 IPM 为例）：\(T_2 = O(\sqrt{m}\log(1/\epsilon)) \cdot O(n^3) = O(n^3\sqrt{m}\log(1/\epsilon))\)

一阶法更快的条件：\(T_1 < T_2 \Leftrightarrow n^2\sqrt{m} > \sqrt{L/\epsilon}/\log(1/\epsilon)\)。

定性结论：

场景	最优选择	理由
\(n\) 小（\(< 10^3\)）, \(\epsilon\) 小（\(< 10^{-8}\)）	IPM	Newton 步少且每步代价可控
\(n\) 大（\(> 10^4\)）, \(\epsilon\) 中（\(\sim 10^{-4}\)）	FISTA/ADMM	每步代价低，中精度够用
MPC（\(n \sim 10^2\), warm-start）	ADMM (OSQP) 或 IPM (HPIPM)	warm-start 让两者都极快
LASSO（\(n \sim 10^6\), 稀疏）	FISTA/CD	IPM 的 \(n^3\) 完全不可行
SDP（\(n \sim 10^2\)）	IPM (Mosek)	SDP 没有好用的一阶法
在线/流式数据	SGD/SVRG	只能用一阶法

14.3 结构利用——超越一阶 vs 二阶的二分法

真正重要的不是"一阶 vs 二阶"，而是"你利用了多少问题结构"。 以 MPC-QP 为例：

方法	是否利用 OCP 结构	复杂度	求解时间（\(n_x=12, T=20\)）
通用 IPM（ECOS）	否	\(O((n_x T)^3 \sqrt{m}\log(1/\epsilon))\)	10 ms
结构化 IPM（HPIPM）	是（Riccati）	\(O(n_x^3 T \sqrt{m}\log(1/\epsilon))\)	0.1 ms
通用 ADMM（SCS）	否	\(O(n_x T \cdot K)\)	5 ms
结构化 ADMM（OSQP）	部分（因子缓存）	\(O(n_x^2 T \cdot K)\)	0.5 ms

HPIPM 比 ECOS 快 100 倍，不是因为"IPM 更好"，而是因为利用了 Riccati 结构（把 \(T^2\) 因子消除了）。

本质洞察：算法选择的黄金法则不是"选一阶还是二阶"，而是"选最能利用你的问题结构的方法"。结构利用带来的加速远大于算法本身的理论改进。HPIPM 利用 OCP 结构从 \(O(n^3 T^3)\) 降到 \(O(n^3 T)\)——这是 \(T^2\) 倍加速（\(T=20\) 时是 400 倍），远超 Nesterov 加速的 \(\sqrt{\kappa}\) 倍。

14.4 混合策略

实践中最好的求解器往往混合使用一阶和二阶思想：

acados 的 SQP + HPIPM：外层用 SQP（一阶，解非线性 MPC），每个 SQP 步内用 HPIPM（二阶，解 QP 子问题）。这结合了一阶法处理非凸性和二阶法的高精度。

L-BFGS：不存储完整 Hessian（\(O(n^2)\) 内存），只用最近 \(m\) 步的梯度差 \(\{s_k, y_k\}\) 隐式近似 Hessian 逆（\(O(mn)\) 内存）。每步代价 \(O(mn)\)——介于 GD 的 \(O(n)\) 和 Newton 的 \(O(n^3)\) 之间，但收敛速度接近 Newton（超线性）。

预条件一阶法：用不完全 Cholesky 或 AMG 预处理后的 CG 求解 Newton 系统。每步代价 \(O(n \cdot \text{nnz})\)，兼顾稀疏性和二阶信息。

⚠️ 常见陷阱

🧠 思维陷阱：认为"大问题一定用一阶法" - 对稀疏大规模 QP（如 SLAM 的 Hessian 是 block-sparse 的），稀疏 Cholesky（二阶）可能比 ADMM（一阶）更快 - Ceres 求解百万变量的 BA 问题用的是 Gauss-Newton + 稀疏 Cholesky——本质是二阶法 - 关键：稀疏性改变了复杂度分析——\(O(n^3)\) 变成 \(O(n \cdot \text{nnz})\)

练习

1. 对 \(n = 500\), \(m = 1000\), \(\epsilon = 10^{-6}\) 的 SOCP，估算 FISTA（忽略投影代价）和 IPM 的总浮点运算次数，判断哪个更快。
2. 解释为什么 L-BFGS 在深度学习中不如 Adam 流行，但在凸优化（如 scikit-learn 的 logistic regression）中是默认选择。
3. 设计一个混合策略：对一个大规模 LASSO 问题，先用 FISTA 到中精度 \(\epsilon_1 = 10^{-3}\)，识别 active set（非零分量），再对 active 变量用 Newton 法精求解到 \(\epsilon_2 = 10^{-10}\)。分析这个策略的总代价。

---

15. 算法选择决策树的深化 ⭐⭐

15.1 按问题类型的详细推荐

无约束光滑凸：

问题：min f(x)，f L-光滑凸
├── n <= 1e3, f 二次 → CG（n 步精确）
├── n <= 1e3, f 非二次 → L-BFGS（超线性收敛）
├── n > 1e4, 已知 mu → Nesterov AGD（O(sqrt(kappa) log(1/eps))）
├── n > 1e4, mu 未知 → AGD + 自适应重启
└── 有限和 f = sum f_i → SVRG/SAGA（if n 大）or GD（if n 小）

复合非光滑 \(f + g\)：

问题：min f(x) + g(x)，f L-光滑，g 有 prox
├── g = lam*||x||_1 → FISTA-LASSO（标准）
│   ├── n > 1e4 → 坐标下降（glmnet）
│   └── 需要高精度 → FISTA + active-set Newton
├── g = delta_C (约束) → 投影梯度 or 加速投影
│   ├── 投影便宜（Box, Ball）→ 加速 PGD
│   └── 投影贵（核范数球）→ Frank-Wolfe
├── g = lam*||X||_* (核范数) → FISTA-SVT or FW
└── g 复杂，无闭式 prox → ADMM 分裂

约束凸：

问题：min f0(x) s.t. f_i(x) <= 0
├── LP → HiGHS/Gurobi（单纯形 or IPM）
├── QP, n <= 1e3 → qpOASES（active-set, 精确）
├── QP, MPC 结构 → HPIPM（Riccati IPM）
├── QP, n > 1e3 → OSQP（ADMM）
├── SOCP → ECOS（IPM）or Clarabel
├── SDP, n <= 1e2 → Mosek（IPM, 高精度）
├── SDP, n > 1e2 → SCS（ADMM, 低精度）
└── NLP → Ipopt（IPM）or SNOPT（SQP）

15.2 按应用场景的选择指南

应用场景	实时性要求	精度要求	推荐求解器	理由
足式 MPC (kHz)	< 1 ms	\(10^{-4}\)	HPIPM/acados	OCP 结构 + Riccati
四旋翼 NMPC	< 10 ms	\(10^{-3}\)	acados (SQP+HPIPM)	非线性 + warm-start
全身控制 WBC	< 2 ms	\(10^{-4}\)	ProxQP/OSQP	中等规模 QP
CBF-QP 安全控制	< 0.5 ms	\(10^{-6}\)	qpOASES	小规模 + 精确
SLAM 后端	< 100 ms	\(10^{-6}\)	Ceres (GN+稀疏)	稀疏结构利用
离线轨迹优化	< 10 s	\(10^{-8}\)	Ipopt / Drake	非凸 + 高精度
稀疏恢复 LASSO	< 1 s	\(10^{-4}\)	FISTA / glmnet	大规模 + prox
实验设计 SDP	< 60 s	\(10^{-6}\)	Mosek	SDP 最可靠

15.3 求解器生态系统地图

机器人常用凸优化求解器
├── 一阶法
│   ├── OSQP (QP, ADMM) ——— MPC 通用
│   ├── SCS (Conic, ADMM) ——— SDP/SOCP 大规模
│   ├── ProxQP (QP, 增广 Lagrangian) ——— WBC/TSID
│   └── Clarabel (Conic, IPM+一阶) ——— 新一代 Rust 实现
├── 二阶法
│   ├── HPIPM (QP, Riccati IPM) ——— MPC 极速
│   ├── ECOS (SOCP, IPM) ——— 嵌入式友好
│   ├── Mosek (LP/QP/SOCP/SDP, IPM) ——— 商业级精度
│   └── Gurobi (LP/QP/MILP, IPM+Simplex) ——— 运筹学标准
├── 混合
│   ├── qpOASES (QP, Active-set) ——— 参数化 QP
│   └── acados (NMPC, SQP+HPIPM) ——— 非线性 MPC
└── 框架
    ├── CVXPY (Python 建模) → 调用上述求解器
    ├── CasADi (C++ 符号微分) → 生成 C 代码
    └── Drake (C++ 全栈) → 机器人仿真+优化

⚠️ 常见陷阱

💡 概念误区：认为"越新的求解器越好" - Clarabel (2023) 比 ECOS (2013) 新，但 ECOS 在 SOCP 上仍然可靠且成熟 - 成熟度和社区支持有时比算法先进性更重要 - 选择求解器的正确流程：先确定问题类型 → 查该类型的 benchmark → 选 benchmark 表现最好且社区活跃的

🧠 思维陷阱：把建模工具和求解器混为一谈 - CVXPY 是建模工具（把问题转换为标准形式），不是求解器 - CasADi 是微分+代码生成工具，不直接求解 - Drake 是全栈框架（仿真+优化+控制），内部调用多个求解器 - 理解这个分层：Problem → Modeling Tool → Standard Form → Solver

练习

1. 你需要为一个 24 自由度人形机器人设计 WBC。全身动力学约束 \(M\ddot{q} + h = S^\top\tau + J_c^\top\lambda\)，加上摩擦锥约束 \(\|f_t\| \le \mu f_n\)。识别问题类型并选择求解器。
2. 比较 OSQP 和 ProxQP 在以下场景中的表现：(a) 10Hz MPC（大 horizon）；(b) 1kHz WBC（小规模、频繁变化）。

---

16. 跨章综合：从理论到实战 ⭐⭐

16.1 从凸分析到算法选择的完整链路

回顾凸分析基础（§1.10）中建立的概念：\(\mu\)-强凸性和 \(L\)-光滑性通过条件数 \(\kappa = L/\mu\) 决定了梯度下降的收敛率。共轭函数与proximal算子章节（§2.3）揭示了强凸-光滑的对偶不变性——\(f\) 的条件数等于 \(f^*\) 的条件数。凸优化问题与对偶理论（§3.5）的 KKT 条件是所有约束优化求解器的终止判据。

完整链路：

章节	提供的工具	在算法选择中的角色
10_凸分析	\(L\), \(\mu\), \(\kappa\), 次梯度 \(\partial f\)	判断问题类型（光滑/非光滑/强凸），确定收敛率
20_共轭/Proximal	\(f^*\), \(\text{prox}_f\), Moreau 包络	判断能否用 proximal 方法，prox 是否有闭式
30_对偶理论	KKT, Slater, 锥规划层级	识别问题类型（LP/QP/SOCP/SDP），选择求解器
40_算法路线图（本章）	收敛率表, 决策树, 实战求解器	落地执行——把理论映射到代码

16.2 三个完整案例

案例 1：LASSO 的算法选择全过程

问题：\(\min \frac{1}{2}\|Ax - b\|^2 + \lambda\|x\|_1\)，\(A \in \mathbb{R}^{1000 \times 5000}\)。

Step 1（凸分析）：光滑部分 \(f(x) = \frac{1}{2}\|Ax-b\|^2\) 的 \(L = \lambda_{\max}(A^\top A) = \sigma_{\max}^2(A)\)。非光滑部分 \(g(x) = \lambda\|x\|_1\) 不可微但有闭式 prox（软阈值）。整体是凸的（不强凸，\(\mu = 0\)）。

Step 2（Proximal）：\(g\) 的 \(\text{prox}_{\alpha g}(v) = \text{soft}_{\alpha\lambda}(v)\)，代价 \(O(n)\)——极便宜。适合 proximal gradient。

Step 3（算法选择）：\(f + g\) 结构，\(f\) 光滑，\(g\) 有廉价 prox → FISTA。收敛率 \(O(L/\epsilon \cdot 1/k^2)\)。

要达到 \(\epsilon = 10^{-6}\)：\(k \ge \sqrt{L/\epsilon}\)。若 \(L = 10^4\)，需要 \(k \ge 10^5\) 步，每步 \(O(mn) = O(5 \times 10^6)\)，总计约 \(5 \times 10^{11}\) 次运算。

替代方案：若加 elastic net \(\frac{\mu_0}{2}\|x\|^2\)（\(\mu_0 = 0.01L\)），条件数 \(\kappa = L/(L\cdot 0.01) = 100\)，FISTA 只需 \(O(\sqrt{\kappa}\log(1/\epsilon)) \approx 10 \cdot 14 = 140\) 步——快了 1000 倍。

案例 2：四足 MPC 的求解器选择

问题：\(\min \frac{1}{2}x^\top Px + q^\top x\) s.t. \(Gx \le h\), \(Ax = b\)，\(n = 200\), \(m = 500\)，每 20 ms 必须解出。

分析： - 问题类型：QP（\(P \succeq 0\)，线性约束）→ 有高效求解器 - 规模：\(n = 200\) 中等，\(m = 500\) 中等 - 实时性：50 Hz → 每步 \(\le 20\) ms，实际留给 QP 的时间 \(\le 5\) ms - 结构：MPC 的 KKT 矩阵是 block-tridiagonal（OCP 结构） - Warm-start：时序求解，上一步的解是极好的初始点

选择：HPIPM（Riccati IPM），利用 OCP 结构，预计 0.1-0.5 ms。备选 OSQP（ADMM + warm-start），预计 0.3-1 ms。

案例 3：Certifiable 旋转同步

问题：\(\min \sum_{(i,j)} \|R_i - R_{ij}R_j\|_F^2\) s.t. \(R_i \in SO(3)\)。

分析： - 非凸（\(SO(3)\) 不凸）→ 不能直接用凸求解器 - SDP 松弛：\(X \succeq 0\)，\(X_{ii} = I_3\)，\(n = 3N\) 维 SDP - 对偶证书：若 SDP 最优解秩为 3，找到了全局最优

策略： 1. 解 SDP 松弛（Mosek, \(N \le 10^2\) 时可行） 2. 检查 \(\text{rank}(X^*)\)——若为 3，恢复旋转矩阵 3. 若松弛不紧，用 Riemannian Staircase 在流形上搜索

实际求解器：SE-Sync（Rosen 等 2019）把 SDP 松弛和 Riemannian 搜索结合在一起，对低噪声 SLAM 问题能在线性时间内给出带证书的全局最优解。

跨章综合练习

📝 综合练习（需要 10/20/30/40 四章知识）

1. 从定义到算法：对 \(\min f(x) = \frac{1}{2}x^\top Ax + b^\top x + \lambda\|x\|_1\)（\(A \succ 0\)）： (a) 计算 \(\mu\), \(L\), \(\kappa\)（Ch1：强凸性/光滑性）。 (b) 计算 \(f^*\)（Ch2：共轭函数）。 (c) 写出 KKT 条件并解释稀疏性（Ch3：Fermat 规则 + KKT）。 (d) 选择算法并估算迭代次数（Ch4：FISTA 收敛率）。 (e) 实现 FISTA 求解并画收敛曲线。

2. 从对偶到算法：考虑摩擦锥约束 QP \(\min \frac{1}{2}\|f - f_d\|^2\) s.t. \(\|f_t\| \le \mu f_n\), \(f_n \ge 0\)： (a) 识别问题类型（Ch3：SOCP）。 (b) 写出 Lagrange 对偶和 KKT 条件（Ch3：对偶理论）。 (c) 选择求解器并分析实时性（Ch4：ECOS vs OSQP）。 (d) 若需要 1 kHz，讨论哪些近似可以让问题更快求解。

3. Moreau 包络与算法：对非光滑 \(f(x) = \|x\|_1\)： (a) 计算 Moreau 包络 \(M_f^\lambda(v)\) 及其梯度（Ch2）。 (b) 用梯度下降最小化 \(M_f^\lambda\)，解释这等价于什么算法（Ch4：proximal point algorithm）。 (c) 讨论 \(\lambda\) 的选择对收敛率的影响。

---

向后指向：从算法路线图到工程实战

本章建立的算法选择能力将在后续工程章节中直接使用：

本章工具	工程应用	关键对接点
GD/AGD 收敛率	SLAM 迭代求解器参数调优	\(L\) 的估计 → 步长选择
FISTA	稀疏重建、鲁棒估计	prox 闭式 → 实时可行性
ADMM (OSQP)	线性 MPC、全身控制	warm-start → kHz 实时
IPM (HPIPM)	高精度 MPC、轨迹优化	Riccati → OCP 结构利用
Frank-Wolfe	低秩矩阵恢复、核范数优化	LMO → 避免昂贵投影
决策树	任何优化问题的算法选型	问题结构 → 求解器选择

下一步学习建议：

1. 足式控制方向：先学 MPC 章节，用 OSQP/HPIPM 求解 QP → 回顾本章 §5/§6 的理论
2. SLAM 方向：先学图优化章节，用 Ceres 解非线性最小二乘 → 回顾本章 §2 的 GD/Newton 理论
3. RL 方向：先学策略优化章节，理解 SGD/Adam 的理论基础 → 回顾本章 §7.4 的随机方法
4. 规划方向：先学轨迹优化章节，用 acados 解 NMPC → 回顾本章 §6 的 IPM + §9 的决策树

---

17. ADMM 的深层分析——从 DR splitting 推导 ⭐⭐⭐

17.1 Proximal Point Algorithm 的对偶形式

回顾共轭函数与proximal算子章节（§2.16）中的不动点解释：\(x^* = \arg\min f \Leftrightarrow x^* = \text{prox}_{\lambda f}(x^*)\)。Proximal Point Algorithm (PPA) 直接迭代这个不动点映射：

\[x_{k+1} = \text{prox}_{\lambda f}(x_k) = \arg\min_x \left[f(x) + \frac{1}{2\lambda}\|x - x_k\|^2\right]\]

PPA 的收敛率：\(f(x_k) - f^* = O(1/k)\)（Guler 1991），加速版本（APM）为 \(O(1/k^2)\)（Salzo-Villa 2012）。

PPA 的问题：每步需要精确求解一个 proximal 子问题——对复杂 \(f\)，这本身就是一个优化问题。PPA 更多是理论框架，实际算法（proximal gradient、ADMM）是 PPA 的"可行化版本"。

17.2 从 Douglas-Rachford 到 ADMM

问题：\(\min h_1(u) + h_2(u)\)（两个闭凸函数之和）。

Douglas-Rachford splitting (DRS)：

\[\hat{u}^{k+1} = \text{prox}_{\gamma h_2}(u^k)$$ $$u^{k+1} = u^k + \text{prox}_{\gamma h_1}(2\hat{u}^{k+1} - u^k) - \hat{u}^{k+1}\]

等价的反射形式：定义反射算子 \(R_f = 2\text{prox}_f - I\)，则 DRS 等价于：

\[u^{k+1} = \frac{1}{2}(R_{h_1} \circ R_{h_2} + I)(u^k)\]

即对当前点做 \(h_2\) 的反射，再做 \(h_1\) 的反射，最后取平均。

DRS → ADMM 的对偶等价（Eckstein-Bertsekas 1992）：

考虑 ADMM 问题 \(\min f(x) + g(z)\) s.t. \(x = z\)。其对偶问题可以写成 \(\min f^*(-y) + g^*(y)\)。

在对偶空间上应用 DRS： $\(y^{k+1/2} = \text{prox}_{\rho g^*}(y^k)\)$ $\(y^{k+1} = y^k + \text{prox}_{\rho f^*(-\cdot)}(2y^{k+1/2} - y^k) - y^{k+1/2}\)$

利用 Moreau 分解 \(\text{prox}_{f^*}(v) = v - \text{prox}_f(v)\)，可以把对偶 DRS 翻译回原空间——恰好得到 ADMM 的三步迭代。

本质洞察：ADMM 不是"凭空发明"的算法——它是 Douglas-Rachford splitting 作用在**对偶问题**上的等价形式。理解这个等价性有三个好处：(1) 可以直接从 DRS 的收敛理论推出 ADMM 的收敛性；(2) 解释了为什么 ADMM 不需要原问题可微；(3) 为设计新的分裂算法提供了系统框架。

17.3 ADMM 参数调优的理论指导

回顾 §5.6 中 \(\rho\) 选择的讨论。这里给出更精确的理论结果。

最优 \(\rho\) 的理论值（Giselsson-Boyd 2017）：对 \(f\) 是 \(\mu_f\)-强凸且 \(L_f\)-光滑、\(g\) 是 \(\mu_g\)-强凸且 \(L_g\)-光滑的问题：

\[\rho^* = \sqrt{\frac{\mu_f L_f}{\mu_g L_g}} \cdot \sqrt[4]{\frac{L_g}{L_f}}\]

简化版本：当 \(f\) 是 QP 的目标（\(\mu_f = \lambda_{\min}(P)\), \(L_f = \lambda_{\max}(P)\)），\(g\) 是指示函数（\(\mu_g = 0\)，退化）：

\[\rho \approx \sqrt{\lambda_{\min}(P) \cdot \lambda_{\max}(P)} = \sqrt{\det(P)^{1/n}} \quad \text{（几何均值启发）}\]

OSQP 的自适应策略详述：

OSQP 使用 Ruiz equilibration 预处理 + 残差平衡。具体地，它维护原残差 \(r_k = Ax_k + Bz_k - c\) 和对偶残差 \(s_k = \rho A^\top B(z_k - z_{k-1})\) 的比值：

\[\text{if } \|r_k\| > \mu_{\text{ratio}} \|s_k\| : \quad \rho \leftarrow \tau \rho \quad \text{（增大 $\rho$，加速原可行性收敛）}$$ $$\text{if } \|s_k\| > \mu_{\text{ratio}} \|r_k\| : \quad \rho \leftarrow \rho/\tau \quad \text{（减小 $\rho$，加速对偶可行性收敛）}\]

OSQP 默认 \(\mu_{\text{ratio}} = 10\), \(\tau = 2\)。

注意：\(\rho\) 改变时 KKT 矩阵的 \((2,2)\) 块变化，需要重新分解。OSQP 通过限制 \(\rho\) 变化频率（每 25-200 步检查一次）来摊销重分解代价。

17.4 Consensus ADMM 与分布式优化

当问题有如下结构时，ADMM 可以并行化：

\[\min \sum_{i=1}^N f_i(x_i) \quad \text{s.t. } x_i = z, \quad i = 1, \ldots, N\]

Consensus ADMM：

\[x_i^{k+1} = \text{prox}_{f_i/\rho}(z^k - u_i^k) \quad \text{（$N$ 个子问题并行求解）}$$ $$z^{k+1} = \frac{1}{N}\sum_{i=1}^N (x_i^{k+1} + u_i^k) \quad \text{（全局平均）}$$ $$u_i^{k+1} = u_i^k + x_i^{k+1} - z^{k+1} \quad \text{（对偶更新）}\]

通信模式：每步只需各 agent 发送 \(x_i + u_i\) 到中心节点，中心计算均值后广播 \(z\)。每步通信量 \(O(Nn)\)。

在多机器人系统中的应用：

场景	\(f_i\)	通信拓扑
分布式 SLAM	第 \(i\) 个机器人的 BA 子问题	星形（中心服务器）
编队控制	第 \(i\) 个机器人的局部 MPC	网状（邻居通信）
协同搬运	第 \(i\) 个手臂的力分配	星形

练习

1. 对 \(\min f(x) + g(x)\)（\(f\), \(g\) 都有廉价 prox），写出 Douglas-Rachford splitting 的完整迭代，并验证当 \(f(x) = \frac{1}{2}\|Ax-b\|^2\), \(g(x) = \lambda\|x\|_1\) 时等价于 ADMM。
2. 实现 consensus ADMM 解 \(\min \sum_{i=1}^{10} \frac{1}{2}\|A_i x - b_i\|^2\)（10 个线性回归子问题），与集中式 GD 对比迭代次数和通信量。
3. 对一个 QP 实例，测试 \(\rho = 0.01, 0.1, 1, 10, 100\) 的 ADMM 迭代次数，画出"\(\rho\) vs 迭代次数"的 U 形曲线，验证最优 \(\rho\) 的理论预测。

---

18. Nesterov 加速的 ODE 视角 ⭐⭐⭐

18.1 Su-Boyd-Candes ODE

回顾 §3.6 中的 ODE 解释。这里深入分析 Su-Boyd-Candes (2016) 的连续时间框架。

Nesterov AGD 的连续极限（\(\alpha \to 0\)，离散步长趋零）是二阶 ODE：

\[\ddot{X}(t) + \frac{r}{t}\dot{X}(t) + \nabla f(X(t)) = 0, \quad t > 0\]

其中 \(r = 3\) 对应经典 AGD 的 \(O(1/k^2)\) 收敛率。

物理类比：这是一个在势能场 \(f\) 中运动的粒子，受到时变阻尼 \(\frac{r}{t}\dot{X}\)。起初阻尼很大（\(t\) 小时 \(r/t\) 大），粒子缓慢启动（避免初始不稳定）；后来阻尼减弱（\(t\) 大时 \(r/t\) 小），粒子加速滑向最小值。

Lyapunov 分析：定义能量函数

\[E(t) = t^2(f(X(t)) - f^*) + 2\|X(t) - x^* + \frac{t}{2}\dot{X}(t)\|^2\]

可以证明 \(\dot{E}(t) \le 0\)（当 \(r \ge 3\) 时），因此 \(E(t) \le E(0)\)。

由 \(E(t) \ge t^2(f(X(t)) - f^*)\)，得 \(f(X(t)) - f^* \le E(0)/t^2 = O(1/t^2)\)。

离散化后（步长 \(h\), \(t = kh\)），恢复出 \(f(x_k) - f^* = O(1/k^2)\)。

18.2 阻尼系数 \(r\) 的意义

\(r\) 的临界值：

\(r\) 值	连续收敛率	离散化行为
\(r < 3\)	\(O(1/t^r)\)（次优）	可能不稳定
\(r = 3\)	\(O(1/t^2)\)（最优）	经典 Nesterov AGD
\(r > 3\)	\(O(1/t^2)\)（同 \(r=3\)）	更稳定但不更快

Su-Boyd-Candes 证明了 \(r = 3\) 是使 \(O(1/t^2)\) 首次成立的最小 \(r\)。更大的 \(r\) 不能加速，但增加阻尼使轨迹更平滑（实践中有助于稳定性）。

反事实推理：如果阻尼恒定（\(\frac{r}{t}\) 替换为常数 \(c\)），ODE 变为 \(\ddot{X} + c\dot{X} + \nabla f = 0\)（阻尼谐振子）。这对应 Heavy-ball 方法——收敛率是 \(O(e^{-ct})\)（指数衰减），但需要精确知道 \(c = 2\sqrt{\mu}\)（强凸参数的平方根）。Nesterov 的时变阻尼不需要知道 \(\mu\)——这正是 AGD 的"自适应"优势。

18.3 从 ODE 设计新算法

ODE 框架不仅解释了已有算法，还可以**设计新算法**：

方法 1（Wibisono-Wilson-Jordan 2016）：考虑广义 ODE

\[\ddot{X} + (\beta + e^{\alpha_t}\mu)\dot{X} + e^{2\alpha_t + \beta_t}(\nabla^2 f \cdot \dot{X} + \nabla f) = 0\]

通过选择不同的 \(\alpha_t, \beta_t\)，可以统一推出 AGD、mirror descent、加速 mirror descent 等一系列算法。

方法 2（Shi et al. 2021）：对强凸函数的更优 ODE

\[\ddot{X} + 2\sqrt{\mu}\dot{X} + \sqrt{s}\nabla^2 f(X)\dot{X} + (1 + \sqrt{\mu s})\nabla f(X) = 0\]

其中 \(s = 1/L\)。这个 ODE 的离散化给出了 triple momentum method——收敛率为 \(\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^k\)，与 AGD 相同但常数更好。

18.4 Symplectic 积分器视角

ODE \(\ddot{X} + \frac{3}{t}\dot{X} + \nabla f = 0\) 可以写成 Hamilton 系统的推广。用辛积分器（symplectic integrator）离散化，可以保持能量守恒的离散类比——这解释了为什么 Nesterov AGD 的 Lyapunov 函数"恰好"单调不增。

与控制理论的联系：这个 ODE 可以看作一个线性时变系统的状态反馈控制。\(\nabla f\) 是"驱动力"，\(\frac{3}{t}\dot{X}\) 是"阻尼控制"。Lyapunov 分析与控制论中的稳定性分析完全同构——这就是为什么控制论工具（如 IQC, Lessard-Recht-Packard 2016）可以用于分析优化算法。

练习

1. 用 scipy.integrate.solve_ivp 数值求解 ODE \(\ddot{X} + \frac{3}{t}\dot{X} + \nabla f(X) = 0\)（\(f(x) = \frac{1}{2}(x_1^2 + 100x_2^2)\)），画出 \(X(t)\) 的轨迹，与离散 AGD 对比。
2. 改变阻尼系数 \(r\)（从 1 到 5），观察 \(f(X(t)) - f^*\) 的衰减速率变化，验证 \(r = 3\) 是临界值。
3. 对 Heavy-ball ODE \(\ddot{X} + c\dot{X} + \nabla f = 0\)（\(c\) 为常数），分析当 \(c = 2\sqrt{\mu}\) 时的收敛率，并解释为什么 \(c\) 的选择依赖于未知的 \(\mu\)。

---

19. 凸优化算法的数值实验框架 ⭐⭐

19.1 如何做有意义的算法对比

凸优化算法的实验对比容易犯以下错误：

错误 1：只报告迭代次数，不报告总计算时间。ADMM 可能 50 步收敛，但每步涉及线性系统求解；GD 可能 5000 步，但每步只有矩阵-向量乘。总时间取决于每步代价。

错误 2：只测试一个问题实例。算法性能对问题结构（条件数、稀疏性、约束数量）高度敏感。必须在多个实例上测试，报告统计量（均值、中位数、最大值）。

错误 3：比较不同精度下的算法。ADMM 在 \(\epsilon = 10^{-3}\) 时可能比 IPM 快，但在 \(\epsilon = 10^{-8}\) 时反过来。必须画"精度 vs 时间"曲线（performance profile）。

19.2 标准测试集

测试集	问题类型	来源	规模范围
Maros-Meszaros	QP	经典基准	\(n: 10^1 \sim 10^4\)
Hans Mittelmann	LP/QP/SOCP/SDP	学术竞赛	\(n: 10^2 \sim 10^5\)
OSQP benchmark	MPC-QP	OSQP 论文	\(n: 10^2 \sim 10^4\)
CVXPortfolio	金融 QP	Boyd 课题组	\(n: 10^2 \sim 10^3\)
SuiteSparse	稀疏线性系统	数据收集	\(n: 10^2 \sim 10^6\)

19.3 Python 实验框架

import numpy as np
import time
from dataclasses import dataclass

@dataclass
class BenchResult:
    """算法基准测试结果"""
    name: str                # 算法名称
    obj_history: list        # 目标值历史
    time_history: list       # 时间戳历史
    n_iters: int             # 总迭代次数
    final_obj: float         # 最终目标值
    total_time: float        # 总时间

def benchmark_lasso(A, b, lam, algorithms, x_star=None, max_iter=5000):
    """
    在 LASSO 问题上对比多种算法

    参数：
        A: m x n 观测矩阵
        b: m 维观测向量
        lam: L1 正则化参数
        algorithms: 算法函数列表，每个接受 (A, b, lam, max_iter) 返回 (x, history)
        x_star: 真实最优解（用 CVXPY 预计算）
    """
    results = []
    f_star = 0.5 * np.linalg.norm(A @ x_star - b)**2 + lam * np.linalg.norm(x_star, 1) if x_star is not None else None

    for alg_name, alg_func in algorithms:
        t0 = time.perf_counter()
        x, obj_hist = alg_func(A, b, lam, max_iter=max_iter)
        total_time = time.perf_counter() - t0

        result = BenchResult(
            name=alg_name,
            obj_history=obj_hist,
            time_history=np.linspace(0, total_time, len(obj_hist)).tolist(),
            n_iters=len(obj_hist),
            final_obj=obj_hist[-1],
            total_time=total_time
        )
        results.append(result)
        print(f"{alg_name:20s}: {result.n_iters:5d} iters, "
              f"obj={result.final_obj:.6e}, time={result.total_time:.3f}s")

    return results

def plot_convergence(results, f_star=None):
    """绘制收敛曲线（目标值差 vs 迭代次数）"""
    import matplotlib.pyplot as plt
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

    for r in results:
        if f_star is not None:
            gaps = [max(obj - f_star, 1e-16) for obj in r.obj_history]
            ax1.semilogy(gaps, label=r.name)
            ax2.semilogy(r.time_history, gaps, label=r.name)
        else:
            ax1.semilogy(r.obj_history, label=r.name)
            ax2.semilogy(r.time_history, r.obj_history, label=r.name)

    ax1.set_xlabel('迭代次数')
    ax1.set_ylabel('f(x_k) - f*')
    ax1.set_title('收敛率对比（按迭代次数）')
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    ax2.set_xlabel('时间 (秒)')
    ax2.set_ylabel('f(x_k) - f*')
    ax2.set_title('收敛率对比（按计算时间）')
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

练习

1. 使用上述框架，对 \(m = 500, n = 2000\) 的 LASSO 问题对比 ISTA、FISTA、ADMM、坐标下降四种算法。改变 \(\lambda\)（从 \(0.01\lambda_{\max}\) 到 \(0.5\lambda_{\max}\)，其中 \(\lambda_{\max} = \|A^\top b\|_\infty\)），观察稀疏度对各算法收敛速度的影响。
2. 对 MPC-QP 问题，对比 OSQP（warm-start 开/关）和 qpOASES 的求解时间，画出"问题规模 vs 时间"的曲线。

---

凸优化算法的核心公式速查 ⭐⭐

收敛率速查

方法	凸（\(L\)-光滑）	强凸（\(\mu\), \(L\)）	适用结构
GD	\(O(LR^2/k)\)	\(O((1-\mu/L)^k)\)	光滑无约束
AGD	\(O(LR^2/k^2)\)	\(O((1-\sqrt{\mu/L})^k)\)	光滑无约束
ISTA	\(O(LR^2/k)\)	\(O((1-\mu/L)^k)\)	\(f\)(光滑)\(+g\)(prox)
FISTA	\(O(LR^2/k^2)\)	\(O((1-\sqrt{\mu/L})^k)\)	\(f\)(光滑)\(+g\)(prox)
ADMM	\(O(1/k)\)	\(O((1-c)^k)\)	可分 / 约束
FW	\(O(LD^2/k)\)	\(O(LD^2/k)\)	紧凸集 + LMO
IPM	\(O(\sqrt{m}\log(1/\epsilon))\)	同左	任意约束凸
次梯度	\(O(GR/\sqrt{k})\)	—	非光滑 Lipschitz

步长速查

方法	标准步长	条件
GD	\(\alpha = 1/L\)	\(L\)-光滑
AGD (一般凸)	\(\alpha = 1/L\), \(\beta_k = (k-1)/(k+2)\)	\(L\)-光滑
AGD (强凸)	\(\alpha = 1/L\), \(\beta = (\sqrt{\kappa}-1)/(\sqrt{\kappa}+1)\)	\(\mu\)-强凸 \(L\)-光滑
ADMM	\(\rho \approx \sqrt{\mu_f L_f}\)	自适应更好
次梯度	\(\alpha_k = c/\sqrt{k}\)	递减步长

下界速查（不可改进）

函数类	一阶方法下界	达到方法
\(L\)-光滑凸	\(\Omega(LR^2/k^2)\)	Nesterov AGD
\(\mu\)-强凸 \(L\)-光滑	\(\Omega((\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1})^k)\)	Nesterov AGD
Lipschitz 凸（非光滑）	\(\Omega(GR/\sqrt{k})\)	次梯度法
\(m\) 约束凸（二阶）	\(\Omega(\sqrt{m}\log(1/\epsilon))\)	IPM

常用 prox 算子速查

回顾共轭函数与proximal算子章节（§2.5）中的核心 prox 闭式解，在本章算法中反复使用：

\(g(x)\)	\(\text{prox}_{\lambda g}(v)\)	出现场景
\(\lambda\|x\|_1\)	\(\text{sign}(v)\odot\max(\|v\|-\lambda,0)\)	LASSO / ISTA / FISTA
\(\lambda\|x\|_2\)	\((1-\lambda/\|v\|)_+ \cdot v\)	Group LASSO
\(\delta_C(x)\)	\(\Pi_C(v)\)	PGD / ADMM 约束处理
\(\frac{\lambda}{2}\|x\|^2\)	\(v/(1+\lambda)\)	Ridge / Elastic Net
\(\lambda\|X\|_*\)	\(U\text{diag}(\text{soft}(\sigma,\lambda))V^\top\)	矩阵补全
\(\delta_{\Delta_n}(x)\)	\(\text{proj}_{\Delta_n}(v)\)	Mirror Descent 对比

共轭函数速查

回顾共轭函数与proximal算子章节（§2.1）中的核心共轭，在对偶分析和 ADMM 推导中使用：

\(f(x)\)	\(f^*(y)\)	连接
\(\frac{1}{2}x^\top Qx\)	\(\frac{1}{2}y^\top Q^{-1}y\)	QP 对偶
\(\|x\|_p\)	\(\delta_{\|y\|_q \le 1}\)	LASSO 对偶
\(\delta_C(x)\)	\(\sigma_C(y) = \sup_{x \in C}\langle y,x\rangle\)	支撑函数
\(\log\sum e^{x_i}\)	\(\sum y_i\log y_i\) (\(y \in \Delta_n\))	RL / softmax
\(-\log x\)	\(-1-\log(-y)\) (\(y < 0\))	Barrier 对偶

---

凸优化算法的历史时间线 ⭐

年代	方法/理论	贡献者	核心创新
1847	最速下降法	Cauchy	沿负梯度方向走——一切优化算法的起点
1952	共轭梯度法	Hestenes-Stiefel	利用共轭方向，\(n\) 步精确解二次
1964	Heavy-ball	Polyak	引入动量——加速的第一次尝试
1970	Proximal point	Rockafellar	隐式梯度步——非光滑优化的理论基石
1975	Douglas-Rachford	Lions-Mercier	两个 prox 交替——分裂方法的原型
1983	加速梯度法	Nesterov	\(O(1/k^2)\) 且证明最优——凸优化的"音障突破"
1984	内点法	Karmarkar	LP 的多项式时间算法——改变了计算复杂度理论
1994	Self-concordance	Nesterov-Nemirovsky	IPM 的统一复杂度理论
2004	凸优化教材	Boyd-Vandenberghe	让凸优化从数学走向工程
2009	FISTA	Beck-Teboulle	加速 + prox = 非光滑也能 \(O(1/k^2)\)
2011	ADMM 综述	Boyd-Parikh 等	分裂方法的工程化——5000+ 引用
2013	SVRG	Johnson-Zhang	方差缩减——有限和问题的最优一阶法
2014	PEP	Drori-Teboulle	用 SDP 自动分析算法——"自动化的 Nesterov"
2016	Su-Boyd-Candes ODE	Su-Boyd-Candes	加速法的连续时间理解——物理直觉
2017	OSQP	Stellato 等	ADMM 用于实时 QP——机器人 MPC 的标准工具
2020	HPIPM	Frison-Diehl	Riccati IPM——机器人 kHz 级 MPC 的计算引擎
2023	TinyMPC	多校合作	凸优化嵌入 MCU——无人机/微型机器人

历史启示：凸优化算法的发展遵循"理论突破 → 工程化 → 嵌入式"的三阶段模式。Nesterov 1983 年的加速理论经过 30 年才变成 OSQP 和 HPIPM 这样的嵌入式求解器。当前的前沿是把这些求解器进一步压缩到微控制器（TinyMPC）——让每一个机器人关节都能在本地解优化问题。

两大传统的融合：

传统	代表人物	核心方法	优势
苏联/东欧	Nesterov, Polyak, Nemirovsky	一阶法 + 复杂度理论	大规模、理论最优
美国/西欧	Boyd, Vandenberghe, Wright	二阶法 + 工程实现	高精度、通用软件

现代最优实践（如 acados = SQP(一阶外层) + HPIPM(二阶内层)）正是两大传统的融合。

---