Brent求根方法

我们已经看到二分法一定能收敛，但其收敛速度是线性，比较慢。牛顿法在一般情况下是二阶收敛的，虽然收敛速度很快，但是迭代初值对其收敛特性有很大影响，很有可能不收敛。在数学形式上，这两种算法都十分简单直接，但由于各自的缺陷，工程上很少直接使用这两种法求解。我们总是希望算法能稳定找到方程的根，同时有超线性的收敛速度。

Brent 方法在效率和可靠性之间做了完美的取舍，不需要求导数，可以说是又快又稳，各大主流科学计算平台和工程软件库都在使用它。本文我们从割线法开始，不断改进，直至 Brent 方法。

1. 割线法

牛顿法收敛很快，但在实际计算过程中效率不一定很高，因为它需要计算 $f(x)$ 的导数。但是导函数 $f'(x)$ 通常都比较复杂，需要更多的数学运算。比如我们有一个函数 $h(x) = f(x)/g(x)$，按照求导法则有，$h'(x) = \frac{f'(x)g(x) - f(x)g'(x)}{g^2(x)}$，显然$h'(x)$的形式要复杂很多。为了避免求导运算，我们可以通过割线来近似。根据导数的定义，我们有: $$ f'(x_{n}) = \lim_{x \to x_{n}} \frac{f(x) - f(x_{n})}{x - x_{n}} $$

如果 $x_{n-1}$ 距离 $x_{n}$ 很近，可以有如下近似:

$$ f'(x_{n}) \approx \frac{f(x_{n-1}) - f(x_{n})}{x_{n-1} - x_{n}} $$

代入牛顿法的迭代公式 $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$，有:

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

如果我们使用上式进行迭代，只需要调用一次 $f(x)$ 就可以了，计算量将大幅减少。这就是所谓的割线法，具体实现参见 SecantRoot。人们之所以称之为割线法，是因为它的迭代公式的几何意义是，用点 $(x_{n-1},f(x_{n-1}))$ 和 $(x_n, f(x_n))$ 构造的一条直线(称之为割线)与 $x$ 轴的交点。

2. 试位法

本质上，割线法只是对牛顿法的迭代公式的一种近似，所以它也是初值敏感的，迭代过程可能发散。在工程中计算效率固然重要，稳定性也是同等甚至是更重要的。考虑到二分法总是能够收敛的，人们在割线法的基础上研究出了一种试位法。

试位法仍然采用割线的几何意义计算下一个点，不同的是它会检查新点的函数值符号，始终维护一个包含根的区间。如此来保证算法一定收敛。

如右图所示，当前迭代根据 $(x_{n-1},f(x_{n-1}))$ 和 $(x_n, f(x_n))$ 构造的割线与 $x$ 轴相交在 $x_{n+1}$ 处。按照割线法的逻辑，下一轮要用 $x_n$ 和 $x_{n+1}$ 迭代了。从图中我们可以明显看出 $f(x_{n}), f(x_{n+1})$ 都是负数，我们不能确定再按照割线法迭代下去是否会收敛。但是 $f(x_{n+1}), f(x_{n-1})$ 是异号的，根据二分法的逻辑，在 $[x_{n+1}, x_{n-1}]$ 之间一定存在一个根。试位法将选择 $x_{n+1}, x_{n-1}$ 进行下一迭代。

试位法要始终维持一个包含根的闭区间 $[a,b]$，满足 $f(a) \cdot f(b) < 0$。所以选择初始迭代点的时候，应当保证它们的函数值异号。在迭代过程中，按照割线法的迭代公式计算新的迭代点，然后检查 $f(x_{n+1})$ 与 $f(x_n)$ 的符号，如果异号用 $x_{n+1}, x_{n}$ 进行下一轮迭代，否则用 $x_{n+1},x_{n-1}$。试位法的迭代区间会像二分法那样不断缩小，直到最后收敛到根上。具体实现代码参见 FalsePosition。

结合了两者，人们希望试位法能像割线法那样快，同时还能如二分法那样保证收敛。实际上很多时候它只能保证收敛，收敛速度并不理想。如上面的示例，如果迭代区间的一端长期不更新，就会导致收敛速度退化为线性的。此时如果不能保证大多数时候迭代区间都能缩短一半，它将甚至比二分法还慢。

3. Dekker 算法

在试位法的基础上，Dekker 算法进一步考虑了迭代区间的缩小长度，如果不能像二分法那样至少缩短一半，就采用二分法来计算新的迭代点。相比于试位法，Dekker 算法主要有两个方面的改动。

其一，对于第 $k$ 轮迭代的两个端点 $x_{n-1}, x_n$，如果 $|f(x_n)| < |f(x_{n-1})|$，那么 Dekker 算法就会认为 $x_n$ 是更好的估计。所以在每轮迭代之初， Dekker 算法都会根据比较两个端点的函数值的绝对值，并进行必要的交换操作，保证 $|f(x_n)| < |f(x_{n-1})|$。这只能说是一种启发式规则，大多数情况下是这样的。因为绝对值 $|f(x_n)|$ 更小，并不能保证 $x_n$ 是更接近根。

其二，每轮迭代时 Dekker 算法都会根据割线法和二分法的迭代公式分别计算采样点 $s$ 和 $m$。

$$ m = \frac{x_{n-1} + x_n}{2} \qquad s = \begin{cases} x_n - \frac{f(x_n)(x_n - x_{n-1})}{f(x_n) - f(x_{n-1})}, & f(x_n) \neq f(x_{n-1}) \\ m, & \text{otherwise} \end{cases} $$

然后判定 $s$ 是否在 $m$ 和 $x_{n-1}$ 之间。如果不在，那么意味着 $s$ 要么已经不在迭代区间内了，要么还在迭代区间里但是缩短长度大概率没有原区间的一半。前者 $s$ 很可能已经发散了，后者可能不如二分法的收敛程度大。这两种情况，Dekker 算法都会采用 $m$ 作为新的迭代点 $x_{n+1}$。如果在，那么 $s$ 就是一个比较合理的估计，Dekker 算法将用 $s$ 为 $x_{n+1}$。

确定了新的迭代点之后就需要像试位法那样，判定 $f(x_{n+1})$ 与 $f(x_n)$ 是否异号，构造新的一定包含根的迭代区间。详细的代码实现参见 DekkerRoot。

Dekker 算法尝试通过讨论新的迭代点 $s$ 是否在合理的区间内，在二分法和割线法之间做出选择。但是它的合理区间计算的略微有些粗糙，在一些极端情况下 Dekker 还是会退化为线性收敛，线性系数可能还达不到二分法的 0.5，比二分法还慢。

4. Brent 算法

Brent 算法在 Dekker 算法的基础上做了两个方面的改进。① 引入了逆二次插值(inverse quadratic interpolation)，相比于割线法的线性插值而言，它能更准确地估计根的位置。② 对插值点做了更严格的检查。 Brent 证明了他的方法最多需要 $N^2$ 次迭代，其中 $N$ 是二分法所需的迭代次数。该算法的逻辑也很复杂，Wiki 是我所找到的最简洁的描述了。我们按照该描述实现了 BrentRoot。

Brent 算法的每一轮迭代都涉及到 $a_k, b_k, b_{k-1}, b_{k-2}$ 四个点。其中 $b_k$ 是当前对根的最佳估计，$b_{k-1}, b_{k-2}$ 则是前两次 $b_k$ 取值。 $a_k$ 是一个辅助点，它始终保持 $f(a_k) \cdot f(b_k) < 0$。在 Wiki 和 BrentRoot 中，我们分别记为 $a := a_k, b := b_k, c := b_{k-1}, d := b_{k-2}$。引入逆二次插值之后，新的估计点通过下式计算:

$$ m = \frac{x_{n-1} + x_n}{2} \qquad p = \begin{cases} x_n - \frac{f(x_n)(x_n - x_{n-1})}{f(x_n) - f(x_{n-1})}, & f(x_n) \neq f(x_{n-1}) \\ m, & \text{otherwise} \end{cases} $$ $$ s = \begin{cases} \frac{a_k f(b_k) f(b_{k-1})}{(f(a_k) - f(b_k))(f(a_k) - f(b_{k-1}))} + \frac{b_k f(a_k) f(b_{k-1})}{(f(b_k) - f(a_k))(f(b_k) - f(b_{k-1}))} + \frac{b_{k-1} f(a_k) f(b_k)}{(f(b_{k-1}) - f(a_k)) (f(b_{k-1}) - f(b_k))}, & f(a_k) \neq f(b_k) \text{and} f(b_k) \neq f(b_{k-1}) \\ p, & \text{otherwise} \end{cases} $$

判定是否需要回退到二分法，Brent 用了 5 个条件，如下图所示，详细的条件解析参见 Wiki。

5. 完

Brent 算法结合了二分法的可靠性与插值法（割线法、逆二次插值）的高效性。它既保证了在函数性态不佳时像二分法一样稳定收敛，又能在函数光滑时实现超线性收敛，速度远快于单纯二分法，是求解一元方程根的稳健且高效的选择。