クロスエントロピー法：モンテカルロ最適化の実践的手法

import numpy as np
import matplotlib.pyplot as plt

# --- Rastrigin関数の定義 ---
def rastrigin(x):
    """Rastrigin関数（最小化対象）"""
    A = 10
    return A * len(x) + np.sum(x**2 - A * np.cos(2 * np.pi * x))

# --- クロスエントロピー法の実装 ---
def cross_entropy_method(
    objective_fn,   # 目的関数（最小化）
    dim,            # 探索空間の次元数
    n_samples=100,  # 各イテレーションでのサンプル数
    elite_frac=0.2, # エリートサンプルの割合
    n_iterations=50,# 最大イテレーション数
    initial_mean=None,  # 初期平均
    initial_std=5.0,    # 初期標準偏差
):
    """クロスエントロピー法による最小化"""

    # ステップ1: 初期化
    mean = initial_mean if initial_mean is not None else np.zeros(dim)
    std = np.full(dim, initial_std)
    n_elite = int(n_samples * elite_frac)

    # 記録用リスト
    best_scores = []      # 各イテレーションの最良スコア
    mean_history = []     # 平均の推移

    for iteration in range(n_iterations):
        # ステップ2: サンプリング（正規分布から）
        samples = np.random.normal(
            loc=mean, scale=std, size=(n_samples, dim)
        )

        # ステップ3: 評価（各サンプルの目的関数値を計算）
        scores = np.array([objective_fn(s) for s in samples])

        # ステップ4: エリートサンプルの選択（スコアが小さい上位P%）
        elite_indices = np.argsort(scores)[:n_elite]
        elite_samples = samples[elite_indices]

        # ステップ5: パラメータの更新（エリートサンプルの平均・標準偏差）
        mean = np.mean(elite_samples, axis=0)
        std = np.std(elite_samples, axis=0)

        # 記録
        best_scores.append(np.min(scores))
        mean_history.append(mean.copy())

        # 収束判定（標準偏差が十分小さくなったら終了）
        if np.all(std < 1e-6):
            print(f"反復 {iteration + 1} で収束しました。")
            break

    return mean, best_scores, mean_history

# --- 実行 ---
np.random.seed(42)
dim = 2  # 2次元問題

best_solution, best_scores, mean_history = cross_entropy_method(
    objective_fn=rastrigin,
    dim=dim,
    n_samples=200,
    elite_frac=0.1,
    n_iterations=100,
    initial_mean=np.array([3.0, -3.0]),  # 最適解から離れた初期点
    initial_std=3.0,
)

print(f"最適解: {best_solution}")
print(f"目的関数値: {rastrigin(best_solution):.6f}")

# --- 結果の可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# (a) 目的関数値の収束推移
axes[0].plot(best_scores, linewidth=2)
axes[0].set_xlabel("Iteration")
axes[0].set_ylabel("Best Score")
axes[0].set_title("Convergence of CEM")
axes[0].set_yscale("log")
axes[0].grid(True, alpha=0.3)

# (b) Rastrigin関数の等高線と探索経路
x_grid = np.linspace(-5, 5, 200)
y_grid = np.linspace(-5, 5, 200)
X, Y = np.meshgrid(x_grid, y_grid)
Z = np.array([
    [rastrigin(np.array([xi, yi])) for xi in x_grid]
    for yi in y_grid
])

axes[1].contourf(X, Y, Z, levels=30, cmap="viridis", alpha=0.8)
axes[1].colorbar = plt.colorbar(axes[1].contourf(X, Y, Z, levels=30, cmap="viridis", alpha=0.8), ax=axes[1])

# 平均の推移を矢印付きでプロット
mean_arr = np.array(mean_history)
axes[1].plot(mean_arr[:, 0], mean_arr[:, 1], "r.-", markersize=8, linewidth=1.5, label="Mean trajectory")
axes[1].plot(mean_arr[0, 0], mean_arr[0, 1], "rs", markersize=12, label="Start")
axes[1].plot(mean_arr[-1, 0], mean_arr[-1, 1], "r*", markersize=15, label="Final")
axes[1].plot(0, 0, "wx", markersize=12, markeredgewidth=3, label="Global optimum")
axes[1].set_xlabel("$x_1$")
axes[1].set_ylabel("$x_2$")
axes[1].set_title("CEM Search Trajectory on Rastrigin Function")
axes[1].legend(loc="upper right")

plt.tight_layout()
plt.savefig("cem_result.png", dpi=150, bbox_inches="tight")
plt.show()