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

クロスエントロピー法 (Cross Entropy Method, CEM) は、モンテカルロ法における重点サンプリング (Importance Sampling) の一種であり、最適化問題や希少事象の確率推定に用いられるアルゴリズムです。

モンテカルロサンプリング

モンテカルロ法は、乱数を用いることで、複雑な問題の近似解を求める数値計算手法の総称です。

例えば、確率変数 \(X\) が確率密度関数 \(f(x)\) に従うとき、関数 \(H(X)\) の期待値 \(l = \mathbb{E}[H(X)] = \int H(x)f(x)dx\) を求めたいとします。 モンテカルロサンプリングでは、この期待値を \(N\) 個の独立同分布なサンプル \(X_1, \dots, X_N\) を用いて以下のように近似します。

重点サンプリング (Importance Sampling)

モンテカルロサンプリングでは、期待値を計算したい確率分布 \(f(x)\) から直接サンプルを生成します。しかし、もし \(H(x)\) が大きな値をとる領域が \(f(x)\) の確率が低い領域にある場合、効率的なサンプリングができません。

重点サンプリングでは、目的の分布 \(f(x)\) ではなく、別のサンプリング分布 \(g(x)\) からサンプルを生成します。そして、そのサンプルに重み \(w(x) = f(x)/g(x)\) を掛けることで、期待値を推定します。

ここで、\(X_i \sim g(x)\) です。重点サンプリングの効率は、サンプリング分布 \(g(x)\) の選択に大きく依存します。特に、分散を最小化する最適なサンプリング分布 \(g^*(x)\) は、\(g^*(x) \propto |H(x)|f(x)\) の形をとることが知られています。以下では、この事実を分散の式から厳密に導出し、CEMがどのようにこの \(g^*(x)\) に近づく手続きなのかを見ていきます。

重点サンプリング推定量の分散の導出

\(l_{IS}\) は \(N\) 個の独立同分布なサンプルの平均なので、その分散は1サンプルあたりの分散を \(N\) で割った

で与えられます。期待値 \(\mathbb{E}_g[\cdot]\) を \(g(x)\) に関する積分に書き直すと、

となります。\(l^2\) は \(g\) に依存しない定数なので、分散を最小にする \(g\) を求める問題は、次の制約付き最小化問題に帰着します。

ラグランジュ未定乗数 \(\lambda\) を導入し、被積分関数を \(g(x)\) で汎関数微分してゼロとおくと、

が得られ、これが冒頭で述べた最適サンプリング分布 \(g^*(x)\) の形です。特に \(H(x) \geq 0\) が恒等的に成り立つ場合（\(H\) が確率を表す指示関数のときなど）、正規化定数を \(l = \mathbb{E}_f[H(X)]\) とすると \(g^*(x) = H(x)f(x)/l\) となり、このとき

は \(x\) によらず定数 \(l\) になります。定数の分散は \(0\) なので、\(\mathrm{Var}_{g^*}[l_{IS}] = 0\) が成立します。つまり理想的な重点サンプリング分布 \(g^*\) を使えば分散ゼロで真の期待値が求まりますが、\(g^*\) 自体の定義に未知の \(l\) が含まれるため、実際には直接使うことができません。CEMは、この \(g^*\) に近い分布を、扱いやすいパラメトリックな分布族 \(g(x;v)\) の中から反復的に探し出すアルゴリズムです。

希少事象確率推定から最適化への転用

CEMはもともと、ある関数 \(S(x)\) がしきい値 \(\gamma\) を超える希少事象の確率 \(l = \mathbb{E}_f[\mathbb{1}\{S(X)\geq\gamma\}] = P_f(S(X)\geq\gamma)\) を推定するために提案されました。ここで \(H(x) = \mathbb{1}\{S(x)\geq\gamma\}\) です。\(\gamma\) が大きいほど事象は希少になり、単純モンテカルロサンプリングでは \(N\) 個のサンプル中にほとんど事象が発生せず、推定量の相対誤差（標準偏差 \(/\) 推定値）が発散してしまいます。重点サンプリングにより、事象が起きやすい領域に確率質量を集めた \(g(x;v)\) からサンプリングすることで、少ないサンプル数でも安定に \(l\) を推定できます。

最適化への転用は、しきい値 \(\gamma\) の役割を変えることで実現されます。希少事象推定では \(\gamma\) は固定された定数でしたが、最適化では \(\gamma_t\) を各イテレーションでサンプルの上位 \(\rho\) 分位点として適応的に決定します。\(\gamma_t\) を「その時点で見つかっている良い値」に反復的に更新し続けることで、サンプリング分布 \(g(x;v_t)\) は目的関数の最適値の近傍に押し込まれていきます。つまりCEMの最適化バージョンは、固定した希少事象を狙うのではなく、動く標的（適応的なエリート閾値）に対する重点サンプリング分布の逐次更新とみなすことができます。次節では、この「エリート閾値を使った分布更新」がクロスエントロピー最小化としてどう定式化されるかを導出します。

クロスエントロピー法 (Cross Entropy Method)

クロスエントロピー法は、この重点サンプリングにおける最適なサンプリング分布 \(g^*(x)\) に近い分布を、反復的に見つけるためのアルゴリズムです。

「クロスエントロピー」は、2つの確率分布 \(p\) と \(q\) の間の類似度を測る尺度の一つで、以下のように定義されます。

KLダイバージェンス \(KL(p||q) = \int p(x) \log \frac{p(x)}{q(x)} dx = H(p,q) - H(p)\) の関係から、真の分布 \(p\) が固定されている場合、KLダイバージェンスを最小化することはクロスエントロピーを最小化することと等価になります。

CEMは、パラメータ化されたサンプリング分布 \(g(x;v)\) のパラメータ \(v\) を、以下の手順で最適化します。

1. 初期化: パラメータ \(v\) を初期化し、サンプリング分布 \(g(x;v)\) を設定します。
2. サンプリング: 現在のサンプリング分布 \(g(x;v)\) から \(N\) 個のサンプル \(X_1, \dots, X_N\) を生成します。
3. 評価: 各サンプル \(X_i\) に対して、目的関数 \(H(X_i)\) を計算します。
4. エリートサンプルの選択: 目的関数値が高い（または低い、最適化の方向による）上位 \(P\) パーセントのサンプルを「エリートサンプル」として選択します。
5. パラメータの更新: エリートサンプルを用いて、サンプリング分布 \(g(x;v)\) のパラメータ \(v\) を更新します。この更新は、エリートサンプルの尤度を最大化するように行われます。これは、エリートサンプルの経験分布と \(g(x;v)\) の間のクロスエントロピーを最小化することに相当します。
6. 収束判定: パラメータ \(v\) の変化が十分に小さくなるか、最大反復回数に達するまで、ステップ2に戻って繰り返します。

この反復プロセスにより、サンプリング分布 \(g(x;v)\) は徐々に最適解の領域に集中していき、効率的な探索が可能になります。

エリート分位点による分布更新の導出

上記の手順4・5「エリートサンプルを選び、その平均・分散でパラメータを更新する」は天下り的に見えますが、実際にはクロスエントロピー最小化問題を解いた結果として厳密に導出できます。

目標分布を、現在のサンプリング分布 \(g(x;v_{t-1})\) に、エリート集合を表す指示関数 \(\mathbb{1}\{H(x) \geq \gamma_t\}\) （最小化問題なら \(\mathbb{1}\{H(x) \leq \gamma_t\}\) ）を掛けて正規化したもの

と定義します。CEMは、この \(g^*_{\gamma_t}\) とパラメトリックな分布族 \(g(x;v)\) とのクロスエントロピー（前節で述べた通り、\(p\) が固定されていればKLダイバージェンス最小化と等価）を最小化する \(v\) を求めます。

\(g^*_{\gamma_t}\) は \(g(x;v_{t-1})\) を台とする（それに指示関数を掛けただけの）分布なので、この積分は \(X \sim g(\cdot;v_{t-1})\) についての期待値として書き直せ、経験分布による近似（現在のサンプル集合上でのモンテカルロ近似）で置き換えられます。

指示関数 \(\mathbb{1}\{H(X_i)\geq\gamma_t\}\) が \(1\) になるサンプルこそが「エリートサンプル」なので、この最大化問題はエリートサンプルの対数尤度を最大化すること、つまりエリートサンプルに対する最尤推定（MLE）そのものです。

\(g(x;v) = \mathcal{N}(x;\mu,\Sigma)\) を対角共分散のガウス分布族とすると、エリート集合 \(\mathcal{E} = \{X_i : H(X_i)\geq\gamma_t\}\) （要素数 \(n_{elite}\) ）に対する対数尤度は

であり、これを \(\mu\) 、\(\Sigma\) について微分してゼロとおくと、ガウス分布のMLEの標準的な結果として

が得られます。これはまさに、後述の実装コードの mean = np.mean(elite_samples, axis=0) と std = np.std(elite_samples, axis=0) に対応しています。つまり「エリートサンプルの平均・分散で更新する」という一見素朴な操作は、クロスエントロピー最小化＝エリート集合に対するガウス分布の最尤推定として厳密に正当化されることが分かります。

なお、\(X_i \sim g(\cdot;v_{t-1})\) を提案分布としてそのまま使い、重み \(f(x)/g(x;v_{t-1})\) を陽に掛けない（\(=1\) とみなす）近似は、CEMの標準的な定式化における簡略化です。前節の重点サンプリングの枠組みに厳密に従う一般形（尤度比重みを明示的に使う版）も存在しますが、実用上はこの簡略化で十分機能することが経験的に知られています。

CEMとMPPIの比較

CEMと密接に関連する手法として MPPI（Model Predictive Path Integral） があります。両者は重点サンプリングに基づく最適化手法ですが、サンプルの重み付け方法が異なります。

CEMの「上位\(P\) %を選ぶ」操作は、MPPIの温度パラメータ\(\lambda \to 0\) の極限に対応します。つまり、CEMはMPPIの特殊ケースとして理解できます。

実用的な使い分けの指針：

• 離散的な最適化や組合せ最適化にはCEMが適している
• 連続的な制御問題やロボティクスにはMPPIが適している
• MPPIはサンプル効率が高い（全サンプルの情報を使うため）が、温度パラメータの調整が必要

PythonによるCEMの実装

ここでは、クロスエントロピー法を用いて連続最適化問題を解くPython実装を紹介します。テスト関数として、多数の局所最小値を持つRastrigin関数を使用します。

Rastrigin関数

Rastrigin関数は、最適化アルゴリズムのベンチマークとして広く使われるテスト関数で、以下のように定義されます。

ここで \(A = 10\) 、\(n\) は次元数です。大域的最小値は \(\mathbf{x}^* = \mathbf{0}\) で \(f(\mathbf{0}) = 0\) です。多数の局所最小値があるため、単純な勾配法では大域的最小値に到達することが困難であり、CEMのような確率的探索手法の有効性を確認するのに適しています。

実装コード

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()

上記のコードを実行すると、以下のことが確認できます。

• 収束の様子: 目的関数値が反復ごとに急速に減少し、大域的最小値 \(f(\mathbf{0}) = 0\) に近づく
• 探索経路: サンプリング分布の平均が初期点 \((3, -3)\) から原点 \((0, 0)\) へ移動していく過程
• 局所最小値の回避: Rastrigin関数の多数の局所最小値を避け、大域的最小値に到達できている

CEMの主要なハイパーパラメータの影響：

• n_samples（サンプル数）：多いほど探索が安定するが、計算コストが増加する
• elite_frac（エリート割合）：小さいほど選択圧が強くなり収束が速いが、早期収束のリスクがある
• initial_std（初期標準偏差）：大きいほど広い範囲を探索できるが、収束に時間がかかる

エリート分位点の選び方：収束速度と多様性のトレードオフ

前節で述べた通り、エリート分位点 \(\rho\) （elite_frac）はCEMの収束速度と探索の多様性を直接トレードオフする最重要ハイパーパラメータです。\(\rho\) が小さいほど選択圧が強く、分布は急速に収縮しますが、収縮が速すぎると局所解の周辺で標準偏差が早期にゼロへ近づいてしまい（早期収束, premature convergence）、大域的最適解を見逃すリスクが高まります。逆に \(\rho\) が大きいと劣ったサンプルの情報も平均・分散の計算に混ざるため、分布の収縮が遅くなり、より多くのイテレーションを要します。

この効果を実測するため、10次元のRastrigin関数（局所最小値が非常に多く、早期収束の影響が顕著に現れる設定）を用いました。サンプル数を n_samples=100 に抑え（次元数に対してサンプルがやや少ない、選択圧の影響が出やすい条件）、エリート分位点 \(\rho \in \{0.02, 0.05, 0.1, 0.2, 0.3, 0.5\}\) ごとに30個の乱数シードで実行し、収束までの平均反復回数と、最終的な目的関数値が \(5.0\) を上回る（＝大域的最適解の谷に到達できず局所解で停滞した）「早期収束」の発生率を測定しました。

import numpy as np

def rastrigin(x):
    A = 10
    return A * len(x) + np.sum(x**2 - A * np.cos(2 * np.pi * x))

def cem_diag(dim, n_samples, elite_frac, n_iterations, initial_mean, initial_std, seed):
    """収束判定付きのCEM（対角共分散）。KLダイバージェンス追跡にも流用する。"""
    rng = np.random.default_rng(seed)
    mean = initial_mean.copy()
    std = np.full(dim, initial_std)
    n_elite = max(1, int(n_samples * elite_frac))

    best_scores = []
    for it in range(n_iterations):
        samples = rng.normal(loc=mean, scale=std, size=(n_samples, dim))
        scores = np.array([rastrigin(s) for s in samples])
        elite_idx = np.argsort(scores)[:n_elite]
        elite = samples[elite_idx]

        mean = np.mean(elite, axis=0)
        std = np.maximum(np.std(elite, axis=0), 1e-12)

        best_scores.append(np.min(scores))
        if np.all(std < 1e-6):
            break

    return {"mean": mean, "n_iter_used": len(best_scores), "final_score": rastrigin(mean)}

elite_fracs = [0.02, 0.05, 0.1, 0.2, 0.3, 0.5]
n_seeds = 30
for ef in elite_fracs:
    iters, finals, stuck = [], [], 0
    for seed in range(n_seeds):
        res = cem_diag(
            dim=10, n_samples=100, elite_frac=ef, n_iterations=200,
            initial_mean=np.full(10, 3.0), initial_std=3.0, seed=seed,
        )
        iters.append(res["n_iter_used"])
        finals.append(res["final_score"])
        if res["final_score"] > 5.0:
            stuck += 1
    print(f"elite_frac={ef:.2f}: mean_iters={np.mean(iters):.1f} "
          f"median_final={np.median(finals):.2f} stuck(>5.0)={stuck}/{n_seeds}")

実行結果は次の通りです。

\(\rho=0.02\) では平均18.7反復という最速で収束しますが、30シード全て（100%）が局所解で停滞しました。\(\rho\) を大きくするにつれて早期収束率は下がっていき（\(\rho=0.5\) では3%）、その代わり収束までの反復回数は125.1回まで増加します。これは \(\rho=0.02\) の約6.7倍の計算コストです。

下図(a)は、代表的な1シードにおける収束曲線（対数スケール）です。\(\rho=0.02\) （赤線）は約20反復で目的関数値が横ばいになり、局所解に「ロックイン」してしまうのに対し、\(\rho=0.5\) （青線）は緩やかにではありますが、より低い値まで収束を続けていることが確認できます。

実用上は \(\rho \in [0.1, 0.2]\) 程度が速度と信頼性のバランスとしてよく使われますが、目的関数の多峰性が強い（局所解が多い）ほど \(\rho\) を大きく取る、あるいは複数の初期値からの再起動（multi-start）を併用することが推奨されます。

高次元における次元の呪い

CEMのパラメータ更新は、共分散行列 \(\Sigma\) をフルに推定する場合、その要素数が次元 \(d\) に対して \(O(d^2)\) で増加します（本記事の実装のように対角成分に制限しても \(O(d)\) です）。共分散を精度良く推定するには次元数に応じてサンプル数 \(N\) を増やす必要があり、これが高次元空間でCEMのサンプル効率が悪化する「次元の呪い」の本質です。

実際に、n_samples=200 を固定したまま次元数を変化させ、100反復後の目的関数値（次元で正規化した「次元あたりスコア」）を10シードの平均で測定しました（コードは前節の cem_diag を dim を変えて呼び出すだけなので省略します）。

次元あたりスコアは次元数にほぼ比例して増加しており（2次元で0.05、100次元で1.36と約27倍）、同じサンプル数・反復回数では高次元ほど探索が相対的に不十分になることが確認できます。上図(b)にこの傾向を示しました。

この悪化を補うために必要なサンプル数の増加も実測しました。「次元あたりスコアが0.5未満」という固定の目標水準に到達するために必要な最小サンプル数を、n_samples を {50, 100, 200, 400, 800, 1600, 3200} の中から探索し、5シードの平均で求めた結果です。

次元が2から50へ25倍になると、目標水準を達成するために必要なサンプル数は100から400へと4倍に増加しています。フル共分散行列を毎回最尤推定するCEMの亜種では、良条件の推定に必要なサンプル数が理論上 \(d^2\) のオーダーで増加することが知られており、それと比べれば対角共分散に制限した本実装の悪化速度は緩やかです。しかしそれでも、サンプル効率の低下は無視できません。

実務的な対処法は次の通りです。

• 対角共分散への制限（本記事の実装のように、次元間の相関を無視して各次元を独立に更新する）は、パラメータ数を \(O(d^2)\) から \(O(d)\) に削減し、高次元でも比較的頑健に動作させる標準的なテクニックです。
• サンプル数を次元に応じてスケールさせる（経験則として \(N \gtrsim 4d\) 程度が目安とされることがあります）。
• 共分散の下限（正則化）を設定し、特定の次元の分散が他より先に潰れてしまうことで探索が偏るのを防ぐ。
• 相関が強い問題では対角制限では表現力が不足するため、低ランク共分散近似（例: \(\Sigma = \sigma^2 I + UU^\top\) ）などの中間的な表現を使う拡張も考えられます。

分布更新前後のKLダイバージェンスの実測

「エリート分位点による分布更新の導出」で見た通り、CEMの各反復はクロスエントロピー（KLダイバージェンス）を最小化する操作です。では実際に、更新前後の分布間のKLダイバージェンスはイテレーションを通してどう推移するのでしょうか。対角共分散の多変量正規分布間のKLダイバージェンスは、次の閉じた式で計算できます。

（\(d\) は次元数。対角行列なので \(\mathrm{tr}\) と \(|\cdot|\) は要素ごとの和・積で計算できます。）

import numpy as np

def rastrigin(x):
    A = 10
    return A * len(x) + np.sum(x**2 - A * np.cos(2 * np.pi * x))

def cem_with_kl(dim, n_samples, elite_frac, n_iterations, initial_mean, initial_std, seed):
    rng = np.random.default_rng(seed)
    mean = initial_mean.copy()
    std = np.full(dim, initial_std)
    n_elite = max(1, int(n_samples * elite_frac))
    stds, kl_trace = [], []

    for it in range(n_iterations):
        samples = rng.normal(loc=mean, scale=std, size=(n_samples, dim))
        scores = np.array([rastrigin(s) for s in samples])
        elite_idx = np.argsort(scores)[:n_elite]
        elite = samples[elite_idx]

        new_mean = np.mean(elite, axis=0)
        new_std = np.maximum(np.std(elite, axis=0), 1e-12)

        var0, var1 = std**2, new_std**2
        kl = 0.5 * np.sum(var0 / var1 + (new_mean - mean) ** 2 / var1 - 1 + np.log(var1 / var0))
        kl_trace.append(kl)

        mean, std = new_mean, new_std
        stds.append(np.mean(std))
        if np.all(std < 1e-6):
            break

    return stds, kl_trace

stds, kl_trace = cem_with_kl(
    dim=2, n_samples=200, elite_frac=0.2, n_iterations=100,
    initial_mean=np.array([3.0, -3.0]), initial_std=3.0, seed=42,
)
print(f"収束反復数: {len(kl_trace)}")
print(f"反復1-5のKL: {[round(k, 3) for k in kl_trace[:5]]}")
print(f"KLの最大値: {max(kl_trace):.3f} (反復{np.argmax(kl_trace)+1})")
print(f"最終反復のKL: {kl_trace[-1]:.3f}")

2次元のRastrigin関数、n_samples=200、elite_frac=0.2（本記事の実装コードとほぼ同じ設定。elite_fracのみ変更）で実行すると、26反復で収束し、次のような推移が観測されました。

• 反復1〜5: \(1.431 \to 0.586 \to 0.202 \to 0.033 \to 0.151\) と、ノイズを伴いながらも概ね減少
• 反復13付近から再上昇し始め、反復18で \(11.114\) というピークに到達
• 反復26（収束時点）では \(4.619\)

直感的には「分布が最適解に収束するにつれてKLダイバージェンスは単調に減少していく」と予想されますが、実測結果はそれに反し、収束の終盤でKLダイバージェンスがむしろ急増するという結果になりました。原因は、この終盤で標準偏差 \(\sigma\) が絶対値としては非常に小さくなっている（反復12で0.388、反復18で0.0028、反復22で \(2.8\times10^{-5}\) ）にもかかわらず、反復間の相対的な収縮率（\(\sigma_t/\sigma_{t-1}\) ）が最後まで大きいままであることです。KLダイバージェンスの \(\log(\Sigma_1/\Sigma_0)\) の項は分散の比に依存するため、絶対的な収縮幅が小さくても比が大きければ値は増大します。

これは実装上重要な注意点です。「KLダイバージェンスが小さい ⇒ 収束が近い」という直感は、分布の収縮が加速するCEMの終盤フェーズでは成立しません。 収束判定にはKLダイバージェンスではなく、本記事の実装のように標準偏差の絶対値としきい値（std < 1e-6）を使う方が安定しています。

近年の研究動向（2023年以降）

CEMは長らく組合せ最適化や希少事象シミュレーションの手法として知られてきましたが、近年はモデルベース強化学習やモデル予測制御（MPC）の文脈で再評価される場面が増えています。

• モデルベース強化学習における行動計画: 学習した環境ダイナミクスモデル上で複数の行動系列をロールアウトし、CEMで評価・エリート選択・再サンプリングを繰り返すことで、実行する行動系列を決定するアプローチが、モデルベースRLの行動計画手法の定番の一つとして使われ続けています。
• MPCの初期化・軌道最適化: 本記事で比較した MPPI を含むサンプリングベースMPC手法では、CEMのハードなエリート選択とMPPIのソフトな重み付けを組み合わせたハイブリッドな更新則（例えば上位分位点内でのみソフト重み付けを行う）が検討されており、ロボティクスの実時間軌道最適化への応用が進んでいます。
• オフライン強化学習・安全な方策探索との接続: オフラインで収集済みのデータからのみ方策を学習する設定では、学習済みの価値関数やダイナミクスモデルを目的関数としてCEM的なサンプリング最適化で行動を選ぶことで、分布外（out-of-distribution）行動への逸脱を抑えつつ性能を改善するアプローチが検討されています。

これらの研究の多くは、本記事で導出した「エリート分位点によるクロスエントロピー最小化」という基本原理をそのまま踏襲しつつ、ニューラルネットワークの価値関数やダイナミクスモデルと組み合わせる点が特徴です。ただし、この分野は現在も活発に発展中であり、具体的な手法名や性能数値については一次文献を直接参照することをお勧めします。

関連記事

• ベイズ最適化の基礎とPython実装 - CEMと同じくブラックボックス最適化の手法ですが、ガウス過程による代理モデルを用いるアプローチを解説しています。
• MPPI（Model Predictive Path Integral）の数理 - CEMの重点サンプリングの枠組みを制御問題に拡張したMPPIを、CEMとの統一的な視点から解説しています。
• 粒子フィルタのPython実装：リサンプリング手法の比較 - CEMと同じくモンテカルロサンプリングに基づく手法を、状態推定の文脈で解説しています。
• 焼きなまし法（Simulated Annealing）の理論とPython実装 - CEMと同じく確率的探索に基づく最適化手法ですが、温度パラメータによる探索・活用のバランスが特徴です。
• 遺伝的アルゴリズム（GA）の理論とPython実装 - CEMのエリートサンプル選択と類似する選択・交叉・突然変異による進化的最適化手法を解説しています。
• SGD・Adamの理論と比較 - CEMが勾配不要な最適化であるのに対し、勾配情報を活用する最適化手法を解説しています。
• 粒子群最適化（PSO）の理論とPython実装 - 群知能ベースの最適化手法。CEMの確率分布更新に対し、PSOは速度・位置更新による集団探索を行います。
• ガウス過程回帰の理論とPython実装 - ベイズ最適化のサロゲートとして使われるGPを、CEMと対比してサンプル効率の観点から理解できます。
• MCMC入門：メトロポリス・ヘイスティングス法とギブスサンプリング - CEM の重点サンプリングと対をなす MCMC サンプリングを解説。同じモンテカルロ枠組みでの探索戦略の対比を学べます。
• モンテカルロ最適化メソッド比較（CEM/SA/GA/MPPI/PSO） - CEM を含む5つのモンテカルロ最適化手法を「サンプリング→評価→更新」の共通骨格で比較するハブ記事。手法選択の指針が整理されています。

参考文献

• Rubinstein, R. Y., & Kroese, D. P. (2013). The Cross-Entropy Method: A Unified Approach to Combinatorial Optimization, Monte-Carlo Simulation, and Machine Learning. Springer Science & Business Media.
• 漆原 勉, “モンテカルロシミュレーションにおける重点サンプリング法に対する大偏差理論の適用について”
  • https://sci-tech.ksc.kwansei.ac.jp/~chiyonobu/gseminar/urushihara.pdf

関連ツール

• CalcBox - 暮らしの計算ツール - 統計計算、確率計算など61種類以上の計算ツール
• DevToolBox - 開発者向け無料ツール集 - JSON整形、正規表現テスターなど85種類以上の開発者向けツール