`suggest` does not work as expected in constrained optimization

[yichenfu]

Hi, I recently run into issues where suggest does not work as expected in constrained optimization. Based on Gardner et. al., 2014, in constrained cases, suggest should return the argmax of acquisition*p_constrain. (Assuming acquisition function is always positive, which is satisfied below.) However, in the following case, it seems like suggest only returns the argmax of acquisition but ignores p_constrain. Below is the plot for such a case.

The circles are data points the algorithm currently has, and the star is the next point suggest provides. The white points satisfy the constraint, and the red points do not. The contour of the GP prediction of target, acquisition function, p_constraint, and acquisition*p_constraint is shown in the figure, where we can clearly see that the next point maximize the acquisition function but ignores the constrain.

My code to reproduce the figure is attached below. Did I do anything wrong in the constrained optimization? Thanks in advance for taking your time!

import matplotlib.pyplot as plt
import numpy as np
from bayes_opt import UtilityFunction, BayesianOptimization
from scipy.optimize import NonlinearConstraint

# Synthetic data
params = np.array([[ 0.55, 0.79],[-0.65,-1.5 ],[-1.25,-0.23],[-0.25,-1.02],[-0.01, 0.13],[-1.96,-0.65],[-1.06, 0.43],[ 0.69,-1.86],[ 0.96,-0.47],[-0.98,-1.18],[-1.69, 0.95],[ 0.07,-1.26],[-0.35, 0.57],[-1.55,-1.65],[-0.75,-0.08],[ 0.26,-0.79],[ 0.39,-0.01],[-0.85,-0.7 ],[-1.46, 0.48],[-0.41,-1.72],[ 0.2, 0.84],[-1.79,-1.36],[-0.9,-0.37],[ 0.9,-1.07],[ 0.8, 0.33],[-1.2,-1.98],[-1.9, 0.25],[-0.09,-0.56],[-0.14,-0.14],[-1.38,-0.94],[-0.58, 0.71]])
target = np.array([-1.13,-1.04,-0.02,-0.74,-0.57,-0.45,-0.94,-1.1,-0.06,-0.87,-1.15,-0.88,-1.02,-1.1,-0.15,-0.49,-0.34,-0.49,-0.96,-1.12,-1.09,-0.98,-0.11,-0.6,-0.89,-1.19,-0.73,-0.3,-0.08,-0.71,-1.08])
constraint_value = np.array([0.56,0.01,0.27,0.03,0.53,0.09,0.69,0.,0.15,0.02,0.53,0.01,0.63,0.,0.36,0.06,0.42,0.08,0.67,0.,0.35,0.01,0.19,0.02,0.66,0.,0.61,0.12,0.32,0.04,0.58])

# Defind optimizer 
pbounds = {'x':[-2,1],'y':[-2,1]}
constraint = NonlinearConstraint(None,0.5,np.inf)
optimizer = BayesianOptimization(f=None, constraint=constraint, pbounds=pbounds)
print(optimizer.is_constrained)

# register data
for i in range(len(target)):
    optimizer.register(params = params[i], target = target[i], constraint_value = constraint_value[i])

# define acquisition and suggest next point
acq = UtilityFunction(kind="ei", xi=0.)
s = optimizer.suggest(acq)

# ========== plot the figure ==========
param_bounds = np.array([[-2,1],[-2,1]])
x = np.linspace(param_bounds[0,0], param_bounds[0,1],100)
y = np.linspace(param_bounds[1,0], param_bounds[1,1],100)
xy = np.array([[x_i,y_j] for y_j in y for x_i in x])
X,Y = np.meshgrid(x,y)

fig = plt.figure(figsize=(10,8))
grid = fig.add_gridspec(ncols=2,nrows=2)

fig.subplots_adjust(wspace=0.2,hspace=0.3)

# Estimation of target
ax0 = plt.subplot(grid[0,0])
target_est = optimizer._gp.predict(xy).reshape(X.shape)
plt.contourf(X,Y,target_est,levels=20)
plt.title('Target',fontsize=20)
plt.colorbar();

# Estimation of acquisition function
ax1 = plt.subplot(grid[0,1])
acq_est = acq.utility(xy,optimizer._gp,optimizer.space.target.max()).reshape(X.shape)
plt.contourf(X,Y,acq_est,levels=20)
plt.title('Acquisition',fontsize=20)
plt.colorbar()

# Estimation of p_constraint
ax2 = plt.subplot(grid[1,0])
p_constraint_est = optimizer.constraint.predict(xy).reshape(X.shape)
plt.contourf(X,Y,p_constraint_est,levels=20)
plt.title('p_constraint',fontsize=20)
plt.colorbar();

# print acquisition * p_constraint
ax3 = plt.subplot(grid[1,1])
plt.contourf(X,Y,acq_est * p_constraint_est,levels=20)
plt.title('p_const * acq',fontsize=20)
plt.colorbar();

# show all data points
res = optimizer.res
keys = list(res[0]['params'].keys())
x_ = np.array([r["params"][keys[0]] for r in res])
y_ = np.array([r["params"][keys[1]] for r in res])
a_ = np.array([r["allowed"] for r in res])

for ax in [ax0,ax1,ax2,ax3]:
    ax.scatter(x_[a_],y_[a_],c='white',s=20,edgecolors='black')
    ax.scatter(x_[~a_],y_[~a_],c='red',s=20,edgecolors='black');
    ax.scatter(s['x'],s['y'],c='magenta',s=100,edgecolors='black',marker='*');

[bwheelz36]

@till-m this one is above my paygrade!

[till-m]

Hey @yichenfu,

I had a quick look yesterday and couldn't figure out what's going on. I will have check back when I get the chance, but I'm also currently swamped with my actual job so this might take some time, apologies.

[till-m]

Hi @yichenfu,

I found the problem. When you define the EI utility function, you're using the maximum value of targets found: optimizer._gp,optimizer.space.target.max(). Internally, however it uses the best allowed result, as returned by optimizer._space._target_max().

import matplotlib.pyplot as plt
import numpy as np
from bayes_opt import UtilityFunction, BayesianOptimization
from scipy.optimize import NonlinearConstraint

# Synthetic data
params = np.array([[ 0.55, 0.79],[-0.65,-1.5 ],[-1.25,-0.23],[-0.25,-1.02],[-0.01, 0.13],[-1.96,-0.65],[-1.06, 0.43],[ 0.69,-1.86],[ 0.96,-0.47],[-0.98,-1.18],[-1.69, 0.95],[ 0.07,-1.26],[-0.35, 0.57],[-1.55,-1.65],[-0.75,-0.08],[ 0.26,-0.79],[ 0.39,-0.01],[-0.85,-0.7 ],[-1.46, 0.48],[-0.41,-1.72],[ 0.2, 0.84],[-1.79,-1.36],[-0.9,-0.37],[ 0.9,-1.07],[ 0.8, 0.33],[-1.2,-1.98],[-1.9, 0.25],[-0.09,-0.56],[-0.14,-0.14],[-1.38,-0.94],[-0.58, 0.71]])
target = np.array([-1.13,-1.04,-0.02,-0.74,-0.57,-0.45,-0.94,-1.1,-0.06,-0.87,-1.15,-0.88,-1.02,-1.1,-0.15,-0.49,-0.34,-0.49,-0.96,-1.12,-1.09,-0.98,-0.11,-0.6,-0.89,-1.19,-0.73,-0.3,-0.08,-0.71,-1.08])
constraint_value = np.array([0.56,0.01,0.27,0.03,0.53,0.09,0.69,0.,0.15,0.02,0.53,0.01,0.63,0.,0.36,0.06,0.42,0.08,0.67,0.,0.35,0.01,0.19,0.02,0.66,0.,0.61,0.12,0.32,0.04,0.58])

# Defind optimizer 
pbounds = {'x':[-2,1],'y':[-2,1]}
constraint = NonlinearConstraint(None,0.5,np.inf)
optimizer = BayesianOptimization(f=None, constraint=constraint, pbounds=pbounds)
print(optimizer.is_constrained)

# register data
for i in range(len(target)):
    optimizer.register(params = params[i], target = target[i], constraint_value = constraint_value[i])

# define acquisition and suggest next point
acq = UtilityFunction(kind="ei", xi=0.)
s = optimizer.suggest(acq)

# ========== plot the figure ==========
param_bounds = np.array([[-2,1],[-2,1]])
x = np.linspace(param_bounds[0,0], param_bounds[0,1],100)
y = np.linspace(param_bounds[1,0], param_bounds[1,1],100)
xy = np.array([[x_i,y_j] for y_j in y for x_i in x])
X,Y = np.meshgrid(x,y)

fig = plt.figure(figsize=(10,8))
grid = fig.add_gridspec(ncols=2,nrows=2)

fig.subplots_adjust(wspace=0.2,hspace=0.3)

# Estimation of target
ax0 = plt.subplot(grid[0,0])
target_est = optimizer._gp.predict(xy).reshape(X.shape)
plt.contourf(X,Y,target_est,levels=20)
plt.title('Target',fontsize=20)
plt.colorbar();

# Estimation of acquisition function
ax1 = plt.subplot(grid[0,1])
acq_est = acq.utility(xy,optimizer._gp,optimizer._space._target_max()).reshape(X.shape)
plt.contourf(X,Y,acq_est,levels=20)
plt.title('Acquisition',fontsize=20)
plt.colorbar()

# Estimation of p_constraint
ax2 = plt.subplot(grid[1,0])
p_constraint_est = optimizer.constraint.predict(xy).reshape(X.shape)
plt.contourf(X,Y,p_constraint_est,levels=20)
plt.title('p_constraint',fontsize=20)
plt.colorbar();

# print acquisition * p_constraint
ax3 = plt.subplot(grid[1,1])
plt.contourf(X,Y,acq_est * p_constraint_est,levels=20)
plt.title('p_const * acq',fontsize=20)
plt.colorbar();

# show all data points
res = optimizer.res
keys = list(res[0]['params'].keys())
x_ = np.array([r["params"][keys[0]] for r in res])
y_ = np.array([r["params"][keys[1]] for r in res])
a_ = np.array([r["allowed"] for r in res])

best_allowed = optimizer.space.max()['params']
best_overall = optimizer.space.array_to_params(optimizer.space.params[optimizer.space.target.argmax()])
for ax in [ax0,ax1,ax2,ax3]:
    ax.scatter(x_[a_],y_[a_],c='white',s=20,edgecolors='black')
    ax.scatter(x_[~a_],y_[~a_],c='red',s=20,edgecolors='black');
    ax.scatter(s['x'],s['y'],c='magenta',s=100,edgecolors='black',marker='*', label='suggestion');
    ax.scatter(best_allowed['x'],best_allowed['y'],c='white',s=100,edgecolors='black',marker='P', label='best allowed');
    ax.scatter(best_overall['x'],best_overall['y'],c='red',s=100,edgecolors='black',marker='D', label='best overall');
ax3.legend(loc='lower right')

Hope this clears things up.

[yichenfu]

Hi @till-m ,

Thank you so much for your reply!

Whe running your code, I got an error saying 'TargetSpace' object has no attribute '_target_max'. Then I realize that I was using bayes_opt with version 1.4.3 instead of the latest version. After I upgrading the version to 1.5 (which was just realsed), the code runs well and the "glitch" is gone. Now suggest correctly captures the constrain just as expected, like the figure you showed.

I believe this issue has been solved in the newest version, so I can close it.

Thanks again for your effort!

[till-m]

Glad to hear everything is clear now!

A nicer function for this was added in v1.5.0, but the behaviour was already like this back then, see this function.