这个问题的起因是我在做leetcode上的一道题延伸出来的，原题目是Leetcode 69. Sqrtx。不过这道题目还是做了很大的简化，它只让我们找出符合要求的整数解。在实际运算中，我们显然希望获得一个开平方的精确小数值，比如说像 \(\sqrt{2} = 1.414214, \sqrt{3} = 1.732051\) 这个样子。这里就需要使用牛顿法来进行求解了。这也是目前计算机学科中进行开平方运算的标准方法。下面我就来详细说说对这一算法的分析。
牛顿法的解题思路
牛顿法告诉我们：
假设方程\(f(x) = 0\)在\(c\)处有一个根，那么可以用以下这个迭代的式子： \[ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} \] 依次计算出\(x_1\)，\(x_2\)，\(x_3\)......，如此数列将无限逼近于方程的根(\(c\))。
我们从几何意义上来理解牛顿法，会比较容易。根据点斜式，可以很轻松的写出在\((x_n, f(x_n))\)处的切线方程为\(y = f'(x_n) * (x - x_n) + f(x_n)\)。递推式中的\(\frac{f(x_n)}{f'(x_n)}\)这一项，相当于是这个切线方程与\(x\)轴相交后形成的直角三角形在\(x\)轴上的距离。每次做这样的递推，就会越来越接近原函数\(f(x) = 0\)的那个零点。
下面这张图很好的展示了这种数列逼近的效果:
由此，我们就可以把牛顿法运用到开平方的运算中来了。对于求解\(\sqrt{a}\)这个式子，相当于求解方程\(f(x) = x^2 - a = 0\)。也就是函数\(f(x)\)的零点。\(f(x)\)的一阶导数可以写为\(f'(x) = 2x\)，这样，牛顿迭代式为： \[ x_{n+1} = x_n - \frac{x_n^2 - a}{2x_n} = \frac{1}{2}(x_n + \frac{a}{x_n}) \] 我们假设求解\(\sqrt{2}\)，任意设置迭代的初值（非负数），比如为\(x_0 = 1\)，带入上面的迭代式进行运算，可以得到： \[ \begin{aligned} &x_0 = 1, \\ &x_1 = \frac{1}{2} * (1 + \frac{2}{1}) = 1.5, \\ &x_2 = \frac{1}{2} * (1.5 + \frac{2}{1.5}) = 1.41666, \\ &x_3 = \frac{1}{2} * (1.41666 + \frac{2}{1.41666}) = 1.414216, \\ &...... \end{aligned} \] 由此我们也可以写出相应的代码来快速的得到答案，下面是一个Princeton大学写的一个Java语言的版本：
// 使用牛顿法解决开平方问题 |
对于1e9+7的运算迭代结果如下
a = 1.000000007E9 |
牛顿法背后的数学原理
我们作为后人，当然可以直接使用400年前牛顿发明的这种算法而不假思索。但是，如果我们作为一个算法最初的设计者，就必须要想的更多一些，那就是：为什么这样做一定可以找到一个数的平方根呢？（从数学的角度来说，就是：满足上述递归条件的数列，是否一定收敛于函数\(f(x) = x^2 - a\)的零点呢）只有先在数学上完成这样的证明，我们才有理由说这个算法是严谨有效的。下面我将给出两种对牛顿法正确性的数学证明。
第一种是英文版的维基百科Newton's method词条中给出的证明，我简单用我理解的中文在此复述如下：
根据泰勒定理，任何函数\(f(x)\)当它有连续的二阶导数时，那么靠近函数\(f(x)\)零点（我们假设为\(\alpha\)，\(\alpha = \sqrt{a}\)）的点都可以写为泰勒展开式如下： \[ f(\alpha) = f(x_n) + f'(x_n)(\alpha - x_n) + R_1 \] 在此，泰勒级数展开式的拉格朗日余项为： \[ R1 = \frac{1}{2!}f''(\xi_n)(\alpha - x_n)^2 \] 这里的\(\xi_n\)在\(x_n ～ \alpha\)的范围内，因为\(\alpha\)是函数的一个根，所以上面的泰勒展开式相当于: \[ \begin{aligned} &0 = f(\alpha) = f(x_n) + f'(x_n)(\alpha - x_n) + \frac{1}{2!}f''(\xi_n)(\alpha - x_n)^2 \\ \Leftrightarrow \quad &\frac{f(x_n)}{f'(x_n)} + (\alpha - x_n) = \frac{-f''(\xi_n)}{2f'(x_n)}(\alpha - x_n)^2 \end{aligned} \] 由于我们定义了数列\(\lbrace x_n \rbrace\)的递推式为\(x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\) ，因此上式可以写为： \[ \begin{aligned} \underbrace{\alpha - x_{n+1}}_{\epsilon_{n+1}} &= \frac{-f''(\xi_n)}{2f'(x_n)}(\underbrace{\alpha - x_n}_{\epsilon_n})^2 \\ \Leftrightarrow \quad \epsilon_{n+1} &= \frac{-f''(\xi_n)}{2f'(x_n)} \cdot \epsilon_n^2 \end{aligned} \] 根据数列极限的定义，我们定义\(\epsilon_n = \alpha - x_n\)，对于上式两边取绝对值，有： \[ \begin{aligned} |\epsilon_{n+1}| &= \frac{|f''(\xi_n)|}{2|f'(x_n)|} \cdot |\epsilon_n|^2 \\ \Leftrightarrow \quad \frac{|\epsilon_{n+1}|}{|\epsilon_n|} &= \underbrace{\frac{|f''(\xi_n)|}{2|f'(x_n)|}}_{M} \cdot |\epsilon_n| \end{aligned} \] 对于函数\(f(x) = x^2 - a\)，我们可以得到函数的一阶导数为\(f'(x) = 2x\)，函数的二阶导数为\(f''(x) = 2\)
此时，上式中的\(M\)项可以写为： \[ M = \frac{1}{2} \cdot \frac{|2|}{|2x_n|} = \frac{1}{2|x_n|} \] 由此，可以得到整个右式： \[ \begin{split} M \cdot |\epsilon_n| &= \frac{1}{2|x_n|} \cdot |x_n - \alpha| \\ &= \frac{1}{2} \cdot \Bigg| \frac{x_n - \alpha}{x_n} \Bigg| \\ &= \frac{1}{2} \cdot \bigg| 1 - \frac{\alpha}{x_n} \bigg| \end{split} \] 对于任意的正数\(x_n\)，\(\alpha\)，能够得到整个右式\(M \cdot |\epsilon_n| < 1\)，因此数列一定收敛。
第二种证明是我和高中同学刘浩（不愧是考研的人，数学就是比我好...）讨论的方法，更加简单易懂一些:
需要用到一个简单的推论： \[ \begin{aligned} \forall a, b \in R^+, \quad (a + b)^2 &\ge 4ab \\ \Rightarrow \ \frac{1}{4}(a + b)^2 &\ge ab \\ \Rightarrow \ \frac{1}{2}(a + b) &\ge \sqrt{ab} \end{aligned} \]
此时， \[ \begin{aligned} x_{n+1} &= \frac{1}{2}(x_n + \frac{a}{x_n}) \\ &\ge \sqrt{x_n \cdot \frac{a}{x_n}} = \sqrt{a} \end{aligned} \] 也就是说，我们证明了对数列中的任意一项\(x_n\)(除了初始值\(x_0\)由用户规定外)都有\(x_n \ge \sqrt{a}\)，因此我们得到以下的不等式恒成立： \[ \begin{gather*} \frac{a}{x_n} \le \frac{a}{\sqrt{a}} = \sqrt{a} \le x_n \\ x_{n+1} = \frac{1}{2}(x_n + \frac{a}{x_n}) \le \frac{1}{2}(x_n + x_n) = x_n \end{gather*} \] 即我们证明了数列\(\lbrace x_n \rbrace\)单调递减并且有下界\(\sqrt{a}\)，根据单调有界必有极限的定理（Monotone convergence theorem），可以证明这一递推数列收敛于\(\sqrt{a}\)
令\(x_{n+1} = x_n = \sqrt{a}\)带入检验，发现递推式成立： \[ \begin{aligned} &左式 = x_{n+1} = \sqrt{a} \\ &右式 = \frac{1}{2}(x_n + \frac{a}{x_n}) = \frac{1}{2}(\sqrt{a} + \frac{a}{\sqrt{a}}) = \sqrt{a} \end{aligned} \]
牛顿法的算法复杂度
在第一种证明的时候，我们可以看到牛顿法的收敛速度是平方级别的：\(|\epsilon_{n+1}| \le M \cdot |\epsilon_n|^2\)。因此它的算法复杂度可以这样来表示\(O(\log{n} \cdot F(n))\)，这里的\(F(n)\)是我们在计算迭代式中\(\frac{f(x)}{f'(x)}\)这一项时需要的n位精确度的开销。这里参见了维基百科中Newton's method这一词条的内容。