牛顿迭代法学习笔记

· · 算法·理论

果然不愧是牛顿发明的东西，是真的很有用。

简介

牛顿迭代法通常用于寻找方程 f(x)=0 的根的迭代算法，从猜测值开始不断逼近，迭代次数通常自己决定。显然地，如果你迭代的次数越多，你的结果就越精确。

公式

我们记初始的猜测值为 x_0，之后第 n 次迭代的结果为 x_{n}。

考虑计算在点 x_n 处的一阶泰勒展开，有：

f(x) \approx f(x_n) + f'(x_n)(x-x_n)

令 f(x)=0，我们开心地发现了这个结论：

x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}

欸，这不就可以一直迭代近似了吗？这正是牛顿迭代的公式。

收敛证明

请注意收敛证明仍然是不可或缺的一部分。

设有函数 f(x)，且 f 在根 r 附近二次连续可微，也即存在 f(x) 的二阶导，并且 f'(r) \neq 0，也即单根。

注意到泰勒展开：

0=f(r)=f(x_n)+f'(x_n)(r-x_n)+\frac{f''(\xi)}{2}(r-x_n)^2

整理，可以得到：

x_n-\frac{f(x_n)}{f'(x_n)}-r=x_{n+1}-r=\frac{f''(\xi)}{2f'(x_n)}(x_n-r)^2

换元（或者不换，这里就不换了），可以得到它是平方收敛的。

例题

求平方根

题意：求 \sqrt a 的值。

这玩意儿怎么弄呢？欸，我们是否可以弄牛顿迭代呢？

考虑构造函数 f(x)=x^2-a，这在 f(x)=0 时，x=\sqrt a，而这个 f(x)=0 的解是可以使用牛顿迭代的。

注意到 f'(x)=2x，则接下来就不难了。

具体而言，我们迭代的方法为：

x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}=x_n-\frac{x_n^2-a}{2x_n}

这里初始迭代的你可以任选。这里迭代次数选了 10^5。

代码如下：

#include <bits/stdc++.h>
using namespace std;

int main() {
    long double n, cur;
    cin >> n;
    cur = n;
    for (int i = 1; i <= 1e5; i++)
        cur -= ((cur * cur - n) / (2 * cur));
    cout << setprecision(8) << fixed << cur << '\n';
    return 0;
}

求 n 次方根

这道题其实和上面的十分类似。

题意：求 \sqrt [n]{a}。

假设有函数 f(x)=x^n-a，考虑牛顿迭代。

注意到 f'(x)=nx^{n-1}，则：

x_{t+1}=x_t-\frac{f(x_t)}{f'(x_t)}=x_t-\frac{x_t^n-a}{nx_t^{n-1}}

代码是容易的，这里就不给出了。

求三角函数根

举个例子，求 f(x)=x^2-2\cos x=0 在 [0,\frac{\pi}{2}] 内的解。

同样，使用牛顿迭代法。

注意到 f'(x)=2x+2\sin x，则：

x_{n+1}=x_n-\frac{x^2-2\cos x}{2x+2 \sin x}

则同样牛顿迭代法。

求指数函数的根

好吧其实这里附赠了一些导数的内容。

例如：求 f(x)=10^x+2x=100 的实根。

这里顺带讲一下指数函数的求导。（本人六年级自学，请轻点喷）

f'(x)=\frac{\mathrm{d}}{\mathrm{d}x} 10^x+\frac{\mathrm{d}}{\mathrm{d}x} 2x

指数函数求导我们通常可以用链式求导法则，也即：

\frac{\mathrm{d}y}{\mathrm{d}x}=\frac{\mathrm{d}y}{\mathrm{d}u} \frac{\mathrm{d} u}{\mathrm{d} x}

请注意这里形如 \frac{\mathrm{d}y}{\mathrm{d} x} 一定不能想当然以为是分式。

注意到

10^x=e^{\ln(10^x)}=e^{x \ln 10}

使用链式求导法则：

\frac{\mathrm{d}}{\mathrm{d} x} e^{x \ln 10}=e^{x \ln 10} \cdot \frac{\mathrm{d}}{\mathrm{d} x}(x \ln 10)=e^{x \ln 10} \cdot \ln 10=10^x \ln 10

另外一部分直接应用普通法则即可。

则：

f'(x)=10^x \ln 10+2

这里先预处理 \ln 10 就行了。

这里为什么要换底数呢？原因很简单，(e^x)'=e^x。

求对数函数的根

求 f(x)=\log_2(x)+x^2-10=0 的正实数根。

那么经过上面的风雨交加的题目，这里就讲求导部分了。其它部分比较简单。

考虑换底：

\log_2 x=\frac{\ln x}{\ln 2}+x^2-10

则函数变为：

f(x)=\frac{\ln x}{\ln 2}+x^2-10

换底的原因同样，因为 \frac{\mathrm{d}}{\mathrm{d}x}\ln x=\frac{1}{x}。

考虑对数项求导（这里后两项省略，因为比较简单）中间使用了链式求导法则：

\frac{\mathrm{d}}{\mathrm{d}x}\left(\frac{\ln x}{\ln 2}\right)=\frac{1}{\ln 2} \cdot \frac{\mathrm{d}}{\mathrm{d}x}(\ln x)=\frac{1}{x \ln 2}

这样就可以了。预处理 \ln 2 即可。

拓展应用

高维牛顿法

这个有趣的东西通常用于求解多元方程组。

对于多元方程组 \textbf{F}(\textbf{x})=0，其中 \textbf{F} : \mathbb{R}^n \to \mathbb{R}^n，则其公式为：

\mathbf{x}_{n+1}=\mathbf{x}_n-J^{-1}(\mathbf{x}_n)\mathbf{F}(x_n)

其中，J 为雅可比矩阵：

J_{ij}=\frac{\partial F_i}{\partial x_j}

举个例子：

\begin{cases} x^2+y^2=4\\ e^x+y=1 \end{cases}

则雅可比矩阵为

J=\begin{bmatrix} 2x & 2y\\e^x & 1\end{bmatrix}

那么就简单很多了。利用线性代数的知识就可以了。

Hessian 矩阵和极小值

梯度下降法可以求无约束的极小值问题。即 \min f(\mathbf{x})。公式：

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

这里 H 为二阶导矩阵，也称 Hessian 矩阵。

这个就是梯度下降法，对比牛顿法复杂度少了一个 n，快了很多。

优化方法

这里给出一个优化方法，也就是当我们迭代的时候迭代的结果小于一个 \epsilon 来终止，具体而言，我们可以通过这个终止条件进行终止，在不必要的精度上一定程度省去了时间。即 |x_{n+1}-x_n| \le \epsilon。

总结

牛顿迭代法是一个非常好用的东西，我们在很多地方都将用到。

特别提醒

我们对于初始项的选择一定不能与结果误差过大，可能会导致迭代发散。例如函数成分复杂，存在拐点或者不连续的十（sang）分（xin）友（bing）好（kuang）的函数，我们需要特殊处理。顺便说一嘴：f(x)=\sqrt[3]{x} 是一个很好的检查初始值设置的函数。

或者说，对于导数为 0，迭代公式会失效，这个时候牛顿迭代法就不可做了。