Nelder Mead Method を python で実装

最近会社で Kaggle 本 Kaggleで勝つデータ分析の技術 の本読みをやっています。昨日は第二章でしきい値の最適化の文脈で Nelder Mead Method を使った最適化が記載されていました。ちょっと気になったので実装してみたという話です。

Nelder Mead Method とは

最適化手法の一種です。いくつかのデータ点(単体)での評価値をベースに更新していく方法で、勾配の情報などを使わないため目的関数さえ評価できれば実行できるお手軽さがメリットです。

アルゴリズムの概要は英語版 wikipedia が詳しいです。(自分が見た限りでは日本語版では一部の更新で符号がまちがっているように見受けられたので英語のほうがよさ気)

en.wikipedia.org

前提として $N$ 次元のベクトルの最適化をしたいとすると

1. 適当に選んだ $N$ 次元ベクトル + 各次元に $+1$ した合計で $N+1$ の $N$ 次元ベクトルを計算
2. 1で計算した各点での目的関数値を計算して、その順序を保存。
3. 一番悪かった点以外の重心を計算
4. 一番悪かった点と重心の対称な点(反射点)を計算

この反射点を使って場合分けを計算していきます。

• 反射点が二番目に悪い点より良い場合は採用として一番悪い点を置き換え
• 反射点が一番良い場合、もっとソッチの方に移動しても良いという発想で、最悪点から反対側に更に移動した点(拡大点)を計算。
  拡大点が反射点より良い場合には拡大点を採用として一番悪い点を置き換えます。(拡大点が反射点より悪い時はすなおに反射点で妥協)
• それ以外の時(2番めのものより悪い時)は重心と最悪点の内点(縮小点)を計算
  • 縮小点が最悪よりは良い時は最悪点を縮小点で置き換え。
  • 縮小点も悪い時は最良点以外の点を、最良点との内点に置き換え

python のコードにすると以下のような感じです。Notebook はこちら https://nyker510.gitlab.io/notebooks/posts/nelde-mead/ にあげています。

import numpy as np

init_point = np.random.uniform(-.1, .1, 2) - np.array([1.5, 1.5])

alpha = 1.
gamma = 2.
rho = .5
sigma = .5
scaling_factor = 1.

eval_points = np.vstack([init_point + np.eye(2) * scaling_factor, init_point])
best_points = [init_point]

for i in range(100):
    # 目的関数で評価して並び替え
    objective_values_by_point = objective(eval_points[:, 0], eval_points[:, 1])
    orders = np.argsort(objective_values_by_point)
    f_values = objective_values_by_point[orders]
    eval_points = eval_points[orders]

    # 最悪点以外で重心を計算
    centroid = np.mean(eval_points[:-1], axis=0)
    p_worst = eval_points[-1]

    # 反射点 (最悪点から重心の対象な点) を計算 + そのときの目的関数を評価
    p_reflect = centroid + alpha * (centroid - p_worst)
    f_ref = objective(*p_reflect)

    f_best = f_values[0]
    f_second_worse = f_values[-2]
    f_worst = f_values[-1]

    # 反射点がそこそこ良い場合 (second worse よりは良い) → 最悪点を反射点で置き換え
    if f_best <= f_ref < f_second_worse:
        eval_points[-1] = p_reflect

    # 反射点が一番良い場合
    elif f_ref < f_best:
        # もうちょっと良い点がないか探りに行くイメージ (拡大点)
        # 最悪点が一番悪い → 反射点の方にもっと移動するイメージ (反射点 + centroid - p_worst [i.e. 最悪点からみた重心方向])
        p_expand = centroid + gamma * (centroid - p_worst)
        f_expand = objective(*p_expand)

        # 拡大点が良かったら採用. だめだったら反射点で妥協
        if f_expand < f_ref:
            eval_points[-1] = p_expand
        else:
            eval_points[-1] = p_reflect

    # 反射点が二番目に悪い点より良くない時
    elif f_second_worse <= f_ref:
        # 重心と最悪点の間を計算 (大体すべての点のまんなかあたりに対応する)
        p_cont = centroid + rho * (p_worst - centroid)
        f_cont = objective(*p_cont)

        # 縮小点が最悪よりは良い時
        # 重心方向への移動は良いと判断して最悪点を置き換え
        if f_cont < f_worst:
            eval_points[-1] = p_cont
        else:
            # そうでない時最良点めがけて縮小させる (最良点との内点に移動して三角形ちっちゃくするイメージ)
            shrink_points =  eval_points[0] + sigma * (eval_points - eval_points[0])
            eval_points[1:] = shrink_points[1:]
            
    best_points = np.vstack([best_points, eval_points[0]])

実験

はじめに自分で適当に作った関数でやります。たぶん (0, -.4) あたりが最適解なはず(あやしい)。

def objective(x, y):
    return (x - .1) ** 2 + (y - .3) ** 2 + np.sin((x + .5) * 2) * np.sin(y * 3) + sigmoid(x * 4)

これでやると以下のような感じ。わりとよさ気。

Iteration ごとの座標

ついでに Rosenbrock function でもやってみましょう。 こちらは最適な座標は (1, 1) になるように設定しています(見えにくいですが黒い丸があるところ)。

Rosenbrock function

実行すると以下のような感じ。たどり着いてるっぽいですね。

Rosenbrock function > Iteration ごとの座標

NOTE

これは Nelder Mead に限った話ではないですが、問題のスケールをちゃんと -1, 1 ぐらいに合わせておく必要があります。というのは一番初めに生成する点列の評価値に依存してその後の点列が決定するため、一番初めの点の性質が悪いと上手く更新が進まないためです。 例えば今回で言うと init_point に対して np.eye(1) を足すことで初期の点を生成していますが、これを np.eye(1) * .00001 などにすると上手く収束してくれません。

初期点の探索スケールが小さすぎる時の path

反対に大きすぎても更新幅が大きすぎて発散気味になります。

初期点の更新スケールが大きすぎるときの path

これと同様にそもそもの初期点をどう選ぶかもとても重要です。コード中では np.random.uniform(-.1, .1, 2) - np.array([1.5, 1.5]) としていますが、たとえばこの周りだけ凹んでいる関数だった場合そこから抜け出せずに局所解に陥って終わりになることもあります。ちょっとうまくいかないなーという時はスケールがあっているかを確認+初期点を変えてやってみるのも大切そうです。