基于 Wilcoxon 剪枝器的早期停止独立评估

本教程展示了 Optuna 的 WilcoxonPruner。此剪枝器对于平均多个评估的优化函数非常有效。

我们使用 模拟退火 (SA) 来解决 旅行商问题 (TSP)。

概述：使用模拟退火解决旅行商问题

旅行商问题 (TSP) 是组合优化中的一个经典问题，涉及为需要一次性访问一组城市并返回起点的推销员找到最短的可行路线。TSP 在数学、计算机科学和运筹学等领域得到了广泛研究，并在物流、制造和 DNA 测序等领域有许多实际应用。该问题被归类为 NP 难问题，因此对于较大的实例通常采用近似算法或启发式方法。

适用于 TSP 的一种简单启发式方法是模拟退火 (SA)。SA 从一个初始解开始（可以由更简单的启发式方法（如贪心法）构建），然后随机检查解的邻域（稍后定义）。如果邻域更好，则将解更新为邻域。如果邻域更差，SA 仍以概率 \(e^{-\Delta c / T}\) 将解更新为邻域，其中 \(\Delta c (> 0)\) 是新解和旧解之间成本（距离之和）的差异，而 \(T\) 是称为“温度”的参数。温度控制着为了逃离局部最小值而容忍解的恶化程度（温度高表示容忍度大）。如果温度太低，SA 将很快陷入局部最小值；如果温度太高，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)

注意

为简化实现，我们使用带 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)
        # Reverse the order of vertices between range [i+1, j].

        # cost difference by 2-opt reversal
        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):
            # accept the 2-opt reversal
            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:
            # If the best solution is not updated for "patience" iteratoins,
            # restart from the best solution.
            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。

T0 和 alpha 定义温度计划

在模拟退火中，确定一个好的温度计划很重要，但没有一个“万能计划”适合所有问题，因此我们必须针对这个问题调整计划。此代码将温度参数化为一个单项函数 T0 * (1 - t) ** alpha，其中 t 从 0 增加到 1。我们尝试优化两个参数 T0 和 alpha。

patience

此参数指定了一个阈值，允许退火过程在不更新最佳值的情况下继续多少次迭代。实际上，模拟退火经常将解推离当前最佳解很远，定期回滚到最佳解通常可以大大提高优化效率。但是，如果回滚发生得太频繁，优化可能会陷入局部最优，因此我们必须将阈值调整到一个合理的范围。

注意

一些采样器，包括默认的 TPESampler，目前无法有效利用剪枝试验的信息（尤其是在最后一个中间值不是最终目标函数的最佳近似时）。作为此问题的解决方法，您可以在 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 = []

    # For best results, shuffle the evaluation order in each trial.
    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():
            # Return the current predicted value instead of raising `TrialPruned`.
            # This is a workaround to tell the Optuna about the evaluation
            # results in pruned trials.
            return sum(results) / len(results)

    return sum(results) / len(results)

我们使用带 WilcoxonPruner 的 TPESampler。

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})  # default params
study.optimize(objective, n_trials=N_TRIALS)

/home/docs/checkouts/readthedocs.org/user_builds/optuna/checkouts/stable/tutorial/20_recipes/013_wilcoxon_pruner.py:255: ExperimentalWarning:

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

我们可以这样显示优化结果：

print(f"The number of trials: {len(study.trials)}")
print(f"Best value: {study.best_value} (params: {study.best_params})")
print(f"Number of evaluation: {num_evaluation} / {len(dataset) * N_TRIALS}")

The number of trials: 50
Best value: 8.362043398560909 (params: {'T0': 0.01742184332150064, 'alpha': 5.35111546857422, 'patience': 69})
Number of evaluation: 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="Evaluations"))
fig.add_trace(go.Bar(x=best_trial_x, y=best_trial_y, name="Best Trial", marker_color="red"))
fig.update_layout(
    title="Number of evaluations in each trial",
    xaxis_title="Trial number",
    yaxis_title="Number of evaluations before pruned",
)
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"greedy solution (initial guess),  cost: {tsp_cost(d, result_idxs):.3f}",
    xaxis=dict(scaleanchor="y", scaleratio=1),
)
fig

可视化 tsp_simulated_annealing 为同一个 TSP 问题找到的解。

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}, cost: {tsp_cost(d, result_idxs):.3f}",
    xaxis=dict(scaleanchor="y", scaleratio=1),
)
fig