BFGS公式による準ニュートン（Quasi-Newton）法はいまのところ実用上最も優れているといえる最適化手法だ。超1次収束かつ各反復における計算オーダが変数の数の2乗で済むため、変数の2乗が許容される規模の問題である限り実用的にベストな解法として使われる。
まずは普通のニュートン法から。
ニュートン法
変数\(\boldsymbol{x}\in\mathbb{R}^n\)、目的関数\(f:\mathbb{R}^n\to\mathbb{R}\)を制約無しで最小化する問題を対象とした反復法では、初期値\(\boldsymbol{x}_0\)からスタートして次のルールで解の候補を更新していく。
\[ \boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\alpha_k\boldsymbol{p}_k\]
\(\boldsymbol{p}_k\)は\(k\)ステップにおける探索方向、\(\alpha_k\)はステップサイズで、直線探索によって決定される。\(k\)ステップにおける関数値を\(f_k\equiv f(\boldsymbol{x}_k)\)とする。最急降下法では\(\boldsymbol{p}_k=-\nabla f_k\)を用いる。\(\nabla f_k\in\mathbb{R}^n\)は各要素が
\[ (\nabla f_k)_i =\left.\frac{\partial f(\boldsymbol{x})}{\partial x_i}\right|_{\boldsymbol{x}=\boldsymbol{x}_k} \quad(i=1,\dots,n)\]
で与えられる勾配ベクトル。ニュートン法ではこの勾配に加えて各要素が以下で与えられるヘッセ行列を用いる。
\[ (\nabla^2 f_k)_{ij}=\left. \frac{\partial^2 f(\boldsymbol{x})}{\partial x_i\partial x_j}\right|_{\boldsymbol{x}=\boldsymbol{x}_k} \quad(i,j=1,\dots,n)\]
\(\nabla^2 f_k\)は\(n\times n\)の対称行列（そうなるような関数だとする）。ヘッセ行列は目的関数の2階微分を全通り並べた行列。ニュートン法はこれを用いて\(\boldsymbol{p}_k=-(\nabla^2f_k)^{-1}\nabla f_k\)という探索方向を用いる。明らかなようにこういう純粋なニュートン法が使えるのは\(\nabla ^2 f_k\)が正則で逆行列が存在するときだけ、対称行列においては正定値であれば十分だ。
ニュートン法における探索方向\(\boldsymbol{p}_k=-(\nabla ^2 f_k)^{-1}\nabla f_k\)を導く。これには、まず関数\(f\)を\(\boldsymbol{x}_k\)の周りで2次の項までテイラー展開する。剰余項を省略すると、
\[ \begin{align} f(\boldsymbol{x}_k+\boldsymbol{p}) &=f_k+(\nabla f_k)^\top\boldsymbol{p}+\frac12 \boldsymbol{p}^\top(\nabla^2f_k)\boldsymbol{p}\\ &\equiv m_k(\boldsymbol{p}) \end{align}\]
これは\(\boldsymbol{p}\)の2次形式になっていてこれを\(m_k(\boldsymbol{p})\)とおいた。これを最小にするような\(\boldsymbol{p}\)を探索方向\(\boldsymbol{p}_k\)として使う。2次関数における平方完成とのアナロジーから
\[ f_k+(\nabla f_k)^\top\boldsymbol{p}+\frac12 \boldsymbol{p}^\top(\nabla^2f_k)\boldsymbol{p} =\frac12(\boldsymbol{p}-\boldsymbol{q})^\top(\nabla^2f_k) (\boldsymbol{p}-\boldsymbol{q})+C\]
と置いて両辺が一致するように\(\boldsymbol{q},C\)を決定する。右辺を展開すると
\[ \frac12\boldsymbol{p}^\top(\nabla^2f_k)\boldsymbol{p} -\boldsymbol{q}^\top(\nabla^2 f_k)\boldsymbol{p} +\frac12\boldsymbol{q}^\top(\nabla^2f_k)\boldsymbol{q}+C\]
となるので\(\boldsymbol{p}\)の1次と0次の項の係数を比較して
\[ \begin{align} (\nabla f_k)^\top&=-\boldsymbol{q}^\top(\nabla^2 f_k)\\ f_k&=\frac12\boldsymbol{q}^\top(\nabla^2f_k)\boldsymbol{q}+C \end{align}\]
がわかるので第1式から\(\boldsymbol{q}=-(\nabla^2 f_k)^{-1}\nabla f_k\)がわかり、第2式にこれを代入して移項することで\(C\)がわかる。長くなるので\(\boldsymbol{q}\)はそのままにすると結局
\[ \begin{gather} f_k+(\nabla f_k)^\top\boldsymbol{p}+\frac12 \boldsymbol{p}^\top(\nabla^2f_k)\boldsymbol{p} =\frac12(\boldsymbol{p}-\boldsymbol{q})^\top(\nabla^2f_k) (\boldsymbol{p}-\boldsymbol{q})+ f_k-\frac12\boldsymbol{q}^\top(\nabla^2f_k)\boldsymbol{q}\\ \boldsymbol{q}=-(\nabla^2 f_k)^{-1}\nabla f_k \end{gather}\]
という恒等式が得られた。明らかに\(\boldsymbol{p}=\boldsymbol{q}\)のときに最小値\(f_k-\boldsymbol{q}^\top(\nabla^2f_k)\boldsymbol{q}/2\)をとる。ここでは準ニュートン法のために\(\nabla m_k\)を考えるもう1つの方法も見ておく。
\[ m_k(\boldsymbol{p})=f_k+\nabla f_k^\top\boldsymbol{p}+\frac12\boldsymbol{p}^\top(\nabla^2 f_k)\boldsymbol{p}\]
においては、考えている変数は\(\boldsymbol{p}\)なので\(f_k\)に関するところは全て定数と見なす。\(\partial/\partial x_i\)が\(\partial /\partial p_i\)になるという見かけ上の変数の違いに注意すると
\[ \nabla m_k(\boldsymbol{p})=\nabla f_k+(\nabla^2f_k)\boldsymbol{p}\]
がわかる。右辺をゼロと置くことで直ちに平方完成と同じ結果が得られる。また、\(m_k(\boldsymbol{p})\)に関して
\[ \begin{align}m_k(\boldsymbol{0})&=f_k\\\nabla m_k(\boldsymbol{0})&=\nabla f_k\end{align}\]
が成り立つことも簡単に確認できる。これは準ニュートン法の導出で使う。
このようにして求めた探索方向を用いてステップサイズ\(\alpha_k\)を直線探索によって決定すれば1回分の反復の更新式が得られる。ニュートン法では次のように関数\(f\)を\(\boldsymbol{x}_k\)のまわりで2次近似している。
\[ \begin{align}f(\boldsymbol{x}_{k+1})&=f(\boldsymbol{x}_k+\alpha_k\boldsymbol{p}_k)\\&=f_k+\alpha_k\nabla f_k^\top\boldsymbol{p}_k+\frac{\alpha_k^2}{2}\boldsymbol{p}_k^\top(\nabla^2 f_k)\boldsymbol{p}_k+o(\alpha_k^2)\\&=m_k(\alpha_k\boldsymbol{p}_k)+o(\alpha_k^2)\end{align}\]
\(\boldsymbol{p}_k\)は\(m_k(\boldsymbol{p})\)を最小にするように取っているので、\(\alpha_k=1\)を基準に\(\alpha_k\)を小さくしていけば後半の項は無視できるほど小さくなる。理論上は\(\alpha_k\)が十分小さければ対象とする関数のいかんにかかわらず（微分はできるとするが）この近似はいい精度で成立し、1次近似ほど\(\alpha_k\)を小さくしなくてもよいことが期待できる。このような直観からニュートン法は2次収束することが知られているが以下の問題点がある。
- 目的関数の2階微分まで与える必要がある（ヘッセ行列の計算）
- ヘッセ行列の逆行列を計算するのに\(n^3\)のコストがかかる。
準ニュートン法では目的関数の1階微分までしか使わない。そのなかでもBFGS公式による準ニュートン法では反復ごとに最大\(n^2\)の計算で済む。
BFGS公式による準ニュートン法
準ニュートン法ではヘッセ行列を逐次計算によって近似したものをつかう。ステップ\(k\)におけるヘッセ行列の推定を\(B_k\)として\(m_k(\boldsymbol{p})\)を次のように修正する。
\[ m_k(\boldsymbol{p})=f_k+\nabla f_k^\top\boldsymbol{p}+\frac12\boldsymbol{p}^\top B_k\boldsymbol{p}\]
これを最小にするような\(\boldsymbol{p}\)として\(\boldsymbol{p}_k=-B_k^{-1}\nabla f_k\)を選ぶことはニュートン法と同じ。\(B_k\)から\(B_{k+1}\)を決めるルールを求めるために\(m_{k+1}\)を考える。
\[ m_{k+1}(\boldsymbol{p})=f_{k+1}+\nabla f_{k+1}^\top\boldsymbol{p}+\frac12\boldsymbol{p}^\top B_{k+1}\boldsymbol{p}\]
\(m_{k+1}\)についても以下が成り立つ。
\[ \begin{align}m_{k+1}(\boldsymbol{0})&=f_{k+1}\\\nabla m_{k+1}(\boldsymbol{0})&=\nabla f_{k+1}\end{align}\]
これは当然としてさらに\(m_k,m_{k+1}\)のもともとの定義に立ち返れば次がわかる。
\[ \begin{gather}f(\boldsymbol{x}_{k+1})=f(\boldsymbol{x}_k+\alpha_k\boldsymbol{p}_k)=m_k(\alpha_k\boldsymbol{p}_k)\\f(\boldsymbol{x}_k)=f(\boldsymbol{x}_{k+1}-\alpha_k\boldsymbol{p}_k)=m_{k+1}(-\alpha_k\boldsymbol{p}_k)\end{gather}\]
これは\(\boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\alpha_k\boldsymbol{p}_k\)というルールを素朴に使っただけだ。両辺微分を取って右辺を計算すると、
\[ \begin{align}\nabla f_{k+1}&=\nabla m_k(\alpha_k\boldsymbol{p}_k)=\nabla f_k+\alpha_kB_k\boldsymbol{p}_k\\\nabla f_k&=\nabla m_{k+1}(-\alpha_k\boldsymbol{p}_k)=\nabla f_{k+1}-\alpha_kB_{k+1}\boldsymbol{p}_k\end{align}\]
第2式に注目すると
\[ \begin{align}B_{k+1}\alpha_k\boldsymbol{p}_k&=\nabla f_{k+1}-\nabla f_k\\B_{k+1}(\boldsymbol{x}_{k+1}-\boldsymbol{x}_k)&=\nabla f_{k+1}-\nabla f_k\end{align}\]
となる。\(\boldsymbol{s}_k=\boldsymbol{x}_{k+1}-\boldsymbol{x}_k\)と\(\boldsymbol{y}_k=\nabla f_{k+1}-\nabla f_k\)を使うと
\[ B_{k+1}\boldsymbol{s}_k=\boldsymbol{y}_k\]
とかける。これのことをセカント条件/セカント方程式という。条件とも方程式ともいわれるのには理由がある。というのも、これは一見して\(\boldsymbol{s}_k\)が未知数であるような方程式の形を取っているが未知なのは\(B_{k+1}\)であって、方程式の数\(n\)に対して未知数の数\(n(n-1)/2\)（対称行列であることから）が大きいのである。よって、そのような\(B_{k+1}\)がいくらでも存在してしまう。
そこでセカント方程式を満たしつつ\(B_k\)とある種の距離がもっとも小さいものを選ぶことにする。つまり、
\[ B_{k+1}=\mathop{\rm argmin}_{B=B^\top}\|B - B_k\|\text{ s.t. }B\boldsymbol{s}_k=\boldsymbol{y}_k\]
とする。セカント方程式を条件とした最小化だ。これがセカント方程式ともセカント条件ともいう理由。\(B_{k+1}\)の選び方は、この問題における「ある種の距離」の取り方によっていくらか存在する。準ニュートン法といわれるものの中でいろいろなバリエーションが生じるのはこの部分による（より一般的には\(B_{k+1}\)もしくはその逆行列の取り方）。
そのなかの1つであるBFGS公式では\(B_{k+1}\)を陽に求めない。というのも、探索方向の決定に用いるのは\(B_k\)ではなく\(B_k^{-1}\)なので\(B_k^{-1}=H_k\)を保持しておけば十分だからだ。\(B_k\)に対するセカント条件を\(H_k\)に対するセカント条件にするには単純に両辺\(B_k^{-1}\)をかけてやれば
\[ H_{k+1}\boldsymbol{y}_k=\boldsymbol{s}_k\]
となる。従って考えるべき問題は
\[ H_{k+1}=\mathrm{argmin}_{H=H^\top}\|H-H_k\|\text{ s.t. }H\boldsymbol{y}_k=\boldsymbol{s}_k\]
実はこの問題は適切な行列ノルムのもとで陽に解ける。
\[ \begin{gather}H_{k+1}=(I-\rho_k\boldsymbol{s}_k\boldsymbol{y}_k^\top)H_k(I-\rho_k\boldsymbol{y}_k\boldsymbol{s}_k^\top)+\rho_k\boldsymbol{s}_k\boldsymbol{s}_k^\top\\\rho_k=\frac{1}{\boldsymbol{y}_k^\top\boldsymbol{s}_k}\end{gather}\]
ある種の距離を少しだけ明らかにしておくと、このような解を得るためには行列のフロベニウスノルム
\[ \|A\|_\mathrm{F}^2=\sum_{i=1}^n\sum_{j=1}^na_{ij}^2\]
に重みを付けた重み付きフロベニウスノルム
\[ \|W^{1/2}AW^{1/2}\|_\mathrm{F}^2\]
において重み\(W\)を\(W\boldsymbol{s}_k=\boldsymbol{y}_k\)を満たすようにとったものを用いる。これによって解は上記のとおりに一意に定まる。以上が直線探索の部分を除いたBFGS公式による準ニュートン法の概略になる。一旦アルゴリズムの全貌をまとめておくと、
- 初期値\(\boldsymbol{x}_0\)と\(H_0\)を与える。
- \(k\colon=0\)
- \(\|\nabla f_k\|\gt\epsilon\)なら以下を繰り返す
- \(\boldsymbol{p}_k\colon=-H_k\nabla f_k\)
- \(\alpha_k\colon=\mathrm{LS}_f(\boldsymbol{x}_k,\boldsymbol{p}_k)\)
- \(\boldsymbol{x}_{k+1}\colon=\boldsymbol{x}_k+\alpha_k\boldsymbol{p}_k\)
- \(\boldsymbol{s}_k\colon=\boldsymbol{x}_{k+1}-\boldsymbol{x}_k\)
- \(\boldsymbol{y}_k\colon=\nabla f_{k+1}-\nabla f_k\)
- \(\rho_k\colon=1/\boldsymbol{y}_k^\top\boldsymbol{s}_k\)
- \(H_{k+1}\colon=(I-\rho_k\boldsymbol{s}_k\boldsymbol{y}_k^\top)H_k(I-\rho_k\boldsymbol{y}_k\boldsymbol{s}_k^\top)+\rho_k\boldsymbol{s}_k\boldsymbol{s}_k^\top\)
- \(k\colon=k+1\)
- 繰り返し終わり
- \(\boldsymbol{x}_k, f_k\)がそれぞれ最適解、最適値
\(\mathrm{LS}_f\)は直線探索を示す。この部分をどうするかについてはBFGS準ニュートン法が指し示すところの範疇ではないが、ステップサイズの基準として用いられるWolfeの条件を満たしていれば十分である。また、\(H_0\)をどのようにして与えるかというところも同じくBFGS準ニュートン法が指し示すところの範疇ではないが
- 有限差分近似で推定したヘッセ行列の逆行列を使う
- 単位行列を使う
- 重み付けした対角行列を使う
などのバリエーションがある。
性質などはともかくとして以上がBFGS公式を用いた準ニュートン法の手順。先に述べたようにこのアルゴリズムのなかでは2次の勾配\(\nabla^2f_k\)使っていないしなんの行列の逆行列ももとめていない。
メモ
BFGS = Broyden, Fletcher, Goldfarb and Shanno
純粋な方のニュートン法の収束する理屈や\(B_k,H_k\)の正定値性などの話は不十分。Sherman-Morrison-Woodbury公式や\(B_k\)を近似するDFP公式について触れていない。
Scipyのoptimize.minimize(method='BFGS')が参照する具体的な直線探索や逆ヘッセ行列の初期値の与え方なども同参考文献にのっていた（これがもともとの動機）。
参考文献
Nocedal, Jorge, and Stephen J. Wright. “Numerical optimization” Second Edition (2006). pp.136- 140, 6.1 The BFGS method.