モンテカルロ最適化（CEM/SA/GA/MPPI/PSO）の基礎・比較・Python実装

モンテカルロ最適化（CEM/SA/GA/MPPI/PSO）の基礎・比較・Python実装をハブ的に整理。サンプリング→評価→更新の共通骨格、各手法の使い分け指針、numpyだけで動く共通フレームワークを解説します。

4 min read

はじめに

モンテカルロ最適化（Monte Carlo optimization） とは、勾配 \(\nabla f\) に頼らず 無作為サンプル を使って目的関数 \(f(x)\) の最小値（あるいは最大値）を探索する手法群の総称です。微分不可能な目的関数、ブラックボックスシミュレータ、多峰性のある探索空間など、勾配法が破綻する状況でこそ威力を発揮します。

クロスエントロピー法（CEM）、焼きなまし法（SA）、遺伝的アルゴリズム（GA）、MPPI、粒子群最適化（PSO） は、いずれもこの「サンプリング→評価→更新」を骨格とするモンテカルロ最適化の代表例です。本記事はそれぞれの手法を 横断的に位置づけ、共通フレームワーク・使い分け指針・Python 実装をハブ的に整理することを目的とします。

モンテカルロ最適化の共通骨格

すべての手法は、以下の 3 ステップを反復しているとみなせます。

  1. サンプリング: 現在の確率分布 \(p_t(x)\) または現解 \(x_t\) の近傍から \(N\) 個のサンプル \(\{x^{(i)}\}_{i=1}^{N}\) を生成
  2. 評価: 各サンプルに対し目的関数値 \(f(x^{(i)})\) を計算
  3. 更新: 評価結果から「次の分布 \(p_{t+1}\) あるいは次解 \(x_{t+1}\) 」を構成

数式で書けば、共通の更新は次のように表せます。

\[ p_{t+1}(x) = \mathcal{U}\big(p_t,\, \{(x^{(i)}, f(x^{(i)}))\}_{i=1}^{N}\big) \tag{1} \]

ここで \(\mathcal{U}\) が 手法ごとに異なる更新オペレータ です。CEM ならエリート分位点の経験分布へのフィット、SA なら Boltzmann 比による確率的受理、GA なら選択・交叉・突然変異、MPPI ならコスト指数重み付け平均、PSO なら速度ベクトルの慣性 + 認知 + 社会更新——これらはすべて同じ枠組みの 異なる \(\mathcal{U}\) にすぎません。

ベイズ最適化 もサンプリング型最適化の一種ですが、サロゲートモデル（ガウス過程）と獲得関数を介する点で「分布更新」よりも「サンプリング点選択」を重視する位置づけになります。また パーティクルフィルタ はモンテカルロサンプリングを 状態推定 に用いる兄弟手法であり、リサンプリング機構が CEM のエリート選択や MPPI の重み付けと数学的に等価です。

各手法の位置づけ表

手法探索戦略解の表現主な用途計算コスト微分
CEM分布更新（エリート分位点）パラメトリック分布連続最適化、強化学習方策最適化中（並列化容易）不要
SA確率的受理（温度スケジュール）単一解組合せ最適化、TSP、ハイパーパラメータ低（逐次）不要
GA個体集団の進化（交叉・突然変異）集団大域最適化、構造設計高（評価回数大）不要
MPPI重点サンプリング（指数重み）軌道列の分布モデル予測制御、ロボティクス中〜高（並列化必須）不要
PSO慣性 + 認知 + 社会更新粒子群連続最適化、群知能タスク中不要
ベイズ最適化サロゲート + 獲得関数GP 後分布評価コストが高い問題高（サロゲート更新）不要（GP 内部で勾配）

「微分不要」「並列評価可能」「多峰性に強い」という 3 つの性質はすべての手法に共通します。違いは どの確率モデルを使い、どんな更新で集中させるか にあります。

共通フレームワークの Python 実装

実際にどの手法も次のように、同じスケルトンに sample と update だけを差し替える だけで書けます。

import numpy as np

def monte_carlo_optimize(
    f, sampler, updater, init_state, n_iter=50, n_samples=100, seed=0
):
    """
    汎用モンテカルロ最適化ループ。

    Parameters
    ----------
    f          : callable  目的関数 f(x) を最小化
    sampler    : callable  state -> shape (n_samples, dim) のサンプル配列
    updater    : callable  (state, samples, scores) -> 新しい state
    init_state : 初期状態（手法ごとに分布パラメータ・個体集団など）
    n_iter     : 反復回数
    n_samples  : 1反復あたりのサンプル数
    """
    rng = np.random.default_rng(seed)
    state = init_state
    history = []
    for t in range(n_iter):
        samples = sampler(state, n_samples, rng)
        scores = np.array([f(x) for x in samples])
        state = updater(state, samples, scores)
        best = scores.min()
        history.append(best)
    return state, np.array(history)


# === 例1: CEMライクな更新（エリート分位点で正規分布をフィット） ===
def cem_sampler(state, n, rng):
    mu, sigma = state
    return rng.normal(mu, sigma, size=(n, mu.shape[0]))


def cem_updater(state, samples, scores, elite_frac=0.2):
    mu, sigma = state
    k = max(1, int(len(scores) * elite_frac))
    elite_idx = np.argsort(scores)[:k]
    elite = samples[elite_idx]
    return elite.mean(axis=0), elite.std(axis=0) + 1e-6


# === 例2: PSOライクな更新（速度更新） ===
def make_pso(dim, n_particles, w=0.7, c1=1.5, c2=1.5):
    def sampler(state, n, rng):
        x, v, pbest, pbest_val, gbest = state
        return x  # PSO は state そのものを評価する

    def updater(state, samples, scores):
        x, v, pbest, pbest_val, gbest = state
        # 個体ベスト更新
        improved = scores < pbest_val
        pbest = np.where(improved[:, None], x, pbest)
        pbest_val = np.where(improved, scores, pbest_val)
        # 大域ベスト更新
        g_idx = pbest_val.argmin()
        gbest = pbest[g_idx]
        # 速度・位置更新
        rng = np.random.default_rng()
        r1, r2 = rng.random(x.shape), rng.random(x.shape)
        v = w * v + c1 * r1 * (pbest - x) + c2 * r2 * (gbest - x)
        x = x + v
        return x, v, pbest, pbest_val, gbest

    return sampler, updater


# === 目的関数: Rastrigin (多峰性ベンチマーク) ===
def rastrigin(x, A=10.0):
    n = x.shape[0]
    return A * n + np.sum(x ** 2 - A * np.cos(2 * np.pi * x))


# === CEM で最適化 ===
dim = 5
init_state = (np.zeros(dim), np.ones(dim) * 2.0)
state, hist = monte_carlo_optimize(
    rastrigin, cem_sampler, cem_updater, init_state,
    n_iter=40, n_samples=200
)
print(f"CEM best: {hist.min():.4f}")

上記のスケルトンは CEM・SA・GA・MPPI・PSO すべてに流用 できます。SA なら sampler が「現解の近傍ランダムウォーク」、updater が「Metropolis 受理判定」になります。GA なら sampler が「現世代そのもの」、updater が「選択 + 交叉 + 突然変異」になります。MPPI なら sampler が「制御列ノイズ付き軌道生成」、updater が「コスト指数重み付き平均で名目制御列を更新」になります。

この一般化は理論的にも重要で、CEM と MPPI が 重点サンプリングの最適提案分布 という同じ起源を共有することは MPPI 解説記事 で詳しく述べています。

各手法の深堀り誘導

CEM（クロスエントロピー法）

KL ダイバージェンス最小化に基づき、エリートサンプルから次の分布を最尤推定します。実装が非常にシンプルで、強化学習の方策パラメータ探索でも使われます。

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

SA（焼きなまし法）

温度 \(T\) を徐々に下げながら、悪化方向への遷移を Boltzmann 確率 \(\exp(-\Delta f / T)\) で確率的に受理します。組合せ問題や TSP で長年標準的に使われてきた手法です。

詳細: 焼きなまし法（Simulated Annealing）の仕組みとPython実装

GA（遺伝的アルゴリズム）

個体集団を選択・交叉・突然変異で進化させます。離散・連続のいずれにも対応し、構造設計やパス計画など多様な分野で使われます。

詳細: 遺伝的アルゴリズム（GA）の基礎とPython実装

MPPI（Model Predictive Path Integral）

制御入力列を確率分布として表し、コストの指数重み付け平均で名目軌道を更新します。CEM の連続時間・連続制御版とみなせ、ロボティクスのリアルタイム MPC で標準化が進んでいます。

詳細: MPPI（Model Predictive Path Integral）の数理

PSO（粒子群最適化）

粒子群の速度を「慣性 + 個体ベスト方向 + 大域ベスト方向」で更新します。実装が直感的で、ハイパーパラメータが少ないため広く使われます。

詳細: 粒子群最適化（PSO）の数理とPython実装

ベイズ最適化

評価コストが高い目的関数に対し、ガウス過程でサロゲートを構築し、獲得関数で次の評価点を選択します。サンプル数が少ない局面で他のモンテカルロ系を凌駕します。

詳細: ベイズ最適化の基礎とPython実装

パーティクルフィルタ（兄弟手法）

最適化ではなく状態推定の手法ですが、リサンプリング機構は CEM のエリート選択や MPPI の重み付けと数学的に等価で、モンテカルロ系の 元祖 にあたります。

詳細: 粒子フィルタのPython実装

使い分けの指針

実問題でどの手法を選ぶかは、以下の観点で判断するのが実用的です。

状況推奨手法理由
目的関数評価が 超高コスト（数分以上）ベイズ最適化サンプル効率が桁違いに良い
連続最適化、低〜中次元（〜数十次元）CEM, PSO実装シンプル、収束が速い
連続最適化、高次元（数百次元以上）CEM, MPPI分布パラメトリック化で次元の呪いを緩和
組合せ最適化、TSP、スケジューリングSA, GA離散近傍/個体表現が自然
リアルタイム制御、ロボティクスMPPIGPU で全サンプル並列、毎時刻 100Hz 以上で更新可能
多峰性が強い、大域解が必須GA, CEM（多モード化）集団による多様性維持
微分情報が一部使える勾配法とのハイブリッドモンテカルロで初期化 → 勾配法で精緻化

簡易フローチャート

目的関数の評価コストは？
├─ 超高コスト ──> ベイズ最適化
└─ 中〜低コスト
    │
    解は連続 or 離散？
    ├─ 離散 ──> SA / GA
    └─ 連続
        │
        リアルタイム性が必要？
        ├─ Yes ──> MPPI
        └─ No
            │
            個体ベストの履歴情報を活かしたい？
            ├─ Yes ──> PSO
            └─ No  ──> CEM

このチャートはあくまで第一近似で、実問題ではハイブリッド化（例: CEM で粗く探索した後 GA で多様性を担保）が有効なケースも多くあります。

信号処理との接点

モンテカルロ最適化は信号処理にも応用が広がっています。

  • 適応フィルタの非凸目的関数: 通常の LMS/RLS が極小に陥る状況で、CEM や GA を初期化に使う
  • フィルタ係数の構造最適化: 整数係数の FIR フィルタ設計（量子化制約あり）に SA / GA
  • モデル予測制御の制御列最適化: 信号系制御では MPPI が標準化しつつある
  • 状態推定: 非線形・非ガウス系では パーティクルフィルタ が線形フィルタ（カルマン）を凌駕

特に MPPI は、伝統的な LQR / MPC では扱えない 微分不可能なコスト（衝突回避、二値制約など）を素直に扱える点で、信号処理 + 制御の融合領域で急速に普及しています。

まとめ

  • モンテカルロ最適化は サンプリング → 評価 → 更新 の共通骨格を持つ
  • CEM / SA / GA / MPPI / PSO は「更新オペレータ \(\mathcal{U}\) 」の違いとして統一的に理解できる
  • Python では sampler と updater を差し替えるだけのスケルトンで実装できる
  • 手法選択は 評価コスト → 連続/離散 → リアルタイム性 → 履歴活用度 の順で絞ると実用的
  • ベイズ最適化やパーティクルフィルタも同じモンテカルロ系の枠組みに位置づけられる

各手法の詳細な数理と実装例は、以下の関連記事で深堀りしてください。

関連記事

参考文献

  • Rubinstein, R. Y., & Kroese, D. P. (2004). The Cross-Entropy Method. Springer.
  • Kirkpatrick, S., Gelatt, C. D., & Vecchi, M. P. (1983). Optimization by Simulated Annealing. Science, 220(4598), 671–680.
  • Holland, J. H. (1992). Adaptation in Natural and Artificial Systems. MIT Press.
  • Williams, G., et al. (2017). Model Predictive Path Integral Control: From Theory to Parallel Computation. Journal of Guidance, Control, and Dynamics, 40(2).
  • Kennedy, J., & Eberhart, R. (1995). Particle Swarm Optimization. Proc. IEEE ICNN.
  • Shahriari, B., et al. (2016). Taking the Human Out of the Loop: A Review of Bayesian Optimization. Proc. IEEE, 104(1).

関連ツール

この記事に関連する無料ツール

DevToolBox Moving Average Calculator SMA/EMA/WMA をブラウザで即座に計算 DevToolBox Frequency Converter Hz/kHz/MHz/rad/s 単位を相互変換 DevToolBox JSON Formatter JSONを整形・圧縮・バリデーション

関連記事