数値解析において、ブレント法は二分法、正割法、逆二次補間を組み合わせたハイブリッド求根アルゴリズムです。二分法の信頼性を持ちながら、信頼性の低い方法と同程度に高速化できます。このアルゴリズムは、可能であれば潜在的に高速に収束する正割法または逆二次補間を使用しようとしますが、必要に応じてより堅牢な二分法にフォールバックします。ブレント法はリチャード・ブレント[1]によるもので、テオドラス・デッカー[2]による初期のアルゴリズムに基づいています。そのため、この方法はブレント・デッカー法としても知られています。
ブレント法の現代的な改良には、根の周りで平坦な関数に対してより単純で高速なチャンドルパトラ法、[3] [4] 、2次関数ではなく指数関数的補間を実行して反復のためのより単純な閉じた式を提供するリダーズ法、および最適な最悪ケースと漸近保証を実現する正規偽法と二分法のハイブリッドである ITP法などがあります。
デッカー法
二分法と割線法を組み合わせるというアイデアは、Dekker (1969) に遡ります。
方程式f ( x ) = 0を解きたいとします。二分法と同様に、デッカー法を 2 つの点、たとえばa 0とb 0で初期化し、 f ( a 0 ) とf ( b 0 ) の符号が逆になるようにする必要があります。 fが [ a 0 , b 0 ]で連続している場合、中間値定理により、 a 0とb 0の間に解が存在することが保証されます。
すべての反復には 3 つのポイントが関係します。
- b k は現在の反復、つまりfの根の現在の推定値です。
- a kは「対点」、つまりf ( a k ) とf ( b k ) が反対の符号を持つ点であり、区間 [ a k , b k ] に解が含まれます。さらに、| f ( b k )| は | f ( a k )|以下である必要があるため、未知の解に対する推定値はa kよりもb k のほうが適切です。
- b k −1は前の反復です(最初の反復では、b k −1 = a 0に設定します)。
次の反復のための 2 つの暫定値が計算されます。最初の値は、セカント法とも呼ばれる線形補間によって与えられます。
2番目は二分法で与えられる
セカント法の結果sがb kとmの間に厳密に収まる場合は、それが次の反復 ( b k +1 = s ) になり、それ以外の場合は中点 ( b k +1 = m ) が使用されます。
次に、新しい対点の値は、f ( a k +1 ) とf ( b 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の値は交換されます。
これで、Dekker 法の 1 回の反復の説明は終了です。
デッカー法は、関数fが適度に適切に動作する場合、適切に機能します。ただし、すべての反復でセカント法が使用される状況がありますが、反復b k は非常にゆっくりと収束します (特に、| b k − b k −1 | が任意に小さい場合があります)。この場合、デッカー法では二分法よりもはるかに多くの反復が必要になります。
ブレント法
ブレント (1973) は、デッカー法の問題を回避するために小さな修正を提案しました。彼は、セカント法の結果を次の反復として受け入れる前に満たさなければならない追加のテストを挿入しました。2 つの不等式が同時に満たされる必要があります。
特定の数値許容値を指定すると、前のステップで二分法を使用した場合、補間を実行するには不等式が保持される必要があります。それ以外の場合は二分法が実行され、その結果が次の反復に使用されます。
前のステップで補間が実行された場合は、代わりに不等式を使用して次のアクションを実行し、補間(不等式が真の場合)または二分法(不等式が真でない場合)を選択します。
また、前のステップで二分法を使用した場合は、不等式が 成立する必要があります。そうでない場合は、二分法が実行され、その結果が次の反復に使用されます。前のステップで補間が実行された場合は、代わりに不等式 が使用されます。
この修正により、k 回目の反復で、最大で回の追加の反復で二分ステップが実行されるようになります。これは、上記の条件により、連続する補間ステップ サイズが 2 回の反復ごとに半分にされ、最大で回の反復後に、ステップ サイズが より小さくなり、二分ステップが呼び出されるためです。Brent は、彼の方法では最大でN 2回の反復が必要であることを証明しました。ここで、N は二分法の反復回数を表します。関数fが適切に動作する場合、Brent の方法は通常、逆二次補間または線形補間のいずれかによって進行し、その場合は超線形 に収束します。
さらに、ブレント法では、セカント法で使用される線形補間の代わりに逆二次補間を使用します。 f ( b k )、f ( a k )、およびf ( b k −1 ) が異なる場合、効率がわずかに向上します。その結果、s (線形補間または逆二次補間によって提案された値) を受け入れるための条件を変更する必要があります。s は(3 a k + b k ) / 4 とb kの間になければなりません。
アルゴリズム
a、b 、およびfの関数(へのポインタ)を入力し、 f ( a )
を計算する。f ( b )
を
計算します。f ( a ) f ( b ) ≥ 0の場合
ルートが括弧で囲まれていないため、関数を終了します。
end if
if | f ( a )| < | f ( b )| then
swap ( a , b )
end if
c := a
set mflag
f ( b or s ) = 0または| b − a |になるまで繰り返します が十分小さい(収束) if f ( a ) ≠ f ( c )かつf ( b ) ≠ f ( c )なら(逆二次補間)そうでないなら(セカント法)終了if (条件 1) sがbの間にないor (条件 2) (mflagが設定されていて| s − b | ≥ | b − c |/2) or (条件 3) (mflagがクリアされていて| s − b | ≥ | c − d |/2) or (条件 4) (mflagが設定されていて| b − c | < | δ |) or (条件 5) (mflagがクリアされていて| c − d | < | δ |) then (二分法) mflag
を設定else mflag を
クリアend if f ( s )
d := c
を計算する(d はここで初めて割り当てられます。mflag が設定) c := b if f ( a ) f ( s ) < 0 then b := s else a := s終了 if if | f ( a )| < |
f ( b )| then
swap ( a , b )
end if
end 繰り返し
出力 b または s (ルートを返す)
例
f ( x ) = ( x + 3)( x − 1) 2で定義される関数の零点を探しているとします。
初期区間として [ a0 , b0 ] =[−4,4/3]をとります。
f ( a 0 ) = −25、f ( b 0 ) = 0.48148(このセクションのすべての数値は四捨五入されています)なので、条件f ( a 0 ) f ( b 0 ) < 0および| f ( b 0 )| ≤ | f ( a 0 )|が満たされます。

- 最初の反復では、( b −1、f ( b −1 )) = ( a 0、f ( a 0 )) = (−4、−25) と ( b 0、f ( b 0 )) = (1.33333、0.48148) の間で線形補間を使用し、 s = 1.23256を生成します。これは、(3 a 0 + b 0 ) / 4 とb 0の間にあるため、この値が受け入れられます。さらに、f (1.23256) = 0.22891 であるため、a 1 = a 0およびb 1 = s = 1.23256 に設定します。
- 2 回目の反復では、 ( a 1、f ( a 1 )) = (−4、−25) と ( b 0、f ( b 0 )) = (1.33333、0.48148) と ( b 1、f ( b 1 )) = (1.23256、0.22891) の間で逆二次補間を使用します。これにより 1.14205 が生成され、これは (3 a 1 + b 1 ) / 4 とb 1の間にあります。さらに、不等式 |1.14205 − b 1 | ≤ | b 0 − b −1 | / 2 が満たされるため、この値が受け入れられます。さらに、f (1.14205) = 0.083582 であるため、 a 2 = a 1およびb 2 = 1.14205に設定します。
- 3 回目の反復では、 ( a 2 , f ( a 2 )) = (−4, −25) と ( b 1 , f ( b 1 )) = (1.23256, 0.22891) と ( b 2 , f ( b 2 )) = (1.14205, 0.083582) の間で逆二次補間を使用します。これにより 1.09032 が生成され、これは (3 a 2 + b 2 ) / 4 とb 2の間にあります。しかし、ここでブレントの追加条件が適用されます。不等式 |1.09032 − b 2 | ≤ | b 1 − b 0 | / 2 は満たされないため、この値は拒否されます。代わりに、区間 [ a 2 , b 2 ] の中点m = −1.42897が計算されます。f ( m ) = 9.26891なので、 a 3 = a 2、b 3 = −1.42897と設定します。
- 4 回目の反復では、( a 3、f ( a 3 )) = (−4、−25) と ( b 2、f ( b 2 )) = (1.14205、0.083582) と ( b 3、f ( b 3 )) = (−1.42897、9.26891) の間で逆二次補間を使用します。これにより 1.15448 が生成されますが、これは (3 a 3 + b 3 ) / 4 とb 3の間の区間内にありません。したがって、これは中点m = −2.71449 に置き換えられます。f ( m ) = 3.93934 なので、 a 4 = a 3、b 4 = −2.71449に設定します。
- 5 回目の反復では、逆二次補間によって −3.45500 が得られ、これは必要な間隔内にあります。ただし、前の反復は二分ステップであったため、不等式 |−3.45500 − b 4 | ≤ | b 4 − b 3 | / 2 を満たす必要があります。この不等式は偽であるため、中点m = −3.35724 を使用します。f ( m ) = −6.78239であるため、m が新しい対点 ( a 5 = −3.35724) になり、反復は同じままです ( b 5 = b 4 )。
- 6 回目の反復では、b 5 = b 4であるため、逆二次補間は使用できません。したがって、 ( a 5、f ( a 5 )) = (−3.35724、−6.78239) と ( b 5、f ( b 5 )) = (−2.71449、3.93934) の間で線形補間を使用します。結果はs = −2.95064 となり、すべての条件を満たします。ただし、反復は前のステップで変更されていないため、この結果を拒否して二分法に戻ります。s = -3.03587、f ( s ) = -0.58418 に更新します。
- 7 回目の反復では、再び逆二次補間を使用できます。結果はs = −3.00219 となり、すべての条件を満たします。ここで、f ( s ) = −0.03515 なので、a 7 = b 6およびb 7 = −3.00219 に設定します ( a 7とb 7は、条件 | f ( b 7 )| ≤ | f ( a 7 )| が満たされるように交換されます)。 (正解 : 線形補間 )
- 8 回目の反復では、 a 7 = b 6 であるため、逆二次補間は使用できません。線形補間ではs = −2.99994 となり、これが受け入れられます。 (正解 : )
- 次の反復では、根x = −3 に急速に近づきます: b 9 = −3 + 6·10 −8およびb 10 = −3 − 3·10 −15。 (正解 : 反復 9: f ( s ) = −1.4 × 10 −7、反復 10: f ( s ) = 6.96 × 10 −12 )
実装
- Brent (1973) はAlgol 60 の実装を公開しました。
- Netlib には、若干の変更を加えたこの実装の Fortran 翻訳が含まれています。
- PARI /GPメソッドは
solveメソッドを実装します。 - アルゴリズムの他の実装 (C++、C、Fortran) については、Numerical Recipes の本を参照してください。
- Apache Commons Math ライブラリは、アルゴリズムをJavaで実装します。
- SciPy最適化モジュールは、Python(プログラミング言語)でアルゴリズムを実装します。
- Modelica 標準ライブラリは、Modelicaのアルゴリズムを実装します。
- この関数は、R (ソフトウェア)
unirootでアルゴリズムを実装します。 - この関数はMATLAB
fzeroでアルゴリズムを実装します。 - Boost (C++ ライブラリ) は、 Math ツールキットで
C++の Brent メソッドに基づく 2 つのアルゴリズムを実装します。
- minima.hpp での関数の最小化と、関数の最小値を見つける例。
- 根探索では、Brent のオリジナルよりも現代的で効率的なアルゴリズムである新しい TOMS748 を実装します (TOMS748)。また、例とともに TOMS748 を内部的に使用する Boost.Math 根探索も実装します。
- Optim.jlパッケージは、 Julia(プログラミング言語)でアルゴリズムを実装します。
- エミーコンピュータ代数システム(Clojure(プログラミング言語)で記述)は、単変数関数の最小化のために設計されたアルゴリズムのバリエーションを実装します。
- Code Project でホストされている C# ライブラリでのルート検索。
参考文献
- ^ ブレント 1973
- ^ デッカー 1969
- ^ Chandrupatla, Tirupathi R. (1997). 「導関数を使用せずに非線形関数のゼロを見つけるための新しいハイブリッド2次/二分アルゴリズム」。エンジニアリングソフトウェアの進歩。28 (3): 145–149。doi : 10.1016 /S0965-9978(96)00051-8。
- ^ 「10 個の小さなアルゴリズム、パート 5: 二次極値補間と Chandrupatla 法 - Jason Sachs」。
- ブレント、RP (1973)、「第 4 章: 関数のゼロを見つけるための収束が保証されたアルゴリズム」、導関数なしの最小化アルゴリズム、イングルウッド クリフス、ニュージャージー: プレンティス ホール、ISBN 0-13-022335-2
- Dekker, TJ (1969)、「連続線形補間によるゼロの検出」、Dejon, B.、Henrici, P. (編)、『代数の基本定理の構成的側面』、ロンドン: Wiley-Interscience、ISBN 978-0-471-20300-1
さらに読む
- アトキンソン、ケンドール E. (1989)。「セクション 2.8」。数値解析入門(第 2 版)。John Wiley and Sons。ISBN 0-471-50023-2。
- Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007)。「セクション 9.3. Van Wijngaarden–Dekker–Brent 法」。数値レシピ:科学計算の芸術(第 3 版)。ニューヨーク: Cambridge University Press。ISBN 978-0-521-88068-8. 2011年8月11日時点のオリジナルよりアーカイブ。2012年2月28日閲覧。
- ジョージア州アレフェルド。ポトラ、FA;シー・イーシュン(1995 年 9 月)。 「アルゴリズム 748: 連続関数のゼロを囲む」。数学ソフトウェア上の ACM トランザクション。21 (3): 327–344。土井: 10.1145/210089.210111。S2CID 207192624。
外部リンク
- Netlibの zeroin.f 。
- C++ のモジュール brent (C、Fortran、Matlab も可) 2018-04-05 にWayback Machineでアーカイブされました(John Burkardt 著)
- GSL 実装。
- Boost C++ 実装。
- Python (Scipy) 実装
