最適化
準ニュートン法
BFGS

逆行列補題について

まずブロック行列とその逆行列から話を始めよう。
$I_{n},I_{m},P,Q,R,S$をそれぞれ以下の条件を満たす行列とする。

$・I_{n}:n×n単位行列$ , $I_{m}:m×m単位行列$
$・P:n×n行列$ , $Q:n×m行列$ , $R:m×m行列$ , $S:m×n行列$
$・P,Rはそれぞれ逆行列を持つ$
$・P+QRS$ $及び$ $R^{-1}+SP^{-1}Q$ はそれぞれ逆行列を持つ

その上で以下のブロック行列を考えよう。

\begin{pmatrix}
P & -Q \\
S & R^{-1}
\end{pmatrix}

この行列は次のように分解できる。

\begin{align}
\begin{pmatrix}
P & -Q \\
S & R^{-1}
\end{pmatrix}
&=\begin{pmatrix}
I_{n} & -QR \\
0 & I_{m}
\end{pmatrix}
\begin{pmatrix}
P+QRS & 0 \\
0 &R^{-1}
\end{pmatrix}
\begin{pmatrix}
I_{n} & 0 \\
RS & I_{m}
\end{pmatrix}\\
\begin{pmatrix}
P & -Q \\
S & R^{-1}
\end{pmatrix}
&=\begin{pmatrix}
I_{n} & 0 \\
SP^{-1} & I_{m}
\end{pmatrix}
\begin{pmatrix}
P & 0 \\
0 &R^{-1}+SP^{-1}Q
\end{pmatrix}
\begin{pmatrix}
I_{n} & -P^{-1}Q \\
0 & I_{m}
\end{pmatrix}
\end{align}

ところで$X,Y,Z$がそれぞれ逆行列を持つ場合、行列$XYZ$の逆行列について、
$(XYZ)^{-1}=Z^{-1}Y^{-1}X^{-1}$が成り立つ。
また上三角行列,下三角行列の逆行列について、

\begin{align}
\begin{pmatrix}
I_{n} & -QR \\
0 & I_{m}
\end{pmatrix}^{-1}
&=
\begin{pmatrix}
I_{n} & QR \\
0 & I_{m}
\end{pmatrix}\\

\begin{pmatrix}
I_{n} & 0 \\
RS & I_{m}
\end{pmatrix}^{-1}
&=
\begin{pmatrix}
I_{n} & 0 \\
-RS & I_{m}
\end{pmatrix}\\

\begin{pmatrix}
I_{n} & 0 \\
SP^{-1} & I_{m}
\end{pmatrix}^{-1}
&=
\begin{pmatrix}
I_{n} & 0 \\
-SP^{-1} & I_{m}
\end{pmatrix}\\

\begin{pmatrix}
I_{n} & -P^{-1}Q \\
0 & I_{m}
\end{pmatrix}^{-1}
&=
\begin{pmatrix}
I_{n} & P^{-1}Q \\
0 & I_{m}
\end{pmatrix}\\

\end{align}

がそれぞれ成り立つ。
これらを踏まえて先のブロック行列の逆行列を考えよう。
すると

\begin{align}
\begin{pmatrix}
P & -Q \\
S & R^{-1}
\end{pmatrix}^{-1}
&=
\begin{pmatrix}
I_{n} & 0 \\
RS & I_{m}
\end{pmatrix}^{-1}
\begin{pmatrix}
P+QRS & 0 \\
0 &R^{-1}
\end{pmatrix}^{-1}
\begin{pmatrix}
I_{n} & -QR \\
0 & I_{m}
\end{pmatrix}^{-1}\\
&=
\begin{pmatrix}
I_{n} & 0 \\
-RS & I_{m}
\end{pmatrix}
\begin{pmatrix}
(P+QRS)^{-1} & 0 \\
0 &R
\end{pmatrix}
\begin{pmatrix}
I_{n} & QR \\
0 & I_{m}
\end{pmatrix}\\
&=
\begin{pmatrix}
(P+QRS)^{-1} & (P+QRS)^{-1}QR \\
-RS(P+QRS)^{-1} & R-RS(P+QRS)^{-1}QR
\end{pmatrix}...①
\end{align}
\begin{align}
\begin{pmatrix}
P & -Q \\
S & R^{-1}
\end{pmatrix}^{-1}
&=
\begin{pmatrix}
I_{n} & -P^{-1}Q \\
0 & I_{m}
\end{pmatrix}^{-1}
\begin{pmatrix}
P & 0 \\
0 &R^{-1}+SP^{-1}Q
\end{pmatrix}^{-1}
\begin{pmatrix}
I_{n} & 0 \\
SP^{-1} & I_{m}
\end{pmatrix}^{-1}\\
&=
\begin{pmatrix}
I_{n} & P^{-1}Q \\
0 & I_{m}
\end{pmatrix}
\begin{pmatrix}
P^{-1} & 0 \\
0 & (R^{-1}+SP^{-1}Q )^{-1}
\end{pmatrix}
\begin{pmatrix}
I_{n} & 0 \\
-SP^{-1} & I_{m}
\end{pmatrix}\\
&=
\begin{pmatrix}
P^{-1}-P^{-1}Q(R^{-1}+SP^{-1}Q)^{-1}SP^{-1} & P^{-1}Q(R^{-1}+SP^{-1}Q)^{-1} \\
-(R^{-1}+SP^{-1}Q)^{-1}SP^{-1} & (R^{-1}+SP^{-1}Q)^{-1}
\end{pmatrix}...②
\end{align}

式①,②を比較すれば、

(P+QRS)^{-1}=P^{-1}-P^{-1}Q(R^{-1}+SP^{-1}Q)^{-1}SP^{-1}

が得られる。これを$逆行列補題(Sherman–Morrison-Woodburyの公式)$という。

さらに逆行列補題の特別な場合を考えよう。
$u,vをn次元列ベクトル,wをスカラー,Pを逆行列を持つn×n行列$とする。
このとき$v^{T}P^{-1}u$がスカラーになることを踏まえれば、

(P+\frac{uv^{T}}{w})^{-1}=P^{-1}-\frac{P^{-1}uv^{T}P^{-1}}{w+v^{T}P^{-1}u}...③

が得られる。
これが後述の$BFGS(H公式)$の導出に深く関わってくる。

準ニュートン法について

$B_{k}$を正定値対称行列として、関数$f(x_{k}+d)$を次のように近似する。

f(x_{k}+d) \approx f(x_{k})+▽f(x_{k})^{T}d+\frac{1}{2}d^{T}B_{k}d\\
d=x-x_{k}

ここで、関数$f(x_{k}+d)$をdの関数としてdで微分しその解を$d_{k}$した、

d_{k}=-B^{-1}_{k}▽f(x_{k})

を利用して、関数$f(x)=0$の数値解を求めるアルゴリズムを準ニュートン法という。
※$B_{k}$が正定値対称行列であるために、以下の条件を満たさなくてはいけない。また(3)をセカント条件という。

\begin{align}
(1)s_{k}&=x_{k+1}-x_{k}\\
(2)y_{k}&=▽f(x_{k+1})-▽f(x_{k})\\
(3)y_{k}&=B_{k+1}s_{k}\\
\end{align}

さてセカント条件を満たす$B_{k}$のなかで代表的なものが以下の$BFGS(B公式)$だろう。

B_{k+1}=
(B_{k}+\frac{y_{k}y^{T}_{k}}{s^{T}_{k}y_{k} } )
-(\frac{B_{k}s_{k}s^{T}_{k}B_{k} }{s^{T}_{k}B_{k}s_{k}})

この逆行列を用いて$d_{k}$を求めるというのが準ニュートン法の流れになるが、しかしわざわざ$B_{k}$の逆行列を求めるのは間怠こしい。そこで先述の逆行列補題の出番となる。
つまり逆行列補題を利用して$B_{k}$の逆行列を$H_{k}$とした上で、$d_{k}=-H_{k}▽f(x_{k})$を直接求めようというわけだ。

逆行列補題によるBFGS(H公式)の導出

式③を思い出してほしい。

(P+\frac{uv^{T}}{w})^{-1}=P^{-1}-\frac{P^{-1}uv^{T}P^{-1}}{w+v^{T}P^{-1}u}

これを利用して、行列$B_{k}$の逆行列を求めよう。
すると、

B^{-1}_{k+1}=
(B_{k}+\frac{y_{k}y^{T}_{k}}{s^{T}_{k}y_{k} } )^{-1}
+
\frac{(B_{k}+\frac{y_{k}y^{T}_{k}}{s^{T}_{k}y_{k}})^{-1}B_{k}s_{k}s^{T}_{k}B_{k}(B_{k}+\frac{y_{k}y^{T}_{k}}{s^{T}_{k}y_{k}})^{-1}}{s^{T}_{k}B_{k}s_{k}-s^{T}_{k}B_{k}(B_{k}+\frac{y_{k}y^{T}_{k}}{s^{T}_{k}y_{k}})^{-1}B_{k}s_{k}}・・・④

が得られる。ここで式④右辺第1項について、

(B_{k}+\frac{y_{k}y^{T}_{k}}{s^{T}_{k}y_{k} } )^{-1}
=B^{-1}_{k}-\frac{B^{-1}_{k}y_{k}y^{T}_{k}B^{-1}_{k}}{s^{T}_{k}y_{k}+y^{T}_{k}B^{-1}_{k}y_{k}}・・・⑤

となる。式⑤を踏まえて式④右辺第2項を鑑みると分子・分母はそれぞれ次の通りとなる。

\begin{align}
(B_{k}+\frac{y_{k}y^{T}_{k}}{s^{T}_{k}y_{k}})^{-1}B_{k}s_{k}s^{T}_{k}B_{k}(B_{k}+\frac{y_{k}y^{T}_{k}}{s^{T}_{k}y_{k}})^{-1}
&=(I-\frac{B^{-1}y_{k}y^{T}_{k}}{s^{T}_{k}y_{k}+y^{T}_{k}B^{-1}_{k}y_{k}})s_{k}s^{T}_{k}(I-\frac{y_{k}y^{T}_{k}B^{-1}}{s^{T}_{k}y_{k}+y^{T}_{k}B^{-1}_{k}y_{k}})\\
&=s_{k}s^{T}_{k}-\frac{s^{T}_{k}y_{k}}{s^{T}_{k}y_{k}+y^{T}_{k}B^{-1}_{k}y_{k}}(B^{-1}_{k}y_{k}s^{T}_{k}+s_{k}y^{T}_{k}B^{-1}_{k})+(\frac{s^{T}_{k}y_{k}}{s^{T}_{k}y_{k}+y^{T}_{k}B^{-1}_{k}y_{k}})
^2B^{-1}_{k}y_{k}y^{T}_{k}B^{-1}_{k}・・・⑥
\end{align}

※$I$は単位行列

{s^{T}_{k}B_{k}s_{k}-s^{T}_{k}B_{k}(B_{k}+\frac{y_{k}y^{T}_{k}}{s^{T}_{k}y_{k}})^{-1}B_{k}s_{k}}
=\frac{(s^{T}_{k}y_{k})^{2}}{s^{T}_{k}y_{k}+y^{T}_{k}B^{-1}_{k}y_{k}}・・・⑦

従って式⑤,⑥,⑦を整理することにより、式④について行列$B_{k}$の逆行列を$H_{k}$と置き換えれば、

\begin{align}
H_{k+1}&=H_{k}-\frac{H_{k}y_{k}s^{T}_{k}+s_{k}y^{T}_{k}H_{k}}{s^{T}_{k}y_{k}}+(1+\frac{y^{T}_{k}H_{k}y_{k}}{s^{T}_{k}y_{k}})\frac{s_{k}s^{T}_{k}}{s^{T}_{k}y_{k}}\\
&= (I-\frac{y_{k}s^{T}_{k}}{s^{T}_{k}y_{k}})^{T}H_{k}(I-\frac{y_{k}s^{T}_{k}}{s^{T}_{k}y_{k}})+\frac{s_{k}s^{T}_{k}}{s^{T}_{k}y_{k}}
\end{align}

が得られる。こうして$BFGS(H公式)$が導かれる。

準ニュートン法の流れ

さて$BFGS(B公式,H公式)$による流れは以下の通り。
・$BFGS(B公式)$
$(step:1)$初期値$x_{0}、B_{0}$(通常は単位行列)を選び、k:=0とする
$(step:2)$$▽f(x_{k})=0$ならば計算終了。$▽f(x_{k})\neq 0$ならば$d_{k}=-B^{-1}_{k}▽f(x_{k})$を求める。
$(step:3)α_{k}$をステップ幅として、$x_{k+1}=x_{k}+α_{k}d_{k}$とする$^{※}$
$(step:4)s_{k}$及び$y_{k}$を算出
$(step:5)BFGS(B公式)$により行列$B_{k}$を更新

B_{k+1}=
(B_{k}+\frac{y_{k}y^{T}_{k}}{s^{T}_{k}y_{k} } )
-(\frac{B_{k}s_{k}s^{T}_{k}B_{k} }{s^{T}_{k}B_{k}s_{k}})

$(step:6)k:=k+1$として$(step:1)$に戻る

・$BFGS(H公式)$
$(step:1)$初期値$x_{0}、H_{0}$(通常は単位行列)を選び、k:=0とする
$(step:2)$$▽f(x_{k})=0$ならば計算終了。$▽f(x_{k})\neq 0$ならば$d_{k}=-H_{k}▽f(x_{k})$を求める。
$(step:3)α_{k}$をステップ幅として、$x_{k+1}=x_{k}+α_{k}d_{k}$とする$^{※}$
$(step:4)s_{k}$及び$y_{k}$を算出
$(step:5)BFGS(H公式)$により行列$H_{k}$を更新

\begin{align}
H_{k+1}&= (I-\frac{y_{k}s^{T}_{k}}{s^{T}_{k}y_{k}})^{T}H_{k}(I-\frac{y_{k}s^{T}_{k}}{s^{T}_{k}y_{k}})+\frac{s_{k}s^{T}_{k}}{s^{T}_{k}y_{k}}
\end{align}

$(step:6)k:=k+1$として$(step:1)$に戻る

こうして関数$f(x_{k})=0$の数値解が求められるわけだ。

※ステップ幅計算については$Appendix$参照

Appendix

上述の$(step:3)$ステップ幅$α_{k}$は以下の$Armijo条件,Wolfe条件$に基づいて算出することが一般的。

・$Armijo条件$
$0<ξ<1$となる定数$ξ$について

f(x_{k}+α_{k}d_{k}) \leq f(x_{k}) + ξα_{k}▽f(x_{k})^{T}d_{k} 

を満たす$α_{k}>0$をステップ幅とする。

・$Wolfe条件$
$0<ξ_{1}<ξ_{2}<1$となる定数$ξ_{1},ξ_{2}$について

f(x_{k}+α_{k}d_{k}) \leq f(x_{k}) + ξ_{1}α_{k}▽f(x_{k})^{T}d_{k} \\
ξ_{2}▽f(x_{k})^{T}d_{k} \leq ▽f(x_{k}+α_{k}d_{k})^{T}d_{k}

を満たす$α_{k}>0$をステップ幅とする。

参考文献

・「統計のための行列代数」(D.A.ハーヴィル/丸善出版)
・「工学基礎 最適化とその応用」(矢部博/数理工学社)
・「数理計画入門 新版」(福島雅夫/朝倉書店)