BFGS法

数理最適化において、ブロイデン・フレッチャー・ゴールドファーブ・シャンノ法 (英: Broyden–Fletcher–Goldfarb–Shanno algorithm)、略してBFGS法は、非制限非線形最適化問題に対する反復的解法の一つである。

BFGS法は山登り法の一種である、準ニュートン法に属しており、︵2回微分可能であることが望ましい︶関数の停留点を探索する。この種の問題の最適性の十分条件は勾配が零となることである。ニュートン法とBFGS法は、最適点の近傍で関数が二次までテイラー展開可能でない場合、収束が保証されない。しかし、BFGS法は関数が滑らかでない場合でもよい性能を発揮することが示されている。

準ニュートン法では、二次微分のヘッセ行列は必ずしも直接計算しなければならないわけではなく、勾配を計算︵あるいは近似計算︶した結果を用いて近似することもできる。準ニュートン法は、一次微分の根を探索する割線法を多次元問題へ一般化したものである。多次元問題では、正割方程式は一意解を与えないため、準ニュートン法はどのように解を拘束するかによってことなる。BFGS法はこの種の手法の中で最も普及しているものの一つである。その他にも、BFGS法を限られたメモリで実行できるようにし、非常に多くの変数を扱えるように変更されたL-BFGS法も広く使われている。単純な矩形拘束を扱うBFGS-B法と呼ばれる変種も存在する。

このアルゴリズムの名前は、チャールズ・ジョージ・ブロイデン、ロジャー・フレッチャー、ドナルド・ゴールドファーブ、デイビッド・シャンノに因む。

原理

最適化問題では、$\mathbf {x}$$\mathbf {x}$を $R^{n}$$R^{n}$上のベクトル、 $f(\mathbf {x} )$$f(\mathbf {x} )$は微分可能なスカラー値関数とし、$\mathbf {x}$$\mathbf {x}$ の取り得る値に制限はないものとして、$f(\mathbf {x} )$$f(\mathbf {x} )$を最小化する。

このアルゴリズムでは、初期推定値 $\mathbf {x} _{0}$$\mathbf {x} _{0}$から初め、各ステージ毎に反復的により良い推定値へと更新していく。

ステージ k における探索方向 p_k  はニュートン方程式に類似した次の方程式を解くことにより得られる。

$B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k}),$$B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k}),$
ここで $B_{k}$$B_{k}$ はヘッセ行列の近似であり、各ステージ毎に反復的に更新する。$\nabla f(\mathbf {x} _{k})$$\nabla f(\mathbf {x} _{k})$ は x_k で計算した勾配である。方向 p_k を得たのち、この方向に向けて直線探索を行い、$f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})$$f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})$ を最小とするようなスカラー $\alpha >0$$\alpha >0$を求める。

$B_{k}$$B_{k}$の更新において課される準ニュートン条件は以下の通りである。

$B_{k+1}(\mathbf {x} _{k+1}-\mathbf {x} _{k})=\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})$$B_{k+1}(\mathbf {x} _{k+1}-\mathbf {x} _{k})=\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})$
$\mathbf {y} _{k}=\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})$$\mathbf {y} _{k}=\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})$ および $\mathbf {s} _{k}=\mathbf {x} _{k+1}-\mathbf {x} _{k}$$\mathbf {s} _{k}=\mathbf {x} _{k+1}-\mathbf {x} _{k}$, $B_{k+1}$$B_{k+1}$ は $B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}$$B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}$曲率条件 $\mathbf {s} _{k}^{\top }\mathbf {y} _{k}>0$$\mathbf {s} _{k}^{\top }\mathbf {y} _{k}>0$ が満たされている必要がある。もし関数が強凸性を持たない場合、この条件を明示的に課す必要がある。

点 x_k+1 における完全なヘッセ行列を計算して B_k+1 とする代わりに、二つの行列を次のように加算して近似ヘッセ行列を更新する。



$B_{k+1}=B_{k}+U_{k}+V_{k}$$B_{k+1}=B_{k}+U_{k}+V_{k}$
U_k および V_k はどちらも階数1の対称行列であるが、これらの和を取ることにより階数2の行列を用いて更新することとなる。BFGS法とDFP法は、以前の行列更新法と階数2の行列を更新に用いる点が異なる。階数1行列を用いる、より単純な手法は対称ランク1法と呼ばれ、正定値性が保証されない。$B_{k+1}$$B_{k+1}$ $B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\top }+\beta \mathbf {v} \mathbf {v} ^{\top }$$B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\top }+\beta \mathbf {v} \mathbf {v} ^{\top }$正割条件 $B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}$$B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}$ $\mathbf {u} =\mathbf {y} _{k}$$\mathbf {u} =\mathbf {y} _{k}$ および $\mathbf {v} =B_{k}\mathbf {s} _{k}$$\mathbf {v} =B_{k}\mathbf {s} _{k}$
$\alpha ={\frac {1}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}}$$\alpha ={\frac {1}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}}$
$\beta =-{\frac {1}{\mathbf {s} _{k}^{T}B_{k}\mathbf {s} _{k}}}$$\beta =-{\frac {1}{\mathbf {s} _{k}^{T}B_{k}\mathbf {s} _{k}}}$
上の $\alpha$$\alpha$ および $\beta$$\beta$ を $B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\top }+\beta \mathbf {v} \mathbf {v} ^{\top }$$B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\top }+\beta \mathbf {v} \mathbf {v} ^{\top }$ に代入して書き下すと、次のようになる。

$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}}}$$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}}}$

アルゴリズム

初期推定値 $\mathbf {x} _{0}$$\mathbf {x} _{0}$ と初期近似ヘッセ行列 $B_{0}$$B_{0}$ から始め、次のステップを反復することにより、$\mathbf {x} _{k}$$\mathbf {x} _{k}$ は解に収束する。

(一)方向 $\mathbf {p} _{k}$$\mathbf {p} _{k}$ を $B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k})$$B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k})$を解くことにより求める。

(二)1次最適化︵直線探索︶によりステップサイズ $\alpha _{k}=\arg \min f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})$$\alpha _{k}=\arg \min f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})$を求める。

(三)$\mathbf {s} _{k}=\alpha _{k}\mathbf {p} _{k}$$\mathbf {s} _{k}=\alpha _{k}\mathbf {p} _{k}$ とし、$\mathbf {x} _{k+1}=\mathbf {x} _{k}+\mathbf {s} _{k}$$\mathbf {x} _{k+1}=\mathbf {x} _{k}+\mathbf {s} _{k}$により推定値を更新する。

(四)$\mathbf {y} _{k}={\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}$$\mathbf {y} _{k}={\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}$
(五)$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}}}$$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(\mathbf {x} )$$f(\mathbf {x} )$ は最小化の対象となる関数である。勾配のノルム $||\nabla f(\mathbf {x} _{k})||$$||\nabla f(\mathbf {x} _{k})||$を監視することにより収束判定を行う。 $B_{0}$$B_{0}$ を $B_{0}=I$$B_{0}=I$ のように初期化した場合、最初のステップは最急降下法と同等となるが、それ以降のステップでは近似ヘッシアン $B_{k}$$B_{k}$の改善に伴って改善される。

このアルゴリズムのステップ1に用いられる $B_{k}$$B_{k}$の逆行列は、シャーマン・モリソンの公式をステップ5に適用することにより次式のように計算することができる。

$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}^{-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}$$B_{k}^{-1}$ が対称行列であり、$\mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}\mathbf {y} _{k}$$\mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}\mathbf {y} _{k}$ および $\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}$$\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}$ がスカラーであることを用いると次のように展開でき、行列を一時記憶領域に記憶する必要なく効率的に計算できる。

$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}}}$$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法

参照文献

• Avriel, Mordecai (2003), Nonlinear Programming: Analysis and Methods, Dover Publishing, ISBN 0-486-43227-0 
• Bonnans, J. Frédéric; Gilbert, J. Charles; Lemaréchal, Claude; Sagastizábal, Claudia A. (2006), Numerical optimization: Theoretical and practical aspects, Universitext (Second revised ed. of translation of 1997 French ed.), Berlin: Springer-Verlag, pp. xiv+490, doi:10.1007/978-3-540-35447-5, ISBN 3-540-35445-X, MR2265882, https://www.springer.com/mathematics/applications/book/978-3-540-35445-1  French Original: Optimisation numérique: Aspects théoriques et pratiques ISBN 3-540-63183-6
• Broyden, C. G. (1970), “The convergence of a class of double-rank minimization algorithms”, Journal of the Institute of Mathematics and Its Applications 6: 76–90, doi:10.1093/imamat/6.1.76 
• Fletcher, R. (1970), “A New Approach to Variable Metric Algorithms”, Computer Journal 13 (3): 317–322, doi:10.1093/comjnl/13.3.317 
• Fletcher, Roger (1987), Practical methods of optimization (2nd ed.), New York: John Wiley & Sons, ISBN 978-0-471-91547-8 
• 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 
• 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, MR2423726 
• Nocedal, Jorge; Wright, Stephen J. (2006), Numerical Optimization (2nd ed.), Berlin, New York: Springer-Verlag, ISBN 978-0-387-30303-1 
• 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, MR274029 
• Shanno, David F.; Kettler, Paul C. (July 1970), “Optimal conditioning of quasi-Newton methods”, Mathematics of Computation 24 (111): 657–664, doi:10.1090/S0025-5718-1970-0274030-6, MR274030 

外部リンク

• Source code of high-precision BFGS