第 09 章 大规模拟牛顿与部分可分优化

学习目标

读完本章，你应当能够：

1. 说清楚为什么稠密 BFGS 在 \(n\) 很大时不可行，以及有限内存、稀疏更新、部分可分三条扩展路线各自的思路与适用范围。
2. 写出 L-BFGS 的双循环递归（Algorithm 9.1），知道它的 \(4mn\) 运算量从何而来，会选初始矩阵 \(H_k^0=\gamma_kI\) 和记忆长度 \(m\)。
3. 理解无记忆 BFGS 与 Hestenes–Stiefel / Polak–Ribière 共轭梯度法的等价关系，从而把 L-BFGS 看作非线性 CG 的自然推广。
4. 读懂 BFGS 和 SR1 的紧凑表示（"对角 + 长窄矩阵 × 小矩阵 × 长窄矩阵"），知道它是 L-BFGS-B 和有限内存 SQP 的基础。
5. 识别部分可分结构与不变子空间，理解"按元素维护小 Hessian 近似"为何比整体近似好得多；并在因子模型协方差 \(\Sigma=BFB^T+D\) 中认出同样的结构思想。

读前导读

这一章在解决什么问题

第 08 章的 BFGS 要存一个 \(n\times n\) 的矩阵。\(n\) 是几十个参数时没问题；但如果你在 3000 只股票上做组合优化，或估计一个有上万个参数的截面模型，这个矩阵就有几百万到上亿个元素，存不下，也算不动。本章讲两种应对办法。

第一种是 L-BFGS：不存矩阵，只存最近 \(m\) 步的"位移 \(s\)"和"梯度变化 \(y\)"（\(m\) 通常取 5 到 20），需要时用这几对向量临时"拼出"矩阵与向量的乘积。这有点像只保留最近 \(m\) 期的数据做滚动估计：旧信息与当前的相关性低，丢掉它损失不大，却省下了大量存储。scipy 里常用的 L-BFGS-B 就是它的带上下界版本。

第二种是利用结构。你对这个思想其实很熟悉：因子模型把 3000×3000 的协方差矩阵写成 \(\Sigma=BFB^T+D\)，即"少数几个因子造成的低秩部分 + 对角的特质部分"。存储从 900 万个数降到 3 万个数，计算 \(\Sigma w\) 也快了两个数量级。本章的"紧凑表示"和"部分可分"都是这种"大矩阵其实由小零件拼成"的想法在优化中的版本。

需要先想起来的数学

1. 矩阵乘法的结合律与计算量。 \((AB)v=A(Bv)\)，结果相同，但计算量可能天差地别。若 \(B\) 是 \(n\times K\)、\(v\) 是长度 \(n\) 的向量，先算 \(B^Tv\)（\(nK\) 次乘法）得到长度 \(K\) 的向量，再乘回去，总共 \(O(nK)\)；先组装 \(BB^T\) 则需要 \(O(n^2K)\)。例：\(n=3000\)、\(K=10\) 时，前者约 6 万次，后者约 9000 万次。本章的双循环递归和因子协方差技巧都靠这一点。见 第 00 册第 06 章 线性代数速成。

2. 秩与低秩矩阵。 矩阵的秩是"线性无关的列的最大个数"。\(uv^T\) 秩为 1；\(BFB^T\)（\(B\) 为 \(n\times K\)）秩至多为 \(K\)。低秩矩阵虽然可能每个元素都非零（"稠密"），但信息量很少。例：全 1 矩阵 \(\mathbf1\mathbf1^T\) 每个元素都是 1，秩却只有 1。见 第 00 册第 06 章。

3. 子空间与零空间。 子空间是对加法和数乘封闭的一组向量，如"所有分量之和为零的向量" \(\{w:\mathbf1^Tw=0\}\)，在组合里就是"美元中性的调仓方向"。矩阵 \(U\) 的零空间是 \(\{w:Uw=0\}\)，即被 \(U\) "看不见"的方向。维数是这个子空间里线性无关向量的最大个数。见 第 00 册第 06 章。

4. 条件数。 对称正定矩阵的条件数 = 最大特征值 / 最小特征值，衡量"碗"有多扁。条件数越大，梯度类方法越容易在狭长谷底来回折返，迭代越多。见 第 00 册第 06 章。

5. 多元链式法则。 若 \(f(x)=\phi(Ux)\)，则 \(\nabla f(x)=U^T\nabla\phi(Ux)\)，\(\nabla^2f(x)=U^T\nabla^2\phi(Ux)U\)。一元类比：\(f(x)=\phi(ax)\) 时 \(f'=a\phi'\)、\(f''=a^2\phi''\)。见 第 00 册第 05 章 多元微积分与优化。

怎么读这一章

核心必读：9.1、9.2（L-BFGS，特别是 9.2.2 双循环和 9.2.5 优缺点）、9.6 与 9.7（量化中的位置与实战）。9.2.6 与 CG 的关系可以只记结论。9.3 紧凑表示第一次只需抓住"对角 + 长窄矩阵 × 小矩阵 × 长窄矩阵"这幅图像，公式细节不必记。9.4 原书自己都说"了解即可"，可以跳过。9.5 部分可分理论性较强，建议先读 9.5.1 的例子和 9.5.3 关于 \((x_1+\cdots+x_n)^2\) 的那段，其余等到真遇到超大规模结构化问题时再回来。

建议顺序：9.1 → 9.2.1–9.2.5 → 9.7.1 → 9.3（只看图像）→ 9.6 → 9.7.2 → 9.5 → 其余。

9.1 问题：稠密近似在大规模下不可承受

第 08 章的拟牛顿近似 \(H_k\) 一般是稠密矩阵，存储需要 \(n^2\) 个数，每步更新和矩阵–向量乘法需要 \(O(n^2)\) 运算。\(n=10^4\) 时光存储就是 \(10^8\) 个双精度数（800 MB），\(n=10^5\) 时完全不可行。原书给出三条扩展路线：

1. 有限内存拟牛顿法（limited-memory quasi-Newton）：只保存若干个长度为 \(n\) 的向量，隐式表示 Hessian 近似。稳健、廉价、易实现，但收敛不快。这是实践中最重要的一条路线。
2. 稀疏拟牛顿法：让近似矩阵模仿真 Hessian 的稀疏模式。效果令人失望，只作简介。
3. 利用部分可分性（partial separability）：大多数大规模目标函数都可以写成许多"只依赖少数变量方向"的元素函数之和，按元素维护小矩阵近似。通常收敛快且稳健，但需要目标函数的详细结构信息。

9.2 有限内存 BFGS（L-BFGS）

9.2.1 思想

适用场景：Hessian 难以以合理代价计算，或者稠密到无法存储。主要思想是只用最近 \(m\) 次迭代的曲率信息构造 Hessian 近似，丢弃更早的信息——它们与当前 Hessian 的关系本来就较小。代价是收敛速率一般只有线性，但常常可以接受。

回顾 BFGS 逆形式：

L-BFGS 只隐式保存 \(m\) 对向量 \(\{s_i,y_i\}\)，\(i=k-m,\dots,k-1\)；每产生新的一对，就删去最旧的一对。每步选一个初始矩阵 \(H_k^0\)（可以随迭代变化，这一点不同于标准 BFGS），反复应用 (9.2) 得

实践中 \(m\) 取 3 到 20 通常就令人满意。

9.2.2 双循环递归

我们其实从不需要 \(H_k\) 本身，只需要乘积 \(H_k\nabla f_k\)。把 (9.5) 从右往左展开，就得到著名的双循环递归（two-loop recursion）：

Algorithm 9.1（L-BFGS 双循环递归）

q ← ∇f_k
for i = k−1, k−2, ..., k−m:          # 第一个循环：从新到旧
    α_i ← ρ_i s_iᵀ q
    q ← q − α_i y_i
r ← H_k^0 q
for i = k−m, k−m+1, ..., k−1:        # 第二个循环：从旧到新
    β ← ρ_i y_iᵀ r
    r ← r + s_i (α_i − β)
return r                              # r = H_k ∇f_k

两个循环各做 \(m\) 次，每次两个内积加一次向量更新，合计约 \(4mn\) 次乘法（不计 \(H_k^0q\)；\(H_k^0\) 为对角阵时再加 \(n\) 次）。\(m=10\)、\(n=10^6\) 时，每步约 \(4\times10^7\) 次乘法、存储 \(2\times10^7\) 个数，完全可以承受。

这个结构还有一个好处：与 \(H_k^0\) 的乘法被隔离在两个循环中间，因此 \(H_k^0\) 可以自由选取、每步改变，甚至可以隐式给出（通过解 \(B_k^0r=q\) 得到 \(r\)），这为预条件留下了接口。

推导拆解：双循环为什么等于 \(H_k\nabla f_k\)。以 \(m=1\) 为例最清楚，此时 \(H_k=V^TH^0V+\rho ss^T\)，\(V=I-\rho ys^T\)。

1. 先算 \(Vg=g-\rho y(s^Tg)\)。令 \(\alpha=\rho s^Tg\)，就是 \(q=g-\alpha y\)。这正是第一个循环做的事。
2. 再乘 \(H^0\)：\(r=H^0q\)。
3. 再乘 \(V^T=I-\rho sy^T\)：\(V^Tr=r-\rho s(y^Tr)=r-\beta s\)，其中 \(\beta=\rho y^Tr\)。
4. 加上 \(\rho ss^Tg=\alpha s\)。合起来 \(r-\beta s+\alpha s=r+s(\alpha-\beta)\)，正是第二个循环的那一行。 \(m>1\) 时，(9.5) 是这样的结构一层套一层：第一个循环从最新的一对往旧剥开右边的 \(V\) 们，并顺手记下每层的 \(\alpha_i\)；第二个循环从最旧的一对往新把左边的 \(V^T\) 和 \(\rho ss^T\) 项依次加回来。全程只做"向量内积 + 向量加减"，从不出现 \(n\times n\) 矩阵。

9.2.3 初始矩阵

实践中最有效的选择是 \(H_k^0=\gamma_kI\)，

也就是第 08 章 (8.21) 的缩放因子：它估计真实 Hessian 沿最近搜索方向的尺度的倒数，使 \(p_k\) 的长度合适，大多数迭代都能接受 \(\alpha_k=1\)。和 BFGS 一样，线搜索必须满足 Wolfe 或强 Wolfe 条件，以保证 \(s^Ty>0\)、更新稳定。

Algorithm 9.2（L-BFGS）：选 \(x_0\) 和整数 \(m>0\)；\(k\leftarrow0\)。重复：选 \(H_k^0\)（如 (9.6)）；用 Algorithm 9.1 算 \(p_k=-H_k\nabla f_k\)；\(x_{k+1}=x_k+\alpha_kp_k\)，\(\alpha_k\) 满足 Wolfe 条件；若 \(k>m\)，丢弃 \(\{s_{k-m},y_{k-m}\}\)；计算并保存 \(s_k,y_k\)；\(k\leftarrow k+1\)；直到收敛。

在前 \(m-1\) 次迭代中，若 \(H_k^0\) 始终等于同一个 \(H_0\)，L-BFGS 与 BFGS 完全等价。理论上把 \(m\) 设得很大就能重现 BFGS，但当 \(m>n/2\) 时这比直接做 BFGS 还贵。"保留最近 \(m\) 对"的策略在实践中效果很好，还没有找到一致更好的策略。

9.2.4 记忆长度 \(m\) 的取舍

原书表 9.1 在 CUTE 测试集上比较了 \(m=3,5,17,29\)（终止准则 \(\|\nabla f_k\|\le10^{-5}\)），节选如下（nfg 为函数/梯度求值次数，时间为 CPU 秒）：

结论：\(m\) 太小时稳健性差（FREUROTH 失败）；\(m\) 增大则函数求值次数下降，但每步代价上升，最佳 CPU 时间常出现在较小的 \(m\)；最优 \(m\) 依问题而定。

9.2.5 优缺点

• 真 Hessian 不稀疏的大问题常首选 L-BFGS（计算并分解真 Hessian 不可行）。它甚至可能胜过用有限差分或自动微分计算 Hessian–向量积的 Newton–CG；经验上比非线性共轭梯度法更快、更稳健。
• 若 Hessian 稠密但目标部分可分（9.5 节），利用该结构的方法在函数求值次数上远胜 L-BFGS；但按总计算时间，L-BFGS 可能因为单步便宜而更有效。
• 主要弱点：常常收敛慢、函数求值多；在高度病态的问题（Hessian 特征值分布很宽）上效率低。量化中常见的病态来源是变量量级悬殊、因子模型中"因子方向很硬、特质方向很软"，见 9.7 节。

金融直觉："硬方向、软方向"可以这样理解。在均值–方差目标里，把权重沿市场因子方向（比如所有股票同时加仓）移动一点，组合方差会急剧上升，这个方向的曲率很大，叫"硬"；而在两只特质风险都很小、因子暴露又相同的股票之间做多空对调，方差几乎不变，曲率很小，叫"软"。L-BFGS 只记得最近 \(m\) 个方向的曲率，面对几千个软硬程度相差上万倍的方向，它记不全，只能一边走一边重新试探，所以迭代多。

9.2.6 与共轭梯度法的关系

有限内存方法最初就是为改进非线性共轭梯度法（第 05 章 5.3 节）而提出的。Hestenes–Stiefel 非线性 CG 的方向可以写成

\(\hat H_{k+1}\) 既不对称也不正定（实际上是奇异的）。既对称正定、又满足割线方程的"最近修正"是

正好是对单位阵做一次 BFGS 更新。以此为方向的方法叫无记忆 BFGS（memoryless BFGS）：每次更新前把近似重置为 \(I\)，只保留最近一对向量；等价于 Algorithm 9.2 取 \(m=1\)、\(H_k^0=I\)。若配合精确线搜索（\(\nabla f_{k+1}^Tp_k=0\)），

这正是 Hestenes–Stiefel CG；而在 \(\nabla f_{k+1}^Tp_k=0\) 时 HS 公式又退化为 Polak–Ribière 公式。

推导拆解：(9.7) 的第二个等号与 (9.9) 的由来。记 \(g=\nabla f_{k+1}\)，并利用 \(s_k=\alpha_kp_k\)。

• (9.7)：\(\dfrac{g^Ty_k}{y_k^Tp_k}p_k=\dfrac{g^Ty_k}{y_k^Ts_k}s_k=\dfrac{s_ky_k^T}{y_k^Ts_k}g\)（分子分母同乘 \(\alpha_k\)；\(s_k(y_k^Tg)\) 可以改写成矩阵 \(s_ky_k^T\) 乘 \(g\)）。于是 \(p_{k+1}=-(I-\frac{s_ky_k^T}{y_k^Ts_k})g\)。
• (9.9)：把 (9.8) 乘到 \(g\) 上。精确线搜索给出 \(s_k^Tg=0\)，于是 \((I-\rho y_ks_k^T)g=g\)，最后一项 \(\rho s_ks_k^Tg=0\)，只剩 \((I-\rho s_ky_k^T)g\)，与 (9.7) 完全一样。 白话：无记忆 BFGS 就是"每步从单位矩阵重新开始，只吸收最近一步的曲率"。在精确线搜索下，它和共轭梯度法走出完全相同的路径。L-BFGS 多记几步，所以比 CG 更稳。

于是：最有效的拟牛顿更新（BFGS）与最有效的非线性 CG（PR、HS）是同一件事的两种看法，L-BFGS 用额外的内存保存更多曲率信息，是 CG 的自然推广。

各方法的存储需求对比：Fletcher–Reeves \(3n\)；Polak–Ribière \(4n\)；Harwell VA14 \(6n\)；CONMIN \(7n\)；L-BFGS \(2mn+4n\)。大体上存储越多，函数求值效率与稳健性越高。

9.3 一般有限内存更新：紧凑表示

L-BFGS 是更新逆矩阵 \(H_k\) 的线搜索方法。信赖域方法需要的是 \(B_k\)；约束优化（L-BFGS-B、SQP）需要反复把 \(B_k\) 投影到约束定义的子空间；我们还希望有基于 SR1 的有限内存方法。这些需求都可以用紧凑表示（compact representation，也叫外积表示）统一满足。

定理 9.1（BFGS 的紧凑表示）：\(B_0\) 对称正定，\(k\) 对向量 \(\{s_i,y_i\}_{i=0}^{k-1}\) 满足 \(s_i^Ty_i>0\)，\(B_k\) 由 \(k\) 次 BFGS 更新得到，则

其中 \(S_k=[s_0,\dots,s_{k-1}]\)、\(Y_k=[y_0,\dots,y_{k-1}]\) 是 \(n\times k\) 矩阵，

可用归纳法证明；\(s_i^Ty_i>0\) 保证中间矩阵非奇异。

有限内存版本：取 \(B_k^0=\delta_kI\)（\(\delta_k=1/\gamma_k\)），\(S_k,Y_k\) 只保留最近 \(m\) 列，

核心图像是：对角基础矩阵 + 两个"长而窄"的 \(n\times2m\) 矩阵夹着一个 \(2m\times2m\) 小矩阵。更新表示约需 \(2mn+O(m^3)\) 运算，矩阵–向量积 \(B_kv\) 约需 \((4m+1)n+O(m^2)\) 次乘法。小矩阵的分解代价可以忽略。

金融直觉：(9.15) 和因子模型协方差是同一个形状。对照一下：

• 因子模型：\(\Sigma=D+BFB^T\)，\(D\) 对角（特质方差），\(B\) 是 \(n\times K\)（因子暴露），\(F\) 是 \(K\times K\)（因子协方差）。
• 紧凑表示：\(B_k=\delta I-WM^{-1}W^T\)，\(\delta I\) 对角，\(W=[\delta S_k\ Y_k]\) 是 \(n\times2m\)，\(M\) 是 \(2m\times2m\)。 区别只是中间那项是减号、\(M\) 不一定正定。计算 \(B_kv\) 的方法也和计算 \(\Sigma w\) 一样：先算 \(W^Tv\)（\(2m\) 个内积），再乘小矩阵 \(M^{-1}\)，再乘回 \(W\)，最后加上 \(\delta v\)。 求这类矩阵的逆，靠的是 Woodbury 公式：\((D+UCU^T)^{-1}=D^{-1}-D^{-1}U(C^{-1}+U^TD^{-1}U)^{-1}U^TD^{-1}\)。风险模型里算 \(\Sigma^{-1}\mu\)（均值–方差最优权重）时也用它，只需求一个 \(K\times K\) 的逆。

定理 9.2（SR1 的紧凑表示）：对称 \(B_0\) 经 \(k\) 次良定的 SR1 更新，

SR1 自对偶，逆公式只需把 \(B,s,y\) 换成 \(H,y,s\)。

为什么不直接"展开"？ BFGS 也可写成 \(B_{k+1}=B_k-a_ka_k^T+b_kb_k^T\)，\(a_k=B_ks_k/(s_k^TB_ks_k)^{1/2}\)，\(b_k=y_k/(y_k^Ts_k)^{1/2}\)，于是 \(B_k=B_k^0+\sum_i(b_ib_i^T-a_ia_i^T)\)。但 \(a_i\) 依赖于被删除的最旧那一对，每次都得重算，代价约 \(\tfrac32m^2n\)，而紧凑表示只要 \(2mn\)。

应用：信赖域无约束优化；更重要的是约束优化——L-BFGS-B（scipy.optimize.minimize(method='L-BFGS-B') 的底层 Fortran 代码）大量使用紧凑表示求解大规模界约束问题；SQP 可以用它近似 Lagrange 函数的 Hessian。逆矩阵 \(H_k\) 也有类似的紧凑表示，基于它的 L-BFGS 与双循环计算量相当，但能用 BLAS-2 矩阵–向量运算，在分层存储的计算机上更快。

9.4 稀疏拟牛顿更新（简介）

思路：要求 \(B_k\) 与真 Hessian 有相同的稀疏模式 \(\Omega\)，在满足割线方程的前提下离 \(B_k\) 最近：

其解可通过解一个稀疏模式为 \(\Omega\) 的线性方程组得到，但不保证正定、不具备尺度不变性，更重要的是实际表现令人失望：所需函数求值至少与 L-BFGS 一样多，单步却更贵。根本原因是这个变分模型不充分，产生的近似质量差。放松为"沿最近几步近似满足割线方程"（\(\min\|BS_k-Y_k\|_F\)）常有改善，但大规模表现仍不突出。结论：了解即可，实务中不用。

9.5 部分可分函数

9.5.1 从可分到部分可分

可分函数如 \(f(x)=f_1(x_1,x_3)+f_2(x_2,x_4,x_6)+f_3(x_5)\)，每个变量只出现在一个函数中，可以拆成独立的小问题。部分可分函数不可拆，但可以写成若干元素函数（element functions）之和，

每个 \(f_i\) 沿大量线性无关方向不变。于是 \(\nabla f=\sum\nabla f_i\)，\(\nabla^2f=\sum\nabla^2f_i\)。关键问题是：分别维护各元素 Hessian 的拟牛顿近似，是否比整体近似更好？答案是肯定的，前提是充分利用每个元素 Hessian 的结构。

例子：

\(f_1\) 只依赖元素变量 \(x_{[1]}=(x_1,x_3)^T=U_1x\)，\(U_1=\begin{bmatrix}1&0&0&0\\0&0&1&0\end{bmatrix}\)。令 \(\phi_1(z_1,z_2)=(z_1-z_2^2)^2\)，则 \(f_1(x)=\phi_1(U_1x)\)，

（精读笔记记录原书此处右下元素印为 \(12x_3^2\)；按 \(\phi_1=(z_1-z_2^2)^2\) 直接求导应为 \(12z_2^2-4z_1\)，即 \(12x_3^2-4x_1\)。）

关键思想：不去近似 \(4\times4\) 的稀疏奇异矩阵 \(\nabla^2f_1\)，而是维护 \(2\times2\) 的 \(B_{[1]}\approx\nabla^2\phi_1\)，每步记录元素级的

用 BFGS 或 SR1 更新。总近似是 \(B=\sum_iU_i^TB_{[i]}U_i\)。\(2\times2\) 矩阵几步之内就被采样充分，因此这个近似比忽略结构的整体近似好得多。

推导拆解：(9.23) 和 \(U_i^TB_{[i]}U_i\) 在做什么。

1. \(U_1x\) 就是"从 \(x\) 中挑出第 1 和第 3 个分量"：\(U_1x=(x_1,x_3)^T\)。\(U_1\) 是一个"选择矩阵"。
2. 链式法则：\(f_1(x)=\phi_1(U_1x)\) 对 \(x\) 求梯度，先对 \(\phi_1\) 的两个参数求导得 \(\nabla\phi_1\)（长度 2），再乘以"内层对 \(x\) 的导数" \(U_1^T\)，就是把这两个数放回 4 维向量的第 1、3 位置，其余补零。
3. Hessian 同理：\(U_1^T\nabla^2\phi_1U_1\) 是把 \(2\times2\) 小矩阵"嵌入"到 \(4\times4\) 矩阵的 \((1,3)\) 行列交叉位置，其余为零。
4. \(\nabla^2\phi_1\) 本身：\(\phi_1=(z_1-z_2^2)^2\)，\(\partial\phi_1/\partial z_1=2(z_1-z_2^2)\)，\(\partial\phi_1/\partial z_2=-4z_2(z_1-z_2^2)\)；再求一次导得 \(\partial^2/\partial z_1^2=2\)，\(\partial^2/\partial z_1\partial z_2=-4z_2\)，\(\partial^2/\partial z_2^2=-4(z_1-z_2^2)+8z_2^2=12z_2^2-4z_1\)，与正文的更正一致。 白话解释：这就像风险系统按交易台汇总风险。每个交易台只管自己那几个风险因子的小矩阵，总部把各台的小矩阵按因子编号"摆进"全局大表再相加。每个小矩阵只有几个参数，几次观测就能估准；直接估一张全局大表则要多得多的数据。

9.5.2 内部变量与最细分解

分解方式往往不唯一，应选最细分解。原书的例子是最小曲面问题，元素函数

依赖 4 个分量，但它在子空间

上不变（\(\dim N_i=n-2\)），实际只沿两个方向变化。定义内部变量 \(u_j=x_j-x_{j+q+1}\)、\(u_{j+1}=x_{j+1}-x_{j+q}\) 和内部函数 \(\phi_i(u_j,u_{j+1})=\frac1{q^2}[1+\frac{q^2}2(u_j^2+u_{j+1}^2)]^{1/2}\)，就得到一般形式

其中 \(U_i\) 是 \(n_i\times n\) 矩阵（\(n_i\ll n\)），其零空间正是不变子空间 \(N_i\)。

9.5.3 不变子空间与定义

定义 9.1：函数 \(f\) 的不变子空间是使 \(f(x+w)=f(x)\) 对所有 \(w\) 成立的最大子空间。

对比两个例子：\(f_i=x_{50}^2\) 的不变子空间是 \(\{w:w_{50}=0\}\)，\(U_i=e_{50}^T\)；\(f_i=(x_1+\dots+x_n)^2\) 的梯度和 Hessian 完全稠密，但不变子空间是 \(\{w:e^Tw=0\}\)，维数也是 \(n-1\)，\(U_i=(1,\dots,1)\)。两者的内部函数都是 \(\phi(z)=z^2\)，只是 \(U_i\) 不同。**稠密 Hessian 也可以有极好的部分可分结构。**这一点在量化里非常有用：预算约束罚项 \((\mathbf 1^Tw-1)^2\)、净敞口罚项 \((\mathbf 1^Tw)^2\) 的 Hessian 是全 1 矩阵，但它只是秩一的。

白话解释：不变子空间就是"怎么动都不影响这个函数的方向集合"。\((\mathbf1^Tw)^2\) 只关心净敞口，所以任何"多一只、空一只、金额相等"的调仓（\(\mathbf1^Tw\) 不变）都不改变它的值，这些调仓方向构成一个 \(n-1\) 维的子空间。剩下"唯一真正起作用"的方向只有 1 个（整体加减仓），所以内部变量只有一个：\(z=\mathbf1^Tw\)。 "稀疏"和"部分可分"的区别在于：稀疏看的是"依赖哪几个坐标"，部分可分看的是"依赖几个方向（可以是坐标的组合）"。后者更一般，能把稀疏抓不住的低秩结构也抓住。 另：9.6 节提到的"(9.35) 型"结构，沿用原书编号，指的就是这里 \((x_1+\cdots+x_n)^2\) 这类"Hessian 稠密但内部只有一个变量"的元素函数；本章正文没有单独列出该式号。

定义 9.2（部分可分）：\(f\) 可写为 (9.32) 的形式，每个 \(U_i\) 的行数远小于 \(n\)。

定理 9.3：二阶连续可微、Hessian 稀疏的函数都是部分可分的。具体地，若 \(\partial^2f/\partial x_i\partial x_j\equiv0\)，则

（\(n=2\) 时的证明：对 \(\partial^2f/\partial x\partial y=0\) 在矩形上积分，得 \(f(x,y)-f(0,y)-f(x,0)+f(0,0)=0\)。）所以部分可分比稀疏更一般。

推导拆解：\(n=2\) 的积分分两步做。

1. 先对 \(y\) 积分：\(\int_0^y\frac{\partial^2f}{\partial x\partial y}(x,t)\,dt=\frac{\partial f}{\partial x}(x,y)-\frac{\partial f}{\partial x}(x,0)\)（微积分基本定理：导数的积分等于原函数之差）。左边被积函数为零，所以 \(\partial f/\partial x\) 不依赖 \(y\)。
2. 再对 \(x\) 从 0 积到 \(x\)：\([f(x,y)-f(0,y)]-[f(x,0)-f(0,0)]=0\)。 移项得 \(f(x,y)=f(0,y)+f(x,0)-f(0,0)\)：一个只依赖 \(y\) 的函数加一个只依赖 \(x\) 的函数（再加常数）。也就是说，"交叉二阶导为零"⇔"没有交互作用"⇔"可以拆开"。这和回归里"没有交互项时，两个解释变量的效应可以分开相加"是一回事。

群部分可分：非线性最小二乘 \(\sum_k(f_k+f_{k+1}+c_k)^2\) 这类形式，用 \(f(x)=\sum_k\psi_k(h_k(x))\) 描述，\(h_k\) 部分可分、\(\psi_k\) 是标量群函数。链式法则 \(\nabla^2[\psi(h)]=\psi''(h)\nabla h\nabla h^T+\psi'(h)\nabla^2h\) 使前述思想全部可以推广，适用于最小二乘、罚函数和价值函数。

9.5.4 部分可分函数的算法

在牛顿法中：截断牛顿法用 CG 近似解 \(\nabla^2f(x_k)p=-\nabla f(x_k)\)，CG 只需 Hessian–向量积，可以按 (9.33) 逐元素计算，不必组装 Hessian。最小曲面例中每个 \(2\times2\) 内部 Hessian 存 3 个数，总共约 \(3n\)，比存全 Hessian 下三角的约 \(5n\) 节省 40%。代价是 \(U_i\) 映射带来的内存访问开销，有时显著。LANCELOT 软件包允许在 CG 与多波前直接分解之间选择。

在拟牛顿法中：存储并更新各内部 Hessian 的近似 \(B_{[i]}\)，组装 \(B=\sum U_i^TB_{[i]}U_i\)，在信赖域框架中用 CG 近似求解。每个元素满足自己的割线方程 \(B_{[i]}s_{[i]}=y_{[i]}\)。小矩阵几步就被采样充分；而忽略结构的 BFGS/L-BFGS 要用一个稠密的 \(n\times n\) 矩阵估计所有元素曲率之和的"总平均曲率"，\(n\) 大时需要多得多的迭代。

曲率条件问题：即使 \(\nabla^2f(x^*)\) 正定，单个元素 Hessian 也可能不定，\(s_{[i]}^Ty_{[i]}>0\) 不保证，BFGS 不一定可用。对策是对每个元素用 SR1（加跳过规则）——这正是第 08 章说 SR1 适合部分可分优化的原因。若找到了最细分解，基于 SR1 的部分可分拟牛顿法性能常与牛顿法相当。

金融直觉：整体是凸的，局部却可以不凸，就像一个整体风险可控的账簿里，单个交易可以是卖出期权（负 Gamma）。例如 \(\phi_1=(z_1-z_2^2)^2\) 在 \(z_1>3z_2^2\) 的区域，\(\nabla^2\phi_1\) 的右下元 \(12z_2^2-4z_1<0\)，这个元素函数本身就是不定的；但加上其他元素函数之后，总 Hessian 仍可能正定。BFGS 要求每个被更新的矩阵都正定，相当于要求"每个交易都是正 Gamma"，太苛刻；SR1 允许单个元素不定，只要总和合理即可。

部分可分拟牛顿法不具有一般仿射不变性（一般仿射变换会破坏可分结构），但对保持可分结构的变换不变，这不算缺陷。

9.6 本章方法在量化中的位置

• 大规模参数估计：数千到数万个参数的模型——大截面上的非线性因子模型估计、带大量哑变量的面板 MLE、logistic/softmax 形式的机器学习因子模型——标准选择就是 L-BFGS。带上下界的参数（波动率参数非负、权重上下限、相关系数在 \([-1,1]\) 内）用 L-BFGS-B。
• 病态是 L-BFGS 变慢的头号原因：未标准化的因子暴露、量级悬殊的参数、因子模型协方差的"硬方向"都会让迭代次数暴涨。先做变量缩放（标准化特征、换单位）往往比调 \(m\) 更有效。
• 组合优化中的结构：因子模型协方差 \(\Sigma=BFB^T+D\) 恰好是"对角 + 低秩"——与紧凑表示 (9.15) 同构。计算 \(\Sigma w=B(F(B^Tw))+Dw\) 只需 \(O(nK)\)，存储只需 \(O(nK)\)，永远不要组装稠密的 \(n\times n\) 矩阵。\((\mathbf 1^Tw)^2\) 这类罚项是 (9.35) 型的"稠密但秩一"结构；按资产独立的交易成本是可分项；多期优化的目标通常是群部分可分的。
• 部分可分拟牛顿法（LANCELOT 一类）在量化日常工作中较少直接使用，属于"知道它存在、遇到超大规模结构化问题时再回原书"的内容。

9.7 量化实战

9.7.1 代码一：手写 L-BFGS 双循环，观察 \(m\) 的影响

下面实现 Algorithm 9.1–9.2（初始矩阵用 (9.6)，强 Wolfe 线搜索），在两个 \(n=1000\) 的问题上比较不同 \(m\)：原书习题 9.1 的扩展 Rosenbrock 函数 \(f(x)=\sum_{i=1}^{n/2}[100(x_{2i}-x_{2i-1}^2)^2+(1-x_{2i-1})^2]\)，以及一个特征值在 \([1,10^4]\) 上对数均匀分布的病态凸二次函数。

import numpy as np, warnings
warnings.simplefilter('ignore')
from scipy.optimize import line_search

def two_loop(g, S, Y, gamma):
    """Algorithm 9.1：返回 H_k g，其中 H_k^0 = gamma*I"""
    q = g.copy(); hist = []
    for s, y in zip(reversed(S), reversed(Y)):            # i = k-1, ..., k-m
        rho = 1.0/(y @ s); a = rho*(s @ q); q -= a*y; hist.append((rho, a))
    r = gamma*q
    for (s, y), (rho, a) in zip(zip(S, Y), reversed(hist)):   # i = k-m, ..., k-1
        b = rho*(y @ r); r += s*(a - b)
    return r

def lbfgs(fun, x0, m=5, tol=1e-5, maxit=20000):
    cnt = [0]; last = {}
    def fg(z):                         # 缓存：同一点的 f 与 g 只算一次
        if 'z' in last and np.array_equal(z, last['z']): return last['v']
        cnt[0] += 1; last['z'] = z.copy(); last['v'] = fun(z); return last['v']
    x = x0.copy(); f, g = fg(x); S, Y = [], []; k = 0
    while np.linalg.norm(g) > tol and k < maxit:
        gamma = (S[-1] @ Y[-1])/(Y[-1] @ Y[-1]) if S else 1.0/np.linalg.norm(g)   # (9.6)
        p = -two_loop(g, S, Y, gamma)
        a, *_ = line_search(lambda z: fg(z)[0], lambda z: fg(z)[1], x, p,
                            gfk=g, old_fval=f, c1=1e-4, c2=0.9)
        if a is None:                  # 线搜索失败：丢弃记忆，改走最速下降方向重试一次
            S.clear(); Y.clear(); p = -g/np.linalg.norm(g)
            a, *_ = line_search(lambda z: fg(z)[0], lambda z: fg(z)[1], x, p, gfk=g, old_fval=f)
            if a is None: break
        xn = x + a*p; fn, gn = fg(xn)
        S.append(xn - x); Y.append(gn - g)
        if len(S) > m: S.pop(0); Y.pop(0)                    # 丢弃最旧的一对
        x, f, g = xn, fn, gn; k += 1
    return x, k, cnt[0]

# 测试 1：扩展 Rosenbrock（原书习题 9.1），n = 1000
def ext_rosen(x, a=100.0):
    xo, xe = x[0::2], x[1::2]; t1 = xe - xo**2; t2 = 1 - xo
    g = np.empty_like(x); g[0::2] = -4*a*xo*t1 - 2*t2; g[1::2] = 2*a*t1
    return np.sum(a*t1**2 + t2**2), g

# 测试 2：病态凸二次函数，特征值在 [1, 1e4] 上对数均匀分布，n = 1000
rng = np.random.default_rng(0)
n = 1000
Q, _ = np.linalg.qr(rng.standard_normal((n, n)))
A = (Q * np.logspace(0, 4, n)) @ Q.T
b = rng.standard_normal(n)
def quad(x): return 0.5*x @ A @ x - b @ x, A @ x - b

print("扩展 Rosenbrock，起点 (-1,...,-1)：")
for m in [1, 3, 5, 17, 29]:
    x, k, nfg = lbfgs(ext_rosen, -np.ones(n), m=m)
    print(f"  m={m:2d}: 迭代 {k:5d}, 函数/梯度求值 {nfg:5d}, ||x-1||={np.linalg.norm(x-1):.1e}")
print("病态二次函数（条件数 1e4）：")
for m in [1, 3, 5, 17, 29]:
    x, k, nfg = lbfgs(quad, np.zeros(n), m=m)
    print(f"  m={m:2d}: 迭代 {k:5d}, 函数/梯度求值 {nfg:5d}")

输出：

扩展 Rosenbrock，起点 (-1,...,-1)：
  m= 1: 迭代    43, 函数/梯度求值    69, ||x-1||=2.8e-10
  m= 3: 迭代    33, 函数/梯度求值    41, ||x-1||=2.7e-09
  m= 5: 迭代    31, 函数/梯度求值    44, ||x-1||=1.8e-08
  m=17: 迭代    31, 函数/梯度求值    41, ||x-1||=4.9e-09
  m=29: 迭代    31, 函数/梯度求值    41, ||x-1||=4.9e-09
病态二次函数（条件数 1e4）：
  m= 1: 迭代   919, 函数/梯度求值  1042
  m= 3: 迭代   949, 函数/梯度求值  1020
  m= 5: 迭代   792, 函数/梯度求值   854
  m=17: 迭代   744, 函数/梯度求值   791
  m=29: 迭代   754, 函数/梯度求值   802

解读：扩展 Rosenbrock 由 500 个互不相关的二维块组成（可分！），曲率结构简单，\(m\ge3\) 后迭代次数就不再下降，1000 维问题只需约 40 次函数求值。病态二次函数上 L-BFGS 需要数百次迭代，增加 \(m\) 只能带来约 20% 的改善——这正是 9.2.5 节说的"高度病态问题上效率低"。在这种问题上，有效的手段是预条件（通过 \(H_k^0\) 接口传入）或变量缩放，而不是一味增大 \(m\)。

9.7.2 代码二：用因子结构做 3000 只股票的带界组合优化

问题：3000 只股票、10 个因子，求解多空组合

\(\Sigma=BFB^T+D\)。净敞口用罚项近似控制（严格的等式约束要到原书第 16–18 章的方法）。对比两种梯度实现：利用因子结构的 \(O(nK)\) 版本与组装稠密 \(\Sigma\) 的 \(O(n^2)\) 版本，都交给 L-BFGS-B。

import numpy as np, time
from scipy.optimize import minimize

rng = np.random.default_rng(1)
n, K = 3000, 10                                   # 3000 只股票，10 个因子
B = rng.standard_normal((n, K)) * 0.8             # 因子暴露
L = rng.standard_normal((K, K)) * 0.02
F = L @ L.T + np.diag(np.full(K, 4e-4))           # 因子协方差（日频量级）
d = rng.uniform(1e-4, 9e-4, n)                    # 特质方差
mu = rng.standard_normal(n) * 2e-4                # alpha 预测（日频）
gamma, kappa, ub = 50.0, 1.0, 0.005               # 风险厌恶、净敞口罚系数、单票上限
S = 1e4                                           # 目标乘 1e4（以基点计），避免相对容差失效

def obj_factor(w):                                # 利用 Σ = B F B' + D：每次 O(nK)
    Sw = B @ (F @ (B.T @ w)) + d * w
    net = w.sum()
    f = 0.5*gamma*(w @ Sw) - mu @ w + 0.5*kappa*net**2
    g = gamma*Sw - mu + kappa*net                 # (1'w)^2 项：Hessian 稠密但只是秩一
    return S*f, S*g

Sigma = B @ F @ B.T + np.diag(d)                  # 稠密 n×n，仅作对照：每次 O(n^2)
def obj_dense(w):
    Sw = Sigma @ w; net = w.sum()
    f = 0.5*gamma*(w @ Sw) - mu @ w + 0.5*kappa*net**2
    return S*f, S*(gamma*Sw - mu + kappa*net)

ev = np.linalg.eigvalsh(gamma*Sigma + kappa*np.ones((n, n)))
print("Hessian 条件数: %.1e" % (ev[-1]/ev[0]))
bounds = [(-ub, ub)] * n
sol = {}
for name, fun in [("因子结构", obj_factor), ("稠密协方差", obj_dense)]:
    t0 = time.perf_counter()
    res = minimize(fun, np.zeros(n), jac=True, method='L-BFGS-B', bounds=bounds,
                   options={'maxcor': 10, 'ftol': 1e-14, 'gtol': 1e-9, 'maxiter': 20000})
    t = time.perf_counter() - t0
    sol[name] = res.x
    print(f"{name}: 迭代 {res.nit}, 求值 {res.nfev}, 用时 {t:.2f}s, "
          f"每次求值 {1e3*t/res.nfev:.2f}ms, 目标 {res.fun/S:.8e}")
w = sol["因子结构"]
print("两种实现的解之差 max|Δw| = %.1e" % np.max(np.abs(w - sol["稠密协方差"])))
print("顶到上下限的股票数:", int(np.sum(np.abs(w) > ub*(1-1e-6))),
      "| 净敞口 %.1e | 年化波动 %.2f%%" % (w.sum(), 100*np.sqrt(252*w @ Sigma @ w)))
print("存储: 因子形式 %.1f 万个数, 稠密形式 %.0f 万个数" % ((n*K + K*K + n)/1e4, n*n/1e4))

输出：

Hessian 条件数: 6.0e+05
因子结构: 迭代 1605, 求值 1640, 用时 0.31s, 每次求值 0.19ms, 目标 -1.64359211e-03
稠密协方差: 迭代 1613, 求值 1652, 用时 0.78s, 每次求值 0.47ms, 目标 -1.64359211e-03
两种实现的解之差 max|Δw| = 1.1e-06
顶到上下限的股票数: 1651 | 净敞口 -3.7e-06 | 年化波动 7.75%
存储: 因子形式 3.3 万个数, 稠密形式 900 万个数

三点观察。第一，两种实现得到同一个解，但因子形式的存储只有稠密形式的 0.4%；\(n=3000\) 时稠密矩阵–向量乘法借助 BLAS 仍然很快，但 \(n=3\times10^4\) 时稠密 \(\Sigma\) 就要 7.2 GB 内存，而因子形式只增长到 33 万个数。第二，Hessian 条件数高达 \(6\times10^5\)：因子方向的特征值（量级约为 \(\gamma\cdot n\bar\beta^2\sigma_f^2\)，本例为几十到约 2000）和净敞口罚项方向的特征值（\(\kappa n=3000\)）都很大，而特质方向的特征值（\(\gamma\sigma_\varepsilon^2\)，仅 0.005–0.045）很小，两者相差悬殊，所以 L-BFGS-B 需要约 1600 次迭代，这就是 9.2.5 节所说的弱点在组合优化中的具体表现。第三，最初版本的代码没有把目标乘 \(10^4\)，日频目标值只有 \(10^{-3}\) 量级，L-BFGS-B 的相对下降准则（ftol）提前触发，两种实现停在了不同的点——目标函数的量级本身也是"缩放"问题的一部分。对这种规模和结构的问题，更好的工具是二次规划内点法或专门的组合优化器（原书第 14、16 章）；L-BFGS-B 的优势在于目标非二次、只能拿到梯度的场合。

本章小结

大规模问题中稠密拟牛顿矩阵不可承受。L-BFGS 只保存最近 \(m\) 对 \((s,y)\)，用双循环递归以 \(4mn\) 次乘法算出 \(H_k\nabla f_k\)，配合 \(H_k^0=\gamma_kI\) 和 Wolfe 线搜索，是大规模无约束优化的默认选择；它是无记忆 BFGS（≈ 非线性 CG）的推广，主要弱点是病态问题上收敛慢。紧凑表示把有限内存 BFGS/SR1 写成"对角 + 低秩"的形式，是 L-BFGS-B 等约束优化器的基础。稀疏拟牛顿更新效果差；部分可分结构更一般、更有用，按元素维护小矩阵（常用 SR1）可接近牛顿法的性能。量化中，因子模型协方差的"对角 + 低秩"结构是同一思想最直接的应用。

练习

基础

1. 证明 (9.7) 中的 \(\hat H_{k+1}=I-\dfrac{s_ky_k^T}{y_k^Ts_k}\) 是奇异矩阵。 提示：找一个非零向量 \(v\) 使 \(\hat H_{k+1}^Tv=0\)，试 \(v=y_k\)。
2. 在精确线搜索（\(\nabla f_{k+1}^Ts_k=0\)）下，从 (9.8) 推导 (9.9)。
3. 把 \(f(x)=x_2x_3e^{x_1+x_3-x_4}+(x_2x_3)^2+(x_3-x_4)\) 写成 \(\sum\phi_i(U_ix)\) 的形式，给出各 \(U_i\)，并尽量减少内部变量个数。 提示：第一项可用内部变量 \((x_2,x_3,x_1+x_3-x_4)\)；第二项只依赖 \(x_2x_3\) 但作为函数依赖 \((x_2,x_3)\)；第三项依赖 \(x_3-x_4\) 一个内部变量。
4. 证明 (9.28) 的维数为 \(n-2\)，\(\{w:e^Tw=0\}\) 的维数为 \(n-1\)。
5. 在 9.7.1 的代码中，把 \(H_k^0\) 固定为 \(I\)（不用 (9.6)），比较扩展 Rosenbrock 上的函数求值次数。

进阶

6. 计算紧凑表示 (9.10) 的 \(B_kv\) 时，哪些步骤可以并行？与双循环递归的并行度比较（原书习题 9.4）。 提示：\([B_0S\ Y]^Tv\) 是 \(2m\) 个独立内积；双循环的每一步依赖上一步的 \(q\) 或 \(r\)，本质上串行。
7. 部分可分拟牛顿近似 \(B=\sum U_i^TB_{[i]}U_i\)，各元素满足 \(B_{[i]}s_{[i]}=y_{[i]}\)，问 \(B\) 是否满足整体割线方程 \(Bs=y\)？（原书习题 9.11） 提示：\(s_{[i]}=U_is\)，\(\sum U_i^Ty_{[i]}=\nabla f(x^+)-\nabla f(x)\)。
8. 在 9.7.2 中，把单票上限放宽到 2% 和收紧到 0.2%，观察迭代次数与顶格股票数的变化，解释界约束活跃集对 L-BFGS-B 的影响。
9. 用 \(D^{-1/2}\) 对 9.7.2 的变量做对角缩放（令 \(w=D^{-1/2}v\)），重写目标和界约束，观察条件数和迭代次数的变化。为什么对角缩放对"因子方向过硬"的病态帮助有限？ 提示：对角缩放只能均衡特质方差的差异，因子方向的大特征值来自低秩项 \(BFB^T\)，不是对角元的问题。

原书推荐习题：9.1（亲手实现 L-BFGS 并观察 \(m\) 的影响）、9.3（无记忆 BFGS 与 CG 的等价）、9.4、9.5（紧凑表示的计算与存储技巧）、9.7、9.9（识别不变子空间与内部变量）。原书第 9 章共有 9.1–9.13 题。

原书对照

页码换算：本段原书正文页码约等于 PDF 页码减 20。