最优化学习笔记：无约束优化算法（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 =