about optimize_acqf()

[LiuYang-1997]

Hello.
When I use optimize_ When acqf function optimizes other acquisition functions I have implemented, this error occurs. How can I avoid it?

Traceback (most recent call last):
File "D:/personalfiles/python/botorch_boa/run_acq/run_wb2s.py", line 48, in
wb2s_new_candidate = wb2s_get_next_point(wb2s_init_x, wb2s_init_y, bounds,wb2s_best_y)
File "D:/personalfiles/python/botorch_boa/run_acq/run_wb2s.py", line 28, in wb2s_get_next_point
candidates, acq_value = optimize_acqf(
File "D:\software\anaconda\lib\site-packages\botorch\optim\optimize.py", line 225, in optimize_acqf
batch_candidates_curr, batch_acq_values_curr = gen_candidates_scipy(
File "D:\software\anaconda\lib\site-packages\botorch\generation\gen.py", line 202, in gen_candidates_scipy
res = minimize(
File "D:\software\anaconda\lib\site-packages\scipy\optimize_minimize.py", line 619, in minimize
return _minimize_lbfgsb(fun, x0, args, jac, bounds,
File "D:\software\anaconda\lib\site-packages\scipy\optimize\lbfgsb.py", line 360, in _minimize_lbfgsb
f, g = func_and_grad(x)
File "D:\software\anaconda\lib\site-packages\scipy\optimize_differentiable_functions.py", line 260, in fun_and_grad
self._update_fun()
File "D:\software\anaconda\lib\site-packages\scipy\optimize_differentiable_functions.py", line 226, in _update_fun
self._update_fun_impl()
File "D:\software\anaconda\lib\site-packages\scipy\optimize_differentiable_functions.py", line 133, in update_fun
self.f = fun_wrapped(self.x)
File "D:\software\anaconda\lib\site-packages\scipy\optimize_differentiable_functions.py", line 130, in fun_wrapped
return fun(x, *args)
File "D:\software\anaconda\lib\site-packages\scipy\optimize\optimize.py", line 74, in call
self._compute_if_needed(x, *args)
File "D:\software\anaconda\lib\site-packages\scipy\optimize\optimize.py", line 68, in _compute_if_needed
fg = self.fun(x, *args)
File "D:\software\anaconda\lib\site-packages\botorch\generation\gen.py", line 181, in f_np_wrapper
raise RuntimeError(msg)
RuntimeError: 15 elements of the 15 element gradient array gradf are NaN. This often indicates numerical issues. Consider using dtype=torch.double.

def wb2s_get_next_point(train_x,train_y,bounds,best_y,beta=100.):
    train_x = normalize(train_x, bounds)
    single_model = SingleTaskGP(train_x, train_y, outcome_transform=Standardize(m=1))
    mll = ExactMarginalLogLikelihood(single_model.likelihood, single_model)
    fit_gpytorch_model(mll)
    wb2s = WB2S(
        model=single_model,
        best_f=best_y,
        beta=beta
    )
    standard_bounds = torch.zeros(2, problem.dim)
    standard_bounds[1] = 1
    candidates, acq_value = optimize_acqf(
        acq_function=wb2s,
        bounds=standard_bounds,
        q=1,
        num_restarts=15,
        raw_samples=256
    )
    new_x = unnormalize(candidates.detach(), bounds).double()
    return new_x

wb2s:

from botorch.acquisition.objective import PosteriorTransform
from botorch.acquisition.analytic import AnalyticAcquisitionFunction
from torch import Tensor
from typing import Dict, Optional, Tuple, Union
from botorch.models.model import Model
import torch
from torch.distributions import Normal

class WB2S(AnalyticAcquisitionFunction):
    def __init__(
            self,
            model: Model,
            best_f: Union[float,Tensor],
            beta: float = 100.,
            posterior_transform: Optional[PosteriorTransform] = None,
            maximize: bool = True,
            **kwargs,
    ):
        super(WB2S, self).__init__(model=model,posterior_transform=posterior_transform,**kwargs)
        self.model = model
        if not torch.is_tensor(best_f):
            best_f = torch.tensor(best_f)
        self.register_buffer("best_f", best_f)
        self.maximize = maximize
        self.beta = beta

    def forward(self, X: Tensor) -> Tensor:
        posterior = self.model.posterior(
            X=X, posterior_transform=self.posterior_transform
        )
        mean, sigma = posterior.mean, posterior.variance.sqrt().clamp_min(1e-9)
        # print(sigma)
        view_shape = mean.shape[:-2] if mean.shape[-2] == 1 else mean.shape[:-1]
        mean = mean.view(view_shape)
        sigma = sigma.view(view_shape)
        u = (mean - self.best_f.expand_as(mean)) / sigma
        if not self.maximize:
            u = -u
        normal = Normal(torch.zeros_like(u), torch.ones_like(u))
        ucdf = normal.cdf(u)
        updf = torch.exp(normal.log_prob(u))
        ei = sigma * (u * ucdf + updf)
        s = self.beta * torch.abs(mean[torch.argmax(ei)])/torch.max(ei)
        wb2s = s*ei - mean
        return wb2s.double()

[saitcakmak]

Hi @LiuYang-1997 . This seems to work for me. Maybe it is just a numerical issue due to the specific problem you're trying to optimize? I'd try using dtype=torch.double for all the arguments you pass into SingleTaskGP, the acquisition function & optimize_acqf. The improved numerical precision typically helps with these issues.

Here's what I ran to replicate the issue:

from typing import Dict, Optional, Tuple, Union

import torch
from botorch.acquisition.analytic import AnalyticAcquisitionFunction
from botorch.acquisition.objective import PosteriorTransform
from botorch.fit import fit_gpytorch_model
from botorch.models import SingleTaskGP
from botorch.models.model import Model
from botorch.models.transforms.outcome import Standardize
from botorch.optim.optimize import optimize_acqf
from botorch.test_functions.synthetic import Branin
from botorch.utils.transforms import normalize, unnormalize
from gpytorch.mlls.exact_marginal_log_likelihood import ExactMarginalLogLikelihood
from torch import Tensor
from torch.distributions import Normal


problem = Branin(negate=True)


class WB2S(AnalyticAcquisitionFunction):
    def __init__(
        self,
        model: Model,
        best_f: Union[float, Tensor],
        beta: float = 100.0,
        posterior_transform: Optional[PosteriorTransform] = None,
        maximize: bool = True,
        **kwargs,
    ):
        super(WB2S, self).__init__(model=model, posterior_transform=posterior_transform, **kwargs)
        self.model = model
        if not torch.is_tensor(best_f):
            best_f = torch.tensor(best_f)
        self.register_buffer("best_f", best_f)
        self.maximize = maximize
        self.beta = beta

    def forward(self, X: Tensor) -> Tensor:
        posterior = self.model.posterior(X=X, posterior_transform=self.posterior_transform)
        mean, sigma = posterior.mean, posterior.variance.sqrt().clamp_min(1e-9)
        # print(sigma)
        view_shape = mean.shape[:-2] if mean.shape[-2] == 1 else mean.shape[:-1]
        mean = mean.view(view_shape)
        sigma = sigma.view(view_shape)
        u = (mean - self.best_f.expand_as(mean)) / sigma
        if not self.maximize:
            u = -u
        normal = Normal(torch.zeros_like(u), torch.ones_like(u))
        ucdf = normal.cdf(u)
        updf = torch.exp(normal.log_prob(u))
        ei = sigma * (u * ucdf + updf)
        s = self.beta * torch.abs(mean[torch.argmax(ei)]) / torch.max(ei)
        wb2s = s * ei - mean
        return wb2s.double()


def wb2s_get_next_point(train_x, train_y, bounds, best_y, beta=100.0):
    train_x = normalize(train_x, bounds)
    single_model = SingleTaskGP(train_x, train_y, outcome_transform=Standardize(m=1))
    mll = ExactMarginalLogLikelihood(single_model.likelihood, single_model)
    fit_gpytorch_model(mll)
    wb2s = WB2S(model=single_model, best_f=best_y, beta=beta)
    standard_bounds = torch.zeros(2, problem.dim)
    standard_bounds[1] = 1
    candidates, acq_value = optimize_acqf(
        acq_function=wb2s, bounds=standard_bounds, q=1, num_restarts=15, raw_samples=256
    )
    new_x = unnormalize(candidates.detach(), bounds).double()
    return new_x


train_x = unnormalize(torch.rand(10, 2), problem.bounds)
train_y = problem(train_x).view(10, 1)

wb2s_get_next_point(train_x, train_y, problem.bounds, train_y.max())

[LiuYang-1997]

Sorry, I forgot to post my optimization problem, my optimization problem is a one-dimensional, single-objective optimization problem.

class TwoSine(SyntheticTestFunction):
    dim = 1
    _bounds = [(2.7,7.5)]
    _optimal_value = -1.899599
    _optimizers = [(5.145735,)]
    def evaluate_true(self, X: Tensor) -> Tensor:
        return torch.sin(X[:,0]) + torch.sin(10/3 * X[:,0])

When I use to optimize this objective function, I get the above error.

[saitcakmak]

So, this is what your function looks like. It is rather simple for BO.

You're trying to maximize this function, and its maximizer is at the boundary of 2.7. If I run a BO loop, after a while I get the error you posted. If I look at the points evaluated up to that point, I see this:

train_x
tensor([[4.7696],
        [6.6976],
        [7.1065],
        [3.4896],
        [4.4484],
        [4.4442],
        [5.5154],
        [4.7971],
        [6.6134],
        [3.0765],
        [7.5000],
        [7.5000],
        [6.1880],
        [6.2702],
        [2.7000],
        [2.7000],
        [6.2186],
        [2.7439],
        [2.7000],
        [2.7000],
        [2.7000],
        [2.7000],
        [7.5000],
        [2.7000],
        [2.7000],
        [2.7000],
        [7.5000],
        [2.7000],
        [2.7000],
        [2.7000],
        [5.1894],
        [7.5000],
        [2.7000],
        [6.2471],
        [6.2252],
        [2.7000],
        [2.7000],
        [7.5000],
        [2.7000],
        [6.2137],
        [2.7000]], dtype=torch.float64)

You may notice that there are a lot of 2.7s in there, so we already found the optimizer. The acquisition function doesn't mind evaluating the same point again, so we have multiple evaluations of 2.7. What ends up happening is that when we try to calculate the GP posterior, we need to invert / compute the Cholesky decomposition of a matrix GP_prior(train_x, train_x), which, due to the repeated occurrence of 2.7 ends up being a singular matrix (numerically singular in this case since there's a tiny bit of inferred observation noise). A singular matrix does not have an inverse or a Cholesky decomposition, so you get a bunch NaN's instead. This is the error you're seeing.

Things that may help: Optimizing a more difficult problem (so that you can't find the optimum so easily), or adding a logic into the acquisition function (or the BO loop) that discourages repeatedly evaluating the same point.

Alternatively, you could force your way around it by preventing the inferred noise level from getting too small (by passing in a custom likelihood to SingleTaskGP) or using a FixedNoiseGP with large enough train_Yvar, but these would just degrade the GP predictions by assuming large observation noise rather than solving a real problem.

[LiuYang-1997]

Thanks, I will try it. Thank you very much.