[Bug] ParEGO fails for >512 data points

[DavidWalz]

Hi, I noticed that ParEGO fails once the dataset has >= 512 entries.

To reproduce

import torch
from botorch.acquisition.monte_carlo import qExpectedImprovement
from botorch.acquisition.objective import GenericMCObjective
from botorch.fit import fit_gpytorch_model
from botorch.models import SingleTaskGP
from botorch.optim.optimize import optimize_acqf
from botorch.sampling.samplers import SobolQMCNormalSampler
from botorch.utils.multi_objective.scalarization import \
    get_chebyshev_scalarization
from botorch.utils.sampling import sample_simplex
from gpytorch.mlls import ExactMarginalLogLikelihood


m = 2
d = 5
n = 512

X = torch.rand((n, d), dtype=torch.double)
Y = torch.rand((n, m), dtype=torch.double)

model = SingleTaskGP(X, Y)
mll = ExactMarginalLogLikelihood(model.likelihood, model)
fit_gpytorch_model(mll)

sampler = SobolQMCNormalSampler(num_samples=1024)
weights = sample_simplex(d=m, n=1).squeeze()
objective = GenericMCObjective(get_chebyshev_scalarization(weights=weights, Y=Y))
bounds = torch.tensor([[0,] * d, [1,] * d], dtype=torch.double)

acq_func = qExpectedImprovement(
    model=model,
    objective=objective,
    best_f=objective(Y).max(),
    sampler=sampler,
)

candidates, _ = optimize_acqf(
    acq_function=acq_func,
    q=1,
    bounds=bounds,
    num_restarts=1,
    raw_samples=1024,
)

raises

IndexError                                Traceback (most recent call last)
2022-07-08-bug.py in <module>
     43 )
     44 
---> 45 candidates, _ = optimize_acqf(
     46     acq_function=acq_func,
     47     q=1,

\botorch\optim\optimize.py in optimize_acqf(acq_function, bounds, q, num_restarts, raw_samples, options, inequality_constraints, equality_constraints, nonlinear_inequality_constraints, fixed_features, post_processing_func, batch_initial_conditions, return_best_only, sequential, **kwargs)
    203             else gen_batch_initial_conditions
    204         )
--> 205         batch_initial_conditions = ic_gen(
    206             acq_function=acq_function,
    207             bounds=bounds,

\botorch\optim\initializers.py in gen_batch_initial_conditions(acq_function, bounds, q, num_restarts, raw_samples, fixed_features, options, inequality_constraints, equality_constraints)
    181                 while start_idx < X_rnd.shape[0]:
    182                     end_idx = min(start_idx + batch_limit, X_rnd.shape[0])
--> 183                     Y_rnd_curr = acq_function(
    184                         X_rnd[start_idx:end_idx].to(device=device)
    185                     ).cpu()

\torch\nn\modules\module.py in _call_impl(self, *input, **kwargs)
   1100         if not (self._backward_hooks or self._forward_hooks or self._forward_pre_hooks or _global_backward_hooks
...
--> 434             new_kernel._parameters[param_name].data = param.__getitem__(index)
    435             ndim_removed = len(param.shape) - len(new_kernel._parameters[param_name].shape)
    436             new_batch_shape_len = len(self.batch_shape) - ndim_removed

IndexError: too many indices for tensor of dimension 1

System information

botorch 0.6.4 gpytorch 1.6.0 torch 1.10.2+cpu

[Balandat]

This looks very related to what's happening in #1287. Could you try the workaround proposed there to make sure that this is indeed the same issue?

[DavidWalz]

Yep, I can confirm that increasing the limit for non-lazy evaluation fixes the issue. Thanks!

[DavidWalz]

Duplicate of #1287.

[DavidWalz]

Mmh I think either this ticket or #1287 should stay open until the issue is fixed without changing gpytorch's defaults.

[esantorella]

This is now not erroring, probably due to #1547 (thanks @saitcakmak ). Here's the code I ran, with slightly updated syntax:

import torch
from botorch.acquisition.monte_carlo import qExpectedImprovement
from botorch.acquisition.objective import GenericMCObjective
from botorch.fit import fit_gpytorch_model
from botorch.models import SingleTaskGP
from botorch.optim.optimize import optimize_acqf
from botorch.sampling.normal import SobolQMCNormalSampler
from botorch.utils.multi_objective.scalarization import get_chebyshev_scalarization
from botorch.utils.sampling import sample_simplex
from gpytorch.mlls import ExactMarginalLogLikelihood

m = 2
d = 5
n = 1024

X = torch.linspace(0, 1, n * d, dtype=torch.double).reshape((n, d))
Y = torch.linspace(0, 1, n * m, dtype=torch.double).reshape((n, m))

model = SingleTaskGP(X, Y)
mll = ExactMarginalLogLikelihood(model.likelihood, model)
fit_gpytorch_model(mll)

sampler = SobolQMCNormalSampler(sample_shape=torch.Size([1024]))
weights = sample_simplex(d=m, n=1).squeeze()
objective = GenericMCObjective(get_chebyshev_scalarization(weights=weights, Y=Y))
bounds = torch.tensor([[0,] * d, [1,] * d], dtype=torch.double)

acq_func = qExpectedImprovement(
    model=model,
    objective=objective,
    best_f=objective(Y).max(),
    sampler=sampler,
)

candidates, acq_values = optimize_acqf(
    acq_function=acq_func,
    q=1,
    bounds=bounds,
    num_restarts=1,
    raw_samples=1024,
)