通过 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.351115468574221, '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