BFGS法
修正、加筆に協力し、現在の表現をより自然な表現にして下さる方を求めています。ノートページや履歴も参照してください。(2024年8月)
数理最適化において、ブロイデン・フレッチャー・ゴールドファーブ・シャンノ法(ブロイデン・フレッチャー・ゴールドファーブ・シャンノほう、英: Broyden–Fletcher–Goldfarb–Shanno algorithm)、略してBFGS法とは、無制約非線形最適化問題に対する反復的解法の一つである[1] 。関連の深いDFP法と同様、BFGS法は勾配のプレコンディショニング[訳語疑問点 ]を曲率の情報を用いて行うことにより降下方向を決定する。その際、損失関数のヘッセ行列の推定値を勾配(またはその推定値)のみを用いて(一般化)割線法により漸進的に改善する[2] 。
BFGS法における曲率行列の更新には逆行列の評価を要さないため、計算複雑度 (英語版)は{\displaystyle {\mathcal {O}}(n^{2})}に留まり、ニュートン法の{\displaystyle {\mathcal {O}}(n^{3})}よりも高速である。L-BFGS法もよく用いられ、メモリ使用量を限定できるため、多変数(e.g. >1000)問題に対する解法に適している。BFGS-B法はシンプルなボックス拘束を扱える[3] 。
このアルゴリズムの名前は、チャールズ・ジョージ・ブロイデン (英語版)、ロジャー・フレッチャー、ドナルド・ゴールドファーブ (英語版)、デイビッド・シャンノ (英語版)に因む[4] [5] [6] 。
理論的根拠
[編集 ]{\displaystyle {\boldsymbol {x}}}を{\displaystyle \mathbb {R} ^{n}}上のベクトル、{\displaystyle f({\boldsymbol {x}})}を微分可能なスカラー値関数とし、{\displaystyle {\boldsymbol {x}}}の取り得る値に制限はないものとして、{\displaystyle f({\boldsymbol {x}})}を最小化する最適化問題を考える。
BFGS法は初期推定値{\displaystyle {\boldsymbol {x}}_{0}}から始め、各ステージ毎に反復的により良い推定値へと更新していく。
ステージkにおける降下方向 (英語版) pkはニュートン方程式に類似した次の方程式を解くことにより得られる。
- {\displaystyle B_{k}{\boldsymbol {p}}_{k}=-\nabla f({\boldsymbol {x}}_{k})}
ここでBkはxkにおけるヘッセ行列の推定値であり、各ステージごとにxkにおける目的関数の勾配{\displaystyle \nabla f({\boldsymbol {x}}_{k})}を用いて反復的に更新される。降下方向pkを得たのち、この方向に向けて直線探索を行い、{\displaystyle f({\boldsymbol {x}}_{k}+\gamma {\boldsymbol {p}}_{k})}を最小とするようなスカラーγ > 0を求め、次の点xk+1を決定する。
Bkの更新においては、以下の式であらわされる準ニュートン条件が課せられる。
- {\displaystyle B_{k+1}({\boldsymbol {x}}_{k+1}-{\boldsymbol {x}}_{k})=\nabla f({\boldsymbol {x}}_{k+1})-\nabla f({\boldsymbol {x}}_{k})}
ここで{\displaystyle {\boldsymbol {y}}_{k}=\nabla f({\boldsymbol {x}}_{k+1})-\nabla f({\boldsymbol {x}}_{k})}および{\displaystyle {\boldsymbol {s}}_{k}={\boldsymbol {x}}_{k+1}-{\boldsymbol {x}}_{k}}とおくと、Bk+1は以下の正割方程式を満たす。
- {\displaystyle B_{k+1}{\boldsymbol {s}}_{k}={\boldsymbol {y}}_{k}}
Bk+1が正定値行列であるためには曲率条件sk⊤yk>0が満たされる必要がある。この条件は正割方程式に左からsk⊤をかけることにより検証できる。目的関数が強凸関数でない場合、この条件は明示的に課す必要があり、これはたとえばxk+1を決定する際にウルフ条件を満たす点を選べばよい。
点xk+1におけるヘッセ行列を全て計算するかわりに、ステージkにおける推定値に次のように2つの行列を足すことによりBk+1を計算する。
- {\displaystyle B_{k+1}=B_{k}+U_{k}+V_{k}}
UkおよびVkはどちらも階数1の対称行列であるが、これらの和を取ることにより階数2の対称行列を用いて更新することとなる。対称ランクワン法と比べ、BFGS法とDFP法はどちらも階数2の行列を更新に用いる点が異なる。より単純な手法である対称ランクワン法は階数1の行列を用いて更新を行うが、正定値性が保証されない。Bkの対称性と正定値性を維持するため、更新式は{\displaystyle B_{k+1}=B_{k}+\alpha {\boldsymbol {u}}{\boldsymbol {u}}^{\top }+\beta {\boldsymbol {v}}{\boldsymbol {v}}^{\top }}のように選ぶ。正割条件{\displaystyle B_{k+1}{\boldsymbol {s}}_{k}={\boldsymbol {y}}_{k}}を課すと、{\displaystyle {\boldsymbol {u}}={\boldsymbol {y}}_{k}}および{\displaystyle {\boldsymbol {v}}=B_{k}{\boldsymbol {s}}_{k}}として以下を得る[7] 。
- {\displaystyle \alpha ={\frac {1}{{\boldsymbol {y}}_{k}^{\top }{\boldsymbol {s}}_{k}}}}
- {\displaystyle \beta =-{\frac {1}{{\boldsymbol {s}}_{k}^{\top }B_{k}{\boldsymbol {s}}_{k}}}}
最後に、αおよびβを{\displaystyle B_{k+1}=B_{k}+\alpha {\boldsymbol {u}}{\boldsymbol {u}}^{\top }+\beta {\boldsymbol {v}}{\boldsymbol {v}}^{\top }}に代入するとBk+1の更新式は以下のように書ける。
- {\displaystyle B_{k+1}=B_{k}+{\frac {{\boldsymbol {y}}_{k}{\boldsymbol {y}}_{k}^{\top }}{{\boldsymbol {y}}_{k}^{\top }{\boldsymbol {s}}_{k}}}-{\frac {B_{k}{\boldsymbol {s}}_{k}{\boldsymbol {s}}_{k}^{\top }B_{k}^{\top }}{{\boldsymbol {s}}_{k}^{\top }B_{k}{\boldsymbol {s}}_{k}}}}
アルゴリズム
[編集 ]非線形関数{\displaystyle f:\mathbb {R} ^{n}\to \mathbb {R} }を対象とした無制約最適化問題{\displaystyle {\begin{aligned}{\underset {{\boldsymbol {x}}\in \mathbb {R} ^{n}}{\text{minimize}}}\quad &f({\boldsymbol {x}})\end{aligned}}}を考える。
初期推定解{\displaystyle {\boldsymbol {x}}_{0}\in \mathbb {R} ^{n}}および初期推定ヘッセ行列{\displaystyle B_{0}\in \mathbb {R} ^{n\times n}}から始め、次の各ステップを反復することによりxkは解に収束する。
- 降下方向pkを{\displaystyle B_{k}{\boldsymbol {p}}_{k}=-\nabla f({\boldsymbol {x}}_{k})}を解くことにより求める。
- 1次元最適化(直線探索)を行い、前ステップで求めた降下方向に向う許容しうるステップサイズαkを求める。厳密な直線探索が行われた場合、{\displaystyle \alpha _{k}=\arg \min f({\boldsymbol {x}}_{k}+\alpha {\boldsymbol {p}}_{k})} となる。実用上はαkがウルフ条件を満たすことをもって許容する、非厳密な直線探索で十分なことが多い。
- {\displaystyle {\boldsymbol {s}}_{k}=\alpha _{k}{\boldsymbol {p}}_{k}}とし、{\displaystyle {\boldsymbol {x}}_{k+1}={\boldsymbol {x}}_{k}+{\boldsymbol {s}}_{k}}により推定解を更新する。
- {\displaystyle {\boldsymbol {y}}_{k}={\nabla f({\boldsymbol {x}}_{k+1})-\nabla f({\boldsymbol {x}}_{k})}}を計算する。
- {\displaystyle B_{k+1}=B_{k}+{\frac {{\boldsymbol {y}}_{k}{\boldsymbol {y}}_{k}^{\top }}{{\boldsymbol {y}}_{k}^{\top }{\boldsymbol {s}}_{k}}}-{\frac {B_{k}{\boldsymbol {s}}_{k}{\boldsymbol {s}}_{k}^{\top }B_{k}^{\top }}{{\boldsymbol {s}}_{k}^{\top }B_{k}{\boldsymbol {s}}_{k}}}}により推定ヘッセ行列を更新する。
何らかの基準値ε > 0のもと、勾配のノルムが{\displaystyle ||\nabla f({\boldsymbol {x}}_{k})||\leq \varepsilon }を満たしたとき解が収束したものとみなしアルゴリズムを終了する。
{\displaystyle B_{0}=I}のように選んだ場合、最初のステップは最急降下法と等価となるが、以降のステップはBkがヘッセ行列を推定することにより徐々に改善される。
このアルゴリズムのステップ1はBkの逆行列を用いて実行されるが、この逆行列はステップ5でSherman–Morrisonの公式 (英語版) を用いることにより次のように効率的に求めることができる。
- {\displaystyle B_{k+1}^{-1}=\left(I-{\frac {{\boldsymbol {s}}_{k}{\boldsymbol {y}}_{k}^{\top }}{{\boldsymbol {y}}_{k}^{\top }{\boldsymbol {s}}_{k}}}\right)B_{k}^{-1}\left(I-{\frac {{\boldsymbol {y}}_{k}{\boldsymbol {s}}_{k}^{\top }}{{\boldsymbol {y}}_{k}^{\top }{\boldsymbol {s}}_{k}}}\right)+{\frac {{\boldsymbol {s}}_{k}{\boldsymbol {s}}_{k}^{\top }}{{\boldsymbol {y}}_{k}^{\top }{\boldsymbol {s}}_{k}}}}
この計算は {\displaystyle B_{k}^{-1}}が対称行列であり、 {\displaystyle {\boldsymbol {y}}_{k}^{\top }B_{k}^{-1}{\boldsymbol {y}}_{k}}および sk⊤ykがスカラーであることを用いて次のように展開でき、一時行列を要せず実行することができる。
- {\displaystyle B_{k+1}^{-1}=B_{k}^{-1}+{\frac {({\boldsymbol {s}}_{k}^{\top }{\boldsymbol {y}}_{k}+{\boldsymbol {y}}_{k}^{\top }B_{k}^{-1}{\boldsymbol {y}}_{k})({\boldsymbol {s}}_{k}{\boldsymbol {s}}_{k}^{\top })}{({\boldsymbol {s}}_{k}^{\top }{\boldsymbol {y}}_{k})^{2}}}-{\frac {B_{k}^{-1}{\boldsymbol {y}}_{k}{\boldsymbol {s}}_{k}^{\top }+{\boldsymbol {s}}_{k}{\boldsymbol {y}}_{k}^{\top }B_{k}^{-1}}{{\boldsymbol {s}}_{k}^{\top }{\boldsymbol {y}}_{k}}}.}
したがって、逆行列を求めるための計算を一切することなく、ヘッセ行列の逆行列{\displaystyle H_{k}{\overset {\operatorname {def} }{=}}B_{k}^{-1}}そのものを推定することが可能である[8] 。
初期推定解x0、ヘッセ行列の逆行列の推定値H0から始め、次の各ステップを反復することによりxkは解へと収束する。
- 降下方向pkを{\displaystyle {\boldsymbol {p}}_{k}=-H_{k}\nabla f({\boldsymbol {x}}_{k})}により得る。
- 1次元最適化(直線探索)を行い、前ステップで求めた降下方向に向う許容しうるステップサイズαkを求める。厳密な直線探索が行われた場合、{\displaystyle \alpha _{k}=\arg \min f({\boldsymbol {x}}_{k}+\alpha {\boldsymbol {p}}_{k})} となる。実用上はαkがウルフ条件を満たすことをもって許容する、非厳密な直線探索で十分なことが多い。
- {\displaystyle {\boldsymbol {s}}_{k}=\alpha _{k}{\boldsymbol {p}}_{k}}とし、{\displaystyle {\boldsymbol {x}}_{k+1}={\boldsymbol {x}}_{k}+{\boldsymbol {s}}_{k}}により推定解を更新する。
- {\displaystyle {\boldsymbol {y}}_{k}={\nabla f({\boldsymbol {x}}_{k+1})-\nabla f({\boldsymbol {x}}_{k})}}を計算する。
- {\displaystyle H_{k+1}=H_{k}+{\frac {({\boldsymbol {s}}_{k}^{\top }{\boldsymbol {y}}_{k}+{\boldsymbol {y}}_{k}^{\top }H_{k}{\boldsymbol {y}}_{k})({\boldsymbol {s}}_{k}{\boldsymbol {s}}_{k}^{\top })}{({\boldsymbol {s}}_{k}^{\top }{\boldsymbol {y}}_{k})^{2}}}-{\frac {H_{k}{\boldsymbol {y}}_{k}{\boldsymbol {s}}_{k}^{\top }+{\boldsymbol {s}}_{k}{\boldsymbol {y}}_{k}^{\top }H_{k}}{{\boldsymbol {s}}_{k}^{\top }{\boldsymbol {y}}_{k}}}}によりヘッセ行列の逆行列の推定値を計算する。
最尤推定やベイズ推定などの統計推定問題においては、最終的なヘッセ行列の逆行列を用いて解の信頼区間もしくは確信区間を推定することができる [要出典 ]。しかし、これらの量は正確には真のヘッセ行列により定義されるものであり、BFGS近似は真のヘッセ行列に収束しない場合がある[9] 。
発展
[編集 ]BFGS更新公式は曲率sk⊤ykが常に正であり、ゼロから離れた下界があることに強く依拠している。この条件は凸な対称関数においてウルフ条件を用いた直線探索を用いる場合は満たされるが、実際の問題(たとえば逐次二次計画法)では負やほぼゼロの曲率があらわれることがしばしば発生する。このようなことは非凸関数を対象とする場合や直線探索ではなく信頼領域アプローチをとった場合に生じるおそれがある。この場合、BFGS法は誤った値をあたえることがある。
このような場合には、減衰BFGS更新[10] などと呼ばれる、skおよび/またはykを修正して頑健にした更新式が用いられることがある。
実装
[編集 ]オープンソースの実装として有名なものは以下のようなものがあげられる。
- ALGLIBはC++およびC#用のBFGSおよびL-BFGS法を実装する。
- GNU Octaveの
fsolve関数は信頼領域を用いた一種のBFGS法を用いる。 - GSLはgsl_multimin_fdfminimizer_vector_bfgs2関数としてBFGSを実装している[11] 。
- R言語では、、BFGS法(および矩形拘束を扱えるL-BFGS-B法)が基本関数optim()のオプションとして実装されている[12] 。
- SciPyでは、scipy.optimize.fmin_bfgs関数がBFGS法を実装している[13] 。パラメータLにとても大きな数を指定することにより、なんらかのL-BFGS法を実行することもできる。
- Juliaでは、Optim.jlパッケージにBFGSおよびL-BFGSが実装されている[14] 。
プロプライエタリな実装としては以下のようなものがあげられる。
- 大規模非線形最適化ソフトウェアArtelys KnitroはBFGS法およびL-BFGS法の両方を実装する。
- MATLAB Optimization Toolboxでは、fminunc関数がBFGS法を3次直線探索と組み合わせたアルゴリズムを「中規模スケール」の問題向けに実装している。
- MathematicaにはBFGS法が含まれる。
- LS-DYNAもBFGS法を用いて陰解を求めている。
脚注
[編集 ]- ^ Fletcher, Roger (1987), Practical Methods of Optimization (2nd ed.), New York: John Wiley & Sons, ISBN 978-0-471-91547-8 , https://archive.org/details/practicalmethods0000flet
- ^ Dennis, J. E. Jr.; Schnabel, Robert B. (1983), "Secant Methods for Unconstrained Minimization", Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Englewood Cliffs, NJ: Prentice-Hall, pp. 194–215, ISBN 0-13-627216-9 , https://books.google.com/books?id=ksvJTtJCx9cC&pg=PA194
- ^ Byrd, Richard H.; Lu, Peihuang; Nocedal, Jorge; Zhu, Ciyou (1995), "A Limited Memory Algorithm for Bound Constrained Optimization", SIAM Journal on Scientific Computing 16 (5): 1190–1208, doi:10.1137/0916069 , http://www.ece.northwestern.edu/~nocedal/PSfiles/limited.ps.gz
- ^ Fletcher, R. (1970), "A New Approach to Variable Metric Algorithms", Computer Journal 13 (3): 317–322, doi:10.1093/comjnl/13.3.317
- ^ Goldfarb, D. (1970), "A Family of Variable Metric Updates Derived by Variational Means", Mathematics of Computation 24 (109): 23–26, doi:10.1090/S0025-5718-1970-0258249-6
- ^ Shanno, David F. (July 1970), "Conditioning of quasi-Newton methods for function minimization", Mathematics of Computation 24 (111): 647–656, doi:10.1090/S0025-5718-1970-0274029-X, MR 274029
- ^ Fletcher, Roger (1987), Practical methods of optimization (2nd ed.), New York: John Wiley & Sons, ISBN 978-0-471-91547-8 , https://archive.org/details/practicalmethods0000flet
- ^ Nocedal, Jorge; Wright, Stephen J. (2006), Numerical Optimization (2nd ed.), Berlin, New York: Springer-Verlag, ISBN 978-0-387-30303-1
- ^ Ge, Ren-pu; Powell, M. J. D. (1983). "The Convergence of Variable Metric Matrices in Unconstrained Optimization". Mathematical Programming 27 (2). doi:10.1007/BF02591941.
- ^ Jorge Nocedal; Stephen J. Wright (2006), Numerical Optimization
- ^ "GNU Scientific Library — GSL 2.6 documentation". www.gnu.org. 2020年11月22日閲覧。
- ^ "R: General-purpose Optimization". stat.ethz.ch. 2020年11月22日閲覧。
- ^ "scipy.optimize.fmin_bfgs — SciPy v1.5.4 Reference Guide". docs.scipy.org. 2020年11月22日閲覧。
- ^ "(L-)BFGS · Optim" (英語). julianlsolvers.github.io. 2024年8月17日閲覧。
関連文献
[編集 ]- Avriel, Mordecai (2003), Nonlinear Programming: Analysis and Methods, Dover Publishing, ISBN 978-0-486-43227-4
- Bonnans, J. Frédéric; Gilbert, J. Charles; Lemaréchal, Claude; Sagastizábal, Claudia A. (2006), "Newtonian Methods", Numerical Optimization: Theoretical and Practical Aspects (Second ed.), Berlin: Springer, pp. 51–66, ISBN 3-540-35445-X
- Fletcher, Roger (1987), Practical Methods of Optimization (2nd ed.), New York: John Wiley & Sons, ISBN 978-0-471-91547-8 , https://archive.org/details/practicalmethods0000flet
- Luenberger, David G.; Ye, Yinyu (2008), Linear and nonlinear programming, International Series in Operations Research & Management Science, 116 (Third ed.), New York: Springer, pp. xiv+546, ISBN 978-0-387-74502-2, MR 2423726
- Kelley, C. T. (1999), Iterative Methods for Optimization, Philadelphia: Society for Industrial and Applied Mathematics, pp. 71–86, ISBN 0-89871-433-8