Getting stuck in local minimum with noisy external objective function evaluations

[vsariola]

Hi,

I'm trying to use GPyOpt to optimize physical experiments, so I started following the example "5.10 External objective function evaluation". My problem is that as soon as I add noise to the objective function when using GP & EI, I tend to get stuck in a local minimum. Here's my code, almost identical to the example, but with noise:

import GPyOpt
from numpy.random import seed
import numpy as np
import matplotlib.pyplot as plt

seed(123)

func = GPyOpt.objective_examples.experiments1d.forrester(0.1) # add some noise

domain =[{'name': 'var1', 'type': 'continuous', 'domain': (0,1)}]

X_init = np.reshape(np.linspace(0,1,3),(-1,1))
Y_init = func.f(X_init)

iter_count = 40
current_iter = 0
X_step = X_init
Y_step = Y_init

while current_iter < iter_count:
    bo_step = GPyOpt.methods.BayesianOptimization(f = None, domain = domain, X = X_step, Y = Y_step, acquisition_type = 'EI')
    x_next = bo_step.suggest_next_locations()
    print(x_next)
    y_next = func.f(x_next)

    X_step = np.vstack((X_step, x_next))
    Y_step = np.vstack((Y_step, y_next))

    current_iter += 1

bo_step._compute_results()
print(bo_step.x_opt)
print(bo_step.fx_opt)

x = np.arange(0.0, 1.0, 0.01)
y = func.f(x)

plt.figure()
plt.plot(x, y)
for i, (xs, ys) in enumerate(zip(X_step, Y_step)):
    plt.plot(xs, ys, 'rD', markersize=2 + 4 * (i+1)/len(X_step))
plt.show()

I tend to get stuck at the local minima around 0.18 - 0.32, never reaching the true minimum near 0.75. What parameters are there in the model or acquisition function that could make the optimization behave better? Obviously I can make more initial experiments, but my real experiments will be really costly (~ 1 week per experiment) so it's worthwhile to get this right in the first place.

Thank you for your great work!

[emc2-2022]

You can try this:

    bo_step = GPyOpt.methods.BayesianOptimization(f=None, domain=domain, X=X_step, Y=Y_step,
                                                  exact_feval=True,
                                                  )

[apaleyes]

You are using Expected Improvement as an acquisition function, which tends to exploit more (see explore-exploit dilemma. Try a different acquisition function, for example entropy search (code for this one might be unstable though) or LCB, these are more explorative.

[matthew-hsr]

Is it true that the parameter "acquisition_jitter" also changes the exploration-exploitation preference?

[raff7]

E.I. is not ment to deal with noise, you need N.E.I. for that, which I believe is not implemented yet.

[ekalosak]

Bayesian optimization is built to deal with measurement noise @raff7 - I'm not sure I understand your comment. @vsariola - What's the signal-to-noise ratio? Do you have a plot of the test function available?

[raff7]

@ekalosak Bayesian optimization can deal with noisy observations, but E.I. cannot, that is because it assumes a "bestY" when computing the expected improvement for a certain set of parameters, and if the observations are noisy you don't have a well defined "bestY"

here a paper by facebook where they go in details on why E.I. will perform very poorly when the observations are noisy, and how they implemented a Noisy version, they also explain that E.I. will result in the model getting stuck in local optioma as in the question title: https://research.fb.com/wp-content/uploads/2018/08/Constrained-Bayesian-Optimization-with-Noisy-Experiments.pdf

[ekalosak]

@raff7 I see - that was an interesting read. I wonder what work would be required to get from EI_mcmc.py to Facebook's NEI? If their formulation were reduced to an un-constrained setting, it looks like equation (6) simplifies to

\alpha_NEI(x | D) = \integral_{f_n} \alpha_EI(x | f) * p(f | D) df_n

and the associated MCMC inference method further simplifies their Algorithm 1. Might be a good place to start an NEI implementation?

[raff7]

@ekalosak Yes in an unconstrained version of it it's basically just EI, repeated in a for loop N times computed on a GP with observations drawn from the mean and sd of every observation according to your original GP. As they say in the paper using classic Montecarlo simulation (with random samples) is slow, so i also added a low discrepancy sequence (Sobol) instead of random numbers. I implemented a NEI acquisition function based on your EI.py file The main function is this one,

def _compute_acq(self,x):
       observations = self.model.model.X.tolist()#Concatenate any pending observation to include them in calculation of Expected Improvement.
        for e in self.batch_elements:
            if(len(e)>0):
                observations.append(e.tolist())
        observations = np.array(observations)
        m,s = self.model.predict(observations)
        m = functools.reduce(operator.iconcat, m, [])
        s = functools.reduce(operator.iconcat, s , [])#Flatten m and s
        se = SobolEngine(dimension=len(observations),scramble=True)#Quasi-random number generator
        tks = np.array(se.draw(self.N_QMC))
        NEI = 0
        gpModel = GPModel(self.space,None, exact_feval=True,verbose=False, ARD=self.model.ARD,optimize_restarts=self.opt_restarts)
        for k in range(self.N_QMC):#quasi monte carlo integration (QMC)
                tk = tks[k]
                Fn = (norm.ppf(tk,loc=m,scale=s))
                Fn = np.array([np.array([xn]) for xn in Fn])
                gpModel.updateModel(observations,Fn,None,None)
                NEI += self._compute_EI_acq(gpModel,x)/self.N_QMC
        return NEI

And i replicated the experiments in facebook paper.. it works really well

Also The implementation i have allows for batch predictions, instead of having a single experiment you can have a whole batch of experiments and perform them at once (you can do it just concatenating to the Observations the X you still didn't test)

[ekalosak]

Brilliant @raff7 , any chance you can make that into a PR? I'm sure others would appreciate greatly and would contribute to its maintenance.

[raff7]

@ekalosak ye sure I can do that I'll need to clean the code quite a bit first tho because I just downloaded the repository offline and modified many parts of it, I need to make sure it works in the original version of the code

[raff7]

@ekalosak ok I sent a PR, i just added the class NEI.py and modified the init.py + argument_manager.py to accommodate for a new acquisition function. I quickly tested it and it seem to work. feel free to do more testing. Also, a note about NEI. because it uses MC integration is considerably slower then EI, so it should only be used in cases with enough Noise that the advantages of NEI outweigh the disadvantages