BFGS法

数理最適化において、**ブロイデン・フレッチャー・ゴールドファーブ・シャンノ法**（ブロイデン・フレッチャー・ゴールドファーブ・シャンノほう、）、略して**BFGS法**とは、無制約非線形最適化問題に対する反復的解法の一つである。関連の深いDFP法と同様、BFGS法は勾配のプレコンディショニングを曲率の情報を用いて行うことにより降下方向を決定する。その際、損失関数のヘッセ行列の推定値を勾配（またはその推定値）のみを用いて（一般化）割線法により漸進的に改善する。

BFGS法における曲率行列の更新には逆行列の評価を要さないため、は\mathcal{O}(n^{2})に留まり、ニュートン法の\mathcal{O}(n^{3})よりも高速である。L-BFGS法もよく用いられ、メモリ使用量を限定できるため、多変数（e.g. >1000）問題に対する解法に適している。BFGS-B法はシンプルなボックス拘束を扱える。

このアルゴリズムの名前は、、ロジャー・フレッチャー、、に因む。

## 理論的根拠
\boldsymbol{x}を\mathbb{R}^n上のベクトル、f(\boldsymbol{x})を微分可能なスカラー値関数とし、\boldsymbol{x}の取り得る値に制限はないものとして、f(\boldsymbol{x})を最小化する最適化問題を考える。

BFGS法は初期推定値\boldsymbol{x}_0から始め、各ステージ毎に反復的により良い推定値へと更新していく。

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

: B_k \boldsymbol{p}_k = -\nabla f(\boldsymbol{x}_k)

ここではにおけるヘッセ行列の推定値であり、各ステージごとににおける目的関数の勾配\nabla f(\boldsymbol{x}_k)を用いて反復的に更新される。降下方向を得たのち、この方向に向けて直線探索を行い、f(\boldsymbol{x}_k + \gamma\boldsymbol{p}_k)を最小とするようなスカラーを求め、次の点を決定する。

の更新においては、以下の式であらわされる準ニュートン条件が課せられる。

: B_{k+1} (\boldsymbol{x}_{k+1} - \boldsymbol{x}_k) = \nabla f(\boldsymbol{x}_{k+1}) - \nabla f(\boldsymbol{x}_k)

ここで\boldsymbol{y}_k = \nabla f(\boldsymbol{x}_{k+1}) - \nabla f(\boldsymbol{x}_k)および\boldsymbol{s}_k = \boldsymbol{x}_{k+1} - \boldsymbol{x}_kとおくと、は以下の正割方程式を満たす。

: B_{k+1} \boldsymbol{s}_k = \boldsymbol{y}_k

が正定値行列であるためには曲率条件が満たされる必要がある。この条件は正割方程式に左からをかけることにより検証できる。目的関数が強凸関数でない場合、この条件は明示的に課す必要があり、これはたとえばを決定する際にウルフ条件を満たす点を選べばよい。

点におけるヘッセ行列を全て計算するかわりに、ステージ**'における推定値に次のように2つの行列を足すことによりを計算する。

: <math>B_{k+1} = B_k + U_k + V_k</math>

およびはどちらも[[行列の階数|階数]]1の[[対称行列]]であるが、これらの和を取ることにより階数2の対称行列を用いて更新することとなる。[[対称ランクワン法]]と比べ、BFGS法と[[DFP法]]はどちらも階数2の行列を更新に用いる点が異なる。より単純な手法である対称ランクワン法は階数1の行列を用いて更新を行うが、[[行列の定値性|正定値性]]が保証されない。の対称性と正定値性を維持するため、更新式は<math>B_{k+1} = B_k + \alpha\boldsymbol{u}\boldsymbol{u}^\top + \beta\boldsymbol{v}\boldsymbol{v}^\top</math>のように選ぶ。正割条件<math>B_{k+1} \boldsymbol{s}_k = \boldsymbol{y}_k </math>を課すと、<math>\boldsymbol{u} = \boldsymbol{y}_k</math>および<math>\boldsymbol{v} = B_k \boldsymbol{s}_k</math>として以下を得る。

: <math>\alpha = \frac{1}{\boldsymbol{y}_k^\top \boldsymbol{s}_k}</math>
: <math>\beta = -\frac{1}{\boldsymbol{s}_k^\top  B_k \boldsymbol{s}_k}</math>

最後に、およびを<math>B_{k+1} = B_k + \alpha\boldsymbol{u}\boldsymbol{u}^\top + \beta\boldsymbol{v}\boldsymbol{v}^\top</math>に代入するとの更新式は以下のように書ける。

: <math>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}</math>

== アルゴリズム ==
非線形関数<math>f : \mathbb{R}^n \to \mathbb{R}</math>を対象とした無制約最適化問題<math display="block">
\begin{aligned}
   \underset{\boldsymbol{x} \in \mathbb{R}^n}{\text{minimize}}\quad & f(\boldsymbol{x})
\end{aligned}
</math>を考える。

初期推定解<math>\boldsymbol{x}_0 \in \mathbb{R}^n</math>および初期推定ヘッセ行列<math>B_0 \in \mathbb{R}^{n\times n}</math>から始め、次の各ステップを反復することによりは解に収束する。

# 降下方向を<math>B_k \boldsymbol{p}_k = -\nabla f(\boldsymbol{x}_k)</math>を解くことにより求める。
# 1次元最適化（[[直線探索]]）を行い、前ステップで求めた降下方向に向う許容しうるステップサイズを求める。厳密な直線探索が行われた場合、<math>\alpha_k=\arg \min f(\boldsymbol{x}_k+\alpha\boldsymbol{p}_k)</math> となる。実用上はがウルフ条件を満たすことをもって許容する、非厳密な直線探索で十分なことが多い。
# <math> \boldsymbol{s}_k = \alpha_k \boldsymbol{p}_k</math>とし、<math>\boldsymbol{x}_{k+1} = \boldsymbol{x}_k + \boldsymbol{s}_k</math>により推定解を更新する。
# <math>\boldsymbol{y}_k = {\nabla f(\boldsymbol{x}_{k+1}) - \nabla f(\boldsymbol{x}_k)}</math>を計算する。
# <math>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}</math>により推定ヘッセ行列を更新する。

何らかの基準値のもと、勾配の[[ノルム]]が<math>||\nabla f(\boldsymbol{x}_k)|| \leq \varepsilon</math>を満たしたとき解が収束したものとみなしアルゴリズムを終了する。

<math>B_0 = I</math>のように選んだ場合、最初のステップは[[最急降下法]]と等価となるが、以降のステップはがヘッセ行列を推定することにより徐々に改善される。

このアルゴリズムのステップ1はの逆行列を用いて実行されるが、この逆行列はステップ5で を用いることにより次のように効率的に求めることができる。

: <math>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}</math>

この計算は <math>B_k^{-1}</math>が対称行列であり、 <math>\boldsymbol{y}_k^{\top} B_k^{-1} \boldsymbol{y}_k</math>および がスカラーであることを用いて次のように展開でき、一時行列を要せず実行することができる。

: <math>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}.</math>

したがって、逆行列を求めるための計算を一切することなく、ヘッセ行列の逆行列<math>H_k \overset{\operatorname{def}}{=} B_k^{-1}</math>そのものを推定することが可能である。

初期推定解、ヘッセ行列の**逆行列の'''推定値から始め、次の各ステップを反復することによりは解へと収束する。

# 降下方向を\boldsymbol{p}_k = -H_k \nabla f(\boldsymbol{x}_k)により得る。
# 1次元最適化（直線探索）を行い、前ステップで求めた降下方向に向う許容しうるステップサイズを求める。厳密な直線探索が行われた場合、\alpha_k=\arg \min f(\boldsymbol{x}_k+\alpha\boldsymbol{p}_k) となる。実用上はがウルフ条件を満たすことをもって許容する、非厳密な直線探索で十分なことが多い。
#  \boldsymbol{s}_k = \alpha_k \boldsymbol{p}_kとし、\boldsymbol{x}_{k+1} = \boldsymbol{x}_k + \boldsymbol{s}_kにより推定解を更新する。
# \boldsymbol{y}_k = {\nabla f(\boldsymbol{x}_{k+1}) - \nabla f(\boldsymbol{x}_k)}を計算する。
# 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近似は真のヘッセ行列に収束しない場合がある。

## 発展
BFGS更新公式は曲率が常に正であり、ゼロから離れた下界があることに強く依拠している。この条件は凸な対称関数においてウルフ条件を用いた直線探索を用いる場合は満たされるが、実際の問題（たとえば逐次二次計画法）では負やほぼゼロの曲率があらわれることがしばしば発生する。このようなことは非凸関数を対象とする場合や直線探索ではなく信頼領域アプローチをとった場合に生じるおそれがある。この場合、BFGS法は誤った値をあたえることがある。

このような場合には、減衰BFGS更新などと呼ばれる、および/またはを修正して頑健にした更新式が用いられることがある。

## 実装
オープンソースの実装として有名なものは以下のようなものがあげられる。

* ALGLIBはC++およびC#用のBFGSおよびL-BFGS法を実装する。
* GNU Octaveの`fsolve関数は`信頼領域を用いた一種のBFGS法を用いる。
* GSLはgsl_multimin_fdfminimizer_vector_bfgs2関数としてBFGSを実装している。
* R言語では、、BFGS法（および矩形拘束を扱えるL-BFGS-B法）が基本関数optim()のオプションとして実装されている。
* SciPyでは、scipy.optimize.fmin_bfgs関数がBFGS法を実装している。パラメータLにとても大きな数を指定することにより、なんらかのL-BFGS法を実行することもできる。
* Juliaでは、Optim.jlパッケージにBFGSおよびL-BFGSが実装されている。

プロプライエタリな実装としては以下のようなものがあげられる。

* 大規模非線形最適化ソフトウェアArtelys KnitroはBFGS法およびL-BFGS法の両方を実装する。
* MATLAB Optimization Toolboxでは、fminunc関数がBFGS法を3次直線探索と組み合わせたアルゴリズムを「中規模スケール」の問題向けに実装している。
* MathematicaにはBFGS法が含まれる。
* LS-DYNAもBFGS法を用いて陰解を求めている。

## 脚注

## 関連文献
* 
* 
* 
* 
* 

## 関連項目
* BHHH法
* DFP法
* 最急降下法
* L-BFGS
* レーベンバーグ・マルカート法
* ネルダー–ミード法
* 
* 準ニュートン法
* 対称ランクワン法

## 外部リンク
* Source code of high-precision BFGS