I have a problem where I need to fit a pair of parameters in a model to best fit observed data for hundreds of instances. I need to determine these parameters to high accuracy, rather than just finding a "reasonable optimum". Because of the cost of running the model (~1 to 20s), and the number of instances to fit, and the fact that my error function must be evaluated numerically (so no analytic derivative function) I was looking at Bayesian optimization as a suitable approach. My work is in Python, so I've been using the skopt package.
Although the parameters can vary significantly from one instance to the next, the general shape of the error surface (error as a function of mismatch of parameters) looks quite consistent. I have evaluated this systematically for one instance below:
There is some noise in this surface that introduces local minima, but there is clearly one global minimum. The surface is quite strongly peaked close to this minimum, as can be seen at higher resolution:
In the first image you can see the points sampled when I run gp_minimize on the error function using the default settings
from skopt import gp_minimize
params_init = res = gp_minimize(model_err, [(-10.,10.),(-10.,10.)])
I was surprised that there was little convergence towards the minimum, and I would appreciate any suggestions on how to improve this. I have already tried increasing the number of random samples to initiate the optimizer, and initiating with a 5 by 5 uniform grid of samples (I can't sample much finer without making the optimization process prohibitively long). Is there some way to make use of my knowledge of the general shape of the error surface within the optimizer? Or is my problem poorly suited to a GP-based approach?
gp_minimize(model_err, [(-10.,10.),(-10.,10.)], noise=1.0)or higher. $\endgroup$