ブレント法

ブレント法（英: Brent's method）は、二分法、割線法、逆2次補間を組み合わせた、複雑ではあるが広く用いられている、数値解析における求根アルゴリズムの1つである。二分法の安定さを有し、かつ安定でない他の手法と同程度に高速に解が求められる場合もある。可能な限り、より収束の早い割線法や逆2次補間を用い、必要に応じてより安定な二分法に切り替えて解を求めるという手法である。ブレント法は、1969年のセオドラス・デッカー（英語版）による初期のアルゴリズムを元にして、1973年にリチャード・ブレント（英語版）により考案されたものである^[1] 。

デッカー法

二分法と割線法を組み合わせるという考えは、デッカーによるものである。

ここでは、方程式 $f (x) = 0$ の解を求めたいとする。まず、二分法と同様に、 $f (a 0)$ と $f (b 0)$ が互いに逆符号を持つような $a 0$ と $b 0$ の2点を初期値とする。 $f$ が区間 $[a 0, b 0]$ で連続であるとき、中間値の定理により、 $a 0$ と $b 0$ の間に解が存在する。

反復計算の $k$ ステップ目において、以下の3点が用いられる：

$b k$ は現在の反復での値、つまりその時点で推定される $f$ の解。
$a k$ は "反対点"、つまり $f (a k)$ と $f (b k)$ が逆符号を持つような点であり、したがって区間 $[a k, b k]$ に解が含まれる。また、 $| f (b k)|$ は $| f (a k)|$ と等しいか、またはより小さい値とする。したがって $b k$ は $a k$ よりも求める解に近い。
$b k -1$ は1つ前の反復での値（最初の反復計算 ( $k = 0$ ) では $b -1 = a 0$ とする）。

次の反復値を求めるため、2つの値が計算される。1つは割線法により以下の式で求められる。

s=b_{k}-{\frac {b_{k}-b_{k-1}}{f(b_{k})-f(b_{k-1})}}f(b_{k}),

2つめは二分法により求められる。

m={\frac {a_{k}+b_{k}}{2}}

割線法を用いた場合の $s$ が $b k$ と $m$ の間にある場合は次の反復で $b k +1 = s$ となり、そうでない場合は中間点が使用され $b k +1 = m$ となる。

そして、 $f (a k +1)$ と $f (b k +1)$ とが逆符号を持つような、新しい反対点 $a k +1$ が選ばれる。 $f (a k)$ と $f (b k +1)$ が逆符号である場合、反対点は移動せず $a k +1 = a k$ となる。両者が同符号であれば、 $f (b k +1)$ と $f (b k)$ が逆符号となり、新たな反対点は $a k +1 = b k$ となる。

最終的に、 $| f (a k +1)| < | f (b k +1)|$ となった場合は、 $a k +1$ は $b k +1$ よりよい解の候補値となり、その結果 $a k +1$ と $b k +1$ の値が交換される。

以上がデッカー法による反復計算の1回分である。

ブレント法

デッカー法では、関数 f が滑らかな形状であれば、効率的に解が求められる。しかし、毎回の収束計算で割線法が選択される場合もあるが、このときは b_k の収束が非常に遅い（特に b_k − b_k−1 が任意に小さい）この場合、デッカー法を用いると、二分法よりも非常に多くの収束計算が必要となる。

ブレントは、この問題を解決する小さな修正を提案した。割線法で得られた値が次の収束値として採用される前に、満たすべき判定式を追加した。具体的には、以下の2つの不等式である：

与えられた許容値 $\delta$ に対し、もし1つ前のステップで二分法を用いた場合は、不等式

|\delta |<|b_{k}-b_{k-1}|

による判定が行われ、満たしていなければ二分法が選択され、得た値を次の収束値とする。もし1つ前のステップが補間である場合は、不等式

|\delta |<|b_{k-1}-b_{k-2}|

が使用される。

また、1つ前のステップが二分法の場合は、同じく

|s-b_{k}|<{\begin{matrix}{\frac {1}{2}}\end{matrix}}|b_{k}-b_{k-1}|

を満たす必要がある。偽であれば二分法が使用され、得た値が次の収束値となる。1つ前が補間であった場合は、

|s-b_{k}|<{\begin{matrix}{\frac {1}{2}}\end{matrix}}|b_{k-1}-b_{k-2}|

が使用される。この修正により、k番目の収束計算で二分法が使用される場合、追加で行われる収束計算は多くとも

2log_{2}(|b_{k-1}-b_{k-2}|/\delta )

回に収まる、という点が保証される。なぜなら、上記の条件によって、2回連続して行われる補間計算の回数を半分にでき、その後、多くとも

2log_{2}(|b_{k-1}-b_{k-2}|/\delta )

回の収束計算により、値の変化は

\delta

よりも小さくでき、二分法に切り替えられるためである。ブレントは、二分法を使用する収束計算回数を N 回とするとき、必要な収束計算回数が多くとも合計 N² 回であることを証明した。関数 f が滑らかな形状であれば、ブレント法は逆2次補間または線形補間を使用し、収束速度は超線形となる。

さらにブレント法では、もし f(b_k)、f(a_k)、f(b_k−1) が明確な差があれば、（割線法で使用される）線形補間の代わりに逆2次補間を使用する。これにより収束性がわずかに上がる。結果として、s（線形補間または逆2次補間で得られる値）が採用される条件は、以下のように変更される：s が (3a_k + b_k) / 4 と b_kの間に存在するかどうか

使用例

関数 f(x) = (x + 3)(x − 1)² の解（f(x) = 0 となる x の値）を求める場合を例に挙げる。初期値として [a₀, b₀] = [−4, 4/3] を考える。f(a₀) = −25、f(b₀) = 0.48148 である（数値はすべて概数とする）したがって、条件 f(a₀)f(b₀) < 0 と |f(b₀)| ≤ |f(a₀)| を満たしている。

関数 f(x) = (x + 3)(x − 1)² のグラフ

最初の収束計算では、以下の2点 (b₋₁, f(b₋₁)) = (a₀, f(a₀)) = (−4, −25) と (b₀, f(b₀)) = (1.33333, 0.48148) の間で線形補間を行い、s = 1.23256 を得る。(3a₀ + b₀) / 4 と b₀ の間にあるため、この値は採用される。さらに f(1.23256) = 0.22891 であるため、a₁ = a₀、b₁ = s = 1.23256 となる。
2回目の収束計算：(a₁, f(a₁)) = (−4, −25) と (b₀, f(b₀)) = (1.33333, 0.48148) と (b₁ を用いて逆2次補間を行い、f(b₁)) = (1.23256, 0.22891) を得る。これにより、(3a₁ + b₁) / 4 と b₁ の間にある 1.14205 が求められる。不等式 1.14205 − b₁ ≤ b₀ − b₋₁ / 2 を満たすため、この値を採用する。さらに f(1.14205) = 0.083582 であるため、a₂ = a₁、b₂ = 1.14205 となる。
3回目の収束計算：(a₂, f(a₂)) = (−4, −25) と (b₁, f(b₁)) = (1.23256, 0.22891) と (b₂, f(b₂)) = (1.14205, 0.083582) で逆2次補間を行う。収束値として (3a₂ + b₂) / 4 と b₂ の間にある 1.09032 が得られる。しかし、ここでブレントによる追加条件が加わる。不等式 1.09032 − b₂ ≤ b₁ − b₀ / 2 を満たさないため、この値は採用されない。代わりに、区間 [a₂, b₂] の中間点 m = −1.42897 を計算しこれを採用する。f(m) = 9.26891 であり、a₃ = a₂、b₃ = −1.42897 となる。
4回目の収束計算：(a₃, f(a₃)) = (−4, −25) と (b₂, f(b₂)) = (1.14205, 0.083582) と (b₃, f(b₃)) = (−1.42897, 9.26891) で逆2次補間を行う。収束値として 1.15448 が得られるが、この値は区間 (3a₃ + b₃) / 4 と b₃) には存在しない。したがって、中間点 m = −2.71449 に置き換わる。f(m) = 3.93934 であり、したがって a₄ = a₃、b₄ = −2.71449 となる。
5回目の収束計算：逆2次補間により、区間内にある −3.45500 を得る。しかし、1つ前の計算に二分法を用いたため、不等式 −3.45500 − b₄ ≤ b₄ − b₃ / 2 を満たす必要がある。この不等式は偽であるため、中間点 m = −3.35724 を採用する。f(m) = −6.78239 であるので、m は新たな反対点 (a₅ = −3.35724) となり、収束値は (b₅ = b₄) のままである。
6回目の収束計算：b₅ = b₄ であるため逆2次補間を使用できない。したがって、(a₅, f(a₅)) = (−3.35724, −6.78239) と (b₅, f(b₅)) = (−2.71449, 3.93934) の2点で線形補間を行う。s = −2.95064 が得られ、この値はすべての条件を満たしている。しかし1つ前のステップから収束値が変わっていないので、この結果は採用せず、再度、二分法を用いる。s = -3.03587、f(s) = -0.58418 となる。
7回目の収束計算：再度、逆2次補間を行う。s = −2.99436 であり、すべての条件を満たしている。f(s) = 0.089961 であるので、a₇ = b₆、 b₇ = −3.00219 となる（条件 |f(b₇)| ≤ |f(a₇)| を満たすよう、a₇ と b₇ の値が交換される。）
8回目の収束計算：a₇ = b₆ であるため逆2次補間は使用できない。線形補間により s = −2.9999、f(s) = 0.0016 を得る。条件は満たしている。
以降の収束計算で、b₉ = −3 + 6·10⁻⁸ と b₁₀ = −3 − 3·10⁻¹⁵ となり、得られる解が x = −3 に急速に近づく。

（訂正：9回目の計算 : f(s) = -1.4E-07、10回目の計算 : f(s) = 6.96E-12）

アルゴリズムの詳細

input a, b, and a pointer to a subroutine for f
calculate f(a)
calculate f(b)
if f(a) f(b) >= 0 then error-exit end if
if |f(a)| < |f(b)| then swap (a,b) end if
c := a
set mflag
repeat until f(b or s) = 0 or |b − a| is small enough (convergence)
    if f(a) ≠ f(c) and f(b) ≠ f(c) then
         $s:={\frac {af(b)f(c)}{(f(a)-f(b))(f(a)-f(c))}}+{\frac {bf(a)f(c)}{(f(b)-f(a))(f(b)-f(c))}}+{\frac {cf(a)f(b)}{(f(c)-f(a))(f(c)-f(b))}}$  (inverse quadratic interpolation)
    else
         $s:=b-f(b){\frac {b-a}{f(b)-f(a)}}$  (secant rule)
    end if
    if (condition 1) s is not between (3a + b)/4 and b
      or (condition 2) (mflag is set and |s−b| ≥ |b−c| / 2)
      or (condition 3) (mflag is cleared and |s−b| ≥ |c−d| / 2)
      or (condition 4) (mflag is set and |b−c| < | $\delta$ |)
      or (condition 5) (mflag is cleared and |c−d| < | $\delta$ |)
    then
         $s:={\frac {a+b}{2}}$  (bisection method)
        set mflag
    else
        clear mflag
    end if
    calculate f(s)
    d := c  (d is assigned for the first time here; it won't be used above on the first iteration because mflag is set)
    c := b
    if f(a) f(s) < 0 then b := s else a := s end if
    if |f(a)| < |f(b)| then swap (a,b) end if
end repeat
output b  or s (return the root)

アルゴリズムの実装

ブレントは1973年にALGOL 60による記述でこのアルゴリズムを公開した。Netlibにはわずかな修正を行ったFortran版が公開されている^[2]。MATLAB関数 fzero にもブレント法が組み込まれている^[3]。PARI/GPのsolveについても同様である。その他の形で実装されたもの（C++、C言語、Fortran）は、書籍 Numerical Recipes (ISBN 4874085601) で参照できる。Apache Commons Math には、Java実装版が含まれている^[4]。GNU Scientific Library では gsl_root_fsolver_brent にて実装されている^[5]。

参考文献

^ Brent (1973). Algorithms for Minimization without Derivatives. Prentice-Hall, Englewood Cliffs, NJ. Reprinted by Dover Publications, Mineola, New York, January 2002. ISBN 0-486-41998-3.
^ zeroin.f - netlib
^ 非線形関数の根 - MATLAB fzero - MathWorks 日本
^ BrentSolver (Apache Commons Math 4.0 API)
^ One Dimensional Root-Finding — GSL 2.5 documentation

Atkinson (1989). An Introduction to Numerical Analysis (2nd ed.), Section 2.8. John Wiley and Sons. ISBN 0-471-50023-2.
R.P. Brent (1973). Algorithms for Minimization without Derivatives, Chapter 4. Prentice-Hall, Englewood Cliffs, NJ. ISBN 0-13-022335-2.
T.J. Dekker (1969). "Finding a zero by means of successive linear interpolation." In B. Dejon and P. Henrici (eds), Constructive Aspects of the Fundamental Theorem of Algebra, Wiley-Interscience, London. ISBN 0471203009

外部リンク

Wikipedia Brent algorithm - Implementation in Borland Delphi (console application) by vpapanik - 24 Aug 2010
Wikipedia Brent example - Solution Table by vpapanik - 8 Sep 2009
zeroin.f at Netlib.
cpp_brent for C++.
Module for Brent's Method by John H. Mathews
Open source Excel VBA functions and examples, including step by step evaluation of the example problem

[1] Brent (1973). Algorithms for Minimization without Derivatives. Prentice-Hall, Englewood Cliffs, NJ. Reprinted by Dover Publications, Mineola, New York, January 2002. ISBN 0-486-41998-3.

[2] zeroin.f - netlib

[3] 非線形関数の根 - MATLAB fzero - MathWorks 日本

[4] BrentSolver (Apache Commons Math 4.0 API)

[5] One Dimensional Root-Finding — GSL 2.5 documentation

[1]

[2]

[3]

[4]

[5]