在求最优解时，前面很多地方都用梯度下降(Gradient Descent)的方法，但由于最优步长很难确定，可能会出现总是在最优解附近徘徊的情况，致使最优解的搜索过程很缓慢。牛顿法(Newton's Method)在最优解的搜索方面有了较大改进，它不仅利用了目标函数的一阶导数，还利用了搜索点处的二阶导数，使得搜索算法能更准确地指向最优解。
我们结合下图所示的一个实例来描述牛顿法的思想。假设我们想要求得参数\(\theta\)，使得\(f(\theta)=0\)。算法的描述如下：
- 随机猜测一个解\(\theta^{(0)}\)，并令\(t=0\)；
- 在\(\theta^{(t)}\)处用一根切线来近似\(f(\theta)\);
- 求得切线与横坐标的交点\(\theta^{(t+1)}\)，作为下一个可能的解;
- \(t=t+1\)；
- 重复2-4，直到收敛，即\(f(\theta^{(t)})\approx 0\)。
那么\(\theta^{(t+1)}\)与\(\theta^{(t)}\)之间存在怎样的迭代关系呢？由切线的斜率可知 \begin{equation} f'(\theta)=\frac{f(\theta)}{\vartriangle}\Rightarrow \vartriangle=\frac{f(\theta)}{f'(\theta)} \end{equation}
观察\(\theta^{(t+1)}\)与\(\theta^{(t)}\)在横坐标上的关系，可知 \begin{equation} \theta^{(t+1)}=\theta^{(t)}-\vartriangle=\theta^{(t)}-\frac{f(\theta)}{f'(\theta)} \end{equation}
牛顿法给出了\(f(\theta)=0\)的求解算法，那么怎样将其运用到求使似然函数\(\mathcal{L}(\theta)\)最大化的参数上呢？一般最优参数\(\theta^{\star}\)在\(\mathcal{L}(\theta)\)的极值点出取得，即\(\mathcal{L}'(\theta^{\star})=0\)。那么，令上面的\(f(\theta)=\mathcal{L}'(\theta)\)，我们很容易就得出了下列的迭代法则 \begin{equation} \theta^{(t+1)}=\theta^{(t)}-\frac{\mathcal{L}'(\theta^{(t)})}{\mathcal{L}''(\theta^{(t)})} \end{equation} 最终求得使\(\mathcal{L}'(\theta)=0\)的参数\(\theta^\star\)，也就是令似然函数\(\mathcal{L}(\theta)\)最大的参数。
上面讨论的参数\(\theta\in\mathbb{R}\)，我们现在将牛顿法则推广到\(n\)维向量\(\theta\in\mathbb{R}^n\)，对应的迭代法则形式如下： \begin{equation} \theta^{(t+1)}=\theta^{(t)}-H^{-1}\nabla_{\theta}\mathcal{L} \end{equation} 其中\(H\)为\(\mathcal{L}\)对向量\(\theta^{(t)}\)的二阶偏导，称为Hessian矩阵,\(H_{ij}=\frac{\partial^2\mathcal{L}}{\partial\theta^{(t)}_i\partial\theta^{(t)}_j}\)。
接下来，我们从另外一个角度来考察牛顿法。用似然函数\(\mathcal{L}(\theta)\)的二阶泰勒展开\(\mathcal{F}(\theta)\)来对其进行逼近。 \begin{equation} \mathcal{L}(\theta)\approx\mathcal{F}(\theta)=\mathcal{L}(\theta^{(t})+\nabla_{\theta^{(t)}}\mathcal{L}(\theta-\theta^{(t)})+\frac{1}{2}(\theta-\theta^{(t)})^TH(\theta-\theta^{(t)}) \end{equation} 令\(\theta=\theta^{(t+1)}\)，可得 \begin{equation} \begin{array}{ll} \mathcal{F}(\theta^{(t+1)})=&\mathcal{L}(\theta^{(t)})+\nabla_{\theta^{(t)}}\mathcal{L}^T(\theta^{(t+1)}-\theta^{(t)})\\ &+\frac{1}{2}(\theta^{(t+1)}-\theta^{(t)})^TH(\theta^{(t+1)}-\theta^{(t)}) \end{array} \end{equation} 现在，我们的目的是求得使\(\mathcal{F}(\theta^{(t+1)})\)最小的参数\(\theta^{(t+1)}\)。将上式对\(\theta^{(t+1)}\)求导并令导数为0，可得 \begin{equation} \frac{\partial\mathcal{F}(\theta^{(t+1)})}{\partial\theta^{(t+1)}}=\nabla_{\theta^{(t)}}\mathcal{L}+H(\theta^{(t+1)}-\theta^{(t)})=0 \end{equation} 等式两侧同时左乘\(H^{-1}\)，化简得 \begin{equation} \theta^{(t+1)}=\theta^{(t)}-H^{-1}\nabla_{\theta}\mathcal{L} \end{equation}
我们用的是二阶泰勒展开式\(\mathcal{F}(\theta)\)逼近似然函数\(\mathcal{L}(\theta)\)。如果\(\mathcal{L}(\theta)\)确实为二次函数，那么\(\mathcal{F}(\theta)\)就是\(\mathcal{L}(\theta)\)的准确展开式，利用牛顿法一步就可以直接求得最优解。一般情况下，\(\mathcal{L}(\theta)\)并非二次函数，那么\(\mathcal{F}(\theta)\)也就存在逼近误差，使得一次迭代不能求得最优解，当\(\mathcal{L}(\theta)\)的次数很高时，往往要经历很多次迭代。一般而言，因为牛顿法利用了二阶导数来修正搜索方向和步长，收敛速度很更快。但是这同样也是要付出代价的，相比梯度下降而言，我们需要额外计算Hessian矩阵并求其逆，这两步的计算代价都很大。只要参数\(\theta\)的维度\(n\)不是很大，可以考虑用牛顿迭代。另外还有一点，如果目标函数不是严格的凸函数，Hessian矩阵\(H\)很可能是奇异矩阵，也就是存在特征值为0的情况，那么它的逆矩阵是不存在的，也就无法用牛顿法。
今年有一道面试题是要求我们写出一段程序，求解\(\sqrt{n}\)。如果把牛顿法用上去，问题就迎刃而解了。我们设定目标函数为\(f(x)=x^2-n\)，那么令\(f(x)=0\)的解很显然就是\(\pm\sqrt{n}\)。要注意的是，我们要选择 合理的 迭代起始点，如果我们从正数开始迭代，求得的是\(\sqrt{n}\);如果从负数开始迭代，求得的就是\(-\sqrt{n}\)；如果从0开始迭代，会出现未定义的计算(0作为除数)。我们根据前面讲的牛顿迭代法则，直接给出该题的迭代法则 \begin{equation} x^{(t+1)}=x^{(t)}-\frac{f(x^{(t)})}{f'(x^{(t)})}=x^{(t)}-\frac{(x^{(t)})^2-n}{2x^{(t)}}=\frac{1}{2}\left(x^{(t)}+\frac{n}{x^{(t)}}\right) \end{equation} 下面是由该算法写出的一段精简的code，浓缩了牛顿算法的精髓
1
double
mysqrt1(
double
n)
2
{
3
if
(n<
0
)
return
-
1
;
4
if
(n==
0
)
return
0
;
5
double
eps=1e-
5
;
6
double
x=
0.1
;
//
start from a positive value
7
while
(fabs(x*x-n)>=
eps)
8
x=(x+n/x)/
2
;
//
Newton's method
9
return
x;
10
}
这道题我还想了另外一个算法，算法的启发点来源于\((x-1)(x+1)+1=x^2=n\)。用这个算法，我们的迭代起始点可以是0。算法的基本思想如下：给定一个初始步长step，从起始点开始每次向前走一个步长，直到超过了\(\sqrt{n}\)；一旦超过了\(\sqrt{n}\)，就要开始慢慢向最终解靠近，每次前进或后退的步长都缩减为以前的一半。很明显，这个算法没有牛顿迭代法快。我只用了少数几个测试用例，两段程序的计算结果都和sqrt库函数的计算结果一致。代码如下：
1
double
mysqrt2(
double
n)
2
{
3
if
(n<
0
)
return
-
1
;
4
double
x=
0
;
5
double
step=
10
;
6
int
threshold=
0
;
7
double
eps=1e-
5
;
8
double
res=x*
x;
9
while
(fabs(res-n)>=
eps)
10
{
11
if
(res<
n)
12
{
13
//
once we have passed the solution,we must walk forward slowly
14
if
(threshold) step/=
2
;
15
x+=step;
//
walking forward
16
}
17
else
//
walk
18
{
19
threshold=
1
;
//
indicating we have passed the solution
20
step/=
2
;
//
reducing the step size to its half
21
x-=step;
//
walking back
22
}
23
res=x*x;
//
compute x*x to estimate real n
24
}
//
end while
25
return
x;
26
}
