BFGS法

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

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

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

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

原理

最適化問題では、 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}}  は解に収束する。

  1. 方向  p k {\displaystyle \mathbf {p} _{k}}  を  B k p k = − ∇ f ( x k ) {\displaystyle B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k})} を解くことにより求める。
  2. 1次最適化（直線探索）によりステップサイズ  α k = arg ⁡ min f ( x k + α p k ) {\displaystyle \alpha _{k}=\arg \min f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})} を求める。
  3. 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}} により推定値を更新する。
  4. y k = ∇ f ( x k + 1 ) − ∇ f ( x k ) {\displaystyle \mathbf {y} _{k}={\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}}
  5. 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法を含んでいる。

関連項目

出典

  1. ^ Fletcher, Roger (1987), Practical methods of optimization (2nd ed.), New York: John Wiley & Sons, ISBN 978-0-471-91547-8 
  2. ^ Lewis, Adrian S.; Overton, Michael (2009), Nonsmooth optimization via BFGS, http://www.cs.nyu.edu/~overton/papers/pdffiles/bfgs_inexactLS.pdf 
  3. ^ Nocedal & Wright (2006), page 24
  4. ^ 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 

参照文献

  • 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
  • 表
  • 話
  • 編
  • 歴
非線形(無制約)
… 関数　
勾配法
収束性
準ニュートン法
その他の求解法
… ヘッセ行列
  • 最適化におけるニュートン法（英語版）
The graph of a strictly concave quadratic function is shown in blue, with its unique maximum shown as a red dot. Below the graph appears the contours of the function: The level sets are nested ellipses.
Optimization computes maxima and minima.
非線形(制約付き)
一般
微分可能
凸最適化
凸縮小化
  • 切除平面法（英語版、デンマーク語版、ドイツ語版、スペイン語版）
  • 簡約勾配法
  • 劣勾配法（英語版）
線型 および
二次
内点法
ベイズ-交換
  • 単体法
  • 改訂シンプレックス法（英語版）
  • 十字法（英語版）
  • レムケの主ピボット操作法（英語版）
組合せ最適化
系列範例
(Paradigms)
グラフ理論
最小全域木
最大フロー問題
メタヒューリスティクス
  • カテゴリ（最適化 • アルゴリズム） • ソフトウェア（英語版）