PythonでNelder-Meadアルゴリズムを用いたRosenbrock関数の最適化とそのgifアニメ化の方法

はじめに

Python を用いた Rosenbrock 関数の最適化をNelder-Meadアルゴリズムで行い、その過程をアニメーションで可視化する方法を紹介します。

Rosenbrock関数は、最適化アルゴリズムのベンチマークとしてよく使われる関数の一つです。今回はこの関数の最適化を行い、そのプロセスをアニメーションで可視化して、Nelder-Meadアルゴリズムの理解を深めてみましょう。

Google Colab にもコードを載せています。

ネルダー・ミード（Nelder-Mead）法の簡単な説明

ネルダー・ミード法（Nelder-Meadアルゴリズム）は、数値最適化問題を解くための手法で、特に微分不可能な関数に対して有効です。ジョン・ネルダーとロジャー・ミードによって1965年に開発されました。

このアルゴリズムは、シンプレックスと呼ばれる多面体（通常は三角形または単体）を使用して関数の最小値を探索します。シンプレックスは、探索空間内の頂点集合で構成され、各頂点は関数の評価点を表します。アルゴリズムは、シンプレックスの頂点を評価し、最も関数値が高い（最悪の）点を識別します。次に、最悪の点を除く他の点を基準にして、最悪の点を反射、拡大、収縮、または縮小することで新しい点を生成し、シンプレックスを更新します。

反射は、最悪の点を対称点に移動させる操作です。拡大は、反射がうまくいった場合に、さらにその方向へ進む操作です。収縮は、反射がうまくいかなかった場合に、最悪の点と他の点の中間に新しい点を作る操作です。縮小は、全体的な性能が悪い場合に、全ての点を最良の点に近づける操作です。

ネルダー・ミード法は勾配情報を必要としないため、微分不可能な問題やノイズの多い問題に適しています。しかし、この手法は局所最適解に収束する傾向があり、大規模な問題や複雑な最適化問題には効果が限られる場合があります。また、シンプレックスの初期配置やパラメータの選択が結果に大きく影響するため、注意深い設定が必要です。

このアルゴリズムは幅広い分野で利用されていて、解析的なアプローチが困難な複雑な問題や、導関数が利用できない場合に有効です。その単純さと適用範囲の広さから、多くの数値最適化ソフトウェアやライブラリで実装されています。

コードの流れ

アルゴリズムの流れについて説明します。

1. まず、Rosenbrock関数を定義しています。
2. 次に、Rosenbrock関数の可視化用のグリッドデータを作成しています。
3. コールバック関数を使用して、最適化の各ステップでの点を記録しています。
4. minimize関数を使用して、Nelder-Meadアルゴリズムを用いてRosenbrock関数の最小値を探索しています。初期のシンプレックスとして与えられた3つの点から始めます。
5. minimizer['allvecs']には、最適化の各ステップでのシンプレックスが格納されています。
6. プロットを初期化して、contourプロットを作成します。
7. アニメーションのための初期化関数とアニメーション関数を定義しています。アニメーション関数では、各ステップのシンプレックスの点をプロットし、その時点でのフレーム番号と最後の点の座標をテキストとして表示しています。
8. FuncAnimationを使用してアニメーションを作成し、GIFとして保存しています。

これがアニメーションを生成するまでの手順になります。このコードを実行すれば、Rosenbrock関数上でのNelder-Meadアルゴリズムの動作を示すGIFが生成されます。

ライブラリのインポート

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from matplotlib.animation import PillowWriter, FuncAnimation

Rosenbrock関数の定義

Rosenbrock関数は以下の式で定義されます：

\begin{equation}
f(x, y)=(1-x)^2+100\left(y-x^2\right)^2
\end{equation}

この関数の最小値は　f(x,y)=0 で、この値は x=1 および y=1 のときに達成されます。したがって、Rosenbrock関数の最小値は点 (1,1) であり、そのときの関数の値は0です。

def rosenbrock(x):
    return (1 - x[0])**2 + 100*(x[1] - x[0]**2)**2

Rosenbrock関数は、最適化アルゴリズムのテストやデモンストレーションのためによく使用されるベンチマーク関数の一つです。以下の理由から、この関数が選ばれることが多いです：

• 形状: Rosenbrock関数は狭い、湾曲した渓谷を持っています。この渓谷の中を探索しながら最小値を見つけるのは非常に挑戦的であり、多くの最適化アルゴリズムがこの関数上で効率的に動作するかを試すのに適しています。

• 非凸性: この関数は非凸であり、多くの局所的な最適解を持たないにも関わらず、その形状がアルゴリズムの探索戦略を試すのに適しています。

• 認識度: 由来や形状がよく知られており、多くの研究や文献で言及されています。したがって、この関数を使うことで、他の研究やデモとの比較が容易になります。

• 計算の簡単さ: Rosenbrock関数は計算が簡単で、微分も解析的に容易に計算できます。これにより、勾配ベースの方法と勾配フリーの方法の両方で使用することができます。

これらの理由から、Rosenbrock関数は最適化の分野でのデモンストレーションやベンチマークテストに広く使われています。

可視化のためのデータの準備

次に、関数の可視化のためのグリッドデータを作成します。

X = np.linspace(-2, 2, 400)
Y = np.linspace(-1, 3, 400)
X, Y = np.meshgrid(X, Y)
Z = (1 - X)**2 + 100*(Y - X**2)**2

最適化

points = []

def callback(xk):
    points.append(xk)

initial_simplex = np.array([[-1.0, 1.0], [-0.5, 1.6], [1.3, 1.5]])
minimizer = minimize(rosenbrock, initial_simplex[0], method='Nelder-Mead', callback=callback, options={'initial_simplex': initial_simplex, 'return_all': True})

アニメーションの作成

最後に、最適化の過程をアニメーションで可視化します。

fig, ax = plt.subplots(figsize=(8,6))
...
ani = FuncAnimation(fig, animate, frames=len(simplexes), init_func=init, blit=True, interval=interval)
ani.save('nelder_mead_rosenbrock.gif', writer=PillowWriter(fps=frame_per_sec))

interval と fps は両方ともアニメーションのタイミングに影響を与えます。

• interval
  • FuncAnimation の interval パラメータは、各フレームの間の時間をミリ秒で指定します。例えば、interval=1 は、各フレームの間に1ミリ秒の遅延を持たせます。interval が小さいほどアニメーションが速く、大きいほどアニメーションが遅くなります。
• fps (frames per second):
  • PillowWriter の fps パラメータは、保存されるGIFまたは動画のフレームレートを指定します。fps=1 は、1秒あたり1フレームの速度でアニメーションが進行することを意味します。
    fps が高いほどアニメーションが速く、低いほどアニメーションが遅くなります。

FuncAnimation での interval は、アニメーションが画面上で実際に再生されるときの速度を制御するのに対し、fps はアニメーションがGIFや動画として保存されたときの再生速度を制御します。一貫性を保つために、これらのパラメータが同じになるように設定しています。

FuncAnimation での blit は、アニメーションを描画する際の最適化手法の一つです。blit=True に設定すると、アニメーションの各フレームで変更される部分だけを再描画するようになります。これにより、全体の描画を毎回行うよりも効率的にアニメーションを更新することができます。

具体的には、fblit=True を使用する場合、次の条件が必要です：

• finit_funfc が提供されていること。この関数はアニメーションの背景を設定する役割を持っています。
• 各アニメーション関数（この場合は animate 関数）は更新されるアーティスト（描画オブジェクト）だけを返す必要があります。

一方、blit=False に設定すると、アニメーションの各フレームで全体が再描画されます。これは効率的ではありませんが、すべての環境やバックエンドで動作します。

注意点として、blit=True は一部のmatplotlibのバックエンド（特にベクトルベースのバックエンドや一部の非インタラクティブなバックエンド）ではサポートされていない場合があります。サポートされているバックエンドの中で最も一般的なものは、TkAggやAgg、GTK3Aggなどです。

要するに、blit を使用するとアニメーションの描画が高速化される可能性がありますが、環境やバックエンドによっては正しく動作しない場合があります。

生成されたアニメーション

生成されたアニメーションがこちらになります。

３角形を描きながら、少しづつ収束していることがわかります。

コード全文

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from matplotlib.animation import PillowWriter, FuncAnimation

# Rosenbrock関数の定義
def rosenbrock(x):
    return (1 - x[0])**2 + 100*(x[1] - x[0]**2)**2

# 可視化用のグリッドデータの作成
X = np.linspace(-2, 2, 400)
Y = np.linspace(-1, 3, 400)
X, Y = np.meshgrid(X, Y)
Z = (1 - X)**2 + 100*(Y - X**2)**2

# コールバック関数での各ステップでの点を記録
points = []

def callback(xk):
    points.append(xk)

# Nelder-Meadアルゴリズムを用いた最適化
initial_simplex = np.array([[-1.0, 1.0], [-0.5, 1.6], [1.3, 1.5]])
minimizer = minimize(rosenbrock, initial_simplex[0], method='Nelder-Mead', callback=callback, options={'initial_simplex': initial_simplex, 'return_all': True})

# 最適化結果を表示
print("Check fitting results")
print(minimizer)

# minimizer['allvecs']には全てのシンプレックスが格納されている
simplexes = minimizer['allvecs']

# 可視化のためのプロット設定
fig, ax = plt.subplots(figsize=(8,6))
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Nelder-Mead algorism using Rosenbrock function")

Z_min, Z_max = np.min(Z), np.max(Z)
levels = np.logspace(np.log10(Z_min), np.log10(Z_max), 35)
#ax.contour(X, Y, Z, levels=np.logspace(0, 5, 35), cmap='jet')
ax.contour(X, Y, Z, levels=levels, cmap='cool', linewidths=1, alpha=0.9)
line, = ax.plot([], [], lw=1, color="red")

def init():
    line.set_data([], [])
    return line,

# グローバル変数としてテキストオブジェクトを定義
text = ax.text(0.005, 1.06, '', transform=ax.transAxes)

def animate(i):
    if i < 3:  # 初期のシンプレックスを表示
        simplex = simplexes[:i+1]
    else:
        simplex = simplexes[i-2:i+1]
    
    x_data, y_data = zip(*simplex)
    
    # 最後の点の座標を取得
    last_x, last_y = x_data[-1], y_data[-1]
    
    # シンプレックスを閉じるための開始点を追加
    x_data += x_data[:1]
    y_data += y_data[:1]
    print(str(i) + "-th Iteratrion x_data = ", x_data)
    
    line.set_data(x_data, y_data)
    
    text.set_text(f"Frame: {i} (x,y)=({last_x:.3f},{last_y:.3f})")
    
    return line, text

frame_per_sec = 5 
interval = int(1e3/5) # ms for FuncAnimation

ani = FuncAnimation(fig, animate, frames=len(simplexes), init_func=init, blit=True, interval=interval)
ani.save('nelder_mead_rosenbrock.gif', writer=PillowWriter(fps=frame_per_sec))

Nelder-Meadアルゴリズムについての補足

Nelder-Mead アルゴリズムとは？

Nelder-Mead アルゴリズムは、最適化問題を解くための直接探索法の1つです。特に、このアルゴリズムは微分不要な最適化手法として知られ、非線形最適化問題に広く使用されています。

アルゴリズムの基本的なアイディアは、シンプレックスと呼ばれる点の集合を用いて関数の最小化を行うことです。シンプレックスは、n 次元空間において n+1 個の点からなる幾何的な形状を指します。

Nelder-Mead アルゴリズムは以下のステップで構成されます：

1. 初期シンプレックスの設定
2. シンプレックスの頂点を関数の値に基づいてソート
3. 新しい試行点を計算（これにはいくつかの手法、例えば反射、拡大、縮小、収縮などがあります）
4. 新しい試行点を使ってシンプレックスを更新
5. 終了基準が満たされるまで2-4のステップを繰り返す

このアルゴリズムの利点は、微分情報が不要であるため、非滑らかな関数や離散的な関数にも適用できることです。一方、収束速度が遅いという欠点もあります。

機械学習で使われるか？

機械学習において、Nelder-Mead アルゴリズムは主に以下のような場面で使用されることがあります：

• ハイパーパラメータの最適化: 一部の機械学習モデルは、良い性能を達成するために適切なハイパーパラメータの設定が必要です。Nelder-Mead アルゴリズムは、このようなハイパーパラメータの最適化問題に使われることがあります。特に、目的関数の勾配が利用できない、または計算が困難な場合に適しています。

• モデルの初期化: 一部の機械学習モデル（例：ニューラルネットワーク）は、初期の重みやバイアスの設定に敏感です。Nelder-Mead アルゴリズムを使って、良好な初期設定を見つけることが考えられます。

• 非微分問題: 一部の最適化問題では、目的関数が非滑らかで微分が困難な場合があります。Nelder-Mead のような微分不要のアルゴリズムは、このような問題に適しています。

しかし、Nelder-Mead アルゴリズムは大規模な機械学習の問題や深いニューラルネットワークの学習には通常適していません。その理由として、収束速度の遅さや多次元空間での効率性の問題が挙げられます。そのため、機械学習の文脈では、よりスケーラブルな最適化アルゴリズム（例：確率的勾配降下法やAdam）が一般的に使用されます。

微分困難でも勾配ベースで大丈夫な場合がある

勾配ベースの最適化アルゴリズムが微分が困難な場合や勾配が不正確な場合にもしばしば適切に動作する理由は以下です：

• 局所的近似: たとえ関数が微分不可能な点を持っていたとしても、その点の近傍で関数は局所的に滑らかであることが多い。このような場合、アルゴリズムはその近傍での勾配情報を利用して適切に動作することができる。

• ノイズの平均化: SGDやミニバッチSGDはランダムなサブセットのデータを使用して勾配を計算するため、ノイズが多い勾配もデータの異なるサブセットでの計算を通じて平均化される傾向があります。

• 正則化とモメンタム: 一部の最適化アルゴリズム（例：Adam）はモメンタムや正則化のテクニックを使用して、不正確な勾配情報やノイズの影響を軽減します。

• 適応的学習率: Adamのようなアルゴリズムは適応的な学習率を持っており、更新の振幅を調整することで不正確な勾配情報を緩和することができます。

• 実践的な成功: 実際のデータセットやタスクにおいて、完全に正確な勾配情報が必ずしも必要でないことが多い。関数が大まかな形状を持っていれば、勾配の方向は最適化の方向として十分なことが多い。

• 多様性の導入: ドロップアウト、データの拡張、ノイズの追加など、機械学習のトレーニング中に導入される正則化のテクニックは、モデルを頑健にするとともに、勾配の不正確さに対する耐性を高める効果があります。

これらの理由から、微分が困難な場面や勾配が不正確な状況でも、勾配ベースの最適化アルゴリズムが実際の機械学習のタスクで効果的に動作することが多いです。そのため、この記事で紹介したNelder-Meadアルゴリズムが、微分不能であれば必ず必要、という訳ではないことに注意していください。