前几天我在知乎上看到了一道数列题： \[ a_n=\begin{cases} 2, & n=1\\ \frac{a_{n-1}^2+1}{2a_{n-1}}, & n\ge2 \end{cases} \quad (n\in{\mathbb{N^{\ast}}}) \] 这一道题基于的背景显然是运用牛顿迭代求解平方根。牛顿迭代法的原理无论是从其直观意义亦或是从微分的角度都是很容易理解的，但是其分析就不然了。我记得自己三年以前试着分析过一次，但是因为当时自己太菜所以最后啥都没分析出来（悲）。
通项
对于计算 \(\sqrt k\) 的牛顿迭代法，其递推式如下（基于求取 \(x^2=k\) 的零点） \[ x_{n+1} = x_n - \frac{x_n^2-k}{2x_n} = \frac{x_n^2+k}{2x_n} \] 或者根据巴比伦人当时的更直观的写法： \[ x_{n+1} = \frac{1}{2}\left(x_n + \frac{k}{x_n}\right) \] 递推式的式子很简单，但是求解的方法并不 trivial，在等式两边同时加上以及减去 \(\sqrt{k}\)： \[ \begin{aligned} x_{n+1} - \sqrt k &= \frac{\left(x_n-\sqrt k\right)^2}{2x_n} \\ x_{n+1} + \sqrt k &= \frac{\left(x_n+\sqrt k\right)^2}{2x_n} \\ \end{aligned} \] 上下两式相除（很显然，除非 \(x_1=-\sqrt {2}\) 分母不可能为 \(0\)） \[ \frac{x_{n+1} - \sqrt k}{x_{n+1} + \sqrt k} = \left(\frac{x_{n} - \sqrt k}{x_{n} + \sqrt k}\right)^2 \] 那么非常显然地，可以得出 \[ \frac{x_{n} - \sqrt k}{x_{n} + \sqrt k} = \left(\frac{x_1 - \sqrt k}{x_1 + \sqrt k}\right)^{2^{n-1}} \] 基于上式求解 \(x_n\) \[ x_n=\sqrt k \cdot \frac{1+\left(\frac{x_1 - \sqrt k}{x_1 + \sqrt k}\right)^{2^{n-1}}}{1-\left(\frac{x_1 - \sqrt k}{x_1 + \sqrt k}\right)^{2^{n-1}}} \] 这个递推式的通项就求出来了。最初的两步是我觉得是很难想到的，非常妙。
收敛性
基于这一通项，牛顿迭代法的收敛性也就显然取决于 \(\left|\frac{x_1 - \sqrt k}{x_1 + \sqrt k}\right|\)，比如说当 \[ \left|\frac{x_1 - \sqrt k}{x_1 + \sqrt k}\right| <1 \Rightarrow x_1 > 0 \] \(\{x_n\}\) 迭代收敛至 \(\sqrt k\)，当 \(x_1<0\) 时收敛至 \(-\sqrt k\)，这和直觉是相符的。
收敛速率
不妨令 \(\alpha=\frac{x_1 - \sqrt k}{x_1 + \sqrt k}\)，在 \(|\alpha| < 1\) 时，牛顿迭代法的误差就可以写作： \[ x_n-\sqrt k =\sqrt k \cdot \frac{2\alpha^{2^{n-1}}}{1-\alpha^{2^{n-1}}} \] 相邻项误差之比为： \[ \frac{\left|x_{n+1}-\sqrt k\right|}{\left|x_n - \sqrt k\right|} =\left|\frac{\alpha^{2^n}(1-\alpha^{2^{n-1}})}{\alpha^{2^{n-1}}(1-\alpha^{2^{n}})}\right| = \frac{\alpha^{2^{n-1}} - \alpha^{2^n}}{1-\alpha^{2^{n}}} \] 考虑极限情况： \[ \begin{aligned} \lim_{n\to \infty}\frac{\left|x_{n+1}-\sqrt k\right|}{\left|x_n - \sqrt k\right|} &= \lim_{n\to \infty}\frac{\alpha^{2^{n-1}} - \alpha^{2^n}}{1-\alpha^{2^{n}}} \\ &= \frac{\lim_{n\to \infty}\alpha^{2^{n-1}}-\lim_{n\to \infty}\alpha^{2^{n}}}{1-\lim_{n\to \infty}\alpha^{2^{n}}} \\ &= 0 \end{aligned} \] 这说明牛顿迭代法的收敛速率是超线性的。
当然，有更强的结论： \[ \begin{aligned} \lim_{n\to \infty}\frac{\left|x_{n+1}-\sqrt k\right|}{\left|x_n - \sqrt k\right|^2} &= \lim_{n\to \infty}\frac{\left(1-\alpha^{2^{n-1}}\right)^2}{1-\alpha^{2^{n}}} \\ &= 1 \end{aligned} \] 也就是说牛顿迭代法的收敛速率是平方级（quadratic）的，也就是说，每迭代一次，结果的有效位数就会翻倍。
而与此相对地，二分法的收敛速率显然是线性的，这也就是为什么牛顿迭代法在实践中比二分不知道高到哪里去了。
更快速收敛的迭代法？
根据以上论证，不难观察发现，牛顿迭代法的收敛速率之所以是平方级的，是因为 \(\alpha^{2^n}\) 当中的 \(2\) 的作用。那如果把 \(2\) 换成 \(3\) 呢？
也就是说，如果我们构造 \(x_n\)，使得 \[ \frac{x_{n+1} - \sqrt k}{x_{n+1} + \sqrt k} = \left(\frac{x_{n} - \sqrt k}{x_{n} + \sqrt k}\right)^3 \] 也就是 \[ x_{n+1} = \frac{x_n\left(x_n^2+3k\right)}{3x_n^2+k} \] 那收敛速度不就是三次方了吗？看起来很不错的样子！但是实用吗？
在硬件层面，计算一次牛顿迭代所需要的基本操作是（不计各类访存）:
- 计算 \(x_n^2\) —— 一次乘法
- 计算 \(x_n^2+k\) —— 一次加法
- 计算 \(2x_n\) —— 一次乘法
- 计算 \(x_{n+1}\) —— 一次除法
即总共两次乘法，一次除法，一次加法，其中 \(2x_n\) 的乘法还可以转化为一次加法，对硬件是非常友好的。
而相对地，计算一次上面的三次方收敛的迭代所需要的基本操作是：
- 计算 \(x_n^2\) —— 一次乘法
- 计算 \(x_n^2+3k\) —— 一次加法（如果我们预先算出了 \(3k\)）
- 计算 \(x_n(x_n^2+3k)\) —— 一次乘法
- 计算 \(3x_n^2\) —— 一次乘法（如果我们缓存了 \(x_n^2\)）
- 计算 \(3x_n^2+k\) —— 一次加法
- 计算 \(x_{n+1}\) —— 一次除法
即总共三次乘法，一次除法，两次加法。看起来比牛顿法没有多多少，但是考虑到访存等，估计在硬件上并不划算。当然，在渐进意义上讲，如果一次迭代所需要的时间比一次牛顿迭代所需要时间的 \(\log_23\) 倍要少，那么这个算法还是 “实用” 的。
类似地对于指数的底数是 \(4,5,6,7,8\) 的情况，我们还有如下的递推式： \[ \begin{aligned} x_{n+1} &= \frac{x_n^4+6kx_n^2+k^2}{4x_n\left(x_n^2+k\right)}\\ x_{n+1} &= \frac{x_n\left(x_n^4+10kx_n^2+5k^2\right)}{5x_n^4+10kx_n^2+k^2}\\ x_{n+1} &= \frac{\left(x_n^2+k\right)\left(x_n^4+14kx_n^2+k^2\right)}{2x_n\left(3x_n^2+k\right)\left(x_n^2+3k\right)}\\ x_{n+1} &= \frac{x_n \left(x_n^6+21 k x_n^4+35 k^2 x_n^2+7 k^3\right)}{7 x_n^6+35 k x_n^4+21 k^2 x_n^2+k^3} \\ x_{n+1} &= \frac{x_n^8+28 k x_n^6+70 k^2 x_n^4+28 k^3 x_n^2+k^4}{8 x_n \left(x_n^2+k\right) \left(x_n^4+6 k x_n^2+k^2\right)} \end{aligned} \] 上面的式子一个比一个收敛得快。
比如说用最后一个八次方收敛的，初值为 \(1\) 计算 \(\sqrt 2\)。
第一次迭代：前 \(5\) 位相同。
第二次迭代：前 \(48\) 位相同。
第三次迭代：前 \(391\) 位相同。
第四次迭代：前 \(3135\) 位相同。
每迭代一次正确位数翻八倍！
啊！我刚刚注意到最后的八次方收敛算法就是把三步牛顿迭代缩成一步化简了而已，自己实在被自己蠢到了
但是不管怎么样三次方，五次方的收敛算法还是比较 non-trivial 的？