<jupyter_start><jupyter_text>Optimization Hello, and welcome to our tutorial on optimization. Here, we will explore the four different optimizers which already have `Tequila` interfaces: a native GD optimizer, alongside interfaces for `SciPy`, `GPyOpt`, and `Phoenics`.<jupyter_code>### start at the start: import statements!
import tequila as tq
import numpy as np
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Overview**How to optimize an `Objective`:** In `tequila`, optimizers are accessed in several ways. They may be instantiated as objects directly, which can then be called; they are also accessible through `tq.minimize`. Either of these methods have an obligatory argument: an `Objective`. `tq.minimize` also requires you supply a `method`, which must be a string; the call methods of the GD, `SciPy`, and `GPyOpt` optimizers accept this key as well. As keywords to any optimization call, you can pass down the usual compilation kwargs for quantum simulation in `Tequila`:* `backend`, a string indicating which quantum simulator to use* `samples`, an int indicating how many shots of the circuit to measure (None means full wf simulation)* `device`, (generally) a string indicating which (real or emulated) quantum computer to sample from (requires samples be specified),* `noise`, the NoiseModel object to apply to the circuits simulated (see the tutorial on noise). additional keywords you might use include:* `variables`,  a `list` of the `Variable`s you want to optimize (the default is all of them).* `initial_values`, which gives a start point to optimization. If you supply arguments to `variables`, you also need to supply arguments to `initial_values` so that all non-optimized parameters have a value to use. The default is random initialization (not recomended) * `gradient`, which specifies which type of gradient will be used : analytical gradients (Default): gradient=None, numerical gradients: gradient={'method':'2-point', "stepsize":1.e-4}, custom gradient objectives: gradient={tq.Variable:tq.Objective}, module specific gradients: gradient="2-point" (to use for example `SciPy` finite difference stencils).* `silent`, silence outputsSome of the optimizers take more, or different, keywords from the others, so check the documentation for each one. In case the optimizer has some degree of verbosity (currently, they all do), you can deactivate this with `silent=True`. The following optimization methods are available on your system in `Tequila`:<jupyter_code>tq.optimizers.show_available_optimizers()<jupyter_output>available methods for optimizer modules found on your system:
method               | optimizer module
--------------------------
NELDER-MEAD          | scipy
COBYLA               | scipy
POWELL               | scipy
SLSQP                | scipy
L-BFGS-B             | scipy
BFGS                 | scipy
CG                   | scipy
TNC                  | scipy
TRUST-KRYLOV         | scipy
NEWTON-CG            | scipy
DOGLEG               | scipy
TRUST-NCG            | scipy
TRUST-EXACT          | scipy
TRUST-CONSTR         | scipy
adam                 | gd
adagrad              | gd
adamax               | gd
nadam                | gd
sgd                  | gd
momentum             | gd
nesterov             | gd
rmsprop              | gd
rmsprop-nesterov     | gd
Supported optimizer modules:  ['scipy', 'phoenics', 'gpyopt', 'gd']
Installed optimizer modules:  ['scipy', 'gd']<jupyter_text>We will use two different `Objective`s for optimization in this tutorial. The first of these is a two qubit expectation value with the tractable but non trivial hamiltonian $[Y(0)+Qm(0)]\otimes X(1)$, where $Qm=\frac{1}{2} (I - Z)$, the projector onto the |1> state.<jupyter_code>### optimizing the circuit in terms of pi makes the result of the optimization easier to interpret.

a = tq.Variable(name="a")*tq.numpy.pi
b = tq.Variable(name="b")*tq.numpy.pi
c = tq.Variable(name="c")*tq.numpy.pi
d = tq.Variable(name='d')*tq.numpy.pi

U1 = tq.gates.H(target=[0])
U1 += tq.gates.H(target=1)
U1 += tq.gates.Ry(target=0, angle=a)
U1 += tq.gates.Rz(target=1, angle=b)
U1 += tq.gates.Z(target=1,control=0)
U1 += tq.gates.Rx(target=0, angle=c)
U1 += tq.gates.Rx(target=1,angle=d)
U1 += tq.gates.Z(target=1,control=0)


### once we have a circuit, we pick a hamiltonian to optimize over
H1=(tq.paulis.Y(0)+tq.paulis.Qm(0))*tq.paulis.X(1)
O1=tq.ExpectationValue(U=U1,H=H1)

### we use the .draw function to pretty-print circuits via backend printers.
print('We will optimize the following objective: \n')
tq.draw(O1,backend='qiskit')<jupyter_output>We will optimize the following objective: 

Objective with 1 unique expectation values
variables = [a, b, c, d]
types     = not compiled

Expectation Value 0
Hamiltonian :  (+1.0000Y(0)X(1)+0.5000X(1)-0.5000Z(0)X(1),)
variables :  [a, b, c, d]
circuit:

     ┌───┐┌───────────────┐   ┌───────────────┐   
q_0: ┤ H ├┤ RY(f((a,))_0) ├─■─┤ RX(f((c,))_2) ├─■─
     ├───┤├───────────────┤ │ ├───────────────┤ │ 
q_1: ┤ H ├┤ RZ(f((b,))_1) ├─■─┤ RX(f((d,))_3) ├─■─
     └───┘└───────────────┘   └───────────────┘   
c: 2/═════════════════════════════════════════════<jupyter_text>Our second `Objective`, O2, will measure a 3-qubit circuit with respect to the Hamiltonian $Y(0)\otimes X(1) \otimes Y(2)$<jupyter_code>### this time, don't scale by pi

H2 = tq.paulis.Y(0)*tq.paulis.X(1)*tq.paulis.Y(2)
U2 = tq.gates.Ry(tq.numpy.pi/2,0) +tq.gates.Ry(tq.numpy.pi/3,1)+tq.gates.Ry(tq.numpy.pi/4,2)
U2 += tq.gates.Rz('a',0)+tq.gates.Rz('b',1)
U2 += tq.gates.CNOT(control=0,target=1)+tq.gates.CNOT(control=1,target=2)
U2 += tq.gates.Ry('c',1) +tq.gates.Rx('d',2)
U2 += tq.gates.CNOT(control=0,target=1)+tq.gates.CNOT(control=1,target=2)
O2 = tq.ExpectationValue(H=H2, U=U2)

print('We will optimize the following objective: \n')
tq.draw(O2, backend="qiskit")<jupyter_output>We will optimize the following objective: 

Objective with 1 unique expectation values
variables = [a, b, c, d]
types     = not compiled

Expectation Value 0
Hamiltonian :  (+1.0000Y(0)X(1)Y(2),)
variables :  [a, b, c, d]
circuit:

     ┌──────────┐┌───────────────┐                                     
q_0: ┤ RY(pi/2) ├┤ RZ(f((a,))_0) ├──■──────────────────────────■───────
     ├──────────┤├───────────────┤┌─┴─┐     ┌───────────────┐┌─┴─┐     
q_1: ┤ RY(pi/3) ├┤ RZ(f((b,))_1) ├┤ X ├──■──┤ RY(f((c,))_2) ├┤ X ├──■──
     ├──────────┤└───────────────┘└───┘┌─┴─┐├───────────────┤└───┘┌─┴─┐
q_2: ┤ RY(pi/4) ├──────────────────────┤ X ├┤ RX(f((d,))_3) ├─────┤ X ├
     └──────────┘                      └───┘└───────────────┘     └───┘
c: 3/══════════════════════════════════════════════════════════════════<jupyter_text>Local OptimizersWe will begin this tutorial by focusing on local optimizers. By local optimization, we mean the any optimization schema where the suggested parameters at step t are always a transformation of the parameters suggested at step t-1. This includes a large number of the standard optimization techniques in use today, like gradient descent.     `Tequila` comes with two local optimizers: a native gradient descent optimizer, implementing a number of the most popular gradient descent algorithms used in classical machine learning, as well as a plugin for the `SciPy` package (which is installed alongside `Tequila`), which allows the use of a number of gradient-free, gradient-based, and hessian-based optimization methods. The GD Optimizerwe will start this tutorial by looking at the GD optimizer. Here is an overview over the available optimization methods.<jupyter_code>tq.show_available_optimizers(module="gd")<jupyter_output>available methods for optimizer module gd
method               | optimizer module
--------------------------
adam                 | gd
adagrad              | gd
adamax               | gd
nadam                | gd
sgd                  | gd
momentum             | gd
nesterov             | gd
rmsprop              | gd
rmsprop-nesterov     | gd<jupyter_text>As one sees, a variety of methods are available for optimization. Here, 'sgd' refers to the standard gradient descent algorithm, without momentum. like all tequila optimizers, the GD optimizer has a minimize function and most of the arguments are the same. However, there is one important difference: the GD optimizer takes a learning rate, lr. This parameter mediates step size in all of the GD optimizer methods; it is a positive float which scales the step in the direction of the gradient.  We will now optimize O1, our two-qubit expectation value, choosing starting angles equivalent to $\frac{1}{4}\pi$ for all four variables, and optimizing via the ['Adam'](https://towardsdatascience.com/_adam-latest-trends-in-deep-learning-optimization-6be9a291375c) method.<jupyter_code>init={'a':0.25,'b':0.25,'c':0.25,'d':0.25}
lr=0.1

### For even more fun, try using sampling with the samples keyword, 
### or pick your favorite backend with the 'backend' keyword!

adam_result=tq.minimize(objective=O1,lr=lr,
              method='adam',
              maxiter=80,
              initial_values=init,
              silent=True)<jupyter_output><empty_output><jupyter_text>The plots below show the trajectory of both the value of the objective and the values of the angles as a function of time.<jupyter_code>adam_result.history.plot('energies')
adam_result.history.plot('angles')
print('best energy: ',adam_result.energy)
print('optimal angles: ',adam_result.angles)<jupyter_output><empty_output><jupyter_text>**We see that, minus a few hiccups, all the angles converge to optimimum values.** **Let's repeat what we did above, but with a few of the other methods! Here's RMSprop:**<jupyter_code>init={'a':0.25,'b':0.25,'c':0.25,'d':0.25}
lr=0.01
rms_result=tq.minimize(objective=O1,lr=lr,
              method='rmsprop',
              maxiter=80,
              initial_values=init,
              silent=True)
print('RMSprop optimization results:')
rms_result.history.plot('energies')
rms_result.history.plot('angles')
print('best energy: ',rms_result.energy)
print('optimal angles: ',rms_result.angles)<jupyter_output>RMSprop optimization results:<jupyter_text>**... And here's Momentum:**<jupyter_code>init={'a':0.25,'b':0.25,'c':0.25,'d':0.25}
lr=0.1
mom_result=tq.minimize(objective=O1,lr=lr,
              method='momentum',
              maxiter=80,
              initial_values=init,
              silent=True)

print('momentum optimization results:')
mom_result.history.plot('energies')
mom_result.history.plot('angles')
print('best energy: ',mom_result.energy)
print('optimal angles: ',mom_result.angles)<jupyter_output>momentum optimization results:<jupyter_text>Note that when using the [RMSprop](https://towardsdatascience.com/understanding-rmsprop-faster-neural-network-learning-62e116fcf29a) method, we reduced the learning rate from 0.1 to 0.01. Different methods may be more or less sensitive to choices of initial learning rate. Try going back to the previous examples, and choosing different learning rates, or different initial parameters, to gain a feel for how sensitive different methods are. **Accelerated gradient descent to achieve numerical accuracy:**One issue with gradient descent is that they are slow when it comes to converging to machine precision. Indeed, as one approaches an optimum, the gradient goes to zero, which means that each successive step is smaller and smaller, and convergence slows down. This is not really a problem in machine learning, where convergence to many digits is not needed, but it is an issue in chemistry.The standard method used to converge quantum chemical calculations is an accelerated form of gradient descent alternatively called Direct Inversion of the Iterative Subspace (DIIS) or Pulay mixing (for its inventor Peter Pulay). Basically, DIIS extrapolates from the error on past function evaluation to where the error is zero using a subspace approach. It is closely related to other subpsace and Krylov methods, such as GMRES.DIIS works best once we are pretty close to our solution. In this case, we start from a significantly better initial guess, and try to achieve a tolerance of 1e-9 with just standard gradient descent,<jupyter_code>init={'a':0.18,'b':0.8,'c':-0.1,'d':-0.4}
lr=0.1
no_diis_result=tq.minimize(objective=O1,lr=lr,
              method='sgd',
              tol=1e-10,
              maxiter=80,
              initial_values=init,
              silent=False)<jupyter_output>Optimizer: <class 'tequila.optimizers.optimizer_gd.OptimizerGD'> 
backend         : qulacs
device          : None
samples         : None
save_history    : True
noise           : None

Objective       : 1 expectationvalues
Gradient        : 8 expectationvalues
gradient instr  : None
active variables : 4
iter.        <O>          Δ<O>      max(d<O>)   rms(d<O>)
  1       -1.61467625   -2.87e-01   1.603e-01   1.169e-01    
  2       -1.61755894   -2.88e-03   3.376e-02   2.077e-02    
  3       -1.61770417   -1.45e-04   2.436e-02   1.566e-02    
  4       -1.61779488   -9.07e-05   2.099e-02   1.332e-02    
  5       -1.61786059   -6.57e-05   1.788e-02   1.135e-02    
  6       -1.61790828   -4.77e-05   1.523e-02   9.665e-03    
  7       -1.61794287   -3.46e-05   1.297e-02   8.230e-03    
  8       -1.61796795   -2.51e-05   1.104e-02   7.008e-03    
  9       -1.61798613   -1.82e-05   9.403e-03   5.966e-03    
 10       -1.61799931   -1.32e-05   8.006e-03   5.079e-03    
 11       -1.61800[...]<jupyter_text>We run out of iterations. Looking at the printout (obtained with `silent=False`), we see that convergence is achieved after 46 steps. We note that specifically convergence slows down significantly after 5 steps, as the gradients (and thus step sizes) become smaller the closer we get to the minima. This is the problem DIIS solves.Let's try again with acceleration.<jupyter_code>diis_result=tq.minimize(objective=O1,lr=lr,
              method='sgd',
              diis={'tol':1e-1},
              tol=1e-10,
              maxiter=80,
              initial_values=init,
              silent=False)<jupyter_output>Optimizer: <class 'tequila.optimizers.optimizer_gd.OptimizerGD'> 
backend         : qulacs
device          : None
samples         : None
save_history    : True
noise           : None

Objective       : 1 expectationvalues
Gradient        : 8 expectationvalues
gradient instr  : None
active variables : 4
iter.        <O>          Δ<O>      max(d<O>)   rms(d<O>)
  1       -1.61467625   -2.87e-01   1.603e-01   1.169e-01    
  2       -1.61755894   -2.88e-03   3.376e-02   2.077e-02    
  3       -1.61770417   -1.45e-04   2.436e-02   1.566e-02    DIIS
  4       -1.61783934   -1.35e-04   1.908e-02   1.203e-02    DIIS
  5       -1.61803383   -1.94e-04   1.855e-03   9.289e-04    DIIS
  6       -1.61802839   +5.45e-06   1.098e-02   5.504e-03    DIIS
  7       -1.61803399   -5.60e-06   7.242e-05   3.631e-05    DIIS
  8       -1.61803399   -2.04e-10   2.974e-05   1.492e-05    DIIS
delta f smaller than tolerance 1e-10. Stopping optimization.<jupyter_text>Convergence is obtained in just 8 steps. We can see here that DIIS kicks in at step 3, that is, when max(gradient) is below 1e-1 (the `tol=` argument in the DIIS dictionary) and we have accumulated at least 3 vectors. Once DIIS starts, convergence is very rapid, with a four orders of magnitude drop in the error in only three steps.We can see this most clearly in a log plot of the error,<jupyter_code>print('SGD + DIIS optimization results:')
Emin = diis_result.history.energies[-1]
plt.semilogy(no_diis_result.history.energies[:-1] - Emin)
plt.semilogy(diis_result.history.energies[:-1] - Emin)
plt.ylabel('Error on energy')
plt.legend(['without DIIS', 'with DIIS'])<jupyter_output>SGD + DIIS optimization results:<jupyter_text>We clearly see here that when DIIS kicks in, convergence rapidly accelerates. In comparison, the convergence rate of standard gradient descent actually gets slower the closer we are to a solution. However, it should be noted that DIIS is not entirely stable and only really works close to an optima. Furthermore, the implementation here is somewhat experimental. Use at your own risk! **The GD optimizer, with the Quantum Natural Gradient:**The Quantum Natural Gradient, or QNG, is a novel method of calculating gradients for quantum systems, inspired by the natural gradient sometimes employed in classical machine learning. The usual gradient we employ is with respect to a euclidean manifold, but this is not the only geometry -- nor even, the optimal geometry -- of quantum space. The QNG is, in essence, a method of taking gradients with respect to (an approximation to) the Fubini-Study metric. For information on how (and why) the QNG is used, see [Stokes et.al](https://arxiv.org/abs/1909.02108). Using the qng in Tequila is as simple as passing in the keyword gradient='qng' to optimizers which support it, such as the GD optimizer.  We will use it to optimize O2, our 3 qubit `Objective`, and then compare the results to optimizing the same circuit with the regular gradient.<jupyter_code>### the keyword stop_count, below, stops optimization if no improvement occurs after 50 epochs.
### let's use a random initial starting point:
init={k:np.random.uniform(-2,2) for k in ['a','b','c','d']}

lr=0.01
qng_result = tq.minimize(objective=O2,
                     gradient='qng',
                     method='sgd', maxiter=200,lr=lr,
                     initial_values=init, silent=True)
qng_result.history.plot('energies')
qng_result.history.plot('angles')
print('best energy with qng: ',qng_result.energy)
print('optimal angles with qng: ',qng_result.angles)<jupyter_output><empty_output><jupyter_text>To gain appreciation for why one might use the QNG, let's optimize the same circuit with the same learning rate and the same method, but without QNG.<jupyter_code>lr=0.01
sgd_noqng_result = tq.minimize(objective=O2,
                     gradient=None,
                     method='sgd', maxiter=200,lr=lr,
                     initial_values=init, silent=True)
print('plotting what happens without QNG')
sgd_noqng_result.history.plot('energies')
sgd_noqng_result.history.plot('angles')
print('best energy without qng: ',sgd_noqng_result.energy)
print('optimal angles without qng: ',sgd_noqng_result.angles)<jupyter_output>plotting what happens without QNG<jupyter_text>Though the starting point was random you will most likely see that the QNG run achieved a greater degree of improvement -- it will not perform worse --, and that the trajectories followed by angles there were different from those followed by angles in the sgd-only optimization. Feel free to play around with other methods, learning rates, or circuits in the space below!<jupyter_code>### Use this space to optimize your own circuits!<jupyter_output><empty_output><jupyter_text>The SciPy Optimizer`SciPy` is one of the most popular optimization packages in `Python`. It offers a wide variety of optimization strategies. We will not cover them here; for a full exploration of all the `SciPy` methods, see [their docs](https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html). Here, we will exhibit a few of the more powerful options available. Most `SciPy` keywords like `method_options` can just be passed to `minimize` directly in the same way as when using `SciPy` directly.<jupyter_code>tq.show_available_optimizers(module="scipy")<jupyter_output>available methods for optimizer module scipy
method               | optimizer module
--------------------------
NELDER-MEAD          | scipy
COBYLA               | scipy
POWELL               | scipy
SLSQP                | scipy
L-BFGS-B             | scipy
BFGS                 | scipy
CG                   | scipy
TNC                  | scipy
TRUST-KRYLOV         | scipy
NEWTON-CG            | scipy
DOGLEG               | scipy
TRUST-NCG            | scipy
TRUST-EXACT          | scipy
TRUST-CONSTR         | scipy<jupyter_text>We will try three different optimizers: `COBYLA`, which is gradient-free, `L-BFGS-B`, which employs gradients, and `NEWTON-CG`, which employs the Hessian.<jupyter_code>print('As a reminder: we will optimize:')
tq.draw(O1, backend="qiskit")
init={'a':0.25,'b':0.25,'c':0.25,'d':0.25}

cobyla_result = tq.minimize(objective=O1,
                            method="cobyla", 
                            initial_values=init, 
                            tol=1.e-3, method_options={"gtol":1.e-3},
                            silent=True)

cobyla_result.history.plot('energies')
cobyla_result.history.plot('angles')
print('best energy with cobyla: ',cobyla_result.energy)
print('optimal angles with cobyla: ',cobyla_result.angles)
lb_result = tq.minimize(objective=O1,
                            method="l-bfgs-b", 
                            initial_values=init, 
                            tol=1.e-3, method_options={"gtol":1.e-3},
                            silent=True)

lb_result.history.plot('energies')
lb_result.history.plot('angles')
print('best energy with L-BFGS-B: ',lb_result.energy)
print('optimal angles with L-BFGS-B: ',lb_result.angles)
newton_result = tq.minimize(objective=O1,
                            method="newton-cg", 
                            initial_values=init, 
                            tol=1.e-3, method_options={"gtol":1.e-3},
                            silent=True)

newton_result.history.plot('energies')
newton_result.history.plot('angles')
print('best energy with NEWTON-CG: ',newton_result.energy)
print('optimal angles with NEWTON-CG: ',newton_result.angles)<jupyter_output><empty_output><jupyter_text>All three of the methods converged to the same minimum, but not necessarily to the same angles; the gradient and hessian based methods converged to approximately the same angles in similar time.  **Scipy Extras: numerical gradients and Hessians**Scipy allows for the use of numerical gradients. To use them, pass down keywords to the `gradient` argument, like `'2-point'`. When using the numerical gradients of `SciPy` it is often crucial to determine a feasible stepsize for the procedure. This can be done with the `method_options` entry `finite_diff_rel_step` (for `SciPy` version 1.5 or higher) or `eps` (for `SciPy` version < 1.5).Here is one example. **Please check your SciPy version!**<jupyter_code>lb_result = tq.minimize(objective=O1,
                            method="l-bfgs-b", 
                            initial_values=init,
                            gradient="2-point",
                            tol=1.e-3, method_options={"gtol":1.e-3, "finite_diff_rel_step":1.e-4}, # eps for scipy version < 1.5
                            silent=True)

lb_result.history.plot('energies')
lb_result.history.plot('angles')
print('best energy with L-BFGS-B: ',lb_result.energy)
print('optimal angles with L-BFGS-B: ',lb_result.angles)<jupyter_output><empty_output><jupyter_text>**Scipy Extras: the QNG in SciPy**Scipy is also configured to use the qng, just as the gd optimizer is. All one needs to do is set `gradient=qng`. Let's See how QNG interacts with the `BFGS` optimizer. We will use 02, our 3-qubit expectationvalue, that we used previously.<jupyter_code>init={k:np.random.uniform(-2,2) for k in ['a','b','c','d']}
lr=0.01
bfgs_qng_result = tq.minimize(objective=O2,
                     gradient='qng',
                     method='bfgs', maxiter=200,lr=lr,
                     initial_values=init, silent=True)
print('plotting what happens with QNG')
bfgs_qng_result.history.plot('energies')
bfgs_qng_result.history.plot('angles')
print('best energy with qng: ',bfgs_qng_result.energy)
print('optimal angles with qng: ',bfgs_qng_result.angles)
bfgs_noqng_result = tq.minimize(objective=O2,
                     gradient=None,
                     method='bfgs', maxiter=200,lr=lr,
                     initial_values=init, silent=True)
print('plotting what happens without QNG')
bfgs_noqng_result.history.plot('energies')
bfgs_noqng_result.history.plot('angles')
print('best energy without qng: ',bfgs_noqng_result.energy)
print('optimal angles without qng: ',bfgs_noqng_result.angles)<jupyter_output>plotting what happens without QNG<jupyter_text>Numerical and Customized GradientsBy default `tequila` compiles analytical gradients of the objectives using `jax`, internal recompilation and the parameter shift rule. The default is not setting the `gradient` keyword or setting it to `None`. The keyword can also be set to a dictionary (keys are the variables, values are the `tequila` objectives which are assumed to evaluate to the corresponding gradients of the objective).  For example `gradient=tq.grad(objective)` will results have the same results as `gradient=None` or simply not setting it.`tequila` offers its own way of compiling numerical gradients which can then be used troughout all gradient based optimizers. It can be activated by setting `gradient` to a dictionary holding the finite difference stencil as `method` as well as the `stepsize`.   Numerical gradients of that type come with the cost of 2*`len(variables)` and can lead to significantly cheaper gradients, especially if many expectation values are involved in the objective and/or if heavy recompilation of parametrized gates is necessary. Here is a small example using our `O2` objective, here the numerical 2-point procedure leads to 4 expectation values in the gradients (while anayltial gradients would lead to 8, set silent to False in the upper example or remove the gradient statement here).<jupyter_code>lr=0.01
num_result = tq.minimize(objective=O2,
                     gradient={"method":"2-point", "stepsize":1.e-4},
                     method='sgd', maxiter=200,lr=lr,
                     initial_values=0.1, silent=False)
num_result.history.plot('energies')<jupyter_output><empty_output><jupyter_text>`tequila` currently offers `2-point` as well as `2-point-forward` and `2-point-backward` stencils as `method`.The method can also be set to a python function performing the task. Here is an example which implements the same as `2-point`. The function can be replaced by any function with the same signature.<jupyter_code>import copy
def my_finite_difference_stencil(obj, var_vals, key, step, *args, **kwargs):
        """
        calculate objective gradient by symmetric shifts about a point.
        Parameters
        ----------
        obj: Objective:
            objective to call.
        var_vals:
            variables to feed to the objective.
        key:
            which variable to shift, i.e, which variable's gradient is being called.
        step:
            the size of the shift; a small float.
        args
        kwargs

        Returns
        -------
        float:
            the approximated gradient of obj w.r.t variable key at point var_vals[key] as a float.

        """
        left = copy.deepcopy(var_vals)
        left[key] += step / 2
        right = copy.deepcopy(var_vals)
        right[key] -= step / 2
        return 1.0 / step * (obj(left, *args, **kwargs) - obj(right, *args, **kwargs))

num_result = tq.minimize(objective=O2,
                     gradient={"method":my_finite_difference_stencil, "stepsize":1.e-4},
                     method='sgd', maxiter=200,lr=lr,
                     initial_values=0.1, silent=True)
num_result.history.plot('energies')<jupyter_output><empty_output><jupyter_text>The `gradient` keyword can also be replaced by a dictionary of `tequila` objectives which evaluate to gradients approximations of it.  Bayesian optimization[Bayesian optimization](https://arxiv.org/abs/1807.02811) is a method of global optimization, often used to tune hyperparameters in classical learning. It has also seen use in the optimization of [quantum circuits](https://arxiv.org/pdf/1812.08862.pdf). Tequila currently supports 2 different bayesian optimization algorithms: [Phoenics](https://github.com/aspuru-guzik-group/phoenics) and [GPyOpt](https://github.com/SheffieldML/GPyOpt), optimizers originally developed for optimizing expensive experimental procedures in chemistry. Click the links to get to the respective github pages, and download the optimizers before continuing this tutorial.  GPyOptGPyOpt can be used like any of our other optimizers. Like the GD and SciPy optimizers, it also takes a 'method' keyword. 3 methods are supported: `'lbfgs'`,`'DIRECT'`, and `'CMA'`. See the `GPyOpt` github for more info.<jupyter_code>print('As a reminder, we will optimize')
tq.draw(O1,backend='qiskit')
### let's use the lbfgs method.
init={'a':0.25,'b':0.25,'c':0.25,'d':0.25}
### note: no lr is passed here! there are fewer tunable keywords for this optimizer.
result=tq.minimize(objective=O1,
              method='lbfgs',
              maxiter=80,
              initial_values=init)

print('GPyOpt optimization results:')
result.history.plot('energies')
result.history.plot('angles')
print('best energy: ',result.energy)
print('optimal angles: ',result.angles)<jupyter_output>Optimizer: <class 'tequila.optimizers.optimizer_gpyopt.OptimizerGPyOpt'> 
backend         : qulacs
device          : None
samples         : None
save_history    : True
noise           : None

method          : lbfgs
Objective       : 1 expectationvalues
num acquisition: 1, time elapsed: 0.22s
num acquisition: 2, time elapsed: 0.47s
num acquisition: 3, time elapsed: 0.73s
num acquisition: 4, time elapsed: 1.46s
num acquisition: 5, time elapsed: 1.67s
num acquisition: 6, time elapsed: 1.96s
num acquisition: 7, time elapsed: 2.39s
num acquisition: 8, time elapsed: 2.78s
num acquisition: 9, time elapsed: 3.20s
num acquisition: 10, time elapsed: 3.68s
num acquisition: 11, time elapsed: 3.83s
num acquisition: 12, time elapsed: 4.06s
num acquisition: 13, time elapsed: 4.24s
num acquisition: 14, time elapsed: 4.41s
num acquisition: 15, time elapsed: 4.59s
num acquisition: 16, time elapsed: 5.07s
num acquisition: 17, time elapsed: 5.25s
num acquisition: 18, time elapsed: 5.50s
num acquisition: [...]<jupyter_text>**Don't worry, the plot's not broken!**Perhaps you are looking at the plots above in horror. But, do take note: bayesian optimization is a global, exploratory optimization method, designed to explore large portions of parameter space while still seeking out optimality. Look at the optimal energy again, and one sees that the best performance of this optimization method matched that of all the gradient descent methods. We will apply gpyopt, next, to the QNG example circuit above, and see how bayesian optimization compares to QNG and SGD.<jupyter_code>print('Hey, remember me?')
tq.draw(O2)
### the keyword stop_count, below, stops optimization if no improvement occurs after 50 epochs.
### let's use a random initial starting point:
init={k:np.random.uniform(-2,2) for k in ['a','b','c','d']}

gpy_result = tq.minimize(objective=O2,maxiter=25,method='lbfgs',
                     initial_values=init)

gpy_result.history.plot('energies')
print('best energy: ',gpy_result.energy)
print('optimal angles: ',gpy_result.angles)<jupyter_output>Hey, remember me?
Objective with 1 unique expectation values
variables = [a, b, c, d]
types     = not compiled

Expectation Value 0
Hamiltonian :  (+1.0000Y(0)X(1)Y(2),)
variables :  [a, b, c, d]
circuit:

0: ───Ry(0.5π)─────Rz(0.318309886183791*pi*f((a,))_0)───@────────────────────────────────────────────@───────
                                                        │                                            │
1: ───Ry(0.333π)───Rz(0.318309886183791*pi*f((b,))_1)───X───@───Ry(0.318309886183791*pi*f((c,))_2)───X───@───
                                                            │                                            │
2: ───Ry(0.25π)─────────────────────────────────────────────X───Rx(0.318309886183791*pi*f((d,))_3)───────X───
Optimizer: <class 'tequila.optimizers.optimizer_gpyopt.OptimizerGPyOpt'> 
backend         : qulacs
device          : None
samples         : None
save_history    : True
noise           : None

method          : lbfgs
Objective       : 1 expectationvalues
[...]<jupyter_text>**In a very, very small number of steps, GPyOpt is able to match the performance of SGD with the QNG.** **There's a few extras you can access if you are well-familiar with GPyOpt.** We return as part of `result` an attribute `result.gpyopt_instance`, which is an instance of the native `GPyOpt` `BayesianOptimization` object -- the one built and run during your optimization. It has some plotting features you can use.<jupyter_code>obj=gpy_result.gpyopt_instance
obj.plot_convergence()<jupyter_output><empty_output><jupyter_text>If your function has 1 or 2 parameters (but no more) you can also see a plot of its acquisition function! see [here](https://www.blopig.com/blog/2019/10/a-gentle-introduction-to-the-gpyopt-module/) for more info!You can also extract the acquisition function of the model itself, and play with it (it takes ones object, an np array, as input), using:`acq=result.gpyopt_instance.acquisition.acquisition_function`  PhoenicsFinally, we turn to [`Phoenics`](https://github.com/aspuru-guzik-group/phoenics). This algorithm, originally developed within the Aspuru-Guzik group, can be accessed in the usual fashion. It's performance on the two-qubit optimization circuit is shown below. Note that the number of datapoints exceeds the provided __maxiter__; __maxiter__ here controls the number of parameter __batches__ suggested by phoenics. phoenics suggests a number of parameter sets to try out, per batch, that scales with the number of parameters (in a nonlinear fashion), so you may want to set maxiter lower if you are only playing around.<jupyter_code>init={'a':0.25,'b':0.25,'c':0.25,'d':0.25}
print('With phoenics we will optimize:')
print(O1)
### to see what you can pass down to phoenics, see the tequila documentation for that module.
p_result=tq.minimize(objective=O1,
            method='Phoenics',
            maxiter=5,
            initial_values=init,
            silent=False)

print('Phoenics optimization results on 2 qubit circuit:')
p_result.history.plot('energies')
p_result.history.plot('angles')
print('best energy: ',p_result.energy)
print('optimal angles: ',p_result.angles)<jupyter_output>With phoenics we will optimize:
Objective with 1 unique expectation values
variables = [a, b, c, d]
types     = not compiled
Phoenics config:

{'general': {'auto_desc_gen': 'False', 'batches': 5, 'boosted': 'False', 'parallel': 'False'}, 'parameters': [{'name': a, 'periodic': 'True', 'type': 'continuous', 'size': 1, 'low': 0, 'high': 6.283185307179586}, {'name': b, 'periodic': 'True', 'type': 'continuous', 'size': 1, 'low': 0, 'high': 6.283185307179586}, {'name': c, 'periodic': 'True', 'type': 'continuous', 'size': 1, 'low': 0, 'high': 6.283185307179586}, {'name': d, 'periodic': 'True', 'type': 'continuous', 'size': 1, 'low': 0, 'high': 6.283185307179586}], 'objectives': [{'name': 'Energy', 'goal': 'minimize'}]}
phoenics has recieved
objective: 

Objective with 1 unique expectation values
variables = [a, b, c, d]
types     = not compiled
noise model : None
samples     : None
maxiter     : 5
variables   : [a, b, c, d]
passive var : {}
now lets begin
[TIME]:   0:00:00.001013   (overall)
[...]<jupyter_text>We also have returned to you the phoenics object. One interesting object we can extract from this is the acquisition function. You can obtain this indirectly, using resut.object.bayesian_network.kernel_contribution. This function takes a numpy array  ( a point in your parameter space) and returns 2 numbers, x and y; the acquisition function then has the value x*y. Note: this is often zero.<jupyter_code>kc=p_result.phoenics_instance.bayesian_network.kernel_contribution
random_point=np.random.uniform(0,2*np.pi,4)
f,s=kc(random_point)
random_ac=f*s
print('random point ', random_point, ' has acquisition function value ',random_ac)<jupyter_output>random point  [3.506517   3.41045123 1.98397105 5.59827826]  has acquisition function value  0.0
{'fraction_non_alphanumeric': 0.12124502352515382, 'fraction_numerical': 0.034215874606754085, 'mean_word_length': 3.5102963335007535, 'pattern_counts': {'":': 10, '<': 103, '<?xml version=': 0, '>': 102, 'https://': 10, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '26', 'dir': 'jupyter-structured-clean-dedup', 'id': '199475', 'license': 'MIT', 'n_tokens_mistral': 12185, 'n_tokens_neox': 11170, 'n_words': 3700, 'path': 'tutorials/Optimizer_Tutorial.ipynb', 'repo_name': 'georgios-ts/tequila', 'size': '35919'}