BFGS法

この項目「BFGS法」は翻訳されたばかりのものです。不自然あるいは曖昧な表現などが含まれる可能性があり、このままでは読みづらいかもしれません。（原文：en: Broyden–Fletcher–Goldfarb–Shanno algorithm） 修正、加筆に協力し、現在の表現をより自然な表現にして下さる方を求めています。ノートページや履歴も参照してください。（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における降下方向（英語版） p_kはニュートン方程式に類似した次の方程式を解くことにより得られる。

${\displaystyle B_{k}{\boldsymbol {p}}_{k}=-\nabla f({\boldsymbol {x}}_{k})}$

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

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

${\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}}$とおくと、B_k+1は以下の正割方程式を満たす。

${\displaystyle B_{k+1}{\boldsymbol {s}}_{k}={\boldsymbol {y}}_{k}}$

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

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

${\displaystyle B_{k+1}=B_{k}+U_{k}+V_{k}}$

U_kおよびV_kはどちらも階数1の対称行列であるが、これらの和を取ることにより階数2の対称行列を用いて更新することとなる。対称ランク1（英語版）法と比べ、BFGS法とDFP法はどちらも階数2の行列を更新に用いる点が異なる。より単純な手法である対称ランク1法は階数1の行列を用いて更新を行うが、正定値性が保証されない。B_kの対称性と正定値性を維持するため、更新式は${\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 }}$に代入するとB_k+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}}$から始め、次の各ステップを反復することによりx_kは解に収束する。

1. 降下方向p_kを${\displaystyle B_{k}{\boldsymbol {p}}_{k}=-\nabla f({\boldsymbol {x}}_{k})}$を解くことにより求める。
2. 1次元最適化(直線探索)を行い、前ステップで求めた降下方向に向う許容しうるステップサイズα_kを求める。厳密な直線探索が行われた場合、${\displaystyle \alpha _{k}=\arg \min f({\boldsymbol {x}}_{k}+\alpha {\boldsymbol {p}}_{k})}$ となる。実用上はα_kがウルフ条件を満たすことをもって許容する、非厳密な直線探索で十分なことが多い。
3. ${\displaystyle {\boldsymbol {s}}_{k}=\alpha _{k}{\boldsymbol {p}}_{k}}$とし、${\displaystyle {\boldsymbol {x}}_{k+1}={\boldsymbol {x}}_{k}+{\boldsymbol {s}}_{k}}$により推定解を更新する。
4. ${\displaystyle {\boldsymbol {y}}_{k}={\nabla f({\boldsymbol {x}}_{k+1})-\nabla f({\boldsymbol {x}}_{k})}}$を計算する。
5. ${\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}$のように選んだ場合、最初のステップは最急降下法と等価となるが、以降のステップはB_kがヘッセ行列を推定することにより徐々に改善される。

このアルゴリズムのステップ1はB_kの逆行列を用いて実行されるが、この逆行列はステップ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}}$および s_k^⊤y_kがスカラーであることをもちいて次のように展開でき、一時行列を要せず実行することができる。

${\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]。

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

1. 降下方向p_kを${\displaystyle {\boldsymbol {p}}_{k}=-H_{k}\nabla f({\boldsymbol {x}}_{k})}$により得る。
2. 1次元最適化(直線探索)を行い、前ステップで求めた降下方向に向う許容しうるステップサイズα_kを求める。厳密な直線探索が行われた場合、${\displaystyle \alpha _{k}=\arg \min f({\boldsymbol {x}}_{k}+\alpha {\boldsymbol {p}}_{k})}$ となる。実用上はα_kがウルフ条件を満たすことをもって許容する、非厳密な直線探索で十分なことが多い。
3. ${\displaystyle {\boldsymbol {s}}_{k}=\alpha _{k}{\boldsymbol {p}}_{k}}$とし、${\displaystyle {\boldsymbol {x}}_{k+1}={\boldsymbol {x}}_{k}+{\boldsymbol {s}}_{k}}$により推定解を更新する。
4. ${\displaystyle {\boldsymbol {y}}_{k}={\nabla f({\boldsymbol {x}}_{k+1})-\nabla f({\boldsymbol {x}}_{k})}}$を計算する。
5. ${\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更新公式は曲率s_k^⊤y_kが常に正であり、ゼロから離れた下界があることに強く依拠している。この条件は凸な対称関数においてウルフ条件を用いた直線探索を用いる場合は満たされるが、実際の問題(たとえば逐次二次計画法)では負やほぼゼロの曲率があらわれることがしばしばある。このようなことは非凸関数を対象とする場合や直線探索ではなく信頼領域アプローチをとった場合におこることがある。この場合、BFGS法は誤った値をあたえることがある。

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

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

• 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法を用いて陰解を求めている。

• BHHH法（英語版）
• DFP法
• 最急降下法
• L-BFGS
• レーベンバーグ・マルカート法
• ネルダー–ミード法
• パターン探索法（英語版）
• 準ニュートン法
• 対称ランク1法（英語版）

• 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, MR2423726 
• Kelley, C. T. (1999), Iterative Methods for Optimization, Philadelphia: Society for Industrial and Applied Mathematics, pp. 71–86, ISBN 0-89871-433-8 

• Source code of high-precision BFGS

表 話 編 歴 数理最適化 • 最適化問題 : メソッド • ヒューリスティック
非線形(無制約)	… 関数　 黄金分割探索 直線探索 ネルダー–ミード法 連続放物線補間（英語版） 勾配法 収束性 信頼領域 ウルフ条件 準ニュートン法 BFGS法 ブロイデン法 L-BFGS（英語版） DFP法 対称ランク1法（英語版） その他の求解法 ガウス・ニュートン法 最急降下法 レーベンバーグ・マルカート法 共役勾配法（非線形共役勾配法） 切り捨てニュートン法（英語版） ドッグレッグ法 … ヘッセ行列 最適化におけるニュートン法（英語版）
非線形(制約付き)	一般 バリア関数 ペナルティ関数法（英語版） 微分可能 ラグランジュの未定乗数法 逐次二次計画法 連続線形計画（英語版）
凸最適化	凸縮小化 切除平面法（英語版、デンマーク語版、ドイツ語版、スペイン語版） 簡約勾配法 劣勾配法（英語版） 線型 および 二次 内点法 カチヤンの楕円体法 カーマーカーの投影アルゴリズム ベイズ-交換 単体法 改訂単体法（英語版） 十文字法（英語版） レムケの主ピボット操作法（英語版） 有効制約法（英語版）
組合せ最適化	系列範例 (Paradigms) 近似アルゴリズム 動的計画法 貪欲法 整数計画問題(分枝限定法 若しくは 分枝カット法) グラフ理論 最小全域木 ベルマン–フォード法 ブルーフカ法 ダイクストラ法 ワーシャル–フロイド法 ジョンソン法（英語版） クラスカル法 最大フロー問題 Dinic法（英語版） エドモンズ・カープ フォード・ファルカーソン プッシュリラベル最大流アルゴリズム（英語版）
メタヒューリスティクス	進化的アルゴリズム（進化戦略） 山登り法 局所探索法 焼きなまし法 タブーサーチ ベイスンホッピング法
カテゴリ（最適化 • アルゴリズム） • ソフトウェア（英語版）