Extend RobustGPSampler to handle environmental disturbance

[msakai]

Contributor Agreements

Please read the contributor agreements and if you agree, please click the checkbox below.

• I agree to the contributor agreements.

Tip

Please follow the Quick TODO list to smoothly merge your PR.

Motivation

RobustGPSampler already handles cases where noise is added to the decision variables (for example, product dimensions typically do not match the specified values exactly, but vary within the machining accuracy of the manufacturing machines).

However, such variation are not always incidental to the decision variables. For example, the behavior of a chemical engineering process may depend on the outside temperature, which cannot be freely determined.

This PR extends RobustGPSampler to handle such cases, by introducing the notion of constant noisy parameter.

For example, if the standard outdoor temperature is C but you want to find a robust solution for variations within the range of ±ε, the user can write it as follows:

def objective(trial: optuna.Trial) -> float:
   x = trial.suggest_float("x", ...)
   outside_temperature = trial.suggest_float("outside_temperature", C-ε, C+ε)
   return f(x, outside_temperature)

sampler = RobustGPSampler(..., const_noisy_param_names=["outside_temperature"])
study = optuna.create_study(direction="minimize", sampler=sampler)
study.optimize(objective, ...)

robust_params = sampler.get_robust_params(study)

Here,

• trial.suggest_float returns a randomly sampled value within the range [C-ε, C+ε] instead of a value that aims to optimize the objective function.
• sampler.get_robust_params returns the optimal robust parameter under outside_temperature=C. (Fluctuations in outside_temperature are accounted for by the acquisition functions which consider nearby points on the GP models).

@nabenabe0928 designed this API and made the initial implementation, and I took over the rest of the implementation.

This API design lacks consistency in the following aspects between decision variables (e.g. dimentions of products) and non-decision variables (e.g. outside temperature), but we plan to address this issue separately in the future.

• For decision variables, we need to specify a nominal range in suggest_float, and it returns nominal values unaffected by disturbances.
• For non-decision variables, we have to specify a range accounting for variation, and it returns values after applying disturbances.

Description of the changes

TODO List towards PR Merge

Please remove this section if this PR is not an addition of a new package.
Otherwise, please check the following TODO list:

• Copy ./template/ to create your package
• Replace <COPYRIGHT HOLDER> in LICENSE of your package with your name
• Fill out README.md in your package
• Add import statements of your function or class names to be used in __init__.py
• (Optional) Add from __future__ import annotations at the head of any Python files that include typing to support older Python versions
• Apply the formatter based on the tips in README.md
• Check whether your module works as intended based on the tips in README.md

[msakai]

Here is the code I used to test its functionality.

This is based on an example from the paper Constrained robust Bayesian optimization of expensive noisy black-box functions with guaranteed regret bounds, though this is merely an artificial example and does not correspond to the intended problem setting.

import matplotlib.pyplot as plt
import numpy as np
import optuna
import optunahub


def g_c1(x, y):
    c1 = (x - 1.5) ** 4 + (y - 1.5) ** 4 - 10.125
    return c1


def g_c2(x, y):
    c2 = -((2.5 - x) ** 3) - (y + 1.5) ** 3 + 15.75
    return c2


def f(x, y):
    t1 = (
        2 * (x**6)
        - 12.2 * (x**5)
        + 21.2 * (x**4)
        - 6.4 * (x**3)
        - 4.7 * (x**2)
        + 6.2 * x
    )
    t2 = (y**6) - 11 * (y**5) + 43.3 * (y**4) - 74.8 * (y**3) + 56.9 * (y**2) - 10 * y
    t3 = -4.1 * x * y - 0.1 * (x**2) * (y**2) + 0.4 * (y**2) * x + 0.4 * (x**2) * y
    return t1 + t2 + t3


def objective(trial: optuna.Trial) -> float:
    x = trial.suggest_float("x", -1, 4)
    y = trial.suggest_float(
        "y", 0.5, 1.5
    )  # The nominal value is 1, but noise of ±0.5 is possible.

    c1 = float(g_c1(x, y))
    c2 = float(g_c2(x, y))
    trial.set_user_attr("c", (c1, c2))

    return float(f(x, y))


def constraints_func(trial: optuna.trial.FrozenTrial) -> tuple[float, float]:
    return trial.user_attrs["c"]


def plot_func():
    xmin, xmax, ymin, ymax = (-1, 4, -1, 4)
    n_bins = 100

    xx = np.linspace(xmin, xmax, n_bins)
    yy = np.linspace(ymin, ymax, n_bins)
    grids = np.meshgrid(xx, yy)

    c1 = g_c1(grids[0].flatten(), grids[1].flatten())
    mask = c1 < 0
    c1[mask] = None
    plt.contour(
        xx, yy, c1.reshape(n_bins, n_bins), colors="red", levels=20, linewidths=0.5
    )

    c2 = g_c2(grids[0].flatten(), grids[1].flatten())
    mask = c2 < 0
    c2[mask] = None
    plt.contour(
        xx, yy, c2.reshape(n_bins, n_bins), colors="blue", levels=20, linewidths=0.5
    )

    fs = f(grids[0].flatten(), grids[1].flatten())
    plt.contour(
        xx, yy, fs.reshape(n_bins, n_bins), colors="black", levels=200, linewidths=0.5
    )


if __name__ == "__main__":
    RobustGPSampler = optunahub.load_local_module(
        package="samplers/value_at_risk", registry_root="optunahub-registry/package/"
    ).RobustGPSampler
    sampler = RobustGPSampler(
        seed=3,
        n_startup_trials=10,
        constraints_func=constraints_func,
        uniform_input_noise_rads={
            "x": 0.5,
        },
        const_noisy_param_names=["y"],
    )
    study = optuna.create_study(direction="minimize", sampler=sampler)
    study.optimize(objective, n_trials=100)

    plot_func()

    plt.scatter(
        [trial.params["x"] for trial in study.trials],
        [trial.params["y"] for trial in study.trials],
    )

    robust_params = sampler.get_robust_params(study)
    print(robust_params)
    plt.scatter([robust_params["x"]], [robust_params["y"]], marker="*", color="red")

    plt.savefig("f.png")

[c-bata]

@HideakiImamura Could you review this PR?

[HideakiImamura]

Thanks for the PR. Basically, LGTM. I have a comment about the initial value of the acquisition function optimization. PTAL.

[HideakiImamura]

This modification likely serves to prevent getting stuck in local optima. The change essentially involves setting the initial values for iterative optimization to random values without altering the acquisition function itself - the target of optimization. However, based on my personal verification, even with simple objective functions, the current GPSampler's acquisition function becomes extremely jagged and becomes ill-defined for gradient-based optimization. Therefore, rather than randomly shifting initial values, wouldn't it be better to smooth the acquisition function itself to make it more amenable to optimization?

Specifically, consider setting ConstrainedLogValueAtRisk's stabilizing_noise to a larger value than currently used (e.g., 1e-5). Note that this value requires adjustment based on the specific application.

[msakai]

This change was made by @nabenabe0928 -san in 6dae507 .

@nabenabe0928 -san removed best_params from CARBO for that reason (#301), but in this PR, it might be simply because None was passed as best_params at every call site of _optimize_acqf.

How about merging this PR as is, and re-add best_params in different PR if necessary?

Anyway, the fact that samples are overly concentrated in specific areas is a big issue, so I'd like to try modifying the stabilizing_noise parameter.

[HideakiImamura]

How about merging this PR as is, and re-add best_params in different PR if necessary?

Sounds good! Sorry for the confusion. I misunderstood the motivation of setting best_params=None.

[HideakiImamura]

LGTM.