Wilcoxon prunerによる早期停止型独立評価

本チュートリアルでは、Optunaの WilcoxonPruner を紹介します。 このプルーナーは、複数の評価値を平均化する目的関数に対して有効です。

巡回セールスマン問題 (TSP) を シミュレーテッドアニーリング (SA) で解きます。

概要: シミュレーテッドアニーリングによる巡回セールスマン問題の解決

巡回セールスマン問題 (TSP) は、組み合わせ最適化における古典的な問題で、 各都市を一度だけ訪れ、出発地に戻る最短経路を求める問題です。 数学、コンピュータサイエンス、オペレーションズリサーチなどの分野で 広く研究されており、物流、製造業、DNA配列解析など多岐にわたる 実用的な応用があります。 この問題はNP困難に分類されるため、大規模なインスタンスでは 近似アルゴリズムやヒューリスティック手法が一般的に用いられます。

TSPに適用可能な単純なヒューリスティック手法の一つがシミュレーテッドアニーリング (SA) です。 SAは初期解（貪欲法などのより単純なヒューリスティックで構築可能）から開始し、 解の近傍（後述）をランダムに探索します。より良い近傍が見つかった場合、 その近傍に解を更新します。より悪い近傍が見つかった場合でも、 SAは確率 \(e^{-\Delta c / T}\) で解を更新します。ここで \(\Delta c (> 0)\) は新解と旧解のコスト（距離の総和）の差、 \(T\) は「温度」と呼ばれるパラメータです。温度は局所最適解からの 脱出の程度を制御します（高いほど許容度が高い）。温度が低すぎると SAは局所最適解に早く陥り、高すぎるとランダムウォーク状態になり 最適化が非効率になります。通常、高温から徐々にゼロへ減少する 「温度スケジュール」を設定します。

TSPの近傍を定義する方法はいくつかありますが、ここでは 2-opt と呼ばれる単純な近傍を使用します。2-opt近傍は現在の解から経路を選択し、 その経路の訪問順序を反転させます。例えば、初期解が a→b→c→d→e→a の場合、 a→d→c→b→e→a は2-opt近傍です（b から d への経路が反転）。 この近傍の利点は、コスト差の計算が定数時間で可能なことです （選択経路の始点と終点のみを考慮すればよい）。

メインチュートリアル: TSPのためのSAパラメータ調整

まず、必要なパッケージをインポートし、SAのパラメータ設定と TSPのコスト関数を定義します。

from dataclasses import dataclass
import math

import numpy as np
import optuna
import plotly.graph_objects as go
from numpy.linalg import norm


@dataclass
class SAOptions:
    max_iter: int = 10000
    T0: float = 1.0
    alpha: float = 2.0
    patience: int = 50


def tsp_cost(vertices: np.ndarray, idxs: np.ndarray) -> float:
    return norm(vertices[idxs] - vertices[np.roll(idxs, 1)], axis=-1).sum()

初期解として貪欲法を使用。

def tsp_greedy(vertices: np.ndarray) -> np.ndarray:
    idxs = [0]
    for _ in range(len(vertices) - 1):
        dists_from_last = norm(vertices[idxs[-1], None] - vertices, axis=-1)
        dists_from_last[idxs] = np.inf
        idxs.append(np.argmin(dists_from_last))
    return np.array(idxs)

Note

実装の簡略化のため、2-opt近傍を用いたSAでTSPを解いていますが、 これはTSPを解く「最良の」方法ではありません。この手法よりも はるかに高度な解法が存在します。

2-opt近傍を用いたSAの実装は以下の通りです。

def tsp_simulated_annealing(vertices: np.ndarray, options: SAOptions) -> np.ndarray:

    def temperature(t: float):
        assert 0.0 <= t and t <= 1.0
        return options.T0 * (1 - t) ** options.alpha

    N = len(vertices)

    idxs = tsp_greedy(vertices)
    cost = tsp_cost(vertices, idxs)
    best_idxs = idxs.copy()
    best_cost = cost
    remaining_patience = options.patience

    for iter in range(options.max_iter):

        i = np.random.randint(0, N)
        j = (i + 2 + np.random.randint(0, N - 3)) % N
        i, j = min(i, j), max(i, j)
        # [i+1, j]の範囲の頂点の順序を反転。

        # 2-opt反転によるコスト変化
        delta_cost = (
            -norm(vertices[idxs[(i + 1) % N]] - vertices[idxs[i]])
            - norm(vertices[idxs[j]] - vertices[idxs[(j + 1) % N]])
            + norm(vertices[idxs[i]] - vertices[idxs[j]])
            + norm(vertices[idxs[(i + 1) % N]] - vertices[idxs[(j + 1) % N]])
        )
        temp = temperature(iter / options.max_iter)
        if delta_cost <= 0.0 or np.random.random() < math.exp(-delta_cost / temp):
            # 2-opt反転を受容
            cost += delta_cost
            idxs[i + 1 : j + 1] = idxs[i + 1 : j + 1][::-1]
            if cost < best_cost:
                best_idxs[:] = idxs
                best_cost = cost
                remaining_patience = options.patience

        if cost > best_cost:
            # 「忍耐」回数分最良の解が更新されない場合、
            # 最良の解から再スタート。
            remaining_patience -= 1
            if remaining_patience == 0:
                idxs[:] = best_idxs
                cost = best_cost
                remaining_patience = options.patience

    return best_idxs

TSPのランダムデータセットを作成。

def make_dataset(num_vertex: int, num_problem: int, seed: int = 0) -> np.ndarray:
    rng = np.random.default_rng(seed=seed)
    return rng.random((num_problem, num_vertex, 2))


dataset = make_dataset(
    num_vertex=100,
    num_problem=50,
)

N_TRIALS = 50

デモ用にSAの反復回数を非常に少なく設定。 実際の使用では1000000程度の大きな値を推奨。

N_SA_ITER = 10000

評価回数をカウントし、剪定された問題数を把握。

num_evaluation = 0

本チュートリアルでは、T0、alpha、 patience の3つのパラメータを最適化。

T0 と alpha は温度スケジュールを定義

シミュレーテッドアニーリングでは適切な温度スケジュールの設定が重要だが、 全ての問題に適用できる「最適スケジュール」は存在しないため、 問題ごとに調整が必要。 本コードでは温度を単項式関数 T0 * (1 - t) ** alpha で表現し、 T0 と alpha の2つのパラメータを最適化。

patience

このパラメータは、最良の値が更新されない場合の許容反復回数を指定。 シミュレーテッドアニーリングでは最適解から大きく離れることが多く、 定期的に最良解に戻ることで最適化効率が向上するが、 頻繁にロールバックすると局所最適解に留まる可能性があるため、 適切な閾値設定が必要。

Note

デフォルトの TPESampler を含む一部のサンプラーは、 剪定された試行の情報（特に最終中間値が最適解でない場合） を効果的に活用できない（optuna.TrialPruned()をraiseする場合）。 この問題の回避策として、 trial.should_prune() がTrueを返す場合、 最終値の推定値（例：全評価値の平均）を返却することで、 サンプラーの性能向上が期待できる。 （raise optuna.TrialPruned()の代わりに）

最適化対象の目的関数を以下のように定義。 プルーナーを使用して評価を早期終了。

def objective(trial: optuna.Trial) -> float:
    global num_evaluation
    options = SAOptions(
        max_iter=N_SA_ITER,
        T0=trial.suggest_float("T0", 0.01, 10.0, log=True),
        alpha=trial.suggest_float("alpha", 1.0, 10.0, log=True),
        patience=trial.suggest_int("patience", 10, 1000, log=True),
    )
    results = []

    # 最良の結果を得るため、各試行で評価順序をシャッフル。
    instance_ids = np.random.permutation(len(dataset))
    for instance_id in instance_ids:
        num_evaluation += 1
        result_idxs = tsp_simulated_annealing(vertices=dataset[instance_id], options=options)
        result_cost = tsp_cost(dataset[instance_id], result_idxs)
        results.append(result_cost)

        trial.report(result_cost, instance_id)
        if trial.should_prune():
            # `TrialPruned` 例外を発生させず、現在の予測値を返す。
            # プルーニングされた試行の評価結果を Optuna に伝えるための回避策。
            return sum(results) / len(results)

    return sum(results) / len(results)

TPESampler と WilcoxonPruner を使用。

np.random.seed(0)
sampler = optuna.samplers.TPESampler(seed=1)
pruner = optuna.pruners.WilcoxonPruner(p_threshold=0.1)
study = optuna.create_study(direction="minimize", sampler=sampler, pruner=pruner)
study.enqueue_trial({"T0": 1.0, "alpha": 2.0, "patience": 50})  # デフォルトパラメータ
study.optimize(objective, n_trials=N_TRIALS)

/mnt/nfs-mnj-hot-99-home/mshibata/sandbox/optuna-documentation-plamo-ja/optuna-doc-plamo-translation/tmp-optuna/tutorial/20_recipes/013_wilcoxon_pruner.py:248: ExperimentalWarning:

WilcoxonPruner is experimental (supported from v3.6.0). The interface can change in the future.

最適化結果を表示。

print(f"試行数: {len(study.trials)}")
print(f"最適値: {study.best_value} (パラメータ: {study.best_params})")
print(f"評価回数: {num_evaluation} / {len(dataset) * N_TRIALS}")

試行数: 50
最適値: 8.362043398560909 (パラメータ: {'T0': 0.01742184332150064, 'alpha': 5.351115468574221, 'patience': 69})
評価回数: 1023 / 2500

最適化履歴を可視化。 完了試行とプルーニング試行を同じ方法で表示。

optuna.visualization.plot_optimization_history(study)

各試行の評価回数を可視化。

x_values = [x for x in range(len(study.trials)) if x != study.best_trial.number]
y_values = [
    len(t.intermediate_values) for t in study.trials if t.number != study.best_trial.number
]
best_trial_y = [len(study.best_trial.intermediate_values)]
best_trial_x = [study.best_trial.number]
fig = go.Figure()
fig.add_trace(go.Bar(x=x_values, y=y_values, name="評価回数"))
fig.add_trace(go.Bar(x=best_trial_x, y=best_trial_y, name="最適試行", marker_color="red"))
fig.update_layout(
    title="各試行の評価回数",
    xaxis_title="試行番号",
    yaxis_title="プルーニング前の評価回数",
)
fig

TSP問題の初期推測として使用される貪欲解を可視化。

d = dataset[0]
result_idxs = tsp_greedy(d)
result_idxs = np.append(result_idxs, result_idxs[0])
fig = go.Figure()
fig.add_trace(go.Scatter(x=d[result_idxs, 0], y=d[result_idxs, 1], mode="lines+markers"))
fig.update_layout(
    title=f"貪欲解（初期推測）, コスト: {tsp_cost(d, result_idxs):.3f}",
    xaxis=dict(scaleanchor="y", scaleratio=1),
)
fig

tsp_simulated_annealing が見つけた解を可視化。

params = study.best_params
options = SAOptions(
    max_iter=N_SA_ITER,
    patience=params["patience"],
    T0=params["T0"],
    alpha=params["alpha"],
)
result_idxs = tsp_simulated_annealing(d, options)
result_idxs = np.append(result_idxs, result_idxs[0])
fig = go.Figure()
fig.add_trace(go.Scatter(x=d[result_idxs, 0], y=d[result_idxs, 1], mode="lines+markers"))
fig.update_layout(
    title=f"n_iter: {options.max_iter}, コスト: {tsp_cost(d, result_idxs):.3f}",
    xaxis=dict(scaleanchor="y", scaleratio=1),
)
fig