牛顿法（Newton’s method）
假如已知函数 f ( x ) f(x) f(x)，想要求 f ( x ) = 0 f(x)=0 f(x)=0 的解（或者叫根）。
牛顿法（Newton’s method）大致的思想是：
（1）选一个初始位置
x
0
x_0
x0?（这个位置最好是在根的附近）；
（2）在这个位置上找一个
f
(
x
)
f(x)
f(x) 的近似函数（通常用泰勒展开）；
（3）令近似函数为
0
0
0 ，求解；
（4）以这个解为新的位置
x
1
x_1
x1?；
（5）重复上述迭代，到第
n
n
n 次迭代得到
x
n
x_n
xn?，当
∣
x
n
?
x
n
?
1
∣
|x_n - x_{n-1}|
∣xn??xn?1?∣ 足够小，结束。
x
n
x_n
xn? 就是
f
(
x
)
=
0
f(x)=0
f(x)=0 的近似解。
下面的解释1是一个比较直观的过程，不过问题的设计不太合理。解释2会好一点。
例子： 求 f ( x ) f(x) f(x) 的极小值。
解释1
假如 f ( x ) f(x) f(x) 长这样：
首先，选一个初始值估计值
x
0
x_0
x0?，由于函数已知，带入得到对应的函数值为
f
(
x
0
)
f(x_0)
f(x0?)。
由于函数在该点处可导，可以求得该点的导数为
f
′
(
x
0
)
f'(x_0)
f′(x0?)。
该点的导数值，就是该点的切线的斜率。
已知一个点和斜率，可以求得直线方程：
y
?
f
(
x
0
)
x
?
x
0
=
f
′
(
x
0
)
?
y
=
f
′
(
x
0
)
(
x
?
x
0
)
+
f
(
x
0
)
(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?x0?y?f(x0?)?=f′(x0?)?y=f′(x0?)(x?x0?)+f(x0?)(1)
公式(1)即为图中红色直线的方程，其中
x
0
x_0
x0?，
f
(
x
0
)
f(x_0)
f(x0?)，
f
′
(
x
0
)
f'(x_0)
f′(x0?) 都是已知量。
现在，求红色直线与
x
x
x 轴的交点：
x
1
x_1
x1? 的值。
往公式(1)带入
y
=
0
y=0
y=0 即可。
求得 x 1 = x 0 ? f ( x 0 ) f ′ ( x 0 ) x_1 = x_0 - \frac{f(x_0)}{f'(x_0)} x1?=x0??f′(x0?)f(x0?)?
把 x 1 x_1 x1? 代回原方程，得到对应的函数值 f ( x 1 ) f(x_1) f(x1?)。在该点处继续求切线方程，重复上述步骤。
求直线与 x x x 轴的交点，得到 x 2 x_2 x2?：
从图中可以看到， f ( x 2 ) f(x_2) f(x2?) 已经几乎是极小值了。
实际操作中会判断两次迭代的差值
f
(
x
n
)
?
f
(
x
n
?
1
)
f(x_n) - f(x_{n-1})
f(xn?)?f(xn?1?)，如果足够小，就可以结束迭代了。
（这貌似有点像梯度下降）（对于这个问题直接求导数，令导数等于
0
0
0 也可以。）
我这个问题得设计好像有点 bug 。继续迭代会出问题， f ( x 2 ) f(x_2) f(x2?) 处可能无法求导， x 3 x_3 x3? 也不知道去到哪儿了。如果导数接近于 0 0 0 也可以结束迭代吧。
总之牛顿法的过程大概是这样。具体可以看百度百科或者其它资料。
贴一个维基百科的图：
像这个例子，是对这种 f ( x ) + a = 0 f(x)+a=0 f(x)+a=0 的问题用牛顿法。
解释2
来自知乎：https://zhuanlan.zhihu.com/p/103724149
一维情况
迭代公式：
x n + 1 = x n ? f ′ ( x n ) f ′ ′ ( x n ) x_{n+1}=x_n - \dfrac{f'(x_n)}{f''(x_n)} xn+1?=xn??f′′(xn?)f′(xn?)?
牛顿法的推导基于二阶可微函数的泰勒展开，以一维函数为例，在
x
0
x_0
x0? 处
f
(
x
)
f(x)
f(x) 的泰勒展开：
f
(
x
)
≈
f
(
x
0
)
+
f
′
(
x
0
)
(
x
?
x
0
)
+
1
2
f
′
′
(
x
0
)
(
x
?
x
0
)
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?)+21?f′′(x0?)(x?x0?)2
即在 x 0 x_0 x0? 附近可以用 f ( x 0 ) + f ′ ( x 0 ) ( x ? x 0 ) + 1 2 f ′ ′ ( x 0 ) ( x ? x 0 ) 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?)+21?f′′(x0?)(x?x0?)2(2) 近似替代 f ( x ) \textcolor{#417991}{f(x)} f(x)。
式子(2)是对
f
(
x
)
\textcolor{#417991}{f(x)}
f(x) 在
x
0
x_0
x0? 处的二阶近似，如下图橙色曲线：
蓝色曲线代表原函数 f ( x ) \textcolor{#417991}{f(x)} f(x)。绿色点代表当前点 x 0 \textcolor{#007500}{x_0} x0?。
对橙色曲线求导，求倒数的零点，得到下一次迭代的位置
x
1
x_1
x1?。
也就是对式子(2)求导，得到：
f
′
(
x
0
)
+
f
′
′
(
x
0
)
(
x
?
x
0
)
(3)
f'(x_0) + f''(x_0)(x-x_0) \tag{3}
f′(x0?)+f′′(x0?)(x?x0?)(3)
然后令式子(3)等于
0
0
0，求得
x
1
=
x
0
?
f
′
(
x
0
)
f
′
′
(
x
0
)
(4)
x_1=x_0 - \dfrac{f'(x_0)}{f''(x_0)} \tag{4}
x1?=x0??f′′(x0?)f′(x0?)?(4)
式子(4)即为一维函数的牛顿法迭代公式。
高维情况
迭代公式：
x n + 1 = x n ? H ( x n ) ? 1 ? f ( x n ) \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 ( x 0 ) + ? f ( x 0 ) T ? ( x ? x 0 ) + 1 2 ( x ? x 0 ) T ? H ( x 0 ) ? ( x ? x 0 ) (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)
上面公式用高维二次曲面在
x
0
\boldsymbol{x}_0
x0? 处逼近原函数。
用下面的图片表示类似的意思。（
x
0
\boldsymbol{x}_0
x0? 就是当前点，相当于图里面的
x
(
k
)
\textbf{x}^{(\text{k})}
x(k) ）
接下来，和一维的情况一样，令式子(5)右边的那部分对
x
\boldsymbol{x}
x 求导，令导数等于
0
0
0：
H
(
x
0
)
(
x
?
x
0
)
+
?
f
(
x
0
)
=
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)
解得 x 1 = x 0 ? H ( x 0 ) ? 1 ? f ( x 0 ) (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)
其中
H
H
H 是 Hessian Matrix （海森矩阵），其实就是个二阶导的矩阵。
对于
n
n
n 元变量：
H
(
f
)
=
[
?
2
f
?
x
1
2
?
2
f
?
x
1
?
x
2
?
?
2
f
?
x
1
?
x
n
?
2
f
?
x
2
?
x
1
?
2
f
?
x
2
2
?
?
2
f
?
x
2
?
x
n
?
?
?
?
?
2
f
?
x
n
?
x
1
?
2
f
?
x
n
?
x
2
?
?
2
f
?
x
n
2
]
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???????????????????
牛顿法的优点是收敛速度快，缺点是需要求矩阵的逆，计算量比较大。
此外，如果矩阵非正定（在一维情况下表现为泰勒展开的二阶导数小于
0
0
0 ），极值点为极大值，而非极小值。
如果初始位置离最优点太远，也会导致迭代过程中目标函数不严格递减的情况。
解决这个问题的方法之一是拟牛顿法（Quasi-Newton Methods）。
高斯-牛顿法（Gauss–Newton algorithm）
介绍：
高斯-牛顿法（Gauss–Newton algorithm）是牛顿法的特例，它是牛顿法的修改版，用于寻找函数的最小值。
和牛顿法不一样，它只能用于解决最小二乘问题。
优点是，不需要二阶导数（二阶导数可能很难计算）。
例子：
设我们有
m
m
m 组样本（数据集）：
{
(
x
(
1
)
,
y
(
1
)
)
,
(
x
(
2
)
,
y
(
2
)
)
,
…
,
(
x
(
m
)
,
y
(
m
)
)
}
，
????
x
∈
R
n
\{ (\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)}
(?) 表示第几组样本（学吴恩达的表示法）。
我们希望找到包含
n
n
n 个参数的非线性函数
f
(
x
,
θ
)
f(\boldsymbol{x}, \boldsymbol{\theta})
f(x,θ)，拟合上面
m
m
m 组数据。
其中
θ
\boldsymbol{\theta}
θ 是我们要找的参数，
θ
∈
R
n
\boldsymbol{\theta} \in \mathbb{R}^n
θ∈Rn。
为了方便描述，设代入第 i i i 个样本后，函数值为 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 = 1 m ∥ 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)
（一般用 J J J 表示代价函数，但是这里用了 ? \epsilon ?，原因是稍后要用 J J J 表示雅可比矩阵）
我们需要找到一组参数 θ = [ θ 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)对第
j
j
j 个参数
θ
j
\theta_j
θj? 求导：
?
?
(
θ
)
θ
j
=
∑
i
=
1
m
2
?
(
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∑m?2?(f(i)(θ)?y(i))??θj??f(i)(θ)?(9)
令雅可比矩阵（jacobi matrix）
J
J
J 为：
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)(θ)????????????????????
所以
J
J
J 是一个
m
m
m 行
n
n
n 列的矩阵，
m
×
n
m\times n
m×n。
第
1
1
1 行的含义为，第
1
1
1 个样本对于每个参数的偏导。
令残差（residual，样本与拟合值之间的差）
r
\textbf{r}
r 为：
r
=
[
f
(
1
)
(
θ
)
?
y
1
f
(
2
)
(
θ
)
?
y
2
?
f
(
m
)
(
θ
)
?
y
m
]
\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)(θ)?y1?f(2)(θ)?y2??f(m)(θ)?ym???????????
则式子(9)可以写成矩阵形式：
?
?
(
θ
)
=
2
J
T
r
(10)
\nabla \epsilon(\boldsymbol{\theta}) = 2 J^T \textbf{r} \tag{10}
??(θ)=2JTr(10)
接下来求海森矩阵第
k
k
k 行
j
j
j 列的元素：
?
2
?
(
θ
)
?
θ
k
?
θ
j
=
?
?
θ
k
?
?
(
θ
)
?
θ
j
=
?
?
θ
k
(
∑
i
=
1
m
2
?
(
f
(
i
)
(
θ
)
?
y
(
i
)
)
?
?
f
(
i
)
(
θ
)
?
θ
j
)
=
2
∑
i
=
1
m
(
f
(
i
)
(
θ
)
?
θ
k
f
(
i
)
(
θ
)
?
θ
j
+
(
f
(
i
)
(
θ
)
?
y
(
i
)
)
?
?
2
f
(
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∑m?2?(f(i)(θ)?y(i))??θj??f(i)(θ)?)=2i=1∑m?(?θk?f(i)(θ)??θj?f(i)(θ)?+(f(i)(θ)?y(i))??θk??θj??2f(i)(θ)?)?(11)
令 H H H 为 ? ( θ ) \epsilon(\boldsymbol{\theta}) ?(θ) 的海森矩阵，由(11)可知： H = 2 ? ( J T J + S ) H= 2 \cdot (J^T J + S) H=2?(JTJ+S) 其中：
S k , j = ∑ i = 1 m ( f ( i ) ( θ ) ? y ( i ) ) ? ? 2 f ( i ) ( θ ) ? θ k ? θ j S_{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
?
(
J
T
J
+
S
)
?
1
?
?
J
T
r
(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) 中的
S
S
S 可以忽略，最终高斯-牛顿法的迭代公式为：
θ
n
+
1
=
θ
n
?
(
J
T
J
)
?
1
?
?
J
T
r
(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)