牛顿法

在实际的算法中，很少需要达到最大可能的精度。在现实世界的数据中，建模和测量误差通常比浮点数舍入导致的误差大了几个数量级，我们通常很愿意选择一种近似方法，以牺牲精度为代价换取速度。

接下来，我们将介绍这种近似的数值算法中最重要的组成部分之一：牛顿法。

Newton’s Method

牛顿方法是一种简单但非常强大的算法，用于寻找实值函数的近似根，即以下通用方程的解：

$f(x) = 0$

对于函数 $f$ ，我们只假设它至少有一个根，并且在搜索区间上 $f(x)$ 是连续且可微的。虽然也有一些令人厌烦的异常情况，但在实际应用中几乎从未出现过这种情况，所以我们可以非正式地说这个函数是"好”的。

牛顿法的主要思想是以某个初始点 $x_0$ 开始，然后不断通过在函数图像上各点处切线与x轴交点的x坐标作为新的近似点来迭代优化。这个算法的直观理解是，如果f是“良好”的（光滑且单调），并且当前的近似点 $x_i$ 已经足够接近根了，那么下一个近似点 $x_{i+1}$ 会更接近这个根。

为了获取交点 $x_n$ ，我们需要另其切线函数等于零：

$0 = f(x_i) + (x_{i+1} - x_i)f'(x_i)$

由此可以得到：

$x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}$

牛顿方法非常重要：它是科学领域和工程领域中各种优化求解器的基础。

平方根

作为一个简单的例子，让我们导出求平方根问题的算法：

$x = \sqrt{n} \iff x^2 = n \iff f(x) = x^2 - n = 0$

如果我们将 $f(x) = x^2 - n$ 替换为上述通用公式，可以获得如下更新的规则：

$x_{i+1} = x_i - \frac{x_i^2 - n}{2x_i} = \frac{x_i + n / x_i}{2}$

在实践中，我们还希望在它足够接近正确答案时立即停止它，我们可以在每次迭代后简单地检查：

const double EPS = 1e-9;

double sqrt(double n) {
    double x = 1;
    while (abs(x * x - n) > eps)
        x = (x + n / x) / 2;
    return x;
}

该算法虽然在许多函数上都能收敛，但只有在某些特定的函数集合（比如凸函数）上，才能保证算法是可靠且可证明的。另外一个问题是如果算法能够收敛，那么它的收敛速度是多快。

收敛速度

让我们运行几次牛顿法来寻找2的平方根。初始值 $x_0$ 设定为1，然后在每次迭代之后，我们检查算法得到的结果有多少位数是准确的：

$\pmb{1}.0000000000000000000000000000000000000000000000000000000000000\\ \pmb{1}.5000000000000000000000000000000000000000000000000000000000000\\ \pmb{1.41}66666666666666666666666666666666666666666666666666666666675\\ \pmb{1.41421}56862745098039215686274509803921568627450980392156862745\\ \pmb{1.41421356237}46899106262955788901349101165596221157440445849057\\ \pmb{1.41421356237309504880168}96235025302436149819257761974284982890\\ \pmb{1.41421356237309504880168872420969807856967187537}72340015610125\\ \pmb{1.4142135623730950488016887242096980785696718753769480731766796}$

仔细观察，我们能够发现每次迭代准确的数字位数会大约翻倍。这种惊人的收敛速度并非巧合。

要定量分析这种收敛速度，我们需要考虑在第 $i$ 次迭代时的相对误差 $δ_i$ ，以及在下一次迭代时的误差 $δ_{i+1}$ 。我们需要确定 $δ_{i+1}$ 相对于 $δ_i$ 变小了多少，以此来衡量我们的算法收敛速度提升了多少：

$|δ_i| = \frac{x_n - x}{x}$

我们可以将 $x_i$ 表示为 $x\cdot(1 + δ_i)$ 。将其带入对应的牛顿迭代公式并约掉两边的 $x$ 可以得到：

$1 + \delta_{i+1} = \frac{1}{2}(1 + \delta_i + \frac{1}{1+\delta_i}) = \frac{1}{2}(1 + \delta_i + 1 - \delta_i + \delta^2_i+o(\delta_i^2)) = 1 + \frac{\delta_i^2}{2} + o(\delta_i^2)$