最近的工作中,在求算子softamx时需要使用牛顿迭代法,记录下。

基本思想

牛顿迭代法的具体内容可以参考 牛顿迭代法的维基百科页面。

几何直觉

观察本文上面的图片,凭借我们的直觉可以发现,**如果在函数$f(x)$的根附近的点$x_n$上画一条切线,这条切线与$x$轴的交点$x_{n+1}$比$x_n$更加接近方程的根**。如果在$x_{n+1}$这个点继续使用上一次的方法,再画一条切线,可以想见新的切线与$x$轴的交点肯定比$x_{n+1}$更接近根,如此迭代就会越来越逼近方程的根。下面这幅图表示的更清晰

所以,据此可以推导出如下的方程, $$ \frac{0 - f(x_n)}{x_{n+1} - x_n} = f’(x_n) $$ 进一步化简可以得到, $$ x_{n+1}=x_{n}-{\frac {f(x_{n})}{f’(x_{n})}} $$ 这就是牛顿迭代法的基本公式。

但是牛顿迭代法不一定总是有效,已有证明牛顿迭代法的二次收敛必须满足以下条件:

  • $f’(x)\neq 0$;
  • 对于所有$x\in I$,其中$I$为区间$[α − r, α + r]$,且$x_{0}$在区间其中$I$内,即$ r\geqslant \left|a-x_{0}\right|$的,对于所有$x\in I$,$f’'(x)$是连续的;
  • $x_{0}$足够接近根 α。

所以使用牛顿迭代法,首先需要选择离方程的根足够近的起点,而且这个起点的切线斜率不能为0。

公式推导

牛顿迭代法的另一种推导方式是使用泰勒展开式 $$ f(x)=f(x_0)+f^\prime(x_0)(x-x_0)+\frac{1}{2}f^{\prime\prime}(x_0)(x-x_0)^2+\dots + \frac{1}{n!}f^{(n)}(x_0)(x-x_0)^n + o(x-x_0)^n $$ 使用一阶展开近似可以得到 $$ f(x)=f(x_0)+f^\prime(x_0)(x-x_0) $$ 化简就可以得到之前的方程(2)。

牛顿迭代法求极值

使用牛顿迭代法可以求函数的极值,通过迭代的方法求方程$f(x)$的极值。根据微积分原理,令$f’(x) = 0$的$x$就是函数的极值所在,同样利用泰勒公式展开到二阶,有 $$ f(x)=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^\prime(x_0)+f^{\prime\prime}(x_0)(x-x_0) = 0 $$ 同样可以得到 $$ x=x_0-{\frac {f’'(x_{0})}{f’(x_{0})}} $$ 这就是牛顿迭代法求极值的理论依据。

指令迭代

假设计算机中有求倒数的指令$y = rec(x) = 1/x$,但是精度不高,如何通过牛顿迭代法提高精度?

可以这么想,假设我们的输入是$a$,那么我们对输入求倒数就等价于求方程$a = 1/x$的根,也就是求方程$f(x) = 1/x -a$的根,那么根据牛顿迭代法,如果我们找到一个初值$x_0$,就可以按照如下的方程来迭代 $$ x_{n+1} = x_n - \frac{f(x_n)}{f’(x_n)}=2x_n - ax_{n}^2 $$ 而刚好这个初值就是指令使用一次之后的结果(相比随意找一个数字,这个值离根更近),即$x_0 = 1/a = rec(a)$。

方程举例

下面我们使用牛顿迭代法用C++的代码求解一个数的立方根,假定这个数是$a$,该问题等价于求方程$x^3 = a$的根,也就是求方程$f(x) = x^3 - a$的根。根据牛顿迭代法,可以按照如下的步骤求根

  1. 确定迭代的终止条件,我们假设假定$|x_n^3 - a |< 0.00001$即停止该迭代;

  2. 确定初始点,即选择合适的$x_0$。很明显如果$a = 0$,方程的根就是0,我们选取1作为初始点;

  3. 确认迭代方程,根据方程(2),我们的迭代方程是

    $$ x_{n+1} = \frac{2x_n}{3}+\frac{a}{3x_n^2} $$

于是,我们的程序如下所示

#include <iostream>
#include <math.h>
using namespace std;

int main()
{
    double a,x0,x1,rsl;
    int times = 0;
    cin >> a;  //输入需要求解的数字

    x0 = 1.0;
    rsl = fabs(x0*x0*x0 - a);

    while(rsl > 0.00001)
    {
        x1 = (2/3.0)*x0 + a/(3.0*x0*x0);
        x0 = x1;
        rsl = fabs(x1*x1*x1 - a);
        times++;
        cout << times << " : " << "actual data -- " << x1*x1*x1 << ", result -- " \
            << rsl << endl;
    }

    cout << "Final x is "<< x1 << " and result is "<< x1*x1*x1 << endl;
    return 0;
}

使用这个程序求解-34.5的立方根结果如下

1 : actual data -- -1271.41, result -- 1236.91
2 : actual data -- -392.257, result -- 357.757
3 : actual data -- -132.242, result -- 97.7418
4 : actual data -- -56.6032, result -- 22.1032
5 : actual data -- -37.2522, result -- 2.75222
6 : actual data -- -34.5672, result -- 0.0672223
7 : actual data -- -34.5, result -- 4.35659e-05
8 : actual data -- -34.5, result -- 1.83391e-11
Final x is -3.25542 and result is -34.5

可以看出通过8轮迭代就找到了-34.5的近似根-3.25542。