<jupyter_start><jupyter_text>Tutorial: Optimization Hello, and welcome to our tutorial on optimization. Here, we will explore three of Tequila's built in optimizers. Chiefly, we will cover the gradient descent (GD) optimizer, we will also discuss the Phoenics and GPyOpt bayesian optimizers that can be accessed through Tequila.  1: The GD optimizer.<jupyter_code>### start at the start: import statements!
import tequila as tq
import numpy as np
from tequila.optimizers.optimizer_gd import minimize as gd_min<jupyter_output>running build_ext
skipping '/Users/sumneralperin-lea/.pyxbld/temp.macosx-10.9-x86_64-3.6/pyrex/BayesianNetwork/kernel_prob_reshaping.cpp' Cython extension (up-to-date)
skipping 'BayesianNetwork.kernel_prob_reshaping' extension (up-to-date)
running build_ext
skipping '/Users/sumneralperin-lea/.pyxbld/temp.macosx-10.9-x86_64-3.6/pyrex/BayesianNetwork/kernel_evaluations.cpp' Cython extension (up-to-date)
skipping 'BayesianNetwork.kernel_evaluations' extension (up-to-date)<jupyter_text>We start by selecting an objective to optimize. We will begin with a fairly simple, 2-qubit expectationvalue. We will optimize our 2-qubit circuit with the simple, but non trivial hamiltonian $[Y(0)+Qm(0)]\otimes X(1)$, where $Qm=\frac{1}{2} (I + Z)$, the projector onto the 0 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

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


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

### we use the .draw function to pretty-print circuits via backend printers.
tq.draw(U,backend='qiskit')
print(O)<jupyter_output>┌───┐┌───────────────┐   ┌───────────────┐   
q_0: |0>┤ H ├┤ Ry(f((a,))_0) ├─■─┤ Rx(f((c,))_2) ├─■─
        ├───┤├───────────────┤ │ ├───────────────┤ │ 
q_1: |0>┤ H ├┤ Rz(f((b,))_1) ├─■─┤ Rx(f((d,))_3) ├─■─
        └───┘└───────────────┘   └───────────────┘   
 c_0: 0 ═════════════════════════════════════════════
                                                     
 c_1: 0 ═════════════════════════════════════════════
                                                     
Objective with 1 unique expectation values
variables = [a, b, c, d]
types     = not compiled<jupyter_text>We are ready to optimize, now! 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.  There are several available optimization methods available to the GD optimizer, including basic SGD, SGD with momentum, and more advanced optimization strategies like Adam or RMS-prop.<jupyter_code>print('the following methods are available for Gradient Descent optimization:\n')
print(tq.optimizers.optimizer_gd.OptimizerGD.available_methods())<jupyter_output>the following methods are available for Gradient Descent optimization:

['adam', 'adagrad', 'adamax', 'nadam', 'sgd', 'momentum', 'nesterov', 'rmsprop', 'rmsprop-nesterov']<jupyter_text>We will now optimize our chosen expectationvalue, chosing 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=gd_min(O,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. Excercise: is this truly the best performance possible, or are we stuck in a local minimum?  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=gd_min(O,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=gd_min(O,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.  1.1: 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 a more complicated circuit below, and then compare the results to optimizing the same circuit with the regular gradient.<jupyter_code>### this time, don't scale by pi

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

print('drawing a more complicated circuit. Hope you like it!')
tq.draw(U)
### 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=E,
                     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=E,
                     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 (and so I, your humble tutorial writer, do not know what your graphs look like), you will most likely see that the QNG run achieved a greater degree of improvement, and that the trajectories followed by angles there was different from that followed by angles in the straight-gd optimization. Feel free to play around with other methods, learning rates, or circuits in the space below!<jupyter_code>### have fun!<jupyter_output><empty_output><jupyter_text>2. 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.  2.1: GPyOpt GPyOpt 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>from tequila.optimizers.optimizer_gpyopt import minimize as gpy_min<jupyter_output><empty_output><jupyter_text>we will use GPyOpt to optimize the same circuits as seen above.<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

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


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

### we use the .draw function to pretty-print circuits via backend printers.
tq.draw(U,backend='qiskit')
print(O)
### 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=gpy_min(O,
              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
backend_options : {}
samples         : None
save_history    : True
noise           : None

method          : lbfgs
Objective       : 1 expectationvalues
num acquisition: 1, time elapsed: 0.51s
num acquisition: 2, time elapsed: 0.96s
num acquisition: 3, time elapsed: 1.67s
num acquisition: 4, time elapsed: 2.48s
num acquisition: 5, time elapsed: 3.10s
num acquisition: 6, time elapsed: 4.26s
num acquisition: 7, time elapsed: 5.01s
num acquisition: 8, time elapsed: 5.88s
num acquisition: 9, time elapsed: 6.88s
num acquisition: 10, time elapsed: 7.18s
num acquisition: 11, time elapsed: 7.69s
num acquisition: 12, time elapsed: 8.07s
num acquisition: 13, time elapsed: 8.59s
num acquisition: 14, time elapsed: 9.35s
num acquisition: 15, time elapsed: 10.14s
num acquisition: 16, time elapsed: 10.54s
num acquisition: 17, time elapsed: 11.08s
num acquisition: 18, time elapsed: 11.48s
num acquisition[...]<jupyter_text>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>### this time, don't scale by pi

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

print('Hey, remember me?')
tq.draw(U)
### 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 = gpy_min(objective=E,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?
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
backend_options : {}
samples         : None
save_history    : True
noise           : None

method          : lbfgs
Objective       : 1 expectationvalues
num acquisition: 1, time elapsed: 0.47s
num acquisition: 2, time elapsed: 0.96s
num acquisition: 3, time elapsed: 1.57s
num acquisition: 4, time elapsed: 2.11s
num acquisition: 5, time elaps[...]<jupyter_text>In a very, very small number of step, GPyOpt is able to match the performance of SGD with the QNG, and discovers the hidden truth: the optimil circuit, here, is one where all angles are zero (modulo 2 $\pi$)  There's a few extras you can access if you are well-familiar with GPyOpt. We return a part of result, result.object, which is a 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.object
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.object.acquisition.acquisition_function  Feel free to play around more with other circuits in the space below!  2.2 Phoenics Finally, we turn to Phoenics. This algorithm, originally developed within the Aspuru-Guzik group (Hey, just like Tequila!), 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>from tequila.optimizers.optimizer_phoenics import minimize as p_min
### 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

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

H=(tq.paulis.Y(0)+tq.paulis.Qm(0))*tq.paulis.X(1)
O=tq.ExpectationValue(U=U,H=H)

init={'a':0.25,'b':0.25,'c':0.25,'d':0.25}
### geez! even fewer keywords! 
### to see what you can pass down to phoenics, see the tequila documentation for that module.
p_result=p_min(O,
            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>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 : {}
backend options qulacs  {}
now lets begin
[TIME]:   0:00:00.001051   (overall)
***********************************************
energy = -0.29609650 , angles= {a: 6.2376809120178[...]<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.object.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  [2.02107126 0.53605093 1.77683872 4.47834312]  has acquisition function value  0.0
StarcoderdataJupyterStructuredCleanDedup