さて, 求根アルゴリズム最終回はデッカー法とそれを改良したブレント法を紹介します。
ブレント法はMATLABの非線形関数の根をみつける関数fzeroにも用いられている有用なアルゴリズムになっています。
デッカー法(Dekker method)
デッカー法は二分法と割線方を単純に組み合わせた手法です。 現在の探索状況から, 二分法と割線法によりそれぞれ探索点を計算し, 次の探索点として, そのどちらかを採用することで, 二分法と割線法のいいとこ取りと欠点の補い合いを行います。[1][Dekker, 1969, Finding a zero by means of successive linear interpolation, in Dejon and Henrici]
ただし割線法による新しい探索点が採用された場合, 探索点がほとんど進まない場合があります。 それを避けるために, 最低ステップ幅 を予め決めておき, 探索中最適
は進むようにします。
input: function f input: point a, b such that f(a)f(b) < 0 and |f(b)| < |f(a)| constant: e 許容誤差 c = a while |f(b)| > e s = (a*f(b)-b*f(a))/(f(b)-f(a)) # 割線法による探索点 m = (a + b) / 2 # 二分法による探索点 c = b if |b - m| < |b - s|: # mがsよりもbに近い b = m else: # sがmよりもbに近い b = s if |b - c| < delta: # 例外処理 b = b + sign(a-b)*delta if f(a)f(b) > 0: a = c if |f(a)| < |f(b)|: a, b = b, a return b
a: 探索区間の端点のうちbではないほう b: 各ループ終了時における最適値 m: 二分法による新しい探索点 s: 割線法による新しい探索点
各ループ終了時にはbが最適値になっていることに留意すると, mとsのどちらを選ぶかをbに近い方を選択するのは自然のように思えますね。
これは, という関数で,
として, 二分法を実行した結果です。
また, を満たす,
が見つかっていない状況でも, 最初はデッカー法の割線法による探索のみを行い, 探索中に
を満たす,
が見つけられれば, 二分法による探索も追加するという変更を加えれば大丈夫です。
利点: 必ず収束し(二分法), 根の近くが線形であれば最後の数回の収束が早い(割線法)
欠点: それぞれの手法, 特に割線法での探索点が常に採用されると, 融合したうまみがなくなってしまう
Coming Soon 欠点を表わす図
このような問題を回避するために, ブレントは次のようにデッカーのアルゴリズムを改良しました。
ブレント法(Brent method)
補間による探索を続けて行うための条件を追加したのが, このブレント法です。 条件とは, 現在の探索点と補間による新しい探索点の幅が前の探索時の更新幅の を上回ってなければ補間による探索を行わないというものです。 つまり, 二分法を用いた探索幅の更新よりもより早いスピードで探索幅を小さくすることを要求するものであり, これを満たさない場合は二分法による探索点が採用されます。[2][Brent, 1973, Algorighms for Minimization without Derivatives]
なお, 補間は2点を用いる線形補間だけでなく, 3点を用いる逆二次補間3を導入しています。 3点を用いる補間としては放物線補間もありますが, これは任意の2点間で関数が単一の値を持たない場合の精度が悪いため推奨されていません。ちなみに逆補間とは, 標本点 に対して,
を満たす関数から,
を満たす関数を見つけることです。
標本点 [tex: (x{n-2}, f{n-2}), (x{n-1}, f{n-1}), (x_n, f_n)] が与えられたときに, を満たす
は逆二次補間を用いると, 次のように予測されます。
input: function f input: point a, b such that f(a)f(b) < 0 and |f(b)| < |f(a)| constant: epsilon 許容誤差 constant: delta 計算機イプシロン c = a e = d = b - a while |f(b)| > epsilon: alpha = 0.5*(a - b) beta = f(b) / f(c) if |e| < delta or |f(c)| < |f(b)|: # 二分法 d = e = alpha else: if a == c: # 線形補間 p = 2 * alpha * beta q = 1 - beta else: # 逆2次補間 q = f(c) / f(a) r = f(b) / f(a) p = beta * (2 * alpha * q * (q - r) - (b - c) * (r - 1)) q = (q - 1) * (r - 1) * (beta - 1) beta, e = e, d if |2 * p| < |3 * alpha * q| and |p| < |0.5 * beta * q|: d = - p / q else: # 二分法 d = e = alpha if |d| > delta: s = b + d elif alpha > 0: s = b + delta else: s = b - delta c, b = b, s if f(b)*f(a) > 0: a = c d = e = b - c if |f(a)| < |f(b)|: a, b, c = b, c, b return b
擬似コード中 はrelative machine precision(計算機イプシロン)を使用します。 とにかく十分小さい数を入れておけば問題ありません。
擬似コード24行目のif文はオーバーフローを監視するものになっています。
なお, Wikipediaのブレント法のページには, これとは形が違う擬似コードが載せられていますが, ブレント本人のアルゴリズムと同値であるか筆者はまだ確認できていません。
これは, という関数で,
として, 二分法を実行した結果です。
関数がリプシッツ連続性を有していれば, ブレント法による収束次数は1.61...(黄金比)となります。
参考
[1] Dekker, 1969, Finding a zero by means of successive linear interpolation, in Dejon and Henrici
[2] Brent, 1973, Algorighms for Minimization without Derivatives
[3] https://en.wikipedia.org/wiki/Inverse_quadratic_interpolation
まえ
