BFGS法は非線形の最適化の際によく用いられるアルゴリズムですが、準ニュートン法の$1$つと見なすことができます。当記事ではBFGS法・準ニュートン法の概要や数式理解、具体的な例に対して計算を行うPythonプログラムなどの確認を行いました。
「新版 数理計画入門(朝倉書店)」の$4.6$節の「準ニュートン法」などの内容を参考に作成を行いました。
・用語/公式解説
https://www.hello-statisticians.com/explain-terms
Contents
BFGS法の仕組み
ニュートン法と準ニュートン法
多次元ニュートン法は最適化の各反復計算の際にヘッセ行列を活用することで収束を速めることができる手法である。
ニュートン法は収束が速いなど有用である一方で、変数が多い時などはヘッセ行列の計算が容易でないなど、必ずしも良い手法であるとは言えない。
上記を受けて準ニュートン法ではヘッセ行列の計算を定義通りに行わずに漸化式的にヘッセ行列の推定を行う。たとえば$1$次元の変数$x \in \mathbb{R}$と目的関数$f(x)$に対し、下記のように$1 \times 1$のヘッセ行列の推定を行うことができる。
$$
\large
\begin{align}
\nabla^2 f(x^{(k+1)}) \simeq \frac{\nabla f(x^{(k+1)}) – \nabla f(x^{(k)})}{x^{(k+1)}-x^{(k)}} \quad (1)
\end{align}
$$
上記のようにヘッセ行列の近似を行うことで、ニュートン法と同様の計算を行うのが準ニュートン法である。以下、当記事では準ニュートン法の$1$つの手法であるBFGS法について詳しく取り扱った。
ヘッセ行列の近似とBFGS法
$n$次元の変数ベクトル$\mathbf{x} \in \mathbb{R}$に対して目的関数$f(\mathbf{x})$を定義する。このとき$k+1$ステップ目のヘッセ行列$\nabla^2 f(\mathbf{x}^{(k+1)})$の近似行列を$B^{(k+1)}$とおくと、下記が成立する。
$$
\large
\begin{align}
B^{(k+1)} \mathbf{s}^{(k)} &= \mathbf{y}^{(k)} \quad (2) \\
\mathbf{s}^{(k)} &= \mathbf{x}^{(k+1)} – \mathbf{x}^{(k)} \\
\mathbf{y}^{(k)} &= \nabla f(\mathbf{x}^{(k+1)}) – \nabla f(\mathbf{x}^{(k)})
\end{align}
$$
$(2)$式は$(1)$式と対応させて理解すると良い。$(2)$式が成立する近似行列の定め方はいくつか提案されてきたが、よく用いられるBFGS法では下記のような式を元に近似行列$B^{(k+1)}$を定義する。
$$
\large
\begin{align}
B^{(k+1)} &= B^{(k)} + \frac{1}{\beta^{(k)}} \mathbf{y}^{(k)} (\mathbf{y}^{(k)})^{\mathrm{T}} – \frac{1}{\gamma^{(k)}} B^{(k)} \mathbf{s}^{(k)} (\mathbf{s}^{(k)})^{\mathrm{T}} B^{(k)} \\
\beta^{(k)} &= (\mathbf{y}^{(k)})^{\mathrm{T}} \mathbf{s}^{(k)} \\
\gamma^{(k)} &= (\mathbf{s}^{(k)})^{\mathrm{T}} B^{(k)} \mathbf{s}^{(k)}
\end{align}
$$
BFGS法におけるヘッセ行列の近似式
$$
\large
\begin{align}
B^{(k+1)} &= B^{(k)} + \frac{1}{\beta^{(k)}} \mathbf{y}^{(k)} (\mathbf{y}^{(k)})^{\mathrm{T}} – \frac{1}{\gamma^{(k)}} B^{(k)} \mathbf{s}^{(k)} (\mathbf{s}^{(k)})^{\mathrm{T}} B^{(k)} \\
\beta^{(k)} &= (\mathbf{y}^{(k)})^{\mathrm{T}} \mathbf{s}^{(k)} \\
\gamma^{(k)} &= (\mathbf{s}^{(k)})^{\mathrm{T}} B^{(k)} \mathbf{s}^{(k)}
\end{align}
$$
上記のように定めた近似行列$B^{(k+1)}$を用いるとき、下記の下記の$(a)$〜$(c)$が成立する。
$(a) \,$ $B^{(k+1)}$に関して$(2)$式が成立する
$(b) \,$ $B^{(k)}$が対称ならば$B^{(k+1)}$も対称である
$(c) \,$ $B^{(k)}$が正定値かつ$\beta^{(k)}>0$ならば$B^{(k+1)}$も正定値である
以下、$(a)$が成立することを確認する。
$$
\large
\begin{align}
& B^{(k+1)} \mathbf{s}^{(k)} = B^{(k)} \mathbf{s}^{(k)} + \frac{1}{\beta^{(k)}} \mathbf{y}^{(k)} (\mathbf{y}^{(k)})^{\mathrm{T}} \mathbf{s}^{(k)} – \frac{1}{\gamma^{(k)}} B^{(k)} \mathbf{s}^{(k)} (\mathbf{s}^{(k)})^{\mathrm{T}} B^{(k)} \mathbf{s}^{(k)} \\
&= B^{(k)} \mathbf{s}^{(k)} + \frac{1}{\cancel{(\mathbf{y}^{(k)})^{\mathrm{T}} \mathbf{s}^{(k)}}} \mathbf{y}^{(k)} \cancel{(\mathbf{y}^{(k)})^{\mathrm{T}} \mathbf{s}^{(k)}} – \frac{1}{\cancel{(\mathbf{s}^{(k)})^{\mathrm{T}} B^{(k)} \mathbf{s}^{(k)}}} B^{(k)} \mathbf{s}^{(k)} \cancel{(\mathbf{s}^{(k)})^{\mathrm{T}} B^{(k)} \mathbf{s}^{(k)}} \\
&= \mathbf{y}^{(k)}
\end{align}
$$
上記より$(a)$が成立する。また、途中計算では$(\mathbf{y}^{(k)})^{\mathrm{T}} \mathbf{s}^{(k)}$と$(\mathbf{s}^{(k)})^{\mathrm{T}} B^{(k)} \mathbf{s}^{(k)}$がスカラーであることに着目して約分を行なった。
ここで計算した近似行列の$B^{(k+1)}$をニュートン法におけるヘッセ行列に置き換えて下記のような計算を行うことで、準ニュートン法を用いた近似解を計算することができる。
$$
\large
\begin{align}
\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} – (B^{(k)})^{-1} \nabla f(\mathbf{x}^{(k)})
\end{align}
$$
Pythonを用いたBFGS法の実行
問題設定
「数理計画入門(朝倉書店)」の$4.4$節の「最急降下法」と同様に下記の関数の最適化にあたって、最急降下法のプログラムの作成を行う。
$$
\large
\begin{align}
f(\mathbf{x}) = f(x_1,x_2) = (x_1-1)^2 + 10(x_1^2-x_2)^2
\end{align}
$$
勾配ベクトルの計算・ヘッセ行列の推定
勾配ベクトル$\nabla f(\mathbf{x})$は下記のように計算できる。
$$
\large
\begin{align}
\nabla f(\mathbf{x}) &= \left(\begin{array}{c} \displaystyle \frac{\partial f}{\partial x_1} \\ \displaystyle \frac{\partial f}{\partial x_2} \end{array} \right) \\
&= \left(\begin{array}{c} 2(x_1-1) + 20(x_1^2-x_2) \cdot 2x_1 \\ 20(x_1^2-x_2) \cdot (-1) \end{array} \right) \\
&= \left(\begin{array}{c} 40x_1^3-40x_1x_2+2x_1-2 \\ -20(x_1^2-x_2) \end{array} \right)
\end{align}
$$
ヘッセ行列は下記の数式を元に近似行列を推定する。
$$
\large
\begin{align}
B^{(k+1)} = B^{(k)} + \frac{1}{\beta^{(k)}} \mathbf{y}^{(k)} (\mathbf{y}^{(k)})^{\mathrm{T}} – \frac{1}{\gamma^{(k)}} B^{(k)} \mathbf{s}^{(k)} (\mathbf{s}^{(k)})^{\mathrm{T}} B^{(k)}
\end{align}
$$
BFGS法の実行
下記を実行することでBFGS法による最適化を行うことができる。
import numpy as np
def calc_func_value(x1,x2):
return (x1-1)**2 + 10*(x1**2-x2)**2
def calc_gradient(x1,x2):
return np.array([40*x1**3 - 40*x1*x2 + 2*x1 - 2, -20*(x1**2-x2)])
x = np.array([0., 1.])
B = np.eye(2)
alpha = np.arange(-5., 5., 0.0000002)
for k in range(10):
grad = np.dot(np.linalg.inv(B),calc_gradient(x[0],x[1]))
x1, x2 = x[0]+alpha*grad[0], x[1]+alpha*grad[1]
f = calc_func_value(x1,x2)
s = (x - np.array([x1[np.argmin(f)], x2[np.argmin(f)]])).reshape([2,1])
y = (calc_gradient(x[0],x[1]) - calc_gradient(x1[np.argmin(f)],x2[np.argmin(f)])).reshape([2,1])
x[0], x[1] = x1[np.argmin(f)], x2[np.argmin(f)]
beta, gamma = np.dot(y.T,s), np.dot(s.T,np.dot(B,s))
B = B + np.dot(y,y.T)/beta - np.dot(np.dot(B,s),np.dot(s.T,B))/gamma
print("Iter: {:.0f}, x: {}".format(k+1,x))・実行結果
Iter: 1, x: [ 0.09988482 0.00115184]
Iter: 2, x: [ 0.32845959 0.00381447]
Iter: 3, x: [ 0.63413382 0.29092239]
Iter: 4, x: [ 0.64276477 0.41585505]
Iter: 5, x: [ 0.83666141 0.66037833]
Iter: 6, x: [ 0.9954299 0.99483118]
Iter: 7, x: [ 1.00116088 1.0024978 ]
Iter: 8, x: [ 0.99998693 0.99998305]
Iter: 9, x: [ 0.99999995 0.99999989]
Iter: 10, x: [ 1. 1.]
上記の結果は「数理計画入門(朝倉書店)」の表$4.3$の結果に基本的には一致する。「数理計画入門(朝倉書店)」に直線探索の設定がないので、厳密な再現ではないことに注意が必要である。