最优化学习笔记：无约束优化算法（7）

6.7 非线性最小二乘问题算法

本节研究非线性最小二乘问题的算法．非线性最小二乘问题是一类特殊的无约束优化问题，它有非常广泛的实际应用背景．例如在统计中，我们经常建立如下带参数的模型：

$$b=\varphi (a;x)+\epsilon$$

其中 $a$ 为自变量，$b$ 为响应变量，它们之间的关系由函数 $\varphi (\cdot ;x)$ 决定且 $x$ 是参数，$\epsilon$ 是噪声项，即观测都是有误差的. 我们的目的是要根据观测 $(a_i,b_i)$，估计未知参数 $x$ 的值. 若 $\epsilon$ 服从高斯分布，则使用 ${\ell }_2$ 范数平方是处理高斯噪声最好的方式：对 ${\ell }_2$ 范数平方损失函数

$$f(x)=\frac{1}{m}\sum _{i=1}^m\Vert b_i-\varphi (a_i;x){\Vert }^2$$

进行极小化即可求出未知参数 $x$ 的估计. 而对 ${\ell }_2$ 范数平方损失函数求解极小值就是一个最小二乘问题.

最小二乘问题一般属于无约束优化问题，但由于问题特殊性，人们针对其结构设计了许多算法快速求解. 一般地，设 $x^{\ast }$ 为最小二乘问题的解，根据最优解处残量 $\sum _{i=1}^m\Vert b_i-\varphi (a_i;x){\Vert }^2$ 的大小，可以将最小二乘问题分为小残量问题和大残量问题. 本节针对小残量问题介绍两种方法：高斯-牛顿算法和 LM 方法；而针对大残量问题简要地引入带结构的拟牛顿法.

6.7.1 非线性最小二乘问题

考虑非线性最小二乘问题的一般形式
$$\begin{array}{r}f(x)=\frac{1}{2}\sum _{j=1}^m{r}_j^2(x)\end{array}(\mathrm{6.7.1})$$

其中 $r_j:{\mathbb{R}}^n\to \mathbb{R}$ 是光滑函数，并且假设 $m⩾n$. 我们称 $r_j$ 为残差. 为了表述问题的方便，定义残差向量 $r:{\mathbb{R}}^n\to {\mathbb{R}}^m$ 为
$$\begin{array}{r}r(x)=(r_1(x),r_2(x),\cdots ,r_m(x){)}^{\mathrm{T}}\end{array}$$

使用这一定义，函数 $f(x)$ 可以写为 $$f(x)=\frac{1}{2}\Vert r(x){\Vert }_2^2$$.

问题 (6.7.1) 是一个无约束优化问题，可以使用前面讲过的任何一种算法求解. 为此直接给出 $f(x)$ 的梯度和海瑟矩阵：
$$\begin{array}{rl}& \nabla f(x)=\sum _{j=1}^m{r}_j(x)\nabla r_j(x)=J(x{)}^{\mathrm{T}}r(x)(\mathrm{6.7.2}a)\\ & {\nabla }^{2}f(x)=\sum _{j=1}^m\nabla r_j(x)\nabla r_j(x{)}^{\mathrm{T}}+\sum _{i=1}^m{r}_i(x){\nabla }^2{r}_i(x)=J(x{)}^{\mathrm{T}}J(x)+\sum _{i=1}^m{r}_i(x){\nabla }^2{r}_i(x)(\mathrm{6.7.2}b)\end{array}$$

其中 $J(x)\in {\mathbb{R}}^{m\times n}$ 是向量值函数 $r(x)$ 在点 $x$ 处的雅可比矩阵: $$J(x)=\left[\begin{array}{c}\nabla r_1(x{)}^{\mathrm{T}}\\ \nabla r_2(x{)}^{\mathrm{T}}\\ ⋮\\ \nabla r_m(x{)}^{\mathrm{T}}\end{array}\right].$$ 这里指出，${\nabla }^{2}f(x)$ 在形式上分为两部分，分别为 $J(x{)}^{\mathrm{T}}J(x)$ 和 $\sum _{i=1}^m{r}_i(x){\nabla }^2{r}_i(x)$. 处理这两部分的难度是截然不同的：注意到计算 $\nabla f(x)$ 时需要 $r(x)$ 的雅可比矩阵，因此海瑟矩阵的前一项是自然得到的，不需要进行额外计算；而海瑟矩阵的第二项则需要计算每个 ${\nabla }^2{r}_i(x)$, 这会导致额外计算量，因此很多最小二乘算法就是根据这个性质来设计的.

6.7.2 高斯 – 牛顿算法

高斯-牛顿法是求解非线性最小二乘问题的经典方法，它可以看成是结合了线搜索的牛顿法的变形. 既然海瑟矩阵中有关 $r_i(x)$ 的二阶导数项不易求出，高斯-牛顿法不去计算这一部分，直接使用 $J(x{)}^{\mathrm{T}}J(x)$ 作为海瑟矩阵的近似矩阵来求解牛顿方程. 用 $J^k$ 简记 $J(x^k)$. 高斯-牛顿法产生的下降方向 $d^k$ 满足

$$(J^k{)}^{\mathrm{T}}{J}^k{d}^k=-(J^k{)}^{\mathrm{T}}{r}^k(\mathrm{6.7.3})$$

上述方程正是法方程的形式，而由线性代数的知识可知，不管 $J^k$ 是否是满秩矩阵，方程 (6.7.3) 一定存在解，而解可以通过伪逆矩阵或最小二乘方法得到。实际上，该方程是如下线性最小二乘问题的最优性条件：
$$\underset{d}{\min }\frac{1}{2}\Vert J^{k}d+r^k{\Vert }^2$$

在求解线性最小二乘问题时，只需要对 $J^k$ 做 QR 分解，因此矩阵 $(J^k{)}^{\mathrm{T}}{J}^k$ 无需计算出来.

高斯-牛顿法的框架如算法 6.10. 为了方便理解，我们将求解线性最小二乘问题的方法进行了展开. 高斯-牛顿法每一步的运算量是来自计算残差向量 $r^k$ 和雅可比矩阵 $J^k$, 和其他算法相比，它的计算量较小. 若 $J^k$ 是满秩矩阵，则高斯-牛顿法得到的方向 $d^k$ 总是一个下降方向，这是因为
$$(d^k{)}^{\mathrm{T}}\nabla f(x^k)=(d^k{)}^{\mathrm{T}}(J^k{)}^{\mathrm{T}}{r}^k=-\Vert J^k{d}^k{\Vert }^2<0$$

这也是高斯-牛顿法的优点. 在此之前的牛顿法，并不总是保证 $d^k$ 是下降方向. 而高斯-牛顿法使用一个半正定矩阵来近似牛顿矩阵，可以获得较好的下降方向.

# 算法 6.10 高斯-牛顿法
import numpy as np
import matplotlib.pyplot as plt
# 生成随机数据
np.random.seed(0)
x = np.random.randn(100, 1)
noise = np.random.normal(0, 1, (100, 1))
para_true = np.array([3.0, 4.0])
y = para_true[0] * x**2 + para_true[1] * x + noise
def func(x, para):
    return para[0] * x**2 + para[1] * x
def jacobi(x, para):
    J = np.zeros((len(x), len(para)))
    J[:, 0] = x[:, 0]**2
    J[:, 1] = x[:, 0]
    return J
def gauss_newton(x, y, max_iter=10000):
    para = np.ones((2, 1))  # 参数初始化
    for i in range(max_iter):
        r = (y - func(x, para))  # 误差
        J = -jacobi(x, para)
        # QR分解
        Q, R = np.linalg.qr(J)
        # 计算下降方向
        d = -np.linalg.inv(R) @ Q.T @ r
        alpha = 0.01 # 步长
        para_next = para + alpha * d
        if np.linalg.norm(para_next - para) < 1e-6:  # 判断迭代终止条件
            break
        para = para_next

    return para

para_est = gauss_newton(x, y)
print("参数真实值：", para_true)
'''参数真实值： [3. 4.]'''
print("参数估计值：", para_est)
'''参数估计值： [[2.98661373]
 [4.12147949]]'''
# 绘制拟合曲线
plt.figure()
plt.scatter(x, y, color='red', label='Data points', s = 20)
x_range = np.linspace(-3, 3, 100)
y_fit = func(x_range, para_est)
plt.plot(x_range, y_fit, color='blue', label='Fitted curve')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

一个很自然的问题是：高斯-牛顿法使用了近似矩阵来求解牛顿方程，那么在什么情况下这个近似是合理的？直观上看，根据海瑟矩阵 (6.7.2b) 的表达式 ${\nabla }^{2}f(x)=J(x{)}^{\mathrm{T}}J(x)+{\sum }_{i=1}^m{r}_i(x){\nabla }^2{r}_i(x)$，当 $(J^k{)}^{\mathrm{T}}{J}^k$ 这一部分起主导时，所使用的近似是有意义的. 一个充分条件就是在最优点 $x^{\ast }$ 处 $r_i(x^{\ast })$ 的值都很小. 此时高斯-牛顿法和牛顿法相近，它们也有很多相似的性质. 如果残差向量 $r(x^{\ast })$ 模长较大，则仅仅使用 $(J^k{)}^{\mathrm{T}}{J}^k$ 并不能很好地近似 ${\nabla }^{2}f(x^k)$, 此时高斯-牛顿法可能收敛很慢甚至发散.

接下来给出高斯-牛顿法的收敛性质. 通过上面的描述可以注意到，雅可比矩阵 $J^k$ 的非奇异性是一个很关键的因素，因此我们在这个条件下建立收敛性. 假设雅可比矩阵 $J(x)$ 的奇异值一致地大于 0, 即存在 $\gamma >0$ 使得
$$\Vert J(x)z\Vert ⩾\gamma \Vert z\Vert ,\forall x\in \mathcal{N},\forall z\in {\mathbb{R}}^n,(\mathrm{6.7.4})$$

其中 $\mathcal{N}$ 是下水平集
L = { x ∣ f ( x ) ⩽ f ( x 0 ) } ( 6.7.5 ) \mathcal{L}=\{x\mid f(x)\leqslant f(x^{0})\}\qquad(6.7.5) L={x∣f(x)⩽f(x0)}(6.7.5)

的一个邻域，$x^0$ 是算法的初始点，且假设 $\mathcal{L}$ 是有界的.

在前面的假设下，有如下收敛性定理：

定理 6.17 (全局收敛性) 如果每个残差函数 $r_j$ 在有界下水平集 (6.7.5) 的一个邻域 $\mathcal{N}$ 内是利普希茨连续可微的，并且雅可比矩阵 $J(x)$ 在 $\mathcal{N}$ 内满足一致满秩条件 (6.7.4), 而步长满足 Wolfe 准则 (6.1.4), 则对高斯-牛顿法得到的序列 { x k } \{x^k\} {xk}
$$\begin{array}{r}\underset{k\to \infty }{\lim }(J^k{)}^{\mathrm{T}}{r}^k=0.\end{array}$$

证明
在这里直接验证 Zoutendijk 条件成立即可.
首先，选择有界下水平集 $\mathcal{L}$ 的邻域 $\mathcal{N}$ 足够小，从而使得存在 $L>0,\beta >0$, 对于任何 $x,\tilde{x}\in \mathcal{N}$ 以及任意的 $j=1,2,\cdots ,m$, 以下条件被满足：
$$\begin{array}{rl}\Vert r_j(x)\Vert & ⩽\beta ,\Vert \nabla r_j(x)\Vert ⩽\beta ,\\ |r_j(x)-r_j(\tilde{x})|& ⩽L\Vert x-\tilde{x}\Vert ,\Vert \nabla r_j(x)-\nabla r_j(\tilde{x})\Vert ⩽L\Vert x-\tilde{x}\Vert .\end{array}$$

那么对任意的 $x\in \mathcal{L}$ 存在 $\tilde{\beta }$ 使得 $\Vert J(x)\Vert =\Vert J(x{)}^{\mathrm{T}}\Vert ⩽\tilde{\beta }$, 并且 $\nabla f(x)=J(x{)}^{\mathrm{T}}r(x)$ 是利普希茨连续函数.
那么 Zoutendijk 条件满足，成立：$$\sum _{k=0}^{\infty }{\cos }^2{\theta }_k\Vert \nabla f(x_k){\Vert }^2<+\infty$$

其中 ${\theta }_k$ 是高斯-牛顿方向 $d^k$ 与负梯度方向的夹角，则
$$\cos {\theta }_k=-\frac{\nabla f(x^k{)}^{\text{T}}{d}^k}{\Vert d^k\Vert \Vert \nabla f(x^k)\Vert }=\frac{\Vert J^k{d}^k{\Vert }^2}{\Vert d^k\Vert \Vert (J^k{)}^{\text{T}}{J}^k{d}^k\Vert }⩾\frac{{\gamma }^2\Vert d^k{\Vert }^2}{{\tilde{\beta }}^2\Vert d^k{\Vert }^2}=\frac{{\gamma }^2}{{\tilde{\beta }}^2}>0.$$

根据推论 6.1(线搜索算法的收敛性) 即可得 $\nabla f(x^k)\to 0.$

定理 6.17 的关键假设是一致满秩条件 (6.7.4). 实际上，若 $J^k$ 不满秩， 则更新方向 $d_k$ 有无穷多个解. 如果对解的性质不提额外要求，则无法推出 $\cos {\theta }_k$ 一致地大于零. 此时收敛性可能不成立.

当 $(J^k{)}^{\mathrm{T}}{J}^k$ 在海瑟矩阵 (6.7.2b) 中占据主导部分时（当 ${\sum }_{i=1}^m{r}_i(x){\nabla }^2{r}_i(x)$ 较小时），高斯-牛顿算法可能会有更快的收敛速度. 类似于牛顿法，我们给出高斯-牛顿法的局部收敛性.

定理 6.18 (局部收敛性) 设 $r_i(x)$ 二阶连续可微，$x^{\ast }$ 是最小二乘问题 (6.7.1) 的最优解，海瑟矩阵 ${\nabla }^{2}f(x)$ 和其近似矩阵 $J(x{)}^{\mathrm{T}}J(x)$ 均在点 $x^{\ast }$ 的一个邻域内利普希茨连续，则当高斯-牛顿算法步长 ${\alpha }_k$ 恒为 1 时，
$$\Vert x^{k+1}-x^{\ast }\Vert ⩽C\Vert ((J^{\ast }{)}^{\text{T}}{J}^{\ast }{)}^{-1}{H}^{\ast }\Vert \Vert x^k-x^{\ast }\Vert +\mathcal{O}(\Vert x^k-x^{\ast }{\Vert }^2)(\mathrm{6.7.8})$$

其中 $H^{\ast }=\sum _{i=1}^m{r}_i(x^{\ast }){\nabla }^2{r}_i(x^{\ast })$ 为海瑟矩阵 ${\nabla }^{2}f(x^{\ast })$ 去掉 $J(x^{\ast }{)}^{\mathrm{T}}J(x^{\ast })$ 的部分，$C>0$ 为常数.

证明
根据迭代公式，
$$x^{k+1}-x^k=x^k+d^k-x^{\ast }=((J^k{)}^{\mathrm{T}}{J}^k{)}^{-1}((J^k{)}^{\mathrm{T}}{J}^k(x^k-x^{\ast })+\nabla f(x^{\ast })-\nabla f(x^k))(\mathrm{6.7.9})$$

由泰勒展开，
$$\nabla f(x^k)-\nabla f(x^{\ast })={\int }_0^1{J}^{\mathrm{T}}J(x^{\ast }+t(x^k-x^{\ast }))(x^k-x^{\ast })dt+{\int }_0^{1}H(x^{\ast }+t(x^k-x^{\ast }))(x^k-x^{\ast })dt$$

其中 $J^{\mathrm{T}}J(x)$ 是 $J^{\mathrm{T}}(x)J(x)$ 的简写，$H(x)={\nabla }^{2}f(x)-J^{\mathrm{T}}J(x)$ 为海瑟矩阵剩余部分. 将泰勒展开式代入(6.7.9)式右边，取范数进行估计，有
$$\begin{array}{rl}& \Vert (J^k{)}^{\mathrm{T}}{J}^k(x^k-x^{\ast })+\nabla f(x^{\ast })-\nabla f(x^k)\Vert \\ & ⩽{\int }_0^1\Vert (J^{\mathrm{T}}J(x^k)-J^{\mathrm{T}}J(x^{\ast }+t(x^k-x^{\ast })))(x^k-x^{\ast })\Vert \mathrm{d}t+{\int }_0^1\Vert H(x^{\ast }+t(x^k-x^{\ast }))(x^k-x^{\ast })\Vert \mathrm{d}t\\ & ⩽\frac{L}{2}\Vert x^k-x^{\ast }{\Vert }^2+C\Vert H^{\ast }\Vert \Vert x^k-x^{\ast }\Vert ,\end{array}$$

其中 $L$ 是 $J^{\mathrm{T}}J(x)$ 的利普希茨常数. 最后一个不等式是因为我们使用 $H^{\ast }$ 来近似 $H(x^{\ast }+t(x^k-x^{\ast }))$, 由连续性，存在 $C>0$ 以及点 $x^{\ast }$ 的一个邻域 $\mathcal{N}$, 对任意的 $x\in \mathcal{N}$ 有 $\Vert H(x)\Vert ⩽C\Vert H(x^{\ast })\Vert$. 将上述估计代入 (6.7.9) 式即可得到(6.7.8)式.

定理 6.18 指出，如果 $\Vert H(x^{\ast })\Vert$ 充分小，则高斯–牛顿法可以达到 Q-线性收敛速度；而当 $\Vert H(x^{\ast })\Vert =0$ 时，收敛速度是 Q-二次的．如果 $\Vert H(x^{\ast })\Vert$ 较大，则高斯–牛顿法很可能会失效．

6.7.3 Levenberg-Marquardt 方法

1. 信赖域型 LM 方法
Levenberg-Marquardt (LM) 方法是由 Levenberg 在 1944 年提出的求解非线性最小二乘问题的方法．它本质上是一种信赖域型方法，主要应用场合是当矩阵 $(J^k{)}^{\mathrm{T}}{J}^k$ 奇异时，它仍然能给出一个下降方向．LM 方法每一步求解如下子问题：

$$\underset{d}{\min }\frac{1}{2}\Vert J^{k}d+r^k{\Vert }^2,s.t.\Vert d\Vert ⩽{\Delta }_k(\mathrm{6.7.10})$$

事实上，LM 方法将如下近似当作信赖域方法中的 $m_k$:

$$m_k(d)=\frac{1}{2}\Vert r^k{\Vert }^2+d^{\mathrm{T}}(J^k{)}^{\mathrm{T}}{r}^k+\frac{1}{2}{d}^{\mathrm{T}}(J^k{)}^{\mathrm{T}}{J}^{k}d(\mathrm{6.7.11})$$

该方法使用 $B^k=(J^k{)}^{\mathrm{T}}{J}^k$ 来近似海瑟矩阵，这个取法是从高斯-牛顿法推广而来的. LM 方法并不直接使用海瑟矩阵来求解. 以下为了方便，省去迭代指标 $k$. 子问题 (6.7.10) 是信赖域子问题，上一节讨论过这个子问题的一些好的性质. 根据定理 6.12, 可直接得到如下推论：

推论 6.4 向量 $d^{\ast }$ 是信赖域子问题
$$\begin{array}{r}\underset{d}{\min }\frac{1}{2}\Vert Jd+r{\Vert }^2,\text{s.t.}\Vert d\Vert ⩽\Delta \end{array}$$

的解当且仅当 $d^{\ast }$ 是可行解并且存在数 $\lambda ⩾0$ 使得
$$\begin{array}{rlr}& (J^{\mathrm{T}}J+\lambda I)d^{\ast }=-J^{\mathrm{T}}r& (\mathrm{6.7.12})\\ & \lambda (\Delta -\Vert d^{\ast }\Vert )=0& (\mathrm{6.7.13})\end{array}$$

注意到 $J^{\mathrm{T}}J$ 是半正定矩阵，因此条件 (6.6.5c) 是自然成立的.
下面简要说明如何求解 LM 子问题 (6.7.10). 实际上，和信赖域子问题中的迭代法相同，我们先通过求根的方式来确定 $\lambda$ 的选取，然后直接求得 LM 方程的迭代方向. 由于 LM 子问题的特殊性，可以不显式求出矩阵 $J^{\mathrm{T}}J+\lambda I$ 的 Cholesky 分解(Cholesky分解是一种将对称正定矩阵分解为一个下三角矩阵和其转置的乘积的方法)，而仍然是借助 QR 分解(QR分解是一种常用的矩阵分解方法，它将一个矩阵分解为一个正交矩阵(Q)与一个上三角矩阵®的乘积)，进而无需算出 $J^{\mathrm{T}}J+\lambda I$．注意，$(J^{\mathrm{T}}J+\lambda I)d=-J^{\mathrm{T}}r$ 实际上是最小二乘问题
$$\underset{d}{\min }\Vert \left[\begin{array}{c}J\\ \sqrt{\lambda }I\end{array}\right]d+\left[\begin{array}{c}r\\ 0\end{array}\right]\Vert (\mathrm{6.7.14})$$

的最优性条件．此问题的系数矩阵带有一定结构，每次改变 $\lambda$ 进行试探时，有关 $J$ 的块是不变的，因此无需重复计算 $J$ 的 QR 分解. 具体来说，设 $J=QR$ 为 $J$ 的 QR 分解，其中 $Q\in {\mathbb{R}}^{m\times n},R\in {\mathbb{R}}^{n\times n}$. 我们有
$$\left[\begin{array}{c}J\\ \sqrt{\lambda }I\end{array}\right]=\left[\begin{array}{c}QR\\ \sqrt{\lambda }I\end{array}\right]=\left[\begin{array}{cc}Q& 0\\ 0& I\end{array}\right]\left[\begin{array}{c}R\\ \sqrt{\lambda }I\end{array}\right]$$

矩阵 $\left[\begin{array}{c}R\\ \sqrt{\lambda }I\end{array}\right]$ 含有较多零元素，可以使用 Householder 变换或 Givens 变换来完成此矩阵的 QR 分解. 如果矩阵 $J$ 没有显式形式，只能提供矩阵乘法，则仍然可以用截断共轭梯度法，即算法 6.9(截断共轭梯度法) 来求解子问题 (6.7.10).

补充：QR 分解
QR 分解是一种重要的矩阵分解方法，通常用于解决线性方程组、最小二乘问题以及特征值计算等数值计算任务。它将一个矩阵分解为一个正交矩阵(Q)与一个上三角矩阵®的乘积。QR 分解有多种形式。
Householder 变换：
Householder 变换通过一系列反射操作将矩阵转化为上三角形式。每次变换都是通过一个特殊的 Householder 矩阵来实现，该矩阵是一个正交矩阵，并且将一个向量反射到另一个向量上。通过一系列的 Householder 变换，可以将矩阵转化为上三角形式，从而实现 QR 分解。
Givens 变换:
Givens 变换通过一系列的旋转操作将矩阵转化为上三角形式。在 Givens 变换中，每次旋转操作都是基于一个特定的 Givens 矩阵，该矩阵是一个二维旋转矩阵，用于将一个特定的位置上的两个元素变为零。

LM 方法的收敛性也可以直接从信赖域方法的收敛性得出，我们直接给出收敛性定理.

定理 6.19 假设常数 $\eta \in \left(0,\frac{1}{4}\right)$, 下水平集 $\mathcal{L}$ 是有界的且每个 $r_i(x)$ 在下水平集 $\mathcal{L}$ 的一个邻域 $\mathcal{N}$ 内是利普希茨连续可微的. 假设对于任意的 $k$, 子问题 (6.7.10) 的近似解 $d_k$ 满足
$$m_k(0)-m_k(d^k)⩾c_1\Vert (J^k{)}^{\mathrm{T}}{r}^k\Vert \min \{{\Delta }_k,\frac{\Vert (J^k{)}^{\mathrm{T}}{r}^k\Vert }{\Vert (J^k{)}^{\mathrm{T}}{J}^k\Vert }\}$$

其中 $c_1>0$ 且 $\Vert d^k\Vert ⩽\gamma {\Delta }_k,\gamma ⩾1$, 则

$$\underset{k\to \infty }{\lim }\nabla f(x^k)=\underset{k\to \infty }{\lim }(J^k{)}^{\text{T}}{r}^k=0$$

2. LMF 方法
信赖域型 LM 方法本质上是固定信赖域半径 $\Delta$, 通过迭代寻找满足条件 (6.7.12) 的乘子 $\lambda$, 每一步迭代需要求解线性方程组
$$(J^{\mathrm{T}}J+\lambda I)d=-J^{\mathrm{T}}r$$

这个求解过程在 LM 方法中会占据相当大的计算量，计算代价较大.

在 LM 方法中，由于 $(J^k{)}^{\mathrm{T}}{J}^k\succ 0$, 那么有 $-{\lambda }_1<0$，此时有 $\lambda >-{\lambda }_1$, 因此若 $\lambda$ 越大，$d$ 的模长越小. 这就意味着 $\Delta$ 的大小被 $\lambda$ 隐式地决定，直接调整 $\lambda$ 的大小就相当于调整了信赖域半径 $\Delta$ 的大小. 因此，我们可构造基于 $\lambda$ 更新的 LM 方法．由于 LM 方程 (6.7.12) 和信赖域子问题的关系由 Fletcher 在 1981 年给出，因此基于 $\lambda$ 更新的 LM 方法也被称为 LMF 方法，即每一步求解子问题：

$$\underset{d}{\min }\frac{1}{2}\Vert Jd+r{\Vert }_2^2+\lambda \Vert d{\Vert }_2^2$$

在 LMF 方法中，设第 $k$ 步产生的迭代方向为 $d^k$, 根据信赖域算法的思想，我们需要计算目标函数的预估下降量和实际下降量的比值 ${\rho }_k$, 来确定下一步信赖域半径的大小. 这一比值很容易通过公式 (6.6.3) 计算，其中 $f(x)$ 和 $m_k(d)$ 分别取为 (6.7.1) 式和 (6.7.11) 式. 即 ${\rho }_k=\frac{f(x^k)-f(x^k+d^k)}{m_k(0)-m_k(d^k)}$, $f(x)=\frac{1}{2}\sum _{j=1}^m{r}_j^2(x)$, $m_k(d)=\frac{1}{2}\Vert r^k{\Vert }^2+d^{\mathrm{T}}(J^k{)}^{\mathrm{T}}{r}^k+\frac{1}{2}{d}^{\mathrm{T}}(J^k{)}^{\mathrm{T}}{J}^{k}d$. 计算 ${\rho }_k$ 后，我们根据 ${\rho }_k$ 的大小来更新 ${\lambda }_k$. 当乘子 ${\lambda }_k$ 增大时，信赖域半径会变小，反之亦然. 所以 ${\lambda }_k$ 的变化策略应该和信赖域算法中 ${\Delta }_k$ 的恰好相反.
接下来就给出了 LMF 算法的框架. 通过比较得知 LMF 方法 (算法 6.11) 和信赖域算法 6.8 的结构非常相似. 算法 6.11 对 ${\gamma }_1$, ${\gamma }_2$ 等参数并不敏感. 但根据信赖域方法的收敛定理 6.15, 参数 $\eta$ 可以取大于 0 的值来改善收敛结果．

和 LM 方法相比，LMF 方法每一次迭代只需要求解一次 LM 方程，从而极大地减少了计算量，在编程方面也更容易实现．所以 LMF 方法在求解最小二乘问题中是很常见的做法．

# 算法 6.11 LMF 方法
import numpy as np
import matplotlib.pyplot as plt
def lmf(x, y, para0, lambda0, eta, rho1_bar = 0.25, rho2_bar = 0.75, gamma1 = 0.1, gamma2 = 2.0, max_iter=1000):
    para = para0
    lambda_k = lambda0
    k = 0
    para_history = [para]
    while k < max_iter:
        r = y - func(x, para)
        J_k = -jacobi(x, para)
        d = -np.linalg.solve((J_k.T).dot(J_k) + lambda_k * np.eye(len(para)), (J_k.T).dot(r))
        rho_k = (np.linalg.norm(r)**2 - np.linalg.norm(r + J_k.dot(d))**2) / (0.5 * np.linalg.norm(r)**2)
        
        if rho_k < rho1_bar:
            lambda_k = gamma2 * lambda_k
        elif rho_k > rho2_bar:
            lambda_k = gamma1 * lambda_k
        else:
            lambda_k = lambda_k
        
        if rho_k > eta:
            para_next = para + d
        
        if np.linalg.norm(para_next - para) < 1e-6:
            break
        
        para = para_next
        para_history.append(para)
        k = k + 1
    return para_history, k 

# 生成随机数据
np.random.seed(0)
x = np.random.randn(100, 1)
noise = np.random.normal(0, 1, (100, 1))
para_true = np.array([3.0, 4.0])
y = para_true[0] * x**2 + para_true[1] * x + noise
def func(x, para):
    return para[0] * x**2 + para[1] * x
def jacobi(x, para):
    J = np.zeros((len(x), len(para)))
    J[:, 0] = x[:, 0]**2
    J[:, 1] = x[:, 0]
    return J
para0 = np.ones((2, 1))  # 参数初始化
lambda0 = 1.0
eta = 0.1

para_hist, k = lmf(x, y, para0, lambda0, eta)
print("参数真实值：", para_true)
'''参数真实值： [3. 4.]'''
print("参数估计值：", para_hist[-1])
'''参数估计值： [[2.98142607]
 [4.09219555]]'''
print("迭代次数：", k)
'''迭代次数： 1'''
print(para_hist)
'''[array([[1.],
       [1.]]), array([[2.98142607],
       [4.09219555]])]'''
# 绘制拟合曲线
plt.figure()
plt.scatter(x, y, color='red', label='Data points', s=20)
x_range = np.linspace(-3, 3, 100)
y_fit = func(x_range, para_est)
plt.plot(x_range, y_fit, color='blue', label='Fitted curve')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

6.7.4 大残量问题的拟牛顿法

在大残量问题中，海瑟矩阵 ${\nabla }^{2}f(x)$ 的第二部分的作用不可忽视，仅仅考虑 $(J^k{)}^{\mathrm{T}}{J}^k$ 作为第 $k$ 步的海瑟矩阵近似则会带来很大误差. 在这种情况下高斯-牛顿法和 LM 方法很可能会失效. 我们可以将最小二乘问题 (6.7.1) 当成一般的无约束问题，并使用之前讨论过的牛顿法和拟牛顿法求解. 但对于很多问题来说，各个残量分量的海瑟矩阵 ${\nabla }^2{r}_i(x)$ 不易求出，使用牛顿法会有很大开销；而直接使用拟牛顿法对海瑟矩阵 ${\nabla }^{2}f(x)$ 进行近似又似乎忽略了最小二乘问题的特殊结构.
海瑟矩阵表达式 ${\nabla }^{2}f(x)=J(x{)}^{\mathrm{T}}J(x)+\sum _{i=1}^m{r}_i(x){\nabla }^2{r}_i(x)$ (6.7.2b) 说明了 ${\nabla }^{2}f(x)$ 由两部分组成：一部分容易得出，但不精确；另一部分较难求得，但在计算中又必不可少. 对于容易部分可以直接保留高斯-牛顿矩阵 $J^{\mathrm{T}}J$, 而对于较难部分则可以利用拟牛顿法来进行近似. 这就是我们求解大残量问题的基本思路，它同时考虑了最小二乘问题的海瑟矩阵结构和计算量，是一种混合的近似方法.

具体来说，我们使用 $B^k$ 来表示 ${\nabla }^{2}f(x^k)$ 的近似矩阵，即$
B^k=(Jk)^\mathrm{T}Jk+T^k,$ 其中 $T^k$ 是海瑟矩阵第二部分 $\sum _{j=1}^m{r}_j(x^k){\nabla }^2{r}_j(x^k)$ 的近似. 问题的关键在于如何构造矩阵 $T^k$, 建立拟牛顿法时，构造拟牛顿格式主要分为两步：一是找出拟牛顿条件，二是根据拟牛顿条件来构造拟牛顿矩阵的低秩更新. 在这里我们使用相似的过程，但注意 $T^k$ 仅仅是拟牛顿矩阵 $B^k$ 的一部分，它不太可能满足割线条件 (6.5.3)（将其中的 $B^{k+1}$ 替换成 $T^{k+1}$）．我们的目标是让 $T^{k+1}$ 和海瑟矩阵的第二部分尽量相似，即
$$T^{k+1}\approx \sum _{j=1}^m{r}_j(x^{k+1}){\nabla }^2{r}_j(x^{k+1}).$$

记 $s^k=x^{k+1}-x^k$，由一阶泰勒展开得知，$T^{k+1}$ 应该尽量保留原海瑟矩阵的性质，即

$$\begin{array}{rl}{T}^{k+1}{s}^k& \approx \left(\sum _{j=1}^m{r}_j(x^{k+1}){\nabla }^2{r}_j(x^{k+1})\right)s^k\\ & =\sum _{j=1}^m{r}_j(x^{k+1})\left({\nabla }^2{r}_j(x^{k+1})\right)s^k\\ & \approx \sum _{j=1}^m{r}_j(x^{k+1})(\nabla r_j(x^{k+1})-\nabla r_j(x^k))\\ & =(J^{k+1}{)}^{\mathrm{T}}{r}^{k+1}-(J^k{)}^{\mathrm{T}}{r}^{k+1}.\end{array}$$

令 ${\hat{y}}^k=(J^{k+1}{)}^{\mathrm{T}}{r}^{k+1}-(J^k{)}^{\mathrm{T}}{r}^{k+1}$, 则 $T^k$ 满足的拟牛顿条件为

$$T^{k+1}{s}^k={\hat{y}}^k.(\mathrm{6.7.15})$$

在这里注意 ${\hat{y}}^k$ 不是原有的 $y^k.$
有了拟牛顿条件后，可以用之前的方法构造拟牛顿格式了。
Dennis, Gay 和 Welsch 给出的一种更新格式为：$$T^{k+1}=T^k+\frac{(y^{\#}-T^k{s}^k)y^{\mathrm{T}}+y(y^{\#}-T^k{s}^k{)}^{\mathrm{T}}}{y^{\mathrm{T}}{s}^k}-\frac{(y^{\#}-T^k{s}^k{)}^{\mathrm{T}}{s}^k}{(y^{\mathrm{T}}s{)}^2}y{y}^{\mathrm{T}}(\mathrm{6.7.16})$$

其中 $$\begin{array}{rl}& s^k=x^{k+1}-x^k\\ & y=(J^{k+1}{)}^{\mathrm{T}}{r}^{k+1}-(J^k{)}^{\mathrm{T}}{r}^k\\ & y^{\#}=(J^{k+1}{)}^{\mathrm{T}}{r}^{k+1}-(J^k{)}^{\mathrm{T}}{r}^{k+1}\end{array}$$.

参考教材：《最优化：建模、算法与理论》