牛顿法，高斯-牛顿法

这是一篇笔记，不是教程。 仅供参考，对与错自行甄别。

---

牛顿法（Newton’s method）

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

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

---

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

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

解释1

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

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

已知一个点和斜率，可以求得直线方程：$$\begin{gathered} \frac{y-f(x_0)}{x-x_0}=f^{\prime }(x_0) \\ \Downarrow \\ y=f^{\prime }(x_0)(x-x_0)+f(x_0)\text{\hspace{1em}(1)} \end{gathered}$$
公式(1)即为图中红色直线的方程，其中 $x_0$，$f(x_0)$，$f^{\prime }(x_0)$ 都是已知量。

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

求得$$x_1=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$$

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

求直线与 $x$ 轴的交点，得到 $x_2$：

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

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

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

贴一个维基百科的图：

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

解释2

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

一维情况

迭代公式：

$$x_{n+1}=x_n-\frac{f^{\prime }(x_n)}{f^{\prime \prime }(x_n)}$$

牛顿法的推导基于二阶可微函数的泰勒展开，以一维函数为例，在 $x_0$ 处 $f(x)$ 的泰勒展开：
$$f(x)\approx f(x_0)+f^{\prime }(x_0)(x-x_0)+\frac{1}{2}{f}^{\prime \prime }(x_0)(x-x_0{)}^2$$

即在 $x_0$ 附近可以用 $$f(x_0)+f^{\prime }(x_0)(x-x_0)+\frac{1}{2}{f}^{\prime \prime }(x_0)(x-x_0{)}^2\text{\hspace{1em}(2)}$$ 近似替代 $f(x)$。

式子(2)是对 $f(x)$ 在$x_0$ 处的二阶近似，如下图橙色曲线：

蓝色曲线代表原函数 $f(x)$。绿色点代表当前点 $x_0$。

对橙色曲线求导，求倒数的零点，得到下一次迭代的位置 $x_1$。
也就是对式子(2)求导，得到：$$f^{\prime }(x_0)+f^{\prime \prime }(x_0)(x-x_0)\text{\hspace{1em}(3)}$$
然后令式子(3)等于 $0$，求得$$x_1=x_0-\frac{f^{\prime }(x_0)}{f^{\prime \prime }(x_0)}\text{\hspace{1em}(4)}$$

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

高维情况

迭代公式：

$${\mathbf{x}}_{n+1}={\mathbf{x}}_n-H({\mathbf{x}}_n{)}^{-1}\nabla f({\mathbf{x}}_n)$$

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

$$f(\mathbf{x})\approx f({\mathbf{x}}_0)+\nabla f({\mathbf{x}}_0{)}^T\cdot (\mathbf{x}-{\mathbf{x}}_0)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_0{)}^T\cdot H({\mathbf{x}}_0)\cdot (\mathbf{x}-{\mathbf{x}}_0)\text{\hspace{1em}(5)}$$

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

接下来，和一维的情况一样，令式子(5)右边的那部分对 $\mathbf{x}$ 求导，令导数等于 $0$：
$$H({\mathbf{x}}_0)(\mathbf{x}-{\mathbf{x}}_0)+\nabla f({\mathbf{x}}_0)=0\text{\hspace{1em}(6)}$$

解得 $${\mathbf{x}}_1={\mathbf{x}}_0-H({\mathbf{x}}_0{)}^{-1}\nabla f({\mathbf{x}}_0)\text{\hspace{1em}(7)}$$

其中 $H$ 是 Hessian Matrix （海森矩阵），其实就是个二阶导的矩阵。
对于 $n$ 元变量：
$$H(f)=\left[\begin{array}{cccc}\frac{{\partial }^{2}f}{\partial x_1^2}& \frac{{\partial }^{2}f}{\partial x_1\partial x_2}& \cdots & \frac{{\partial }^{2}f}{\partial x_1\partial x_n}\\ \frac{{\partial }^{2}f}{\partial x_2\partial x_1}& \frac{{\partial }^{2}f}{\partial x_2^2}& \cdots & \frac{{\partial }^{2}f}{\partial x_2\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}f}{\partial x_n\partial x_1}& \frac{{\partial }^{2}f}{\partial x_n\partial x_2}& \cdots & \frac{{\partial }^{2}f}{\partial x_n^2}\end{array}\right]$$

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

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

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

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

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

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

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

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

$$ϵ(\mathbf{\theta })=\sum _{i=1}^m\Vert f^{(i)}(\mathbf{\theta })-y^{(i)}{\Vert }^2\text{\hspace{1em}(8)}$$

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

我们需要找到一组参数 $\mathbf{\theta }=[{\theta }_1,{\theta }_2,\dots ,{\theta }_n{]}^T$，使式子(8)最小。也就是：$$\arg \underset{\mathbf{\theta }}{\min }ϵ(\mathbf{\theta })$$

用式子(8)对第 $j$ 个参数 ${\theta }_j$ 求导：
$$\frac{\partial ϵ(\mathbf{\theta })}{{\theta }_j}=\sum _{i=1}^{m}2\cdot \left(f^{(i)}(\mathbf{\theta })-y^{(i)}\right)\cdot \frac{\partial f^{(i)}(\mathbf{\theta })}{\partial {\theta }_j}\text{\hspace{1em}(9)}$$

令雅可比矩阵（jacobi matrix）$J$ 为：
$$J=\left[\begin{array}{cccc}\frac{\partial f^{(1)}(\mathbf{\theta })}{\partial {\theta }_1}& \frac{\partial f^{(1)}(\mathbf{\theta })}{\partial {\theta }_2}& \cdots & \frac{\partial f^{(1)}(\mathbf{\theta })}{\partial {\theta }_n}\\ \frac{\partial f^{(2)}(\mathbf{\theta })}{\partial {\theta }_1}& \frac{\partial f^{(2)}(\mathbf{\theta })}{\partial {\theta }_2}& \cdots & \frac{\partial f^{(2)}(\mathbf{\theta })}{\partial {\theta }_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial f^{(m)}(\mathbf{\theta })}{\partial {\theta }_1}& \frac{\partial f^{(m)}(\mathbf{\theta })}{\partial {\theta }_2}& \cdots & \frac{\partial f^{(m)}(\mathbf{\theta })}{\partial {\theta }_n}\end{array}\right]$$

所以 $J$ 是一个 $m$ 行 $n$ 列的矩阵，$m\times n$。
第 $1$ 行的含义为，第 $1$ 个样本对于每个参数的偏导。

令残差（residual，样本与拟合值之间的差）$\text{r}$ 为：
$$\text{r}=\left[\begin{array}{c}{f}^{(1)}(\mathbf{\theta })-y_1\\ f^{(2)}(\mathbf{\theta })-y_2\\ ⋮\\ f^{(m)}(\mathbf{\theta })-y_m\end{array}\right]$$

则式子(9)可以写成矩阵形式：
$$\nabla ϵ(\mathbf{\theta })=2J^T\text{r}\text{\hspace{1em}(10)}$$

接下来求海森矩阵第 $k$ 行 $j$ 列的元素：
$$\begin{array}{rl}\frac{{\partial }^2ϵ(\mathbf{\theta })}{\partial {\theta }_k\partial {\theta }_j}& =\frac{\partial }{\partial {\theta }_k}\frac{\partial ϵ(\mathbf{\theta })}{\partial {\theta }_j}\\ & =\frac{\partial }{\partial {\theta }_k}\left(\sum _{i=1}^{m}2\cdot \left(f^{(i)}(\mathbf{\theta })-y^{(i)}\right)\cdot \frac{\partial f^{(i)}(\mathbf{\theta })}{\partial {\theta }_j}\right)\\ & =2\sum _{i=1}^m\left(\frac{f^{(i)}(\mathbf{\theta })}{\partial {\theta }_k}\frac{f^{(i)}(\mathbf{\theta })}{\partial {\theta }_j}+\left(f^{(i)}(\mathbf{\theta })-y^{(i)}\right)\cdot \frac{{\partial }^2{f}^{(i)}(\mathbf{\theta })}{\partial {\theta }_k\partial {\theta }_j}\right)\end{array}\text{\hspace{1em}(11)}$$

令 $H$ 为 $ϵ(\mathbf{\theta })$ 的海森矩阵，由(11)可知：$$H=2\cdot (J^{T}J+S)$$ 其中：

$$S_{k,j}=\sum _{i=1}^m\left(f^{(i)}(\mathbf{\theta })-y^{(i)}\right)\cdot \frac{{\partial }^2{f}^{(i)}(\mathbf{\theta })}{\partial {\theta }_k\partial {\theta }_j}$$

把 (10) 和 (12) 带入 (7) 得：
$$\begin{array}{rl}{\mathbf{\theta }}_{n+1}& ={\mathbf{\theta }}_n-H^{-1}\nabla ϵ(\mathbf{\theta })\\ & ={\mathbf{\theta }}_n-{\left(J^{T}J+S\right)}^{-1}\cdot J^T\text{r}\end{array}\text{\hspace{1em}(13)}$$

很多时候 (13) 中的 $S$ 可以忽略，最终高斯-牛顿法的迭代公式为：
$${\mathbf{\theta }}_{n+1}={\mathbf{\theta }}_n-{\left(J^{T}J\right)}^{-1}\cdot J^T\text{r}\text{\hspace{1em}(14)}$$