第 10 章非线性最小二乘

学习目标

读完本章，你应当能够：

写出最小二乘目标的梯度 $J^Tr$ 与 Hessian $J^TJ+\sum r_j\nabla^2r_j$，说明"知道 Jacobian 就免费得到 Hessian 的主要部分"这一结构为什么是所有算法的出发点，并能从正态误差的最大似然推出最小二乘。
比较求解线性最小二乘的正规方程/Cholesky、QR、SVD 三种方法的精度与代价，解释"条件数平方"效应，并用 SVD 理解共线因子回归的不稳定性。
推导 Gauss–Newton 方向，证明它在 Jacobian 满秩时是下降方向，并用 $\|(J^TJ)^{-1}H\|$ 判断它的局部收敛速度（小残差快、零残差二次、大残差可能不收敛）。
把 Levenberg–Marquardt 理解为"Gauss–Newton + 信赖域"，写出 $(J^TJ+\lambda D^2)p=-J^Tr$ 及其增广最小二乘形式，并看出它与岭回归是同一个公式。
知道大残差问题的混合方法（Fletcher–Xu、NL2SOL）、大规模问题的处理方式以及正交距离回归的用途。
能用 scipy.optimize.least_squares 和自己写的 GN/LM 拟合收益率曲线，并识别参数不可识别和多个局部极小的问题。

读前导读

这一章在解决什么问题

你在 CFA 里学过的 OLS 回归，就是"选系数，使残差平方和最小"。那是线性最小二乘：模型对参数是线性的，有现成的公式 $\hat\beta=(X^TX)^{-1}X^Ty$。本章处理更一般的情形：模型对参数是非线性的。Nelson–Siegel 收益率曲线里的衰减参数、SABR/Heston 对期权报价的校准、冲击成本模型 $c=a\cdot(\text{成交量占比})^b$ 中的指数 $b$，都让残差平方和不再是参数的简单二次函数，没有闭式解，只能迭代求解。

本章的核心发现是：最小二乘问题有一个特别好的结构。只要算出每个残差对每个参数的一阶导数（Jacobian $J$），就能"免费"得到 Hessian 的主要部分 $J^TJ$。Gauss–Newton 法每一步把模型在当前点线性化，解一个普通的 OLS 问题来决定往哪走；Levenberg–Marquardt（LM）法在此基础上加一个"别走太远"的约束，结果是 $(J^TJ+\lambda I)p=-J^Tr$，和岭回归是同一个公式。LM 是业界校准模型的标准工具。

本章还会回答一个 CFA 只点到为止的问题：因子高度共线时，回归系数为什么不稳定？用 SVD 可以精确地看出来是哪个方向、被放大了多少倍。

需要先想起来的数学

1. Jacobian 与梯度。 有 $m$ 个残差、$n$ 个参数时，Jacobian $J$ 是 $m\times n$ 矩阵，第 $j$ 行第 $i$ 列是 $\partial r_j/\partial x_i$，即"第 $j$ 个观测的残差对第 $i$ 个参数的敏感度"。线性回归里 $r=X\beta-y$，所以 $J=X$，就是设计矩阵。见第 00 册第 05 章多元微积分与优化。

2. 正交矩阵与范数。 正交矩阵 $Q$ 满足 $Q^TQ=I$，作用在向量上相当于旋转或反射，不改变长度：$\|Qv\|=\|v\|$。例：二维旋转 $\begin{bmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{bmatrix}$。QR 和 SVD 都靠这一性质把问题变简单。$\|v\|_2$ 就是向量的欧氏长度 $\sqrt{\sum v_i^2}$。见第 00 册第 06 章线性代数速成。

3. 奇异值分解（SVD）。 任何矩阵都可以写成 $J=U\Sigma V^T$：$V$ 的列是输入空间的一组正交方向，$U$ 的列是输出空间的一组正交方向，奇异值 $\sigma_i\ge0$ 是 $J$ 沿第 $i$ 个方向的"放大倍数"。它和 PCA 是一回事：对去均值的收益率矩阵做 SVD，$V$ 的列就是主成分方向，$\sigma_i^2/(T-1)$ 就是对应的方差。见第 00 册第 06 章。

4. 条件数。 $\kappa(J)=\sigma_{\max}/\sigma_{\min}$，衡量"输入的小误差最多被放大多少倍"。$\kappa=10^6$ 意味着数据第 16 位的舍入误差可能影响到结果的第 10 位。关键事实：$\kappa(J^TJ)=\kappa(J)^2$。见第 00 册第 06 章。

5. 最大似然。 选参数使"观测到这批数据的概率"最大；通常取对数，把连乘变成求和。CFA 统计里你见过它与 OLS 的关系：误差正态时两者一致。见第 00 册第 07 章概率中的分析工具。

怎么读这一章

核心必读：10.1（梯度与 Hessian 的结构）、10.2.3（线性最小二乘的三种解法，直接关系到你写回归代码的数值稳定性）、10.3.1 与 10.3.3（GN 的方向与收敛速度）、10.4.1–10.4.2（LM 与岭回归），以及 10.7 的两段实战。10.3.2 的全局收敛证明和 10.4.3 的 Householder/Givens 实现细节第一次可以跳过。10.5 大残差与大规模问题可以只看每段第一句。10.6 正交距离回归与"beta 估计误差导致的衰减偏误"有关，有兴趣再读。

建议顺序：10.1 → 10.2 → 10.7.2 → 10.3 → 10.4 → 10.7.3 → 10.5、10.6。

10.1 问题与结构

10.1.1 目标函数

最小二乘问题的目标是

f(x)=\frac12\sum_{j=1}^mr_j^2(x),\tag{10.1}

每个 $r_j:\mathbb R^n\to\mathbb R$ 是光滑函数，叫做残差（residual）。全章假设 $m\ge n$。原书开篇就说：凡是给化学、物理、金融、经济建立参数化模型的人，几乎都用 (10.1) 度量模型与观测之间的差异，并最小化它来选取最匹配数据的参数。量化里的收益率曲线拟合、波动率曲面校准、因子回归、冲击成本模型估计，全都属于这一类。

记残差向量 $r(x)=(r_1(x),\dots,r_m(x))^T$，于是 $f=\tfrac12\|r\|_2^2$；Jacobian $J(x)=[\partial r_j/\partial x_i]$ 是 $m\times n$ 矩阵，第 $j$ 行是 $\nabla r_j^T$。则

\nabla f(x)=\sum_{j=1}^mr_j\nabla r_j=J^Tr,\tag{10.4}

\nabla^2f(x)=J^TJ+\sum_{j=1}^mr_j\nabla^2r_j.\tag{10.5}

推导拆解：两个公式都只用了链式法则和乘积法则。

梯度：$\tfrac12r_j^2$ 对 $x$ 求导，外层 $\tfrac12u^2$ 的导数是 $u$，内层是 $\nabla r_j$，所以得 $r_j\nabla r_j$。对 $j$ 求和；而"$\sum_jr_j\times(J$ 的第 $j$ 行$)$"按矩阵乘法正好是 $J^Tr$。

Hessian：对 $r_j\nabla r_j$ 再求导，用乘积法则：一部分是"$r_j$ 在变"，得 $\nabla r_j\nabla r_j^T$；另一部分是"$\nabla r_j$ 在变"，得 $r_j\nabla^2r_j$。而 $\sum_j\nabla r_j\nabla r_j^T=J^TJ$（$J$ 的各行做外积再相加）。一元例子：$f=\tfrac12r(x)^2$，$f'=rr'$，$f''=r'^2+rr''$。$r'^2$ 只需一阶导，$rr''$ 在残差小或模型近线性时很小。统计上，$J^TJ$ 对应回归里的 $X^TX$；在正态误差下，$J^TJ/\sigma^2$ 就是 Fisher 信息矩阵，它的逆给出参数估计的协方差。

10.1.2 最小二乘的独特之处

知道 Jacobian，就"免费"得到 Hessian 的第一部分 $J^TJ$。 而且这一项往往比第二项重要得多：要么因为模型在解附近近似线性（$\nabla^2r_j$ 小），要么因为残差小（$r_j$ 小）。几乎所有最小二乘算法都利用这一结构，再嵌入线搜索或信赖域框架。

一个例外值得记住：某些大规模应用中，用自动微分直接算 $\nabla f$ 比单独算出 $J$ 更实际（例如训练神经网络时）。这时无法利用最小二乘结构，本章算法不适用，应改用第 06、09 章的一般大规模无约束方法。

10.2 背景：建模、统计与线性最小二乘

10.2.1 例 10.1 与固定回归量模型

例 10.1（药物浓度）：在服药后时刻 $t_j$ 抽血测得浓度 $y_j$，选模型

\phi(x;t)=x_1+tx_2+e^{-x_3t},\tag{10.6}

最小化 $\tfrac12\sum_j[y_j-\phi(x;t_j)]^2$，即 $r_j(x)=y_j-\phi(x;t_j)$。每一项是数据点 $(t_j,y_j)$ 到模型曲线的竖直距离的平方。这是统计学中的固定回归量模型（fixed-regressor model）：假设 $t_j$ 精确已知，误差只出现在 $y_j$ 上。若 $t_j$ 也有误差，需要 10.6 节的正交距离回归。

也可以用其他准则，比如绝对值之和 $\sum|y_j-\phi(x;t_j)|$ 或六次方和，各有统计动机；但最小二乘有最坚实的统计基础。

10.2.2 最大似然解释

设偏差 $\epsilon_j=y_j-\phi(x;t_j)$ 独立同分布，方差 $\sigma^2$，密度 $g_\sigma$。似然为

p(y;x,\sigma)=\prod_{j=1}^mg_\sigma(y_j-\phi(x;t_j)).\tag{10.9}

若 $g_\sigma$ 是正态密度 $\frac1{\sqrt{2\pi\sigma^2}}\exp(-\frac{\epsilon^2}{2\sigma^2})$，则

p=(2\pi\sigma^2)^{-m/2}\exp\Big(-\frac1{2\sigma^2}\sum_j[y_j-\phi(x;t_j)]^2\Big),

对任意固定的 $\sigma^2$，最大化 $p$ 等价于最小化平方和。**结论：误差独立同正态分布时，最大似然估计就是最小二乘解。**若误差相关或异方差，可推广为广义形式 $r^TVr$（$V$ 对称），即加权/广义最小二乘。统计背景见第 03 册第 13a 章（线性回归与最小二乘推断）和第 05 册第 19b 章（广义最小二乘）。

10.2.3 线性最小二乘

当各 $r_j$ 都是线性函数时，$J$ 为常数矩阵，

f(x)=\tfrac12\|Jx+r\|_2^2,\qquad r=r(0),\tag{10.10}

$\nabla f=J^T(Jx+r)$，$\nabla^2f=J^TJ$（第二项消失），$f$ 是凸函数（非线性问题未必）。解满足正规方程（normal equations）

J^TJx^*=-J^Tr.\tag{10.11}

以下设 $J$ 列满秩。三种主要算法：

1. 正规方程 + Cholesky。 计算 $J^TJ$ 与 $-J^Tr$，做 Cholesky 分解 $J^TJ=\bar R^T\bar R$，再两次三角回代。快且常用，但有一个根本缺陷：$J^TJ$ 的条件数是 $J$ 条件数的平方。解的相对误差通常与系数矩阵的条件数成正比，所以精度可能远差于直接处理 $J$ 的方法；$J$ 病态时舍入误差甚至会使分解过程中出现负的对角元而失败。

白话解释：先补两个概念。(10.11) 来自"梯度为零"：$\nabla f=J^T(Jx+r)=0$，移项即得；它就是 OLS 的 $X^TX\hat\beta=X^Ty$（这里 $r=-y$）。Cholesky 分解是把正定矩阵写成"上三角矩阵的转置 × 上三角矩阵"，相当于矩阵版的开平方，之后解方程只需两次简单的逐行回代。为什么条件数会平方：若 $J$ 的奇异值是 $\sigma_i$，则 $J^TJ$ 的特征值是 $\sigma_i^2$，最大与最小之比自然平方。数值上，双精度约 16 位有效数字。$\kappa(J)=10^6$ 时，直接处理 $J$ 的 QR 大约损失 6 位，还剩 10 位；正规方程损失 12 位，只剩 4 位。10.7.2 的实验正是这个结果。实务含义：自己写回归时，别用 inv(X.T @ X) @ X.T @ y，用 np.linalg.lstsq 或 scipy.linalg.lstsq。

2. QR 分解。 正交变换不改变欧氏范数。带列主元的 QR 分解

J\Pi=Q\begin{bmatrix}R\\0\end{bmatrix}=[Q_1\ Q_2]\begin{bmatrix}R\\0\end{bmatrix}=Q_1R,\tag{10.14}

于是

\|Jx+r\|_2^2=\|R(\Pi^Tx)+Q_1^Tr\|_2^2+\|Q_2^Tr\|_2^2.\tag{10.15}

第二项与 $x$ 无关，令第一项为零即得 $x^*=-\Pi R^{-1}Q_1^Tr$（实际是解三角系统 $Rz=-Q_1^Tr$ 后置换）。相对误差通常正比于 $J$ 的条件数而不是其平方，一般可靠，是默认首选。

推导拆解：(10.15) 的由来。

正交矩阵不改变长度，所以 $\|Jx+r\|^2=\|Q^T(Jx+r)\|^2$。

$Q^TJ\Pi=\begin{bmatrix}R\\0\end{bmatrix}$，记 $z=\Pi^Tx$（$\Pi$ 是置换矩阵，只是把变量换个顺序，$\Pi\Pi^T=I$），则 $Q^TJx=Q^TJ\Pi z=\begin{bmatrix}Rz\\0\end{bmatrix}$。

$Q^Tr=\begin{bmatrix}Q_1^Tr\\Q_2^Tr\end{bmatrix}$。把上下两块分开算长度平方：上块 $\|Rz+Q_1^Tr\|^2$，下块 $\|Q_2^Tr\|^2$。下块与 $x$ 无关，就是最优时剩下的残差平方和（回归里的 RSS）；上块可以通过解 $Rz=-Q_1^Tr$ 精确地变成 0。直观上，$Q_1$ 的列张成"模型能解释的空间"，$Q_2$ 的列张成"模型解释不了的空间"，和 $R^2$ 分解中"解释平方和 + 残差平方和"的划分一致。

3. SVD。 $J=U\begin{bmatrix}S\\0\end{bmatrix}V^T=U_1SV^T$，$S=\operatorname{diag}(\sigma_1\ge\dots\ge\sigma_n>0)$，同理

\|Jx+r\|_2^2=\|S(V^Tx)+U_1^Tr\|_2^2+\|U_2^Tr\|_2^2,\tag{10.17}

x^*=-VS^{-1}U_1^Tr=-\sum_{i=1}^n\frac{u_i^Tr}{\sigma_i}v_i.\tag{10.18}

（精读笔记指出抽取文本中 (10.18) 漏了负号；按 (10.17) 推导应带负号，(10.19) 同理。）这个表达式给出了敏感性信息：$\sigma_i$ 小时，$x^*$ 对影响 $u_i^Tr$ 的数据扰动特别敏感，误差被放大 $1/\sigma_i$ 倍，方向是 $v_i$。$J$ 接近秩亏（$\sigma_n/\sigma_1\ll1$）时这一信息尤其有用，QR 一般不提供。

金融直觉：把 (10.18) 读成"解 = 各主成分方向上的贡献之和"。$v_i$ 是参数空间里的第 $i$ 个方向（比如"市场 beta 加 0.7、规模 beta 加 0.3、行业 beta 减 1"这样一个组合），$u_i^Tr$ 是数据在对应方向上的"信号"，除以 $\sigma_i$ 就是这个方向上系数该走多远。若某个 $\sigma_i$ 很小，说明数据几乎无法区分沿 $v_i$ 的不同参数组合：三个因子沿这个方向此消彼长，拟合效果几乎一样。此时噪声在 $u_i^Tr$ 上哪怕有一点点，除以很小的 $\sigma_i$ 后也会让系数大幅摆动。这就是 CFA 里说的"多重共线性使系数标准误变大、符号可能反转"背后的精确机制。

三者比较与秩亏情形。 Cholesky 在 $m\gg n$ 且只能存 $J^TJ$ 而不能存 $J$ 时特别有用。QR 中病态通常表现为 $R$ 右下角元素远小于其余元素。SVD 对病态问题最稳健。$J$ 真正秩亏时部分 $\sigma_i=0$，任意

x^*=-\sum_{\sigma_i\ne0}\frac{u_i^Tr}{\sigma_i}v_i+\sum_{\sigma_i=0}\tau_iv_i\tag{10.19}

都是极小点，通常取 $\tau_i=0$ 的最小范数解。$J$ 满秩但病态时，可以略去小 $\sigma_i$ 对应的项，得到稳健的近似解——截断 SVD，统计上就是主成分回归。

10.3 Gauss–Newton 法

10.3.1 方向与四个优点

把线搜索牛顿法中 Hessian 的二阶项去掉，解

J_k^TJ_k\,p_k^{\text{GN}}=-J_k^Tr_k.\tag{10.20}

Gauss–Newton 法有四个优点：

几乎免费。 近似 $\nabla^2f_k\approx J_k^TJ_k$ 省去了每个残差的 Hessian；算梯度 $J_k^Tr_k$ 时已经有了 $J_k$。
常常接近牛顿法。 当各项 $|r_j|\,\|\nabla^2r_j\|$ 明显小于 $J^TJ$ 的特征值时（小残差或近线性），GN 表现接近牛顿法。实践中许多问题在解处残差小，常观察到快速局部收敛。
总是下降方向（只要 $J_k$ 满秩且 $\nabla f_k\ne0$）： $(p^{\text{GN}})^T\nabla f_k=(p^{\text{GN}})^TJ_k^Tr_k=-(p^{\text{GN}})^TJ_k^TJ_kp^{\text{GN}}=-\|J_kp^{\text{GN}}\|^2\le0,$ 等号仅当 $J_kp^{\text{GN}}=0$，而这等价于 $J_k^Tr_k=\nabla f_k=0$。
本身是线性最小二乘问题。 (10.20) 是 $\min_p\ \|J_kp+r_k\|_2^2\tag{10.22}$ 的正规方程（精读笔记记录原书此处误印为 $f_k$，应为 $r_k$），因此可用 QR 或 SVD 求方向，不必显式形成 $J_k^TJ_k$。

另一个理解角度：GN 不是对 $f$ 建二次模型，而是对向量函数 $r$ 建线性模型 $r(x+p)\approx r(x)+J(x)p$，代入 $\tfrac12\|r\|^2$ 后对 $p$ 最小化。

推导拆解：把线性模型代进去展开：$\tfrac12\|r+Jp\|^2=\tfrac12r^Tr+p^TJ^Tr+\tfrac12p^TJ^TJp$。对 $p$ 求梯度并令其为零：$J^Tr+J^TJp=0$，正是 (10.20)。与牛顿法比较，牛顿法的二次项是完整 Hessian $J^TJ+\sum r_j\nabla^2r_j$，GN 只保留了第一部分。第 3 条"下降方向"的推导用了 (10.20) 本身：把 $J_k^Tr_k$ 换成 $-J_k^TJ_kp^{\text{GN}}$，于是 $p^TJ^Tr=-p^TJ^TJp=-\|Jp\|^2$。

金融直觉：用收益率曲线拟合来理解 GN。手上有一组债券报价，模型价格是曲线参数的非线性函数。GN 的每一步相当于：在当前曲线附近，用每只债券对每个参数的"久期式敏感度"（即 $J$ 的一行）把定价线性化，"价格误差 ≈ 当前误差 + 敏感度 × 参数调整量"，然后用一次 OLS 求出能最大程度消除定价误差的参数调整量。走过去之后重新计算敏感度，再来一次。若模型在解附近接近线性（凸性小）、或者拟合得很好（残差小），这种"用久期近似"的误差很小，所以 GN 很快。

10.3.2 全局收敛

沿 GN 方向做满足 Wolfe 条件的线搜索。设 Jacobian 的奇异值一致有下界：存在 $\gamma>0$ 使

\|J(x)z\|_2\ge\gamma\|z\|_2\tag{10.23}

对水平集 $\mathcal L=\{x:f(x)\le f(x_0)\}$ 邻域内所有 $x$ 成立。

定理 10.1：各 $r_j$ 在水平集邻域内 Lipschitz 连续可微，$J$ 满足 (10.23)，步长满足 Wolfe 条件，则 $\lim_{k\to\infty}J_k^Tr_k=0$。

证明：$J$ 连续，故存在 $\beta$ 使 $\|J(x)\|\le\beta$。于是

\cos\theta_k=-\frac{\nabla f^Tp^{\text{GN}}}{\|p^{\text{GN}}\|\,\|\nabla f\|}=\frac{\|Jp^{\text{GN}}\|^2}{\|p^{\text{GN}}\|\,\|J^TJp^{\text{GN}}\|}\ge\frac{\gamma^2\|p^{\text{GN}}\|^2}{\beta^2\|p^{\text{GN}}\|^2}=\frac{\gamma^2}{\beta^2}>0.

方向与负梯度的夹角一致远离 $90°$，由 Zoutendijk 定理即得结论。若某步 $J_k$ 秩亏，(10.20) 仍有（无穷多）解，但 $\cos\theta_k$ 不再保证远离零，定理不能一般地成立——这正是 LM 方法要解决的问题。

10.3.3 局部收敛速度

记 $H(x)=\sum_jr_j\nabla^2r_j$ 为被忽略的二阶项。类似牛顿法的分析可得

\|x_k+p_k^{\text{GN}}-x^*\|\lesssim\big\|[J^TJ(x^*)]^{-1}H(x^*)\big\|\,\|x_k-x^*\|+O(\|x_k-x^*\|^2).\tag{10.25}

由此三种情形一目了然：

零残差（$H(x^*)=0$）：二次收敛，与牛顿法一样；
小残差（$\|(J^TJ)^{-1}H\|\ll1$）：线性收敛但速率常数很小，实际上很快；
大残差（$\|(J^TJ)^{-1}H\|\ge1$）：单位 GN 步甚至可能不会更接近解。

白话解释：$(J^TJ)^{-1}H$ 是"被忽略的那部分 Hessian"与"保留的那部分"之比。(10.25) 说：每一步之后的误差 ≈ 这个比值 × 之前的误差。比值是 0.01，误差每步缩小到 1%，非常快；比值是 0.5，每步只缩一半；比值超过 1，误差可能不减反增。 $H=\sum_jr_j\nabla^2r_j$ 是"残差大小 × 模型弯曲程度"。所以 GN 好不好用，取决于两件事的乘积：模型在解处拟合得多好（$r_j$ 小不小），以及模型对参数有多非线性（$\nabla^2r_j$ 大不大）。只要其中一个足够小，GN 就快。用 Nelson–Siegel 拟合报价时，残差在几个基点，属于典型的小残差问题。

10.4 Levenberg–Marquardt 法

10.4.1 GN + 信赖域

LM 方法把 GN 的线搜索换成信赖域，消除了 GN 在 $J$ 秩亏或接近秩亏时的弱点；它仍然忽略二阶项，所以局部收敛性质与 GN 相同。历史上 LM 通常被视为一般信赖域方法的鼻祖。子问题是

\min_p\ \tfrac12\|J_kp+r_k\|_2^2\quad\text{s.t.}\quad\|p\|\le\Delta_k,\tag{10.26}

相当于模型 $m_k(p)=\tfrac12\|r_k\|^2+p^TJ_k^Tr_k+\tfrac12p^TJ_k^TJ_kp$。

定理 10.2：在一般信赖域算法框架（原书 Algorithm 4.1，$\eta\in(0,\tfrac14)$）中，$r_j$ 在水平集邻域内二阶连续可微，每步近似解的模型下降量不小于 Cauchy 点的一个固定比例：

m_k(0)-m_k(p_k)\ge c_1\|J_k^Tr_k\|\min\Big(\Delta_k,\frac{\|J_k^Tr_k\|}{\|J_k^TJ_k\|}\Big),\tag{10.28}

则 $\lim\nabla f_k=\lim J_k^Tr_k=0$。注意这里不需要满秩假设 (10.23)。

10.4.2 子问题的刻画

若 GN 步 $\|p^{\text{GN}}\|<\Delta$，它就是 (10.26) 的解；否则存在 $\lambda>0$ 使解满足 $\|p\|=\Delta$ 且

(J^TJ+\lambda I)p=-J^Tr.\tag{10.29}

引理 10.3：$p^{\text{LM}}$ 是 $\min\|Jp+r\|^2$ s.t. $\|p\|\le\Delta$ 的解，当且仅当存在 $\lambda\ge0$ 使

(J^TJ+\lambda I)p^{\text{LM}}=-J^Tr,\qquad \lambda(\Delta-\|p^{\text{LM}}\|)=0.\tag{10.30}

（第二式是互补条件，第 12 章会看到它正是 KKT 条件的一部分。）一般信赖域子问题还需要 $B+\lambda I\succeq0$，这里 $J^TJ\succeq0$、$\lambda\ge0$，自动成立。

(10.29) 又是线性最小二乘问题

\min_p\ \frac12\Big\|\begin{bmatrix}J\\\sqrt\lambda I\end{bmatrix}p+\begin{bmatrix}r\\0\end{bmatrix}\Big\|^2\tag{10.31}

的正规方程，所以求 LM 步同样不必形成 $J^TJ$。

与岭回归的关系。（第 04 章 4.6 节已从信赖域子问题的角度讨论过同一结构，这里从最小二乘角度再看一次。）把 $J$ 换成设计矩阵 $X$、$-r$ 换成 $y$，(10.29) 就是岭回归 $(X^TX+\lambda I)\beta=X^Ty$。用 SVD 写出解：$p=-\sum_i\frac{\sigma_i}{\sigma_i^2+\lambda}(u_i^Tr)v_i$。与 (10.18) 比较，$1/\sigma_i$ 被换成了 $\sigma_i/(\sigma_i^2+\lambda)$：大奇异值方向几乎不变，小奇异值方向被压缩。LM 的信赖域约束与岭回归的 $\ell_2$ 惩罚是同一个数学对象，前者限制步长，后者限制系数大小。

推导拆解：(10.31) 为什么等价于 (10.29)，以及滤波因子怎么来。

增广矩阵 $A=\begin{bmatrix}J\\\sqrt\lambda I\end{bmatrix}$，增广右端 $b=\begin{bmatrix}r\\0\end{bmatrix}$。线性最小二乘 $\min\|Ap+b\|^2$ 的正规方程是 $A^TAp=-A^Tb$。

按分块计算：$A^TA=J^TJ+\sqrt\lambda I\cdot\sqrt\lambda I=J^TJ+\lambda I$，$A^Tb=J^Tr+\sqrt\lambda I\cdot0=J^Tr$。正好是 (10.29)。

直观上，增广的 $n$ 行 $\sqrt\lambda I$ 相当于加了 $n$ 个"虚拟观测"，每个都说"这个参数的调整量应该是 0"，权重为 $\sqrt\lambda$。这和岭回归的贝叶斯解释（参数先验均值为 0）一致。

滤波因子：把 $J=U_1SV^T$ 代入，$J^TJ+\lambda I=V(S^2+\lambda I)V^T$，$J^Tr=VSU_1^Tr$，所以 $p=-V(S^2+\lambda I)^{-1}SU_1^Tr$，逐个方向就是 $-\dfrac{\sigma_i}{\sigma_i^2+\lambda}(u_i^Tr)v_i$。$\sigma_i^2\gg\lambda$ 时约为 $1/\sigma_i$（不变），$\sigma_i^2\ll\lambda$ 时约为 $\sigma_i/\lambda$（几乎被压没）。所以 $\lambda$ 小时 LM ≈ GN（大胆走），$\lambda$ 大时 $p\approx-J^Tr/\lambda$，即沿负梯度走一小步（谨慎走）。LM 根据上一步效果在两者之间自动切换。

历史。 Levenberg 与 Marquardt 的原始方法里没有信赖域的概念，而是直接调整 $\lambda$：上一步有效就减小，否则增大。这一启发式与信赖域半径调整如出一辙，Moré 后来牢固地建立了两者的联系，并由此催生了 MINPACK 中高效稳健的实现。

10.4.3 实现

用信赖域一章的求根算法寻找与给定 $\Delta$ 匹配的 $\lambda$。保护很容易：$J^TJ$ 半正定，只要 $\lambda>0$，Cholesky 因子就一定存在。更好的做法是对 (10.31) 的系数矩阵做 QR，得 $R_\lambda^TR_\lambda=J^TJ+\lambda I$，并用 Householder + Givens 组合提高效率：先用 Householder 只对 $J$ 做一次 QR，$J=Q\begin{bmatrix}R\\0\end{bmatrix}$，于是

\begin{bmatrix}R\\0\\\sqrt\lambda I\end{bmatrix}=\begin{bmatrix}Q^T&\\&I\end{bmatrix}\begin{bmatrix}J\\\sqrt\lambda I\end{bmatrix}.\tag{10.34}

左端只有 $\sqrt\lambda I$ 的 $n$ 个对角元破坏了上三角结构，用 $n(n+1)/2$ 次 Givens 旋转即可消去。求根过程中改变 $\lambda$ 时只需重做 Givens 部分：$m\gg n$ 时，首次 $O(mn^2)$ 之后每个 $\lambda$ 只需 $O(n^3)$。

尺度与椭球信赖域。 最小二乘问题的尺度常常很差（有的变量约 $10^4$，有的约 $10^{-6}$）。改用椭球信赖域 $\|D_kp\|\le\Delta_k$（$D_k$ 为正对角阵），解满足

(J_k^TJ_k+\lambda D_k^2)p_k^{\text{LM}}=-J_k^Tr_k,\tag{10.36}

等价于 $\min_p\Big\|\begin{bmatrix}J_k\\\sqrt\lambda D_k\end{bmatrix}p+\begin{bmatrix}r_k\\0\end{bmatrix}\Big\|$。$D_k$ 可以随迭代变化，在一定范围内收敛理论不受影响。Seber–Wild 建议令 $D_k^2$ 等于 $J_k^TJ_k$ 的对角部分，使算法对变量的对角缩放不变。MINPACK 的 LM 实现默认采用类似做法：以 Jacobian 各列范数（取迭代过程中的历史最大值）作为 $D_k$ 的对角元。

局部行为：与 GN 相同。小残差解附近模型准确，信赖域最终不活跃，算法取单位 GN 步，快速线性收敛。

10.5 大残差与大规模问题

10.5.1 大残差问题

二阶项太大时，(10.26) 的二次模型不再充分。有统计学者认为残差大说明模型不合格，应当重建模型；但大残差也常由人为错误造成的离群值引起（读错仪表、把地震读数归错事件），这时仍需求解，以便识别离群值并删除或降权。

GN 与 LM 在大残差情形下表现较差：渐近只有线性收敛，慢于牛顿/拟牛顿法的超线性收敛。但在远离解的早期迭代中，牛顿/拟牛顿法又可能不如 GN/LM。由于事先不知道残差大小，人们设计了混合算法：

Fletcher–Xu 方法：维护一个正定矩阵 $B_k$。若 GN 步使 $f$ 下降达到某个固定比例（如因子 5），就接受该步并用 $J_k^TJ_k$ 覆盖 $B_k$；否则用 $B_k$ 求方向并线搜索。两种情形都对 $B_k$ 做 BFGS 式更新。零残差时最终总取 GN 步（二次收敛），非零残差时最终退化为 BFGS（超线性收敛）。
只近似二阶项（NL2SOL）：维护 $S_k\approx\sum_jr_j(x_k)\nabla^2r_j(x_k)$，用 $B_k=J_k^TJ_k+S_k$ 作模型。Dennis–Gay–Welsch 算法要求 $S_{k+1}$ 满足割线方程 $S_{k+1}s=y^\sharp,\qquad s=x_{k+1}-x_k,\quad y^\sharp=J_{k+1}^Tr_{k+1}-J_k^Tr_{k+1},$ （精读笔记指出原书此处一式误印为 $J_{k+1}^Tr_{k+1}-J_k^Tr_k$，按推导应为 $J_k^Tr_{k+1}$，与后文 $y^\sharp$ 的定义一致），并在对称、变化最小的意义下更新： $S_{k+1}=S_k+\frac{(y^\sharp-S_ks)y^T+y(y^\sharp-S_ks)^T}{y^Ts}-\frac{(y^\sharp-S_ks)^Ts}{(y^Ts)^2}yy^T,\qquad y=J_{k+1}^Tr_{k+1}-J_k^Tr_k.\tag{10.38}$ 这是 DFP 的轻微变体。为保证零残差时 $S_k$ 能够消失，更新前把 $S_k$ 缩放为 $\tau_kS_k$，$\tau_k=\min(1,|s^Ty^\sharp|/|s^TS_ks|)$；当 GN 模型足够好时干脆略去 $S_k$。

10.5.2 大规模问题

$n$ 小、$m$ 巨大（如 $m\sim10^6$、$n\sim100$）：$J$ 可能存不下，但可以逐个计算 $r_j$ 和 $\nabla r_j$，累加 $J^TJ=\sum_j\nabla r_j\nabla r_j^T$ 和 $J^Tr=\sum_jr_j\nabla r_j$，直接解正规方程；LM 中改变 $\lambda$ 时只需在 $J^TJ$ 上加 $\lambda I$ 重新分解。这是"用条件数平方换内存"的典型权衡。量化中逐笔数据、超长样本的模型估计属于这一类。
$n$、$m$ 都大且 $J$ 稀疏：类似非精确牛顿法，用 CG 类迭代法近似求解 GN/LM 方程，$J^TJ$ 的半正定性简化了算法。除非 $\nabla^2f(x^*)=J^TJ$，否则不能期望超线性收敛。
Wright–Holt 非精确 LM：步 $\bar p$ 只需满足 $\|(J^TJ+\lambda I)\bar p+J^Tr\|\le\eta\|J^Tr\|$，不用信赖域半径，而是直接选 $\lambda$ 使实际下降与模型预测下降之比不小于 $\gamma_1\in(0,1)$，且 $\lambda$ 不比满足该条件的最小值大太多，即可证明全局收敛。（精读笔记记录该比值的分母中正则项印作 $\lambda^2\|\bar p\|^2$；若按 (10.29) 的参数化，LM 模型的正则项应为 $\lambda\|\bar p\|^2$，两者只差 $\lambda$ 的参数化方式，需要时请核对原始论文。）它借助 Paige–Saunders 的 LSQR 算法同时求解多个 $\lambda$ 对应的 (10.31)，每多一个 $\lambda$ 只多存两个向量，代价与解单个最小二乘问题相差无几。

10.6 正交距离回归

例 10.1 假设自变量 $t_j$ 没有误差。若 $t_j$ 本身也有测量误差而被忽略，结果可能严重失真。统计学称这类模型为变量含误差模型（errors-in-variables）；线性情形的优化问题叫总体最小二乘（total least squares），非线性情形叫正交距离回归（orthogonal distance regression, ODR）。

为 $t_j$ 引入扰动 $\delta_j$、为 $y_j$ 引入 $\epsilon_j$：

y_j=\phi(x;t_j+\delta_j)+\epsilon_j,\qquad \min_{x,\delta,\epsilon}\ \tfrac12\sum_{j=1}^mw_j^2\epsilon_j^2+d_j^2\delta_j^2.\tag{10.41–10.42}

权重全相等时，每一项就是数据点到曲线的最短距离的平方，最短路径在交点处与曲线正交——名称由此而来。消去 $\epsilon_j$ 得到无约束最小二乘

\min_{x,\delta}\ \tfrac12\sum_{j=1}^mw_j^2[y_j-\phi(x;t_j+\delta_j)]^2+d_j^2\delta_j^2=\tfrac12\sum_{j=1}^{2m}r_j^2(x,\delta),\tag{10.43}

这是 $2m$ 个残差、$n+m$ 个未知数的标准问题（精读笔记指出原书此处把未知数个数和残差个数写混了）。朴素求解很贵，但 Jacobian 有块结构

J(x,\delta)=\begin{bmatrix}\hat J&V\\0&D\end{bmatrix},\tag{10.45}

$V$、$D$ 是 $m\times m$ 对角阵。Boggs–Byrd–Schnabel 对 (10.43) 用 LM，LM 方程的右下块是对角阵，可以直接消去 $p_\delta$，得到只含 $p_x$ 的 $n\times n$ 系统，总代价只比普通 $m\times n$ 最小二乘略高。ODRPACK（scipy.odr 的底层）实现了这一算法。

金融直觉：两阶段 Fama–MacBeth 回归是典型的变量含误差问题。第一阶段用时间序列回归估出每只股票的 beta，第二阶段把估计出的 beta 当作解释变量回归平均收益。beta 是估计值、带误差，相当于这里的 $t_j$ 有 $\delta_j$。普通最小二乘只量竖直距离、把 $t_j$ 当成精确值，结果是第二阶段的斜率（风险溢价）系统性地偏向 0，这就是"衰减偏误"。实务中常用组合（而非单只股票）的 beta 来降低 $\delta_j$ 的方差，思路与 ODR 不同，但针对的是同一个问题。块结构 (10.45) 的直观：每个 $\delta_j$ 只影响第 $j$ 个观测，所以 $V$、$D$ 都是对角的；对角块可以逐个消去，这是 ODR 便宜的原因。

10.7 量化实战

10.7.1 在哪里用

因子回归与数值稳定性。 时间序列因子暴露估计、Fama–MacBeth 横截面回归都是线性最小二乘。因子高度共线时，教科书公式 $(X^TX)^{-1}X^Ty$ 的条件数平方效应会放大数值误差。numpy.linalg.lstsq 用 SVD，statsmodels 的 OLS 默认用伪逆，scipy.linalg.lstsq 默认用 LAPACK 的 gelsd（SVD），都比手写正规方程可靠。更重要的是 (10.18)：共线方向上的系数对数据扰动极其敏感，这是统计问题而不仅是数值问题，需要截断 SVD（主成分回归）或岭回归（即 LM 的 $J^TJ+\lambda I$）来稳定估计。
曲线拟合与模型校准。 Nelson–Siegel(–Svensson) 收益率曲线、SVI 波动率微笑、SABR/Heston 对市场报价的校准、冲击成本"平方根律"参数估计，都是非线性最小二乘，业界标准工具就是 LM（scipy.optimize.least_squares(method='lm') 调用 MINPACK；'trf' 是带界约束的信赖域反射法）。参数量级悬殊时，$D_k$ 缩放或手工标准化很关键；报价噪声大、有离群报价时，GN/LM 收敛变慢，可以加权或用稳健损失（least_squares 的 loss='huber'/'soft_l1'）。
标准误。 在正态误差下，$\hat\sigma^2(J^TJ)^{-1}$ 是参数的渐近协方差估计——GN 的 Hessian 近似恰好就是 Fisher 信息的估计。一般 MLE 中与之对应的是用梯度外积近似信息矩阵的 BHHH/OPG 估计量。
变量含误差。 用估计出来的 beta 做第二阶段回归时，解释变量本身带误差，会导致系数衰减偏误。ODR/总体最小二乘是一种计算层面的处理思路，统计上更常用工具变量或误差修正（第 05 册第 09 章讨论变量测量误差偏误，第 12 章讲工具变量）。

10.7.2 代码一：共线因子回归——正规方程、QR、SVD

构造一个"行业因子 ≈ 0.7×市场 + 0.3×规模"的近似共线设计矩阵，用无噪声的 $y$ 检验三种算法的纯数值误差。

import numpy as np
from scipy.linalg import cho_factor, cho_solve, qr, solve_triangular

rng = np.random.default_rng(7)
T = 500
mkt = rng.standard_normal(T)
size = rng.standard_normal(T)
# 第三个“因子”几乎是前两个的线性组合（如行业指数 ≈ 市场 + 规模）
ind = 0.7*mkt + 0.3*size + 1e-6*rng.standard_normal(T)
X = np.column_stack([np.ones(T), mkt, size, ind])
beta_true = np.array([0.001, 1.2, -0.4, 0.5])
y = X @ beta_true                                  # 无噪声：任何误差都来自数值计算

print("cond(X) = %.1e, cond(X'X) = %.1e" % (np.linalg.cond(X), np.linalg.cond(X.T @ X)))
# (1) 正规方程 + Cholesky (10.12)
b_ne = cho_solve(cho_factor(X.T @ X), X.T @ y)
# (2) QR 分解 (10.14)：X = Q1 R，解 R b = Q1' y
Q1, R = qr(X, mode='economic')
b_qr = solve_triangular(R, Q1.T @ y)
# (3) SVD (10.18)：b = Σ (u_i' y / σ_i) v_i
U, s, Vt = np.linalg.svd(X, full_matrices=False)
b_svd = Vt.T @ ((U.T @ y) / s)
for name, b in [("正规方程", b_ne), ("QR", b_qr), ("SVD", b_svd)]:
    print(f"{name:6s} 相对误差 {np.linalg.norm(b-beta_true)/np.linalg.norm(beta_true):.1e}")
print("奇异值:", np.array2string(s, precision=3))
print("最小奇异值对应的右奇异向量 v_4:", np.round(Vt[-1], 3))

输出：

cond(X) = 1.4e+06, cond(X'X) = 2.1e+12
正规方程   相对误差 3.9e-04
QR     相对误差 6.0e-12
SVD    相对误差 6.1e-11
奇异值: [2.711e+01 2.194e+01 2.086e+01 1.873e-05]
最小奇异值对应的右奇异向量 v_4: [-0.     0.557  0.239 -0.796]

$\text{cond}(X)\approx1.4\times10^6$，$\text{cond}(X^TX)\approx2\times10^{12}$，正好是平方关系。正规方程的相对误差是 $4\times10^{-4}$，QR 和 SVD 在 $10^{-11}$ 量级——相差七个数量级，与"误差正比于条件数"的预言一致。最小奇异值对应的右奇异向量 $v_4\propto(0,\,0.7,\,0.3,\,-1)$，正是我们埋进去的共线关系：(10.18) 告诉我们，数据沿这个方向的任何扰动都会被放大 $1/\sigma_4\approx5\times10^4$ 倍反映到系数上。在真实数据中（$y$ 有噪声），QR 也救不了这一点——需要的是删除冗余因子、截断 SVD 或岭回归。

10.7.3 代码二：Nelson–Siegel 收益率曲线拟合

Nelson–Siegel 模型

y(\tau)=\beta_0+\beta_1\frac{1-e^{-\tau/\theta}}{\tau/\theta}+\beta_2\Big(\frac{1-e^{-\tau/\theta}}{\tau/\theta}-e^{-\tau/\theta}\Big)

对 $(\beta_0,\beta_1,\beta_2)$ 线性、对衰减参数 $\theta$ 非线性（代码中 NS 的衰减参数记为 th，以免与 LM 的 $\lambda$ 混淆）。用 11 个期限、带 3bp 噪声的模拟收益率，对比手写 GN（QR/SVD 求步 + 回溯）、手写 LM（$D_k^2=\operatorname{diag}(J^TJ)$，按实际/预测下降比调 $\lambda$）和 scipy 的 MINPACK LM，然后检查局部极小。

import numpy as np
from scipy.optimize import least_squares

rng = np.random.default_rng(3)
tau = np.array([0.25, 0.5, 1, 2, 3, 5, 7, 10, 15, 20, 30.0])          # 期限（年）
def ns(x, t):
    b0, b1, b2, th = x                  # th：Nelson–Siegel 衰减参数 θ
    z = t/th; L1 = (1-np.exp(-z))/z
    return b0 + b1*L1 + b2*(L1 - np.exp(-z))
x_true = np.array([0.045, -0.020, 0.015, 2.0])
yobs = ns(x_true, tau) + 0.0003*rng.standard_normal(tau.size)          # 3bp 报价噪声

def resid(x): return ns(x, tau) - yobs                                  # r_j(x)
def jac(x, h=1e-7):                                                     # 前向差分 Jacobian
    r0 = resid(x); J = np.empty((tau.size, 4))
    for i in range(4):
        e = np.zeros(4); e[i] = h*max(1, abs(x[i])); J[:, i] = (resid(x+e) - r0)/e[i]
    return J

def gauss_newton(x, iters=30):
    for k in range(iters):
        r, J = resid(x), jac(x)
        p = np.linalg.lstsq(J, -r, rcond=None)[0]      # (10.22)：用最小二乘求 GN 步，不形成 J'J
        a = 1.0                                         # 简单回溯
        while 0.5*np.sum(resid(x+a*p)**2) > 0.5*np.sum(r**2) + 1e-4*a*(r @ J @ p) and a > 1e-8:
            a *= 0.5
        x = x + a*p
        if np.linalg.norm(J.T @ r) < 1e-12: break
    return x, k+1

def levenberg_marquardt(x, iters=200, lam=1e-2):
    f = 0.5*np.sum(resid(x)**2)
    for k in range(iters):
        r, J = resid(x), jac(x)
        g = J.T @ r
        if np.linalg.norm(g) < 1e-12: break
        D2 = np.diag(np.diag(J.T @ J))                  # D_k^2 = diag(J'J)：对角缩放不变 (10.36)
        A = np.vstack([J, np.sqrt(lam)*np.sqrt(D2)])    # (10.37) 的增广最小二乘形式
        p = np.linalg.lstsq(A, -np.concatenate([r, np.zeros(4)]), rcond=None)[0]
        fn = 0.5*np.sum(resid(x+p)**2)
        pred = f - 0.5*np.sum((r + J @ p)**2)
        rho = (f - fn)/pred if pred > 0 else -1
        if rho > 0.25: x, f, lam = x + p, fn, lam/3     # 好步：接受并减小 λ（相当于放大信赖域）
        else: lam *= 4                                   # 坏步：拒绝并增大 λ
    return x, k+1

x0 = np.array([0.03, 0.0, 0.0, 0.5])                                    # 粗糙初值
for name, solver in [("Gauss-Newton", gauss_newton), ("LM（手写）", levenberg_marquardt)]:
    x, k = solver(x0.copy())
    print(f"{name:12s} 迭代 {k:3d}  参数 {np.round(x, 5)}  RMSE {1e4*np.sqrt(np.mean(resid(x)**2)):.2f}bp")
sol = least_squares(resid, x0, method='lm')
print(f"{'scipy lm':12s} 求值 {sol.nfev:3d}  参数 {np.round(sol.x, 5)}  RMSE {1e4*np.sqrt(np.mean(sol.fun**2)):.2f}bp")

J = jac(sol.x)
print("解处 cond(J) = %.1e，各列范数:" % np.linalg.cond(J), np.round(np.linalg.norm(J, axis=0), 3))
s2 = np.sum(sol.fun**2)/(tau.size - 4)
se = np.sqrt(np.diag(s2*np.linalg.inv(J.T @ J)))       # GN 近似 J'J 给出的渐近标准误
print("标准误（J'J 近似）:", np.round(se, 5))

# 剖面法：给定 θ 时模型对 (β0,β1,β2) 线性，可在 θ 网格上只解线性最小二乘
def profile_rss(th):
    z = tau/th; L1 = (1-np.exp(-z))/z
    X = np.column_stack([np.ones_like(tau), L1, L1 - np.exp(-z)])
    b = np.linalg.lstsq(X, yobs, rcond=None)[0]
    return np.sum((X @ b - yobs)**2), b
grid = np.linspace(0.2, 8, 157)
rss = np.array([profile_rss(t)[0] for t in grid])
loc = [i for i in range(1, len(grid)-1) if rss[i] < rss[i-1] and rss[i] < rss[i+1]]
print("剖面 RSS 的局部极小:", [(float(round(grid[i], 2)), float(round(rss[i]*1e8, 1))) for i in loc], "（θ, bp²）")
sol2 = least_squares(resid, np.r_[profile_rss(2.5)[1], 2.5], method='lm')   # 从另一个盆地出发
print(f"从 θ=2.5 出发 → 参数 {np.round(sol2.x, 5)}  RMSE {1e4*np.sqrt(np.mean(sol2.fun**2)):.2f}bp")
grid_t = np.linspace(0.25, 30, 300)
print("两条拟合曲线的最大差异: %.1f bp" % (1e4*np.max(np.abs(ns(sol.x, grid_t) - ns(sol2.x, grid_t)))))

输出：

Gauss-Newton 迭代  10  参数 [ 0.04576 -0.01888 -0.0186   0.54367]  RMSE 3.79bp
LM（手写）       迭代  11  参数 [ 0.04576 -0.01888 -0.0186   0.54367]  RMSE 3.79bp
scipy lm     求值   8  参数 [ 0.04576 -0.01888 -0.0186   0.54367]  RMSE 3.79bp
解处 cond(J) = 7.4e+02，各列范数: [3.317 1.185 0.543 0.024]
标准误（J'J 近似）: [0.00029 0.00166 0.00636 0.10193]
剖面 RSS 的局部极小: [(0.55, 158.1), (1.8, 197.7)] （θ, bp²）
从 θ=2.5 出发 → 参数 [ 0.04565 -0.02057  0.0105   1.78822]  RMSE 4.24bp
两条拟合曲线的最大差异: 4.1 bp

几点值得细看。

三个求解器收敛到同一点。 小残差（RMSE 约 3.8bp，与 3bp 噪声相当）问题上，GN 与 LM 都只用约 10 次迭代——符合 10.3.3 节"小残差时 GN 很快"的结论。

参数估计与真值相去甚远。 真值是 $\beta_2=0.015$、$\theta=2$，估计却是 $\beta_2=-0.019$、$\theta=0.54$。这不是算法的错：Jacobian 各列范数相差两个数量级（$\theta$ 列只有 0.024），$\theta$ 和 $\beta_2$ 的标准误分别是 0.10 和 0.006，说明在 11 个带噪声的点上，这两个参数几乎不可识别。

多个局部极小。 利用"给定 $\theta$ 后是线性最小二乘"做剖面（profile）扫描，发现 RSS 在 $\theta\approx0.55$ 和 $\theta\approx1.8$ 各有一个局部极小；从 $\theta=2.5$ 出发的 LM 收敛到后者。两条曲线最大只差 4bp，与报价噪声同量级，但参数完全不同。实务启示：(1) 对"部分线性"的模型，先对非线性参数做网格扫描、再用 LM 精修，是简单可靠的全局化手段（这种思想叫变量投影，variable projection）；(2) 若要对参数做时间序列分析（如把 $\beta_0,\beta_1,\beta_2$ 当作水平、斜率、曲率因子），常见做法是固定 $\theta$（Diebold–Li 固定衰减参数），把问题变成线性回归，以换取参数的稳定性。

本章小结

非线性最小二乘的核心结构是 $\nabla f=J^Tr$、$\nabla^2f=J^TJ+\sum r_j\nabla^2r_j$：只要有 Jacobian，Hessian 的主要部分就免费得到。正态误差下最小二乘就是最大似然。线性最小二乘应优先用 QR，病态或秩亏时用 SVD，正规方程会把条件数平方。Gauss–Newton 线性化残差、解线性最小二乘求步，在 Jacobian 满秩时是下降方向，零残差时二次收敛、小残差时很快、大残差时可能失效。Levenberg–Marquardt 是 GN 加信赖域，等价于解 $(J^TJ+\lambda D^2)p=-J^Tr$，与岭回归同构，是业界校准的标准工具。大残差用混合法或二阶项拟牛顿近似；自变量有误差时用正交距离回归。

概念	公式 / 要点
梯度与 Hessian	$\nabla f=J^Tr$；$\nabla^2f=J^TJ+\sum_jr_j\nabla^2r_j$
正规方程	$J^TJx=-J^Tr$；$\kappa(J^TJ)=\kappa(J)^2$
QR 解	$x^*=-\Pi R^{-1}Q_1^Tr$；误差 $\propto\kappa(J)$
SVD 解	$x^*=-\sum_i\frac{u_i^Tr}{\sigma_i}v_i$；秩亏取最小范数解
Gauss–Newton	$J^TJp=-J^Tr$ ⇔ $\min_p\|Jp+r\|$；$\cos\theta\ge\gamma^2/\beta^2$
GN 局部速率	由 $\|(J^TJ)^{-1}H(x^*)\|$ 决定；零残差二次收敛
Levenberg–Marquardt	$(J^TJ+\lambda D^2)p=-J^Tr$；$\lambda(\Delta-\|Dp\|)=0$
增广形式	$\min\|[J;\sqrt\lambda D]p+[r;0]\|$，不形成 $J^TJ$
岭回归联系	滤波因子 $\sigma_i/(\sigma_i^2+\lambda)$ 代替 $1/\sigma_i$
大残差	Fletcher–Xu 混合法；NL2SOL：$B=J^TJ+S_k$
ODR	自变量有误差；Jacobian 块结构使代价接近普通最小二乘

练习

基础

证明：$m\ge n$ 时，$J$ 列满秩当且仅当 $J^TJ$ 非奇异。（原书习题 10.1）
$J$ 秩亏时，证明 (10.19) 中取所有 $\tau_i=0$ 得到的是最小范数解。（原书习题 10.3）提示：$v_i$ 两两正交，$\|x\|^2=\sum(\cdot)^2+\sum\tau_i^2$。
对例 10.1 的模型 $\phi(x;t)=x_1+tx_2+e^{-x_3t}$，写出残差 $r_j$ 的梯度 $\nabla r_j$ 和 Hessian $\nabla^2r_j$。哪一个参数使问题非线性？提示：$\nabla r_j=-(1,t_j,-t_je^{-x_3t_j})^T$；只有 $\partial^2r_j/\partial x_3^2=-t_j^2e^{-x_3t_j}$ 非零。
用 SVD 写出 (10.29) 的解 $p(\lambda)$ 及 $\|p(\lambda)\|^2$，证明 $\lambda\to0$ 时 $p(\lambda)$ 趋于最小范数 GN 解，并说明 $\|p(\lambda)\|$ 关于 $\lambda$ 单调递减。（原书习题 10.6）提示：$p(\lambda)=-\sum_i\frac{\sigma_i(u_i^Tr)}{\sigma_i^2+\lambda}v_i$。
在 10.7.2 的代码中把噪声幅度 1e-6 改为 1e-3 和 1e-8，记录 cond(X) 与三种方法的误差，验证"正规方程误差 ∝ cond(X)²·机器精度"的规律。

进阶

用 SVD 表示 $\nabla f_k^Tp^{\text{GN}}$、$\|p^{\text{GN}}\|$、$\|\nabla f_k\|$ 和 $\cos\theta_k$，说明 $J_k$ 秩亏（或接近秩亏）时为什么不能保证 $\cos\theta_k$ 一致远离零。（原书习题 10.5）
从 ODR 的 LM 方程（原书 (10.46)）中消去 $p_\delta$，说明得到的 $n\times n$ 系统是哪个线性最小二乘问题的正规方程。（原书习题 10.7）
把 10.7.3 的问题改为 Diebold–Li 做法：固定 $\theta=1.37$（年），只做线性最小二乘估计 $(\beta_0,\beta_1,\beta_2)$。对 200 组独立噪声重复估计，比较"固定 $\theta$"与"自由估计 $\theta$"时 $\beta_2$ 估计值的标准差。
设计一个大残差例子：在 10.7.3 的数据中把一个期限的收益率加 50bp 离群值，比较 least_squares 默认损失与 loss='soft_l1' 的拟合结果和迭代次数，并联系 10.5.1 节解释。
对 $m=10^6$、$n=20$ 的线性回归，比较"逐块累加 $X^TX$ 再 Cholesky"与"整体 QR"的内存占用和精度。在什么条件下前者是可以接受的？

原书推荐习题：10.1、10.3（正规方程可解性与最小范数解）、10.5（GN 在秩亏时为何失效）、10.6（LM 步随 $\lambda$ 的变化，对应岭回归收缩）、10.7（ODR 的块消元）。原书第 10 章共有 10.1–10.7 题。

原书对照

本章内容	原书位置	PDF 页码
10.1 目标、梯度、Hessian 结构	第 10 章引言，式 (10.1)–(10.5)	PDF p.270–273
10.2.1–10.2.2 例 10.1、最大似然解释	§10.1 Background: Modeling, Regression, Statistics，式 (10.6)–(10.9)	PDF p.273–275
10.2.3 线性最小二乘：Cholesky、QR、SVD	§10.1 Linear Least-Squares Problems，式 (10.10)–(10.19)	PDF p.275–279
10.3 Gauss–Newton、定理 10.1、(10.25)	§10.2 The Gauss–Newton Method，式 (10.20)–(10.25)	PDF p.279–282
10.4 Levenberg–Marquardt、定理 10.2、引理 10.3、实现	§10.2 The Levenberg–Marquardt Method，式 (10.26)–(10.37)	PDF p.282–287
10.5 大残差与大规模问题	§10.2 Large-Residual Problems；Large-Scale Problems，式 (10.38)–(10.40)	PDF p.287–290
10.6 正交距离回归	§10.3 Orthogonal Distance Regression，式 (10.41)–(10.46)	PDF p.291–293
注释与参考、习题 10.1–10.7	Notes and References；Exercises	PDF p.293–295

页码换算：本段原书正文页码约等于 PDF 页码减 20。§10.2 内部各小节的 PDF 页码为按篇幅估计的大致位置。