[数理最適化][実装]Python で Nelder-Mead法

概要

• Nelder-Mead法は導関数を使わずに凸関数を最適化するアルゴリズムです．
• 本稿ではNelder-Mead法を実装しベンチマークテストを行います．
• 言語はPythonで実装します．
• 本稿では実装とテストに焦点を置いているため数学的な背景についてはあまり触れませんのでご容赦ください．
• 説明には個人的な解釈が多分に含まれますのでご注意ください．

Nelder-Mead法

• Nelder-Mead法は$N$次元空間上で$N+1$個の頂点からなる凸包（二次元の場合は三角形，三次元では四面体であり，ここでは単体：シンプレックスといいます）を動かすことで凸関数の最小値を探索するヒューリスティックアルゴリズムです．
• シンプレックスに対する操作はすべて線形であるため，導関数を用いず，関数の陽的または陰的な表式が得られていない場合（すなわち，最適化する対象がブラックボックス関数である場合）でも適用できます．

Nelder-Mead法の概要

• Nelder-Mead法では関数値が最も悪い（最小値探索の場合は，最も関数値が大きい）頂点（最悪点）を選び，シンプレックスの重心との外分/内分点に移動させる操作またはシンプレックスを相似形を保ったまま縮小する操作を繰り返します．

• 頂点と重心に対する外分/内分操作には，反射（REFLECTION），拡張（EXPANSION），外側収縮（OUTSIDE CONTRACT），内側収縮（INSIDE CONTRACT）の四種類があり，それぞれについて外分比/内分比が決まっています．

• 繰り返しにおいてはあらかじめ設定しておいたハイパーパラメータ$\alpha,\beta,\gamma,\delta$を用います．

外分/内分操作

1. 反射点　　　：最悪点と重心を$1+\alpha:\alpha$ $(0<\alpha)$ の比率で外分した点

2. 拡張点　　　：反射点と重心を$\beta:\beta+1$ （$0<\beta$）の比率で外分した点
   　　　　　　　すなわち，最悪点と重心を$1+\alpha(\beta+1):\alpha(\beta+1)$で外分した点

3. 外側収縮点　：反射点と重心を$1-\gamma:\gamma$ （$0<\gamma<1$）の比率で内分した点
   　　　　　　　すなわち，最悪点と重心を$1+\alpha\gamma:\alpha\gamma$で外分した点

4. 内側収縮点　：反射点と重心を$1+\gamma:\gamma$ （$0<\gamma<1$）の比率で外分した点
   　　　　　　　すなわち，最悪点と重心を$1-\alpha\gamma:\alpha\gamma$で内分（または外分）した点

• $N=2$のときのそれぞれの内分/外分点と最悪点，重心の位置関係は下図の通りです．
  
  
  　

$N>2$のときでもこれらの点はすべて一直線上に並びます．また，$\alpha\gamma>1$の場合，内側収縮点は最悪点と重心を外分し，画像で最悪点より上側に来ます．

　

• 反射点での関数値と最良点，第二最悪点の関数値の大小関係に応じて，拡張，外側収縮，内側収縮のいずれかを行います．
• これらの操作で第二最悪点より良い関数値を得られないときは，縮小操作を行います．

縮小操作

• 縮小操作では，シンプレックスの各頂点を最良点との$1-\delta:\delta$ （$0<\delta<1$）の内分点に移動します．

Nelder-Mead法のアルゴリズムと実装

アルゴリズム概要

• Nelder-Mead法のアルゴリズムは以下の通りです．

---

Step1.

事前にハイパーパラメータとして$\alpha\in(0,\infty)$，$\beta\in(0,\infty)$，$\gamma\in(0,1)$，$\delta\in(0,1)$をきめておきます．

Step2.

初期のシンプレックスの各頂点の位置をきめておきます．

Step3.

以下を停止条件を満たすまで繰り返します．

Step3-1.

シンプレックスの各頂点を関数値の大きい順にソートし，最良点での関数値$F_1$，第二最悪点，最悪点，重心を求めます．

シンプレックスの各頂点での関数値はすべて最悪点での関数値より良いため，最悪点から重心方向にむかって関数値が下がっていると仮定します．（もちろん必ずしもそうとは言えません）

Step3-2.

反射点を求め，反射点での関数値$F_r$を求めます． > 反射点での関数値は，最悪点から重心方向への”下り坂”がどこまで続いているかの指標になります．

Step3-3.

反射点での関数値が，

Case1.

最良点での関数値より良いとき

拡大点とそこでの関数値を求めます．
拡大点での関数値が反射点での関数値より良いとき
　->最悪点を拡大点に移動します．
そうでないとき
　->最悪点を反射点に移動します．

このとき，下り坂は反射点にまで達していると考え，さらに先にある拡大点まで下り坂が続いているか調べています．

Case2.

最良点での関数値より悪いが第二最悪点より良いとき

最悪点を反射点に移動します．

このとき，下り坂は反射点あたりで終わると考え，反射点への移動で手を打ちます．

Case3.

第二最悪点での関数値より悪いが最悪点より良いとき

外側収縮点とそこでの関数値を求めます．
外側収縮点での関数値が反射点での関数値より良いとき
　->最悪点を外側収縮点に移動します．
そうでないとき
　->縮小します．

このとき，重心方向に下り坂になっているが反射点まで続いていないと考え，反射点より手前の外側収縮点を調べています．外側収縮点が悪いときに収縮するのは，最小値がより近場にあると考え，探索範囲を狭めるためです．

Case4.

最悪点での関数値より悪いとき

内側収縮点とそこでの関数値を求めます．
内側収縮点での関数値が反射点での関数値より良いとき
　->最悪点を内側収縮点に移動します．
そうでないとき
　->縮小します．

このとき，重心方向に下り坂になっていないか，なっていても反射点まで続いていないと考え，重心より手前の内側収縮点を調べています．内側収縮点が悪いときに収縮するのは，最小値がより近場にあると考え，探索範囲を狭めるためです．

---

疑似コード

• 頂点を大きさ$N+1$の配列$X$，各頂点での関数値を大きさ$N+1$の配列$F$でもちます．
• 最小値の更新幅や最後に最小値が更新されてからのイテレーション回数などで停止条件$\textrm{condition}$を定めています．

• 関数値の配列$F$のソートにかかる時間は$O(N\log N)$です．
• $N$が大きい場合は，頂点の配列と関数の配列をまとめて平衡二分木でもつことで縮小以外の操作を$O(\log N)$に短縮できます．

実装

Pythonによる実装

def NelderMead(f,X,n,alpha,beta,gamma,delta):
    hst=[]
    F=np.apply_along_axis(f,1,X)
    for i in range(n):
        asF=F.argsort()
        X=X[asF,:]
        F=F[asF]
        hst.append(X.copy())
        x0=X[:-1,:].mean(axis=0)
        xr=(1+alpha)*x0-alpha*X[-1,:]
        fr=f(xr)
        if (fr<F[0]):
            xe=(1+beta)*xr-beta*x0
            fe=f(xe)
            if (fe<fr):
                X[-1]=xe
                F[-1]=fe
            else:
                X[-1]=xr
                F[-1]=fr
        elif (F[0]<=fr<F[-2]):
            X[-1]=xr
            F[-1]=fr
        elif (F[-2]<=fr<F[-1]):
            xoc=(1-gamma)*x0+gamma*xr
            foc=f(xoc)
            if (foc<=fr):
                X[-1]=xoc
                F[-1]=foc
            else:
                for i in range(1,len(X)):
                    X[i]=(1-delta)*X[0]+delta*X[i]
                    F[i]=f(X[i])
        else:
            xic=(1+gamma)*x0-gamma*xr
            fic=f(xic)
            if (fic<=fr):
                X[-1]=xic
                F[-1]=fic
            else:
                for i in range(1,len(X)):
                    X[i]=(1-delta)*X[0]+delta*X[i]
                    F[i]=f(X[i])
    hst.append(X)
    return X[np.argmin(F),:],hst

ベンチマークテスト

• 先ほど実装したNelder-Mead法のPythonコードを，N=2でテストします．
• $x=y=0$に極（小）値を１つだけ持つ関数$f(x,y)=\exp(-x^2/\sigma_1^2-y^2/\sigma_2^2)$でテストします．
• 初期化されたシンプレックスの内側に最小値が含まれている場合と，含まれていない場合の二通りを試します．

X=np.array([[-0.35,0.1],[0.1,-0.55],[0.55,0.35]])
func=lambda x,y:-np.exp(-(x/0.5)**2-(y/0.3)**2)
x,hst=NelderMead(lambda x:func(*x),X,20,1,1,0.5,0.5)
def plot(i,ax):
    if int(i/5)<len(hst):j=int(i/5)
    else:j=-1
    x,y=hst[j].T
    z=func(x,y)
    ax.scatter(x,y,z,c="#0000FF")
    ax.plot([x[0],x[1]],[y[0],y[1]],[z[0],z[1]],color="#000000")
    ax.plot([x[0],x[2]],[y[0],y[2]],[z[0],z[2]],color="#000000")
    ax.plot([x[1],x[2]],[y[1],y[2]],[z[1],z[2]],color="#000000")

• プロットすると以下のようになりました．

• 最小値に向かってシンプレックスが収縮を繰り返していく様子がわかります．
• この関数では最小値から離れるほど関数値が大きくなるので，ほとんどの場合内側収縮か外側収縮が行われます．
• シンプレックスが反射または拡張して重心位置を並進させる様子を見るには，初期シンプレックスの内側に最小値が含まれないようにする必要があります．