BFGS法

数理最適化において、ブロイデン・フレッチャー・ゴールドファーブ・シャンノ法 (英: Broyden–Fletcher–Goldfarb–Shanno algorithm)、略してBFGS法は、非制限非線形最適化問題に対する反復的解法の一つである。 BFGS法は山登り法の一種である、準ニュートン法に属しており、（2回微分可能であることが望ましい）関数の停留点を探索する。この種の問題の最適性の十分条件は勾配が零となることである。ニュートン法とBFGS法は、最適点（英語版）の近傍で関数が二次までテイラー展開可能でない場合、収束が保証されない。しかし、BFGS法は関数が滑らかでない場合でもよい性能を発揮することが示されている。 準ニュートン法では、二次微分のヘッセ行列は必ずしも直接計算しなければならないわけではなく、勾配を計算（あるいは近似計算）した結果を用いて近似することもできる。準ニュートン法は、一次微分の根を探索する割線法を多次元問題へ一般化したものである。多次元問題では、正割方程式は一意解を与えないため、準ニュートン法はどのように解を拘束するかによってことなる。BFGS法はこの種の手法の中で最も普及しているものの一つである。その他にも、BFGS法を限られたメモリで実行できるようにし、非常に多くの変数を扱えるように変更されたL-BFGS法も広く使われている。単純な矩形拘束を扱うBFGS-B法と呼ばれる変種も存在する。 このアルゴリズムの名前は、チャールズ・ジョージ・ブロイデン（英語版）、ロジャー・フレッチャー（英語版）、ドナルド・ゴールドファーブ（英語版）、デイビッド・シャンノ（英語版）に因む。 最適化問題では、 x {\displaystyle \mathbf {x} } を  R n {\displaystyle R^{n}} 上のベクトル、  f ( x ) {\displaystyle f(\mathbf {x} )} は微分可能なスカラー値関数とし、 x {\displaystyle \mathbf {x} }  の取り得る値に制限はないものとして、 f ( x ) {\displaystyle f(\mathbf {x} )} を最小化する。 このアルゴリズムでは、初期推定値  x 0 {\displaystyle \mathbf {x} _{0}} から初め、各ステージ毎に反復的により良い推定値へと更新していく。 ステージ k における探索方向 pk  はニュートン方程式に類似した次の方程式を解くことにより得られる。 B k p k = − ∇ f ( x k ) , {\displaystyle B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k}),} ここで  B k {\displaystyle B_{k}}  はヘッセ行列の近似であり、各ステージ毎に反復的に更新する。 ∇ f ( x k ) {\displaystyle \nabla f(\mathbf {x} _{k})}  は xk で計算した勾配である。方向 pk を得たのち、この方向に向けて直線探索を行い、 f ( x k + α p k ) {\displaystyle f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})}  を最小とするようなスカラー  α > 0 {\displaystyle \alpha >0} を求める。 B k {\displaystyle B_{k}} の更新において課される準ニュートン条件は以下の通りである。 B k + 1 ( x k + 1 − x k ) = ∇ f ( x k + 1 ) − ∇ f ( x k ) {\displaystyle B_{k+1}(\mathbf {x} _{k+1}-\mathbf {x} _{k})=\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})} y k = ∇ f ( x k + 1 ) − ∇ f ( x k ) {\displaystyle \mathbf {y} _{k}=\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}  および  s k = x k + 1 − x k {\displaystyle \mathbf {s} _{k}=\mathbf {x} _{k+1}-\mathbf {x} _{k}} , B k + 1 {\displaystyle B_{k+1}}  は  B k + 1 s k = y k {\displaystyle B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}} 曲率条件  s k ⊤ y k > 0 {\displaystyle \mathbf {s} _{k}^{\top }\mathbf {y} _{k}>0}  が満たされている必要がある。もし関数が強凸性を持たない場合、この条件を明示的に課す必要がある。 点 xk+1 における完全なヘッセ行列を計算して Bk+1 とする代わりに、二つの行列を次のように加算して近似ヘッセ行列を更新する。 ￼ B k + 1 = B k + U k + V k {\displaystyle B_{k+1}=B_{k}+U_{k}+V_{k}} Uk および Vk はどちらも階数1の対称行列であるが、これらの和を取ることにより階数2の行列を用いて更新することとなる。BFGS法とDFP法は、以前の行列更新法と階数2の行列を更新に用いる点が異なる。階数1行列を用いる、より単純な手法は対称ランク1（英語版）法と呼ばれ、正定値性が保証されない。 B k + 1 {\displaystyle B_{k+1}}   B k + 1 = B k + α u u ⊤ + β v v ⊤ {\displaystyle B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\top }+\beta \mathbf {v} \mathbf {v} ^{\top }} 正割条件  B k + 1 s k = y k {\displaystyle B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}}   u = y k {\displaystyle \mathbf {u} =\mathbf {y} _{k}} および  v = B k s k {\displaystyle \mathbf {v} =B_{k}\mathbf {s} _{k}} α = 1 y k T s k {\displaystyle \alpha ={\frac {1}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}}} β = − 1 s k T B k s k {\displaystyle \beta =-{\frac {1}{\mathbf {s} _{k}^{T}B_{k}\mathbf {s} _{k}}}} 上の  α {\displaystyle \alpha } および  β {\displaystyle \beta }  を  B k + 1 = B k + α u u ⊤ + β v v ⊤ {\displaystyle B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\top }+\beta \mathbf {v} \mathbf {v} ^{\top }}  に代入して書き下すと、次のようになる。 B k + 1 = B k + y k y k T y k T s k − B k s k s k T B k T s k T B k s k {\displaystyle B_{k+1}=B_{k}+{\frac {\mathbf {y} _{k}\mathbf {y} _{k}^{\mathrm {T} }}{\mathbf {y} _{k}^{\mathrm {T} }\mathbf {s} _{k}}}-{\frac {B_{k}\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} }B_{k}^{\mathrm {T} }}{\mathbf {s} _{k}^{\mathrm {T} }B_{k}\mathbf {s} _{k}}}} 初期推定値  x 0 {\displaystyle \mathbf {x} _{0}}  と初期近似ヘッセ行列  B 0 {\displaystyle B_{0}}  から始め、次のステップを反復することにより、 x k {\displaystyle \mathbf {x} _{k}}  は解に収束する。 方向  p k {\displaystyle \mathbf {p} _{k}}  を  B k p k = − ∇ f ( x k ) {\displaystyle B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k})} を解くことにより求める。 1次最適化（直線探索）によりステップサイズ  α k = arg ⁡ min f ( x k + α p k ) {\displaystyle \alpha _{k}=\arg \min f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})} を求める。 s k = α k p k {\displaystyle \mathbf {s} _{k}=\alpha _{k}\mathbf {p} _{k}}  とし、 x k + 1 = x k + s k {\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\mathbf {s} _{k}} により推定値を更新する。 y k = ∇ f ( x k + 1 ) − ∇ f ( x k ) {\displaystyle \mathbf {y} _{k}={\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}} B k + 1 = B k + y k y k T y k T s k − B k s k s k T B k s k T B k s k {\displaystyle B_{k+1}=B_{k}+{\frac {\mathbf {y} _{k}\mathbf {y} _{k}^{\mathrm {T} }}{\mathbf {y} _{k}^{\mathrm {T} }\mathbf {s} _{k}}}-{\frac {B_{k}\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} }B_{k}}{\mathbf {s} _{k}^{\mathrm {T} }B_{k}\mathbf {s} _{k}}}} f ( x ) {\displaystyle f(\mathbf {x} )}  は最小化の対象となる関数である。勾配のノルム  | | ∇ f ( x k ) | | {\displaystyle ||\nabla f(\mathbf {x} _{k})||} を監視することにより収束判定を行う。  B 0 {\displaystyle B_{0}}  を  B 0 = I {\displaystyle B_{0}=I}  のように初期化した場合、最初のステップは最急降下法と同等となるが、それ以降のステップでは近似ヘッシアン B k {\displaystyle B_{k}} の改善に伴って改善される。 このアルゴリズムのステップ1に用いられる  B k {\displaystyle B_{k}} の逆行列は、シャーマン・モリソンの公式（英語版）をステップ5に適用することにより次式のように計算することができる。 B k + 1 − 1 = ( I − s k y k T y k T s k ) B k − 1 ( I − y k s k T y k T s k ) + s k s k T y k T s k {\displaystyle B_{k+1}^{-1}=\left(I-{\frac {s_{k}y_{k}^{T}}{y_{k}^{T}s_{k}}}\right)B_{k}^{-1}\left(I-{\frac {y_{k}s_{k}^{T}}{y_{k}^{T}s_{k}}}\right)+{\frac {s_{k}s_{k}^{T}}{y_{k}^{T}s_{k}}}} この計算は、 B k − 1 {\displaystyle B_{k}^{-1}}  が対称行列であり、 y k T B k − 1 y k {\displaystyle \mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}\mathbf {y} _{k}}  および  s k T y k {\displaystyle \mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}}  がスカラーであることを用いると次のように展開でき、行列を一時記憶領域に記憶する必要なく効率的に計算できる。 B k + 1 − 1 = B k − 1 + ( s k T y k + y k T B k − 1 y k ) ( s k s k T ) ( s k T y k ) 2 − B k − 1 y k s k T + s k y k T B k − 1 s k T y k {\displaystyle B_{k+1}^{-1}=B_{k}^{-1}+{\frac {(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}+\mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}\mathbf {y} _{k})(\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} })}{(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k})^{2}}}-{\frac {B_{k}^{-1}\mathbf {y} _{k}\mathbf {s} _{k}^{\mathrm {T} }+\mathbf {s} _{k}\mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}}{\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}}}} 最尤推定やベイズ推定などの統計推定問題においては、最終的なヘッセ行列の逆行列を用いて解の信頼区間もしくは確信区間を推定することができる。しかし、これらの量は正確には真のヘッセ行列により定義されるものであり、BFGS近似は真のヘッセ行列に収束しない場合がある。 GSL は gsl_multimin_fdfminimizer_vector_bfgs2 でBFGS法を実装している。Ceres Solver はBFGS法とL-BFGS法の両方を実装している。 SciPyでは、scipy.optimize.fmin_bfgs 関数がBFGS法を実装している。パラメータ L にとても大きな数を指定することにより、L-BFGS法を実行することもできる。 Octave は double-dogleg 近似を用いたBFGS法を  cubic line search に用いている。 R言語では、BFGS法（および矩形拘束を扱えるL-BFGS-B法）が基本関数 optim() のオプションとして実装されている。 MATLAB Optimization Toolbox（英語版） では、fminunc 関数が問題サイズを「中程度」に指定した場合にBFGS法を利用する。 C++ で実装され、高精度算術パッケージ ARPREC に統合されている高精度算術版BFGS法 (pBFGS) は（丸め誤差などの）数値的不安定性に対して頑強である。 他にも、C++ によるBFGS（およびL-BFGS, L-BFGS-B, CG, ニュートン法）の実装としては、MIT License で公開されている Eigen（英語版） ライブラリが存在する。 オープンソースの C 実装としては、Gnu Regression, Econometrics and Time-series Library (gretl) がBFGSおよび L-BFGS法を含んでいる。 準ニュートン法 DFP法 L-BFGS法 最急降下法 ネルダー–ミード法 レーベンバーグ・マルカート法 パターン探索法（英語版） BHHH法（英語版） 対称ランク1法（英語版）