HGDL Constrained Optimization of Schwefel’s Function
In this script, we show how HGDL is used for constrained optimization. Unconstrained optimization is simpler and automatically included.
First some function to make nice plots
%load_ext autoreload
%autoreload 2
import numpy as np
import plotly.graph_objects as go
def plot(x,y,z,data = None, constr = None):
fig = go.Figure()
fig.add_trace(go.Surface(x = x, y = y,z=z))
if data is not None:
fig.add_trace(go.Scatter3d(x=data[:,0], y=data[:,1], z=data[:,2] + 50,
mode='markers'))
if constr is not None:
fig.add_trace(go.Scatter3d(x=constr[:,0], y=constr[:,1], z=constr[:,2],
mode='markers'))
fig.update_layout(title='Surface Plot', autosize=True,
width=800, height=800, font=dict(
family="Courier New, monospace",
size=18),
margin=dict(l=65, r=50, b=65, t=90))
fig.show()
def make_plot(data = None, constraint1 = None, constraint2 = None):
x1,x2 = np.linspace(-500,500,100),np.linspace(-500,500,100)
x_pred = np.transpose([np.tile(x1, len(x2)), np.repeat(x2, len(x1))])
r1 = np.sqrt(160000.)
r2 = np.sqrt(40000.)
c1,c2 = r1*np.cos(np.linspace(0,2.*np.pi,100)),r1*np.sin(np.linspace(0,2.*np.pi,100))
d1,d2 = r2*np.cos(np.linspace(0,2.*np.pi,100)),r2*np.sin(np.linspace(0,2.*np.pi,100))
c3 = np.zeros((len(c2)))
d3 = np.zeros((len(c2)))
x1,x2 = np.meshgrid(x1,x2)
z = np.zeros((10000))
zc1 = np.zeros((10000))
zc2 = np.zeros((10000))
for i in range(10000):
z[i] = schwefel(x_pred[i], 1, 1)
if constraint1: zc1[i] = constraint1(x_pred[i])
if constraint2: zc2[i] = constraint2(x_pred[i])
for i in range(len(c1)):
c3[i] = schwefel(np.array([c1[i],c2[i]]))
d3[i] = schwefel(np.array([d1[i],d2[i]]))
plot(x1,x2,z.reshape(100,100).T, data = data,
constr = np.row_stack([np.column_stack([c1,c2,c3]),np.column_stack([d1,d2,d3])]))
Defining the Constraints and some Bounds
Keep in mind that not all local optimizers allow any combination of bounds and constraints Visit https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html for more information on that
bounds = np.array([[-500,500],[-500,500]])
def g1(x):
return (np.linalg.norm(x)**2/10.0) - 16000.0
def g2(x):
return (np.linalg.norm(x)**2/10.0) - 4000.0
Now we import HGDL and run a constrained optimization
from hgdl.hgdl import HGDL as hgdl
from hgdl.support_functions import *
import time
#example arguments that will be passed to the function, the gradient and the Hessian function
a = 5
b = 6
#constraint definitions form scipy
from scipy.optimize import NonlinearConstraint
nlc = NonlinearConstraint(g1, -np.inf, 0)
nlc = NonlinearConstraint(g2, 0,np.inf)
a = hgdl(schwefel, schwefel_gradient, bounds,
hess = None, ##if this is None, the Hessian will be approximated if the local optimizer needs it
global_optimizer = "random", #there are a few options to choose from for the global optimizer
#global_optimizer = "genetic",
local_optimizer = "dNewton", #dNewton is an example and will be changed automatically to "SLSQP" because constraints are used
number_of_optima = 30000, #the number fo optima that will be stored and used for deflation
args = (a,b), num_epochs = 1000, #the number of total epochs. Since this is an asynchronous algorithms, this number can be very high
constraints = (nlc,) #the constraints
)
a.optimize(x0=None)
time.sleep(10)
/home/docs/checkouts/readthedocs.org/user_builds/hgdl/envs/latest/lib/python3.9/site-packages/hgdl/hgdl.py:134: UserWarning: Warning: dNewton will not adhere to bounds. It is recommended to formulate your objective function such that it is defined on R^N by simple non-linear transformations.
warnings.warn("Warning: dNewton will not adhere to bounds. It is recommended to formulate your objective function such that it is defined on R^N by simple non-linear transformations.")
/home/docs/checkouts/readthedocs.org/user_builds/hgdl/envs/latest/lib/python3.9/site-packages/hgdl/hgdl.py:137: UserWarning: Constraints provided, local optimizer changed to 'SLSQP'
warnings.warn("Constraints provided, local optimizer changed to 'SLSQP'")
/home/docs/checkouts/readthedocs.org/user_builds/hgdl/envs/latest/lib/python3.9/site-packages/scipy/optimize/_minimize.py:565: RuntimeWarning: Method SLSQP does not use Hessian information (hess).
warn('Method %s does not use Hessian information (hess).' % method,
/home/docs/checkouts/readthedocs.org/user_builds/hgdl/envs/latest/lib/python3.9/site-packages/scipy/optimize/_optimize.py:404: RuntimeWarning: Values in x were outside bounds during a minimize step, clipping to bounds
warnings.warn("Values in x were outside bounds during a "
res = a.get_latest()
for entry in res: print(entry)
{'x': array([420.96869994, 420.96881868]), 'f(x)': 2.5456064349782537e-05, 'classifier': 'minimum', 'Hessian eigvals': array([0.25236707, 0.2523671 ]), 'df/dx': array([-1.17159581e-05, 1.82504597e-05]), '|df/dx|': 2.1687391595558583e-05, 'radius': 3.96248214619082}
{'x': array([ 420.96873994, -302.524948 ]), 'f(x)': 118.43836006959498, 'classifier': 'minimum', 'Hessian eigvals': array([0.25236706, 0.25328927]), 'df/dx': array([-1.61919261e-06, -3.13676379e-06]), '|df/dx|': 3.5300243293411744e-06, 'radius': 3.9624823449145543}
{'x': array([-302.52487368, 420.96867488]), 'f(x)': 118.43836007070081, 'classifier': 'minimum', 'Hessian eigvals': array([0.25328924, 0.25236705]), 'df/dx': array([ 1.56867941e-05, -1.80380122e-05]), '|df/dx|': 2.3904924004299042e-05, 'radius': 3.962482372805606}
{'x': array([203.81430055, 420.96870517]), 'f(x)': 217.13969484621202, 'classifier': 'minimum', 'Hessian eigvals': array([0.25487081, 0.25236703]), 'df/dx': array([ 1.22083086e-05, -1.03956144e-05]), '|df/dx|': 1.6034699789225413e-05, 'radius': 3.962482665661676}
{'x': array([420.968747 , 203.81418887]), 'f(x)': 217.13969484622396, 'classifier': 'minimum', 'Hessian eigvals': array([0.25236706, 0.25487074]), 'df/dx': array([ 1.61561543e-07, -1.62564221e-05]), '|df/dx|': 1.625722490255197e-05, 'radius': 3.9624822821596886}
{'x': array([-302.52493208, -302.52493719]), 'f(x)': 236.8766946840102, 'classifier': 'minimum', 'Hessian eigvals': array([0.25328926, 0.25328926]), 'df/dx': array([ 8.93367525e-07, -4.00129678e-07]), '|df/dx|': 9.788816550153248e-07, 'radius': 3.9480553114919847}
{'x': array([ 420.96874253, -124.8293765 ]), 'f(x)': 296.10673921370426, 'classifier': 'minimum', 'Hessian eigvals': array([0.25236706, 0.2579167 ]), 'df/dx': array([-9.65513559e-07, -5.17831090e-06]), '|df/dx|': 5.2675535323262e-06, 'radius': 3.962482320509884}
{'x': array([-124.8293633 , 420.96887463]), 'f(x)': 296.10673921573266, 'classifier': 'minimum', 'Hessian eigvals': array([0.25791668, 0.2523671 ]), 'df/dx': array([-1.77466883e-06, 3.23720583e-05]), '|df/dx|': 3.24206664216044e-05, 'radius': 3.962481673786233}
{'x': array([-302.5249374 , 203.81425056]), 'f(x)': 335.57802946014425, 'classifier': 'minimum', 'Hessian eigvals': array([0.25328926, 0.25487079]), 'df/dx': array([-9.44954047e-06, 6.36246836e-05]), '|df/dx|': 6.432257910380985e-05, 'radius': 3.948055242271176}
{'x': array([ 203.81419486, -302.52499945]), 'f(x)': 335.57802946108507, 'classifier': 'minimum', 'Hessian eigvals': array([0.25487075, 0.25328928]), 'df/dx': array([-1.47281135e-05, -1.61705648e-05]), '|df/dx|': 2.1872459744635076e-05, 'radius': 3.9480549359691253}
{'x': array([420.96873564, 65.54787437]), 'f(x)': 355.34793077603774, 'classifier': 'minimum', 'Hessian eigvals': array([0.25236705, 0.26492048]), 'df/dx': array([-2.70633505e-06, 2.45819790e-06]), '|df/dx|': 3.65608893933296e-06, 'radius': 3.9624823832647507}
{'x': array([ 65.54784503, 420.96875312]), 'f(x)': 355.34793077607094, 'classifier': 'minimum', 'Hessian eigvals': array([0.26492044, 0.25236706]), 'df/dx': array([-5.31490717e-06, 1.70544021e-06]), '|df/dx|': 5.581824481315708e-06, 'radius': 3.962482224634397}
{'x': array([420.96880213, -25.87749467]), 'f(x)': 394.8999525057472, 'classifier': 'minimum', 'Hessian eigvals': array([0.25236709, 0.28661014]), 'df/dx': array([ 1.40734309e-05, -2.21622625e-05]), '|df/dx|': 2.625313957690689e-05, 'radius': 3.9624817696616934}
{'x': array([-124.82936334, -302.52496024]), 'f(x)': 414.5450738281712, 'classifier': 'minimum', 'Hessian eigvals': array([0.25791668, 0.25328928]), 'df/dx': array([-1.78443476e-06, -6.23707967e-06]), '|df/dx|': 6.487323806889281e-06, 'radius': 3.9480549532743248}
{'x': array([-302.52489901, -124.8292712 ]), 'f(x)': 414.5450738291944, 'classifier': 'minimum', 'Hessian eigvals': array([0.25328923, 0.2579166 ]), 'df/dx': array([9.26993547e-06, 2.19795442e-05]), '|df/dx|': 2.3854393024516587e-05, 'radius': 3.94805572854741}
{'x': array([420.96876934, 5.23920473]), 'f(x)': 415.0376111023545, 'classifier': 'minimum', 'Hessian eigvals': array([0.25236707, 0.40385525]), 'df/dx': array([5.79929101e-06, 2.19413164e-06]), '|df/dx|': 6.200483032745995e-06, 'radius': 3.962482074720008}
{'x': array([203.8141565 , 203.81423792]), 'f(x)': 434.279364237484, 'classifier': 'minimum', 'Hessian eigvals': array([0.25487072, 0.25487078]), 'df/dx': array([-2.45044561e-05, -3.75501400e-06]), '|df/dx|': 2.4790492121139397e-05, 'radius': 3.9235578528389694}
{'x': array([-302.52492771, 65.54786309]), 'f(x)': 473.7862653904581, 'classifier': 'minimum', 'Hessian eigvals': array([0.25328925, 0.26492046]), 'df/dx': array([ 2.00109399e-06, -5.29953476e-07]), '|df/dx|': 2.070079186312567e-06, 'radius': 3.948055366868634}
{'x': array([ 65.54781926, -302.52487938]), 'f(x)': 473.7862653911284, 'classifier': 'minimum', 'Hessian eigvals': array([0.26492039, 0.25328925]), 'df/dx': array([-1.21417441e-05, 1.42425058e-05]), '|df/dx|': 1.8715526165816096e-05, 'radius': 3.948055496657661}
{'x': array([ 203.81425329, -124.8293453 ]), 'f(x)': 513.2464086042393, 'classifier': 'minimum', 'Hessian eigvals': array([0.25487078, 0.25791666]), 'df/dx': array([1.62917781e-07, 2.86813143e-06]), '|df/dx|': 2.872754794991073e-06, 'radius': 3.923556844463549}
{'x': array([-124.82935895, 203.81426394]), 'f(x)': 513.2464086042404, 'classifier': 'minimum', 'Hessian eigvals': array([0.25791668, 0.25487078]), 'df/dx': array([-6.52215661e-07, 2.87731699e-06]), '|df/dx|': 2.950311563570122e-06, 'radius': 3.9235568427544383}
{'x': array([ -25.87738614, -302.52490261]), 'f(x)': 513.3382871192132, 'classifier': 'minimum', 'Hessian eigvals': array([0.28660976, 0.25328923]), 'df/dx': array([8.94577950e-06, 8.35812831e-06]), '|df/dx|': 1.2242764379321665e-05, 'radius': 3.9480556835538745}
{'x': array([ 65.54787753, 203.814248 ]), 'f(x)': 572.487600166608, 'classifier': 'minimum', 'Hessian eigvals': array([0.26492049, 0.25487079]), 'df/dx': array([6.26591541e-05, 1.81193862e-05]), '|df/dx|': 6.522638843814873e-05, 'radius': 3.9235567453351456}
{'x': array([203.81426151, -25.87741215]), 'f(x)': 612.0396218950847, 'classifier': 'minimum', 'Hessian eigvals': array([0.25487079, 0.28660985]), 'df/dx': array([2.2591487e-06, 1.4901038e-06]), '|df/dx|': 2.706318934236745e-06, 'radius': 3.923556692352725}
{'x': array([-25.87743925, 203.81421628]), 'f(x)': 612.0396218953081, 'classifier': 'minimum', 'Hessian eigvals': array([0.28660995, 0.25487076]), 'df/dx': array([-6.27771063e-06, -9.26950275e-06]), '|df/dx|': 1.1195237021420129e-05, 'radius': 3.923557143557899}
{'x': array([203.81425302, 5.23919544]), 'f(x)': 632.1772804928578, 'classifier': 'minimum', 'Hessian eigvals': array([0.25487079, 0.40385528]), 'df/dx': array([ 9.43866808e-08, -1.55905716e-06]), '|df/dx|': 1.561911673513773e-06, 'radius': 3.9235566513340783}
{'x': array([ 5.2392053 , 203.81422478]), 'f(x)': 632.177280492961, 'classifier': 'minimum', 'Hessian eigvals': array([0.40385525, 0.25487077]), 'df/dx': array([ 2.4230036e-06, -7.1024131e-06]), '|df/dx|': 7.504346622145413e-06, 'radius': 3.9235569863197206}
res = a.kill_client()
2023-09-29 22:45:16,572 - distributed.worker.state_machine - WARNING - Async instruction for <Task cancelled name='execute(hgdl-b6d02f5aace9cb386f1fb2d05149742e)' coro=<Worker.execute() done, defined at /home/docs/checkouts/readthedocs.org/user_builds/hgdl/envs/latest/lib/python3.9/site-packages/distributed/worker_state_machine.py:3609>> ended with CancelledError
2023-09-29 22:45:16,579 - distributed.worker.state_machine - WARNING - Async instruction for <Task cancelled name='execute(local_method-f6b6184f36b49568d8d60fc4926da551)' coro=<Worker.execute() done, defined at /home/docs/checkouts/readthedocs.org/user_builds/hgdl/envs/latest/lib/python3.9/site-packages/distributed/worker_state_machine.py:3609>> ended with CancelledError
Making a Plot
You should see the constraints and the foudn optima. If everything worked, the found points are in between the two constraints.
data = [np.append(entry["x"],entry["f(x)"]) for entry in res]
make_plot(data = np.array(data), constraint1 = g1, constraint2 = g2)