牛顿法（Newton’s method）

假如已知函数 f(x)f(x)f(x)，想要求 f(x)=0f(x)=0f(x)=0 的解（或者叫根）。

牛顿法（Newton’s method）大致的思想是：
（1）选一个初始位置 x0x_0x0（这个位置最好是在根的附近）；
（2）在这个位置上找一个 f(x)f(x)f(x) 的近似函数（通常用泰勒展开）；
（3）令近似函数为 000 ，求解；
（4）以这个解为新的位置 x1x_1x1；
（5）重复上述迭代，到第 nnn 次迭代得到 xnx_nxn，当 ∣xn−xn−1∣|x_n - x_{n-1}|∣xn−xn−1∣ 足够小，结束。
xnx_nxn 就是 f(x)=0f(x)=0f(x)=0 的近似解。

下面的解释1是一个比较直观的过程，不过问题的设计不太合理。解释2会好一点。

例子： 求 f(x)f(x)f(x) 的极小值。

解释1

假如 f(x)f(x)f(x) 长这样：

首先，选一个初始值估计值 x0x_0x0，由于函数已知，带入得到对应的函数值为 f(x0)f(x_0)f(x0)。
由于函数在该点处可导，可以求得该点的导数为 f′(x0)f'(x_0)f′(x0)。
该点的导数值，就是该点的切线的斜率。

已知一个点和斜率，可以求得直线方程：y−f(x0)x−x0=f′(x0)⇓y=f′(x0)(x−x0)+f(x0)(1)\dfrac{y-f(x_0)}{x-x_0} = f'(x_0) \\ \Downarrow \\ y = f'(x_0)(x-x_0) + f(x_0) \tag{1} x−x0y−f(x0)=f′(x0)⇓y=f′(x0)(x−x0)+f(x0)(1)
公式(1)即为图中红色直线的方程，其中 x0x_0x0，f(x0)f(x_0)f(x0)，f′(x0)f'(x_0)f′(x0) 都是已知量。

现在，求红色直线与 xxx 轴的交点：x1x_1x1 的值。
往公式(1)带入 y=0y=0y=0 即可。

求得x1=x0−f(x0)f′(x0)x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}x1=x0−f′(x0)f(x0)

把 x1x_1x1 代回原方程，得到对应的函数值 f(x1)f(x_1)f(x1)。在该点处继续求切线方程，重复上述步骤。

求直线与 xxx 轴的交点，得到 x2x_2x2：

从图中可以看到，f(x2)f(x_2)f(x2) 已经几乎是极小值了。

实际操作中会判断两次迭代的差值 f(xn)−f(xn−1)f(x_n) - f(x_{n-1})f(xn)−f(xn−1)，如果足够小，就可以结束迭代了。
（这貌似有点像梯度下降）（对于这个问题直接求导数，令导数等于 000 也可以。）

我这个问题得设计好像有点 bug 。继续迭代会出问题，f(x2)f(x_2)f(x2) 处可能无法求导，x3x_3x3 也不知道去到哪儿了。如果导数接近于 000 也可以结束迭代吧。
总之牛顿法的过程大概是这样。具体可以看百度百科或者其它资料。

贴一个维基百科的图：

像这个例子，是对这种 f(x)+a=0f(x)+a=0f(x)+a=0 的问题用牛顿法。

解释2

来自知乎：https://zhuanlan.zhihu.com/p/103724149

一维情况

迭代公式：

xn+1=xn−f′(xn)f′′(xn)x_{n+1}=x_n - \dfrac{f'(x_n)}{f''(x_n)}xn+1=xn−f′′(xn)f′(xn)

牛顿法的推导基于二阶可微函数的泰勒展开，以一维函数为例，在 x0x_0x0 处 f(x)f(x)f(x) 的泰勒展开：
f(x)≈f(x0)+f′(x0)(x−x0)+12f′′(x0)(x−x0)2\textcolor{#417991}{ f(x)} \approx \textcolor{#D97600}{f(x_0) + f'(x_0)(x-x_0) + \dfrac{1}{2} f''(x_0)(x-x_0)^2} f(x)≈f(x0)+f′(x0)(x−x0)+21f′′(x0)(x−x0)2

即在 x0x_0x0 附近可以用 f(x0)+f′(x0)(x−x0)+12f′′(x0)(x−x0)2(2)\textcolor{#D97600}{f(x_0) + f'(x_0)(x-x_0) + \dfrac{1}{2} f''(x_0)(x-x_0)^2} \tag{2}f(x0)+f′(x0)(x−x0)+21f′′(x0)(x−x0)2(2) 近似替代 f(x)\textcolor{#417991}{f(x)}f(x)。

式子(2)是对 f(x)\textcolor{#417991}{f(x)}f(x) 在x0x_0x0 处的二阶近似，如下图橙色曲线：

蓝色曲线代表原函数 f(x)\textcolor{#417991}{f(x)}f(x)。绿色点代表当前点 x0\textcolor{#007500}{x_0}x0。

对橙色曲线求导，求倒数的零点，得到下一次迭代的位置 x1x_1x1。
也就是对式子(2)求导，得到：f′(x0)+f′′(x0)(x−x0)(3)f'(x_0) + f''(x_0)(x-x_0) \tag{3} f′(x0)+f′′(x0)(x−x0)(3)
然后令式子(3)等于 000，求得x1=x0−f′(x0)f′′(x0)(4)x_1=x_0 - \dfrac{f'(x_0)}{f''(x_0)} \tag{4}x1=x0−f′′(x0)f′(x0)(4)

式子(4)即为一维函数的牛顿法迭代公式。

高维情况

迭代公式：

xn+1=xn−H(xn)−1∇f(xn)\boldsymbol{x}_{n+1} = \boldsymbol{x}_n - H(\boldsymbol{x}_n)^{-1} \nabla f(\boldsymbol{x}_n)xn+1=xn−H(xn)−1∇f(xn)

对于高维函数，推导过程基于多元函数的泰勒展开：

f(x)≈f(x0)+∇f(x0)T⋅(x−x0)+12(x−x0)T⋅H(x0)⋅(x−x0)(5)f(\boldsymbol{x}) \approx f(\boldsymbol{x}_0) + \nabla f(\boldsymbol{x}_0)^T \cdot (\boldsymbol{x} - \boldsymbol{x}_0) + \frac{1}{2} (\boldsymbol{x} - \boldsymbol{x}_0)^T \cdot H(\boldsymbol{x}_0) \cdot (\boldsymbol{x} - \boldsymbol{x}_0) \tag{5} f(x)≈f(x0)+∇f(x0)T⋅(x−x0)+21(x−x0)T⋅H(x0)⋅(x−x0)(5)

上面公式用高维二次曲面在 x0\boldsymbol{x}_0x0 处逼近原函数。
用下面的图片表示类似的意思。（ x0\boldsymbol{x}_0x0 就是当前点，相当于图里面的 x(k)\textbf{x}^{(\text{k})}x(k) ）

接下来，和一维的情况一样，令式子(5)右边的那部分对 x\boldsymbol{x}x 求导，令导数等于 000：
H(x0)(x−x0)+∇f(x0)=0(6)H(\boldsymbol{x}_0) (\boldsymbol{x} - \boldsymbol{x}_0) + \nabla f(\boldsymbol{x}_0) = 0 \tag{6} H(x0)(x−x0)+∇f(x0)=0(6)

解得 x1=x0−H(x0)−1∇f(x0)(7)\boldsymbol{x}_{1} = \boldsymbol{x}_0 - H(\boldsymbol{x}_0)^{-1} \nabla f(\boldsymbol{x}_0) \tag{7} x1=x0−H(x0)−1∇f(x0)(7)

其中 HHH 是 Hessian Matrix （海森矩阵），其实就是个二阶导的矩阵。
对于 nnn 元变量：
H(f)=[∂2f∂x12∂2f∂x1∂x2⋯∂2f∂x1∂xn∂2f∂x2∂x1∂2f∂x22⋯∂2f∂x2∂xn⋮⋮⋱⋮∂2f∂xn∂x1∂2f∂xn∂x2⋯∂2f∂xn2]H(f) = \begin{bmatrix} \dfrac{\partial^2 f}{\partial x_1^2} & \dfrac{\partial^2 f}{\partial x_1 \partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_1 \partial x_n} \\[1.5em] \dfrac{\partial^2 f}{\partial x_2 \partial x_1} & \dfrac{\partial^2 f}{\partial x_2^2} & \cdots & \dfrac{\partial^2 f}{\partial x_2 \partial x_n} \\[1.5em] \vdots & \vdots & \ddots & \vdots \\[1.5em] \dfrac{\partial^2 f}{\partial x_n \partial x_1} & \dfrac{\partial^2 f}{\partial x_n \partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_n^2} \end{bmatrix} \\[1.5em] H(f)=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡∂x12∂2f∂x2∂x1∂2f⋮∂xn∂x1∂2f∂x1∂x2∂2f∂x22∂2f⋮∂xn∂x2∂2f⋯⋯⋱⋯∂x1∂xn∂2f∂x2∂xn∂2f⋮∂xn2∂2f⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤

牛顿法的优点是收敛速度快，缺点是需要求矩阵的逆，计算量比较大。
此外，如果矩阵非正定（在一维情况下表现为泰勒展开的二阶导数小于 000 ），极值点为极大值，而非极小值。
如果初始位置离最优点太远，也会导致迭代过程中目标函数不严格递减的情况。

解决这个问题的方法之一是拟牛顿法（Quasi-Newton Methods）。

高斯-牛顿法（Gauss–Newton algorithm）

介绍：
高斯-牛顿法（Gauss–Newton algorithm）是牛顿法的特例，它是牛顿法的修改版，用于寻找函数的最小值。
和牛顿法不一样，它只能用于解决最小二乘问题。
优点是，不需要二阶导数（二阶导数可能很难计算）。

例子：
设我们有 mmm 组样本（数据集）：{(x(1),y(1)),(x(2),y(2)),…,(x(m),y(m))}，x∈Rn\{ (\boldsymbol{x}^{(1)}, y^{(1)}), (\boldsymbol{x}^{(2)}, y^{(2)}), \dots, (\boldsymbol{x}^{(m)}, y^{(m)}) \} ，\; \; \boldsymbol{x}\in \mathbb{R}^n{(x(1),y(1)),(x(2),y(2)),…,(x(m),y(m))}，x∈Rn 上标括号 (⋅)^{(\cdot)}(⋅) 表示第几组样本（学吴恩达的表示法）。

我们希望找到包含 nnn 个参数的非线性函数 f(x,θ)f(\boldsymbol{x}, \boldsymbol{\theta})f(x,θ)，拟合上面 mmm 组数据。
其中 θ\boldsymbol{\theta}θ 是我们要找的参数，θ∈Rn\boldsymbol{\theta} \in \mathbb{R}^nθ∈Rn。

为了方便描述，设代入第 iii 个样本后，函数值为 f(i)(θ)f^{(i)}(\boldsymbol{\theta})f(i)(θ) 。也就是 f(i)(θ)=f(x(i),θ)f^{(i)}(\boldsymbol{\theta}) = f( \boldsymbol{x}^{(i)}, \boldsymbol{\theta} )f(i)(θ)=f(x(i),θ)。

则最小二乘目标函数（或者叫平方误差函数？）为：

ϵ(θ)=∑i=1m∥f(i)(θ)−y(i)∥2(8)\epsilon (\boldsymbol{\theta}) = \sum_{i=1}^m \| f^{(i)}(\boldsymbol{\theta}) - y^{(i)} \| ^2 \tag{8}ϵ(θ)=i=1∑m∥f(i)(θ)−y(i)∥2(8)

（一般用 JJJ 表示代价函数，但是这里用了 ϵ\epsilonϵ，原因是稍后要用 JJJ 表示雅可比矩阵）

我们需要找到一组参数 θ=[θ1,θ2,…,θn]T\boldsymbol{\theta} = [\theta_1, \theta_2, \dots, \theta_n]^Tθ=[θ1,θ2,…,θn]T，使式子(8)最小。也就是：arg⁡min⁡θϵ(θ)\arg\min_{\boldsymbol{\theta}} \epsilon(\boldsymbol{\theta})argθminϵ(θ)

用式子(8)对第 jjj 个参数 θj\theta_jθj 求导：
∂ϵ(θ)θj=∑i=1m2⋅(f(i)(θ)−y(i))⋅∂f(i)(θ)∂θj(9)\frac{\partial \epsilon(\boldsymbol{\theta}) }{ \theta_j } = \sum_{i=1}^m 2 \cdot \left( f^{(i)}(\boldsymbol{\theta}) - y^{(i)} \right) \cdot \frac{ \partial f^{(i)}(\boldsymbol{\theta}) }{ \partial \theta_j } \tag{9} θj∂ϵ(θ)=i=1∑m2⋅(f(i)(θ)−y(i))⋅∂θj∂f(i)(θ)(9)

令雅可比矩阵（jacobi matrix）JJJ 为：
J=[∂f(1)(θ)∂θ1∂f(1)(θ)∂θ2⋯∂f(1)(θ)∂θn∂f(2)(θ)∂θ1∂f(2)(θ)∂θ2⋯∂f(2)(θ)∂θn⋮⋮⋱⋮∂f(m)(θ)∂θ1∂f(m)(θ)∂θ2⋯∂f(m)(θ)∂θn]J= \begin{bmatrix} \dfrac{\partial f^{(1)}(\boldsymbol{\theta})}{\partial \theta_1} & \dfrac{\partial f^{(1)}(\boldsymbol{\theta}) }{\partial \theta_2} & \cdots & \dfrac{\partial f^{(1)}(\boldsymbol{\theta}) }{\partial \theta_n} \\[1.5em] \dfrac{\partial f^{(2)}(\boldsymbol{\theta})}{\partial \theta_1} & \dfrac{\partial f^{(2)}(\boldsymbol{\theta}) }{\partial \theta_2} & \cdots & \dfrac{\partial f^{(2)}(\boldsymbol{\theta}) }{\partial \theta_n} \\[1.5em] \vdots & \vdots & \ddots & \vdots \\[1.5em] \dfrac{\partial f^{(m)}(\boldsymbol{\theta})}{\partial \theta_1} & \dfrac{\partial f^{(m)}(\boldsymbol{\theta}) }{\partial \theta_2} & \cdots & \dfrac{\partial f^{(m)}(\boldsymbol{\theta}) }{\partial \theta_n} \end{bmatrix}J=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡∂θ1∂f(1)(θ)∂θ1∂f(2)(θ)⋮∂θ1∂f(m)(θ)∂θ2∂f(1)(θ)∂θ2∂f(2)(θ)⋮∂θ2∂f(m)(θ)⋯⋯⋱⋯∂θn∂f(1)(θ)∂θn∂f(2)(θ)⋮∂θn∂f(m)(θ)⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤

所以 JJJ 是一个 mmm 行 nnn 列的矩阵，m×nm\times nm×n。
第 111 行的含义为，第 111 个样本对于每个参数的偏导。

令残差（residual，样本与拟合值之间的差）r\textbf{r}r 为：
r=[f(1)(θ)−y1f(2)(θ)−y2⋮f(m)(θ)−ym]\textbf{r} = \begin{bmatrix} f^{(1)}(\boldsymbol{\theta}) -y_1 \\[0.5em] f^{(2)}(\boldsymbol{\theta}) -y_2 \\[0.5em] \vdots \\[0.5em] f^{(m)}(\boldsymbol{\theta}) -y_m \end{bmatrix} r=⎣⎢⎢⎢⎢⎢⎢⎡f(1)(θ)−y1f(2)(θ)−y2⋮f(m)(θ)−ym⎦⎥⎥⎥⎥⎥⎥⎤

则式子(9)可以写成矩阵形式：
∇ϵ(θ)=2JTr(10)\nabla \epsilon(\boldsymbol{\theta}) = 2 J^T \textbf{r} \tag{10} ∇ϵ(θ)=2JTr(10)

接下来求海森矩阵第 kkk 行 jjj 列的元素：
∂2ϵ(θ)∂θk∂θj=∂∂θk∂ϵ(θ)∂θj=∂∂θk(∑i=1m2⋅(f(i)(θ)−y(i))⋅∂f(i)(θ)∂θj)=2∑i=1m(f(i)(θ)∂θkf(i)(θ)∂θj+(f(i)(θ)−y(i))⋅∂2f(i)(θ)∂θk∂θj)(11)\begin{aligned} \frac{\partial ^2\epsilon(\boldsymbol{\theta})} { \partial \theta_k \partial \theta_j } &= \frac{\partial}{\partial \theta_k} \frac{\partial \epsilon(\boldsymbol{\theta})}{ \partial \theta_j} \\[1.5em] &=\frac{\partial}{\partial \theta_k} \left( \sum_{i=1}^m 2 \cdot \left( f^{(i)}(\boldsymbol{\theta}) - y^{(i)} \right) \cdot \frac{ \partial f^{(i)}(\boldsymbol{\theta}) }{ \partial \theta_j } \right)\\[1.5em] &= 2 \sum_{i=1}^m \left( \frac{ f^{(i)}(\boldsymbol{\theta}) }{ \partial \theta_k } \frac{ f^{(i)}(\boldsymbol{\theta}) }{ \partial \theta_j } + \left( f^{(i)}(\boldsymbol{\theta}) - y^{(i)} \right) \cdot \frac{\partial ^2 f^{(i)} (\boldsymbol{\theta})} { \partial \theta_k \partial \theta_j } \right) \\ \tag{11} \end{aligned} ∂θk∂θj∂2ϵ(θ)=∂θk∂∂θj∂ϵ(θ)=∂θk∂(i=1∑m2⋅(f(i)(θ)−y(i))⋅∂θj∂f(i)(θ))=2i=1∑m(∂θkf(i)(θ)∂θjf(i)(θ)+(f(i)(θ)−y(i))⋅∂θk∂θj∂2f(i)(θ))(11)

令 HHH 为 ϵ(θ)\epsilon(\boldsymbol{\theta})ϵ(θ) 的海森矩阵，由(11)可知：H=2⋅(JTJ+S)H= 2 \cdot (J^T J + S) H=2⋅(JTJ+S) 其中：

Sk,j=∑i=1m(f(i)(θ)−y(i))⋅∂2f(i)(θ)∂θk∂θjS_{k,j} = \sum_{i=1}^m \left( f^{(i)}(\boldsymbol{\theta}) - y^{(i)} \right) \cdot \frac{\partial ^2 f^{(i)} (\boldsymbol{\theta})} { \partial \theta_k \partial \theta_j } Sk,j=i=1∑m(f(i)(θ)−y(i))⋅∂θk∂θj∂2f(i)(θ)

把 (10) 和 (12) 带入 (7) 得：
θn+1=θn−H−1∇ϵ(θ)=θn−(JTJ+S)−1⋅JTr(13)\begin{aligned} \boldsymbol{\theta}_{n+1} & = \boldsymbol{\theta}_{n} - H^{-1} \nabla \epsilon(\boldsymbol{\theta}) \\ & = \boldsymbol{\theta}_n - \left( J^T J + S \right) ^{-1} \cdot \, J^T \textbf{r} \end{aligned} \tag{13} θn+1=θn−H−1∇ϵ(θ)=θn−(JTJ+S)−1⋅JTr(13)

很多时候 (13) 中的 SSS 可以忽略，最终高斯-牛顿法的迭代公式为：
θn+1=θn−(JTJ)−1⋅JTr(14)\boldsymbol{\theta}_{n+1} = \boldsymbol{\theta}_n - \left( J^T J \right) ^{-1} \cdot \, J^T \textbf{r} \tag{14}θn+1=θn−(JTJ)−1⋅JTr(14)

牛顿法，高斯-牛顿法相关推荐

最优化方法总结——梯度下降法、最速下降法、牛顿法、高斯牛顿法、LM法、拟牛顿法
目录 1 最优化方法的结构 2 常用最优化方法对比分析 3 相关计算公式 1 最优化方法的结构最优化问题的一般形式为: 其中为决策变量,是目标函数,为约束集或可行域.特别地,如果,则最优化问题成为无 ...
SLAM从0到1——状态估计之最小二乘问题解法：最速下降法、牛顿法、高斯牛顿法、LM法...
点击上方"3D视觉工坊",选择"星标" 干货第一时间送达学习3D视觉核心技术,扫描查看介绍,3天内无条件退款圈里有高质量教程资料.可答疑解惑.助你高效解决问 ...
高斯牛顿法在具体工程中的应用——C++版
高斯牛顿法在具体工程中的应用--C++版说明:本文章没有用大量篇幅来讲述高斯牛顿的原理和数学中的应用,而是用具体的代码来说明,具体是怎么应用的. 如果对高斯牛顿法的原理比较感兴趣,可以阅读以下链接中 ...
牛顿法、梯度下降法、高斯牛顿法、Levenberg-Marquardt算法
何为梯度? 一般解释: f(x)在x0的梯度:就是f(x)变化最快的方向举个例子,f()是一座山,站在半山腰, 往x方向走1米,高度上升0.4米,也就是说x方向上的偏导是 0.4 往y方向走1米,高 ...
漫步最优化三十四——高斯-牛顿法
你的温柔像羽毛,\textbf{你的温柔像羽毛,} 秘密躺在我怀抱.\textbf{秘密躺在我怀抱.} 你的微笑像拥抱,\textbf{你的微笑像拥抱,} 只有我能看到.\textbf{只有我能看到. ...
高斯牛顿法 matlab代码实现
高斯牛顿法的解释请参见 http://blog.csdn.net/dsbatigol/article/details/12448627 %% Gauss-Newton method clc; clea ...
牛顿法和高斯牛顿法对比
文章目录一.非线性最小二乘一.牛顿法二.高斯牛顿法三.列文伯格-马夸尔特法(LM) 四.ceres求解优化问题一.非线性最小二乘考虑最小二乘函数F(x), 其等于: 通过求F(x)导数为零 ...
高斯牛顿法去畸变（C++实现）
1. 背景在实际三维重建时,通常会遇到如下问题: 已知有畸变的点坐标x′,y′x',y'x′,y′,已知畸变模型f(x,y),g(x,y)f(x,y),g(x,y)f(x,y),g(x,y), 求无 ...
最优化八：高斯牛顿法、LM法
梯度法:,负梯度方向牛顿法:,A为Hession矩阵高斯牛顿法:,为的解 LM法:,为的解 1 高斯牛顿法(Gauss-Newton) 针对优化问题求解x使得f(x)取得最小值,采用高斯牛顿法,步 ...
LIO-SAM：在高斯牛顿法求解过程中用SO3代替欧拉角
LIO-SAM发表于IROS2020,是一个效果非常好的惯性-激光紧耦合里程计我打算给我们的机器人搞一个激光里程计,于是打算把LIO-SAM改一改搞过来修改过程中发现一个问题:在里程计求解(map ...

牛顿法，高斯-牛顿法

牛顿法（Newton’s method）

解释1

解释2

一维情况

高维情况

高斯-牛顿法（Gauss–Newton algorithm）

牛顿法，高斯-牛顿法相关推荐

最新文章

热门文章