はじめに

以前、以下の単体法の記事を書いた。

yamagensakam.hatenablog.com

今回はこのシリーズ第2弾で、以下の制約なし非線形計画問題を準ニュートン法で解くプログラムを実装する。

$\displaystyle \min_{{\bf x} \in \mathbb{R}^n} f({\bf x}) \tag{1}$

手法

制約なし非線形計画問題の解き方

非線形計画問題を解く際よく行われるのは、最適性の1次の必要条件である $\nabla f({\bf x}_{*}) = {\bf 0}$ を満たす ${\bf x}_{*}$ （に十分に近い点）を以下の反復計算で見つける方法である。

$\displaystyle {\bf x}^{(k + 1)} = {\bf x}^{(k)} + \alpha_{k} {\bf d}({\bf x}^{(k)}) \tag{2}$

ここで ${\bf d}({\bf x}^{(k)})$ は探索方向で、 $f({\bf x}^{(k)} + \alpha_{k} {\bf d}({\bf x}^{(k)}))$ が減少する方向に探索方向を選ぶ。探索方向の選び方には様々なバリエーションがあり、その1つが今回実装した準ニュートン法。なお $\alpha_{k}$ はステップ幅で、これについても目的関数が減少するよう適切に選ぶ必要があるが、こちらの選び方の説明は今回省略する*1。

ニュートン法とその問題点

準ニュートン法の前にニュートン法について書く。 $f({\bf x})$ のヘッセ行列を $\nabla^{2} f({\bf x})$ とするとニュートン法では探索方向を

$\displaystyle {\bf d}({\bf x}^{(k)}) = - {\nabla^{2} f({\bf x}^{(k)})}^{-1} \nabla f({\bf x}^{(k)}) \tag{3}$

とする。これは、目的関数を ${\bf x}^{(k)}$ まわりでの2次近似した関数

$\displaystyle f({\bf x}^{(k)} + {\bf d}) \sim f({\bf x}^{(k)}) + \nabla f({\bf x}^{(k)})^T {\bf d} + \frac{1}{2} {\bf d}^T \nabla^2 f({\bf x}^{(k)}) {\bf d} \tag{4}$

の ${\bf d}$ に関する停留点に当たり、ヘッセ行列が正定値なら2次近似の最小解となる。従って、ニュートン法は目的関数を2次近似してその停留点を求めるという操作を収束するまで繰り返すアルゴリズムだといえる。ニュートン法は2次収束するため収束性はよい。

一方問題点としては、目的関数が狭義凸関数でない場合ヘッセ行列が正定値にならないケースがある*2。そのため、 ${\bf d}$ が下降方向になることは保証されない。そもそもヘッセ行列が逆行列を持たない可能性もある。更にはヘッセ行列を計算すること自体が難しい場合もある。

準ニュートン法のアルゴリズム

ニュートン法の上記の問題を解消するため、ヘッセ行列そのものではなくヘッセ行列を近似する正定値対称行列 ${\bf B}_k$ を使って探索方向を

$\displaystyle {\bf d}({\bf x}^{(k)}) = - {{\bf B}_{k}}^{-1} \nabla f({\bf x}^{(k)}) \tag{5}$

とするのが準ニュートン法。 ${\bf B}_k$ が正定値なので探索方向は常に下降方向であり逆行列も存在する。

ヘッセ行列の近似

直前の反復でヘッセ行列を近似した正定値対称行列 ${\bf B}_k$ が得られているとして、これを用いて次の近似行列 ${\bf B}_{k+1}$ を求めることを考える。まず、 $\nabla f({\bf x}^{(k)})$ の1次近似として、

$\displaystyle \nabla f({\bf x}^{(k)}) = \nabla f({\bf x}^{(k + 1)} + ({\bf x}^{(k)} - {\bf x}^{(k + 1)} )) \sim \nabla f({\bf x}^{(k + 1)}) + \nabla^2 f({\bf x}^{(k + 1)}) ({\bf x}^{(k)} - {\bf x}^{(k + 1)} ) \tag{6}$

が得られる。式（6）の関係よりヘッセ行列 $\nabla^{2} f({\bf x}^{(k + 1)})$ の近似行列 ${\bf B}_{k+1}$ としては以下の条件を満たすのが適切だといえる。

$\displaystyle {\bf B}_{k + 1} ({\bf x}^{(k + 1)} - {\bf x}^{(k)} ) =\nabla f({\bf x}^{(k + 1)}) - \nabla f({\bf x}^{(k)})\tag{7}$

この式（7）はセカント条件と呼ばれる。これに加えて ${\bf B}_{k + 1}$ は正定値対称行列でなければならない。これらの条件を満たす行列はいくつか提案されているが、中でもよく知られているのはBFGS公式と呼ばれるもので、これは直前の反復で得られた正定値対称行列 ${\bf B}_{k}$ および ${\bf s}^{(k)} ={\bf x}^{(k + 1)} - {\bf x}^{(k)} , \ {\bf y}^{(k)} = \nabla f({\bf x}^{(k + 1)}) - \nabla f({\bf x}^{(k)})$ を用いて

$\displaystyle {\bf B}_{k + 1} =　{\bf B}_{k} - \frac{{\bf B}_{k} {\bf s}^{(k)} ({\bf B}_{k} {\bf s}^{(k)})^T}{({\bf s}^{(k)})^T {\bf B}_{k} {\bf s}^{(k)}} + \frac{{\bf y}^{(k)} ({\bf y}^{(k)})^T }{({\bf s}^{(k)})^T {\bf y}^{(k)} }\tag{8}$

と計算できる。詳細は省くが式（8）の両辺に右から ${\bf s}^{(k)}$ をかければ ${\bf B}_{k + 1} {\bf s}^{(k)} = {\bf y}^{(k)}$ となりセカント条件を満たすことを示せる。また、ステップ幅 $\alpha$ を適切に決めれば正定値性は満たされる*3。

逆行列の更新

式（5）にあるように探索方向を求めるために ${\bf B}_{k}$ の逆行列が必要になるが、次元が大きくなると更新毎に式（8）で ${\bf B}_{k+1}$ を計算してその逆行列を計算する、なんてやってると計算量的に厳しい。そこで式（8）の右辺第2項、第3項が共にランク1の行列であることを考えると、Sherman-Morrisonの公式*4を2回適用することで逆行列そのものの更新式を導くことができる。具体的な式展開などは以下の記事が参考になる。

qiita.com

結果としては、

$\displaystyle ({\bf B}_{k + 1})^{-1} =　\left({\bf I} - \frac{{\bf s}^{(k)} ({\bf y}^{(k)})^T }{({\bf s}^{(k)})^T {\bf y}^{(k)} } \right)({\bf B}_{k})^{-1} \left({\bf I} - \frac{{\bf y}^{(k)} ({\bf s}^{(k)})^T }{({\bf s}^{(k)})^T {\bf y}^{(k)} } \right) + \frac{{\bf s}^{(k)} ({\bf s}^{(k)})^T }{({\bf s}^{(k)})^T {\bf y}^{(k)} } \tag{9}$

という更新式が得られる。

BFGS公式について

式（8）でやや唐突にBFGS公式を書いたが、ここではBFGS公式の考え方を示す。まずBFGS公式含めて提案されている多くの更新式は、正定値対称性とセカント条件を満たす行列のうち、直前の行列 ${\bf B}_{k}$ にもっとも近い行列を ${\bf B}_{k+1}$ としている。これはヘッセ行列は急激に変化しないという事前知識を基づく。これを式で書くと適当な行列間の非類似度を測る関数 $D(\cdot, \cdot)$ を導入して、

$\displaystyle \begin{aligned} \min_{{\bf B}} \ & D({\bf B}, {\bf B}_{k}) \\ {\rm s.t.} \ & {\bf B} {\bf s}^{(k)} = {\bf y}^{(k)} \\ \ & {\bf B} \succ {\bf O} \\ \ & {\bf B} = {\bf B}^T \end{aligned} \tag{10}$

という最適化問題の解を ${\bf B}_{k+1}$ としている。

BFGS公式は非類似度としてKLダイバージェンスを採用して式（10）を解くことに相当する。KLダイバージェンスは確率分布間の非類似度を測る尺度であるが、ここでは平均は共通で共分散行列として ${\bf B}$ と ${\bf B}_{k}$ を用いた多次元正規分布同士のKLダイバージェンスを $D({\bf B}, {\bf B}_{k})$ とする。多次元正規分布のKLダイバージェンスについては以下が参考になる。

qiita.com

これを踏まえて $D({\bf B}, {\bf B}_{k})$ を式で表すと、

$\displaystyle \begin{aligned} D({\bf B}, {\bf B}_{k}) &= D_{KL} \left[{\mathcal N}({\bf \mu}, {\bf B}) \parallel {\mathcal N}({\bf \mu}, {\bf B}_{k}) \right] \\ &= {\rm tr}({\bf B} {\bf B}_k^{-1}) - \log \det ({\bf B} {\bf B}_k^{-1}) - n \end{aligned} \tag{11}$

となり、この式（11）を式（10）に当てはめて最適化問題を解くと式（8）が得られる*5。

ちなみに、KLダイバージェンスは対称性がないため、式（11）の ${\bf B}$ と ${\bf B}_{k}$ を入れ替えたKLダイバージェンス $D({\bf B}_{k}, {\bf B})$ を使うと別の更新式が得られ、こちらはDFP公式と呼ばれている。実はKLダイバージェンスに基づく非類似度には $D({\bf B}_{k}, {\bf B}) = D({\bf B}^{-1}, {\bf B}_{k}^{-1})$ という性質があり、式（10）の ${\bf B}^{-1}$ に関する最適化問題を解くことと等価になるため、その最適解である式（9）の ${\bf B}_{k}^{-1}$ を ${\bf B}_{k}$ に置き換えたものがDFP公式の更新式となる（従って、その逆行列の更新式は式（8）の ${\bf B}_{k}$ を ${\bf B}_{k}^{-1}$ に置き換えたものとなる）。

アルゴリズムの流れ

BFGS公式に基づく準ニュートン法のアルゴリズムは以下の通り。

近似行列の初期値 $({\bf B}_0)^{-1}$ *6と初期点 ${\bf x}^{(0)}$ を定める。
$\parallel \nabla f(x^{(k)}) \parallel$ が十分に小さければ終了。
式（5）で探索方向を求める。
直線探索で適切なステップ幅 $\alpha_k$ を求める。
式（2）で $x^{(k+1)}$ を求める。
式（9）で $({\bf B}_{k + 1})^{-1}$ を求める。
2.に戻る。

実装

以下が実装した準ニュートン法のpythonコード。

import numpy as np
import copy

# アルミホ条件を満たすαを線形探索
def backtracking_line_search(objective, x, d, tau1, beta):
    alpha = 1
    g = lambda alpha : objective.func(x + alpha * d)
    g_dash = lambda alpha : objective.grad(x + alpha * d) @ d    
    while(1):
        if ((g(alpha) <= g(0) + tau1 * g_dash(0) * alpha)):
                break
        alpha *= beta
    return alpha

# 目的関数のクラスobjectiveは以下のメンバを持っている必要あり
# objective.func(x) : return f(x)
# objective.grad(x) : return ∇f(x)
# objective.dim : 変数の次元
def quasi_newton(objective, eps = 1e-9, tau1 = 0.5, beta = 0.9):
    H = np.eye(objective.dim) # = inv(B)
    x = np.random.randn(objective.dim)
    
    while(1):
        grad = objective.grad(x)
        if(np.linalg.norm(grad) <= eps):
            break
        d = - H @ grad
        alpha = backtracking_line_search(objective, x, d, tau1, beta)
        x_new = x + alpha * d
        s = (x_new - x)[:, None]
        y = (objective.grad(x_new) - grad)[:, None]
        I = np.eye(H.shape[0])
        sy = (s.T @ y)
        H = (I - (s @ y.T) / sy)  @ H @ (I - (y @ s.T) / sy) + (s @ s.T / sy)
        x = x_new
    return x, objective.func(x)

引数objectiveに目的関数を表すクラスのインスタンスを入力するが、これは $f({\bf x})$ の値を返すメソッドfuncと $\nabla f({\bf x})$ を返すメソッドgradを実装してる必要がある。準ニュートン法の処理自体は素直に書いており特別なことはしていない。なお、脚注*3にあるように $({\bf B}_k)^{-1}$ の正定値性を保証するならウルフ条件を満たすように直線探索するべきだが、alphaの初期値次第ではウルフ条件を満たすalphaが見つからず無限ループに陥るなど扱いにくかったのでアルミホ条件のバックトラッキングで妥協した*7。

前回同様、上記のソースコードと検証用コードを以下のgithub リポジトリに上げている。

github.com

最後に

今回は制約なしの非線形計画問題を解く準ニュートン法を実装したが、次は制約ありの非線形計画問題を解く内点法を実装する予定。実は内点法の実装自体はもう出来てるので後は記事を書くだけ。。。

参考文献

今回の内容も主にこちらの書籍を参考にして書いた。
しっかり学ぶ数理最適化モデルからアルゴリズムまで (KS情報科学専門書)
- 作者:梅谷俊治
- 発売日: 2020/10/26
- メディア: 単行本
BFGS公式がKLダイバージェンスの最小化として導出されることは以下の書籍に書かれている。
機械学習のための連続最適化 (機械学習プロフェッショナルシリーズ)
- 作者:金森敬文,鈴木大慈,竹内一郎,佐藤一誠
- 発売日: 2017/09/08
- メディア: Kindle版

*1:実装ではアルミホ条件を満たすように、バックトラック法で直線探索している。

*2:逆に狭義凸関数なら正定値行列、凸関数なら半正定値行列になる。

*3: 参考文献である『しっかり学ぶ数理最適化』では $g(\alpha) = f({\bf x}^{(k)} + \alpha {\bf d}({\bf x}^{(k)}) )$ の $g'(\alpha) = 0$ を達成する $\alpha$ を選択したときの正定値性が示されているが、別の参考文献『機械学習のための連続最適化』ではウルフ条件を満たした $\alpha$ でも正定値になることが示されている。

*4: $\displaystyle({\bf A} + {\bf u}{\bf v}^T)^{-1} = {\bf A}^{-1} - \frac{{\bf A}^{-1} {\bf u}{\bf v}^T {\bf A}^{-1}}{1 + {\bf v}^T{\bf A}^{-1} {\bf u}}$