非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (IV)

Title: 非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (IV)

姊妹博文**

非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (I)

非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (II)

非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (III)

非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (IV) $⟵$ 本篇

$\uparrow$ 理论部分

---

$\downarrow$ 实例部分

非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (实例篇 V)

---

目录

I.引言
II.牛顿迭代法解非线性方程
   1.一元非线性方程形式的牛顿迭代法
   2.高维非线性方程组形式的牛顿迭代法
   3.牛顿迭代法的雅可比矩阵
III.牛顿迭代法解非线性最小二乘问题
   1.从方程问题到最小二乘问题的演化
   2.最小二乘问题 Jacobian 矩阵的推导
   3.最小二乘问题 Jacobian 矩阵的性质
     A.维度
     B.整体 Hessian 矩阵
     C.局部 Hessian 矩阵
     D.Jacobian 矩阵/Hessian 矩阵
     E.Hessian 矩阵的对称性
IV.高斯-牛顿法解非线性最小二乘问题
   1.高斯-牛顿法的获得
   2.高斯-牛顿法的优势
   3.高斯-牛顿法的解读 —— 优化观点
V.最小二乘法与高斯的贡献
VI.总结
参考文献

---

[点击回到上一章节]

IV. 高斯-牛顿法解非线性最小二乘问题

1. 高斯-牛顿法的获得

结合式 (III-2-9) 和式 (III-3-3) 得到
$${\mathbf{x}}_{[i+1]}={\mathbf{x}}_{[i]}-\mathbf{H}({\mathbf{x}}_{[i]}{)}^{+}\nabla g({\mathbf{x}}_{[i]})\text{\hspace{1em}(IV-1-1)}$$
由式 (III-3-3) 可知, 牛顿迭代法的整体 Hessian 矩阵 (Jacobian 矩阵) 分为两个部分.

而高斯-牛顿法直接将 “第二部分” 省略, 而只保留 “第一部分”.

定义简化 Hessian 矩阵
$$\tilde{\mathbf{H}}({\mathbf{x}}_{[i]})\triangleq {{\left[\frac{\partial \mathbf{r}(\mathbf{x})}{\partial \mathbf{x}}\right]}^{\mathrm{T}}\frac{\partial \mathbf{r}(\mathbf{x})}{\partial \mathbf{x}}|}_{{\mathbf{x}}_{[i]}}\text{\hspace{1em}(IV-1-2)}$$
上述简化后的 Hessian 矩阵有一个很好的性质:

在 ${\frac{\partial \mathbf{r}(\mathbf{x})}{\partial \mathbf{x}}|}_{{\mathbf{x}}_{[i]}}$ 列满秩情况下, ${{\left[\frac{\partial \mathbf{r}(\mathbf{x})}{\partial \mathbf{x}}\right]}^{\mathrm{T}}\frac{\partial \mathbf{r}(\mathbf{x})}{\partial \mathbf{x}}|}_{{\mathbf{x}}_{[i]}}$ 正定, 则 $\tilde{\mathbf{H}}({\mathbf{x}}_{[i]})$ 可逆.

(为了保证列满秩衍生出各种高斯-牛顿变形算法, 此处不展开)

将简化 Hessian 矩阵替换式 (IV-1-1) 中的整体 Hessian 矩阵得到高斯-牛顿法的迭代更新公式
$$\begin{array}{rl}{\mathbf{x}}_{[i+1]}& ={\mathbf{x}}_{[i]}-\tilde{\mathbf{H}}({\mathbf{x}}_{[i]}{)}^{-1}\nabla g({\mathbf{x}}_{[i]})\\ & ={\mathbf{x}}_{[i]}-{\left({{\left[\frac{\partial \mathbf{r}(\mathbf{x})}{\partial \mathbf{x}}\right]}^{\mathrm{T}}\frac{\partial \mathbf{r}(\mathbf{x})}{\partial \mathbf{x}}|}_{{\mathbf{x}}_{[i]}}\right)}^{-1}{\left({\left[\frac{\partial \mathbf{r}(\mathbf{x})}{\partial \mathbf{x}}\right]}^{\mathrm{T}}\mathbf{r}(\mathbf{x}))\right|}_{{\mathbf{x}}_{[i]}}\end{array}\text{\hspace{1em}(IV-1-3)}$$
除了 Hessian 矩阵的替换, 另外已用普通的矩阵逆代替矩阵伪逆, 而其他部分和牛顿迭代法保持一致, 这样就得到了高斯-牛顿法.

这样简化的依据/考虑^[4,5]是什么?

[1] 接近极小值 (最优参数) ${\mathbf{x}}^{\ast }$ 时, 预测值与测量值之间的偏差 $r_i(\mathbf{x})$ 很小

[2] 接近极小值 ${\mathbf{x}}^{\ast }$ 时, $\frac{{\partial }^2{r}_i}{\partial x_j\partial x_k}$ 的值很小. 此时 $r_i(\mathbf{x})$ 近似呈仿射函数形式 (一阶泰勒展开)
$$r_i(\mathbf{x})\approx r_i({\mathbf{x}}^{\ast })+\nabla r_i({\mathbf{x}}^{\ast })(\mathbf{x}-{\mathbf{x}}^{\ast })\text{\hspace{1em}(IV-1-4)}$$
对其进行二阶导数得到, $\mathbf{D}(\mathbf{x})\approx 0$.

这样未经简化的整体 Hessian 矩阵中的第二部分由两个很小的量相乘组合而成, 第二部分进一步迅速减小.

可以说在接近极小值 ${\mathbf{x}}^{\ast }$ 时, 式 (III-3-3) 中的第一部分占主导, 故省略第二部分.

(这些依据也是建立在比较 “好/正常” 的函数基础上, 或者残差函数的非线性程度比较 “温和”)

---

2. 高斯-牛顿法的优势

如此简化 Hessian 矩阵的优势^[4,5]如下:

[1] 减少公式复杂度, 减少计算量, 这是肉眼可见

[2] 保证可逆, 避免求伪逆, 进一步减少计算量

在 ${\frac{\partial \mathbf{r}(\mathbf{x})}{\partial \mathbf{x}}|}_{{\mathbf{x}}_{[i]}}$ 列满秩情况下, $\tilde{\mathbf{H}}({\mathbf{x}}_{[i]})$ 正定和可逆, 且 $\tilde{\mathbf{H}}({\mathbf{x}}_{[i]}{)}^{-1}$ 也正定.

而牛顿迭代法中, 无法保证 $\mathbf{D}({\mathbf{x}}_{[i]})$ 正定, 故而无法保证整体 Hessian 矩阵 $\mathbf{H}({\mathbf{x}}_{[i]})$ 正定/可逆.

[3] 保证迭代计算向着下降方向进行

所谓下降方向是指, 代价函数 $g(\mathbf{x})$ 随着 $\mathbf{x}$ 的迭代更新而逐渐下降 (逐渐接近极小值). 下面利用一阶泰勒展开验证高斯-牛顿方法的迭代方向.
$$\begin{array}{rl}g({\mathbf{x}}_{[i+1]})& =g\left({\mathbf{x}}_{[i]}-\tilde{\mathbf{H}}({\mathbf{x}}_{[i]}{)}^{-1}\nabla g({\mathbf{x}}_{[i]})\right)\\ & \approx g\left({\mathbf{x}}_{[i]}\right)-\underset{>0}{\underbrace{{\left(\nabla g({\mathbf{x}}_{[i]})\right)}^{\mathrm{T}}\tilde{\mathbf{H}}({\mathbf{x}}_{[i]}{)}^{-1}\nabla g({\mathbf{x}}_{[i]})}}\end{array}\text{\hspace{1em}(IV-2-1)}$$
因为 $\tilde{\mathbf{H}}({\mathbf{x}}_{[i]}{)}^{-1}$ 正定, 故式 (IV-2-1) 中右手边第二项大于零. 那么可知
$$g({\mathbf{x}}_{[i+1]})<g({\mathbf{x}}_{[i]})\text{\hspace{1em}(IV-2-2)}$$
代价函数 $g(\mathbf{x})$ 确实随着 $\mathbf{x}$ 的迭代更新而下降.

而牛顿迭代法中无法保证 $\mathbf{H}({\mathbf{x}}_{[i]})$ 正定, 故也无法保证代价函数向着下降的方向更新.

---

3. 高斯-牛顿法的解读 —— 优化观点

对代价函数在 ${\mathbf{x}}_{[i]}$ 处取二阶泰勒近似
$$g(\mathbf{x})\approx g({\mathbf{x}}_{[i]})+{\left(\nabla g({\mathbf{x}}_{[i]})\right)}^{\mathrm{T}}(\mathbf{x}-{\mathbf{x}}_{[i]})+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_{[i]}{)}^{\mathrm{T}}\mathbf{H}({\mathbf{x}}_{[i]})(\mathbf{x}-{\mathbf{x}}_{[i]})\text{\hspace{1em}(IV-3-1)}$$
由 “极值的必要条件”, 代价函数取极小值需要满足
$$\begin{array}{r}\nabla g(\mathbf{x})=0\\ \Rightarrow \frac{\partial }{\partial \mathbf{x}}\left(g({\mathbf{x}}_{[i]})+{\left(\nabla g({\mathbf{x}}_{[i]})\right)}^{\mathrm{T}}(\mathbf{x}-{\mathbf{x}}_{[i]})+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_{[i]}{)}^{\mathrm{T}}\mathbf{H}({\mathbf{x}}_{[i]})(\mathbf{x}-{\mathbf{x}}_{[i]})\right)=0\\ \Rightarrow {\left(\nabla g({\mathbf{x}}_{[i]})\right)}^{\mathrm{T}}+(\mathbf{x}-{\mathbf{x}}_{[i]}{)}^{\mathrm{T}}\mathbf{H}({\mathbf{x}}_{[i]})=0\\ (\text{Symmetric Hessian})\Rightarrow \mathbf{x}={\mathbf{x}}_{[i]}-\mathbf{H}({\mathbf{x}}_{[i]}{)}^{+}\nabla g({\mathbf{x}}_{[i]})\end{array}\text{\hspace{1em}(IV-3-2)}$$
这样就又到达了高斯-牛顿法的起点式 (IV-1-1).

式 (IV-3-2) 中最后一步得到的 $\mathbf{x}$ 就是基于第 $i$ 步参数迭代值 ${\mathbf{x}}_{[i]}$ 的第 $i+1$ 步参数迭代值. 每次迭代中都会在上一迭代值附近取代价函数 $g(\mathbf{x})$ 的二阶近似, 寻找使得二阶近似代价函数取得极小值所对应的参数点作为新的参数迭代值.

即新的迭代值是都是局部近似最优参数, 其中 “在上一迭代值附近” 体现局部性, “对代价函数进行二阶泰勒近似” 体现近似性, “求使得近似代价函数取得极小值时的参数点” 体现优化性.

这是对牛顿迭代法和高斯-牛顿法在非线性最小二乘法问题解算过程的另一种解读 —— 优化观点, 区别于从解非线性方程推广得到的牛顿迭代法 —— 梯度观点.

---

V. 最小二乘法与高斯的贡献

这个高斯-牛顿法似乎是牛顿迭代法的简单变形而来 (式 (IV-1-2) 中的 Hessian 矩阵简化), 所以我也曾疑问为什么这个方法不叫 “牛顿-高斯法”？

事实上, 这个方法的核心是最小二乘法的创建及其数值计算方法. 将最小二乘问题变换为解方程问题后, 借鉴和引入牛顿迭代法作为该问题的一个计算步骤. 牛顿时代还没有最小二乘法及其在优化问题中的应用. 即使是解非线性方程的工作, 牛顿实际使用的方法和现在版本的 “牛顿迭代法” 也相去甚远^[6]. 对非线性最小二次问题及求解提出创造性推动的是高斯, 所以这样命名也比较公正.

这里简单回顾一下最小二乘法的历史^[7,8]:

最小二乘法 (Least Squares) 是一种使误差平方和达到最小以寻求最优估计值的方法, 其具有重要的理论和应用价值, “最小二乘法之于数理统计学犹如微积分之于数学” (美国统计学家 S. M. Stigler).

1805 年，法国数学家勒让德 (Adrien-Marie Legendre) 在著作《计算彗星轨道的新方法》中第一次公开提出最小二乘法. 勒让德认为: “赋予误差的平方和为极小, 则意味着在这些误差间建立了一种均衡性, 它阻止了极端情形所施加的过分影响. 这非常好地适用于揭示最接近真实情形的系统状态”.

德国数学家高斯 (Carl Friedrich Gauss) 则声称从 1795 年起就开始使用最小二乘法. 1801 年, 意大利天文学家皮亚齐 (Giuseppe Piazzi) 发现了第一颗小行星谷神星. 经过 40 天的跟踪观测后, 由于谷神星运行至太阳背后, 使得皮亚齐失去了谷神星的位置. 随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星, 但根据大家的计算结果来寻找谷神星都无果. 时年 24 岁的高斯利用自己创立的最小二乘法也计算了谷神星的轨道. 奥地利天文学家奥尔伯斯 (Heinrich Wilhelm Matthias Olbers) 根据高斯计算的轨道重新发现了谷神星. 高斯使用的最小二乘法的方法发表于1809 年他的著作《天体运动论》, 高斯-牛顿法就首现于该著作.

勒让德和高斯之间关于最小二乘法的优先权之争延续了多年. 现在一般认为, 二人各自独立地发明了最小二乘法, 但高斯较之于勒让德把最小二乘法推进得更远. 高斯由误差函数推导出最小二乘法并详尽阐述了该方法的理论依据, 包括 (1) 正态误差理论, (2) 高斯-马尔可夫定理 (最佳线性无偏估计).

以下引用大英百科全书的表述:

Least Squares Method^[9]

In 1805 the French mathematician Adrien-Marie Legendre published the first known recommendation to use the line that minimizes the sum of the squares of these deviations—i.e., the modern least squares method. The German mathematician Carl Friedrich Gauss, who may have used the same method previously, contributed important computational and theoretical advances. The method of least squares is now widely used for fitting lines and curves to scatterplots (discrete sets of data).

---

VI. 总结

本篇博客按照如下先后顺序回顾了, 高斯-牛顿法这一广泛应用于最小二乘问题求解的数值方法的来龙去脉:

- 牛顿迭代法解一元非线性方程;

- 牛顿迭代法推广到高维形式;

- 最小二乘问题的牛顿迭代法求解;

- 由牛顿迭代法推出高斯-牛顿法;

- 高斯-牛顿法的依据、优势、优化观点下的解读;

- 最小二乘法历史与高斯的贡献.

需要强调的是, 高斯-牛顿法的应用基于如下先决条件:

- 误差函数 $r_i$ ($i=1,2,\cdots ,m$) 在包含最优参数 ${\mathbf{x}}^{\ast }$ 的局部凸集合内部是二次连续可微的;

- 迭代初始点 ${\mathbf{x}}_{[0]}$ 接近于局部极小点 ${\mathbf{x}}^{\ast }$;

- ${\frac{\partial \mathbf{r}(\mathbf{x})}{\partial \mathbf{x}}|}_{{\mathbf{x}}^{\ast }}$ 列满秩 (其实每一步都要求列满秩, 只是各迭代值接近于 ${\mathbf{x}}^{\ast }$ 的缘故, 利用连续性, 故只要求 ${\mathbf{x}}^{\ast }$ 上的列满秩);

- 局部极小值 $|g(\mathbf{x})|$ 很小 (足够小).

需要说明:

- 高斯-牛顿法得到的是局部最优值, 非全局最优;

- 高斯-牛顿法的收敛速度、改进版本等此处都没有涉及.

(谢谢! 如有问题请指出.)

---

参考文献

[1] 张池平, 施云慧, 计算方法, 科学出版社, 2002

[2] 王绵森, 马知恩, 工科数学分析基础, 高等教育出版社, 2006

[3] 百度百科, “连续可微”, https://baike.baidu.com/item/%E8%BF%9E%E7%BB%AD%E5%8F%AF%E5%BE%AE/3420868

[4] J. E. Dennis, Jr., and R. B. Schabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, SIAM, 1996

[5] L. Vandenberghe, “16. Gauss–Newton method”, http://www.seas.ucla.edu/~vandenbe/236C/lectures/gn.pdf, ECE236C - Optimization Methods for Large-Scale Systems, UCLA, 2022

[6] P. Deuflhard, “A Short History of Newton’s Method”, Mathematics, 2012

[7] 王文平, 朱春浩, 最小二乘法: 勒让德与高斯, 武汉船舶职业技术学院学报, 2011 年第 6 期

[8] Jan R. Magnus, Gauss on least-squares and maximum-likelihood estimation, Archive for History of Exact Sciences, 2022

[9] Britannica, “least squares method”, https://www.britannica.com/topic/least-squares-approximation

---