The purpose of this notebook is to describe a naive but pedagogical, first baby step, in the implementation of what is called in Qubiter the Mean Hamiltonian Minimization problem.
The qc history of this problem started with quantum chemists planning to use on a qc the phase estimation algo invented by Kitaev? (an algo that is also implemented in Qubiter) to estimate the energy levels ( eigenvalues) of simple molecules, initially H2. Then a bunch of people realized, heck, rather than trying to estimate the eigenvalues of a Hamiltonian by estimating the phase changes it causes, we can estimate those eigenvalues more efficiently by estimating the mean value of that Hamiltonian as measured empirically on a qc. Basically, just the Rayleigh-Ritz method, one of the oldest tricks in the book. One of the first papers to propose this mean idea is https://arxiv.org/abs/1304.3061 Their algo is commonly referred to by the ungainly name VQE (Variational Quantum Eigensolver) VQE was originally applied to do quantum chemistry with a qc. But now Rigetti and others have renamed it hybrid quantum-classical quantum computing and pointed out that it's an algo that has wide applicability, not just to quantum chemistry.
The idea behind hybrid quantum-classical is very simple. One has a classical box CBox and a quantum box QBox. The gates of QBox depend on N gate parameters. QBox sends info to CBox. CBox sends back to QBox N new gate parameters that will lower some cost function. This feedback process between CBox and QBox continues until the cost is minimized. The cost function is the mean value of a Hamiltonian which is estimated empirically from data obtained from the qc which resides inside the QBox.
To minimize a function of N continuous parameters, one can use some methods like simulated annealing and Powell that do not require calculating derivatives, or one can use methods that do use derivatives. Another possible separation is between methods that don't care which local minimum they find, as long as they find one of them, and those methods that try to find the best local minimum of them all, the so called global minimum. Yet another separation is between methods that allow constraints and those that don't.
Among the methods that do use derivatives, the so called gradient based methods only use the 1st derivative, whereas other methods use both first (Jacobian) and second (Hessian) derivatives. The performance of those that use both 1st and 2nd derivatives degrades quickly as N grows. Besides, calculating 2nd derivatives is very expensive. Hence, methods that use the 2nd derivatives are practically useless in the neural network field where N is usually very large. In that field, gradient based methods rule.
A method that uses no derivatives is Powell. A gradient based method that is designed to have a fast convergence rate is the Conjugate Gradient (CG) method. Another gradient based method is back-propagation (BP). BP can be implemented as distributed computing much more easily than other gradient based methods so it is favored by the most popular computer programs for doing distributed AI, such as PyTorch and Tensorflow.
Qubiter can perform minimization using various minlibs (minimization software libraries) such as 'scipy', 'autograd', 'pytorch', 'tflow'. It can also use various devices (aka simulators or backends), either virtual or real, to do the minimization.
Non-scipy minlibs implement backprop.
The 'scipy' minlib is a wrapper for the scipy function
scipy.optimize.minimize
. This scipy umbrella method implements many
minimization methods, including Powell and CG.
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
By a native device, we mean one that uses Qubiter native simulators like SEO_simulator.
So, without further ado, here is an example of the use of
class MeanHamilMinimizer
with a scipy minlib and native device.
$\newcommand{\bra}[1]{\left\langle{#1}\right|}$ $\newcommand{\ket}[1]{\left|{#1}\right\rangle}$ test: $\bra{\psi}M\ket{\phi}$
First change the directory to the Qubiter directory and add it to the path environment variable
import os
import sys
print(os.getcwd())
os.chdir('../../')
print(os.getcwd())
sys.path.insert(0,os.getcwd())
/home/rrtucci/PycharmProjects/qubiter/qubiter/jupyter_notebooks /home/rrtucci/PycharmProjects/qubiter
Next we construct a simple 4 qubit circuit that depends on two placeholder variables
#1
and #2
. These are the continuous variables that we will vary
to minimize a cost function. The cost function will be specified later.
from qubiter.adv_applications.MeanHamil_native import *
from qubiter.adv_applications.MeanHamilMinimizer import *
loaded OneQubitGate, WITHOUT autograd.numpy
num_qbits = 4
file_prefix = 'mean_hamil_native_test'
emb = CktEmbedder(num_qbits, num_qbits)
wr = SEO_writer(file_prefix, emb)
wr.write_Rx(2, rads=np.pi/7)
wr.write_Rx(1, rads='#2*.5')
wr.write_Rn(3, rads_list=['#1', '-#1*3', '#2'])
wr.write_Rx(1, rads='-my_fun#2#1')
wr.write_cnot(2, 3)
wr.close_files()
The code above wrote inside Qubiter's io_folder
, two files,
an English file and a Picture file. We next ask the writer object
wr to print those two files for us
wr.print_eng_file(jup=True)
1 | ROTX 25.714286 AT 2 | 2 | ROTX #2*.5 AT 1 | 3 | ROTN #1 -#1*3 #2 AT 3 | 4 | ROTX -my_fun#2#1 AT 1 | 5 | SIGX AT 3 IF 2T |
wr.print_pic_file(jup=True)
1 | | Rx | | | 2 | | | Rx | | 3 | R | | | | 4 | | | Rx | | 5 | X---@ | | |
The circuit depends on a placeholder function called my_fun
.
This function will remain fixed during the
minimization process but it must be defined
and it must be passed in inside an input dictionary to
the constructor of class MeanHamil
def my_fun(x, y):
return x + .5*y
fun_name_to_fun = {'my_fun': my_fun}
One must also pass into the constructor of MeanHamil
, the value of a Hamiltonian called hamil
.
Qubiter stores hamil
in an object of the class QubitOperator
from the open source Python library OpenFermion
.
Note that the QubitOperator
constructor simplifies hamil
automatically.
Hamiltonians are Hermitian, so after QubitOperator
finishes simplifying hamil
, the coefficient of every "term"
of hamil
must be real.
hamil = QubitOperator('X1 Y3 X1 Y1', .4) + QubitOperator('Y2 X1', .7)
print('hamil=\n', hamil)
hamil= 0.7 [X1 Y2] + 0.4 [Y1 Y3]
So what is the purpose of this Hamiltonian hamil
?
The cost function to be minimized is defined as the mean value of hamil
.
More precisely, if $\ket{\psi}$ is the ouput of the circuit
specified above, then the cost function equals $\bra{\psi} | H | \ket{\psi}$.
The initial values for parameters #1
and #2
to be minimized
must also be passed in inside an input dictionary to
the constructor of class MeanHamilMinimizer
. Variable all_var_nums
contains the keys
of that dictionary.
Note: Results are very heavily dependent on initial x_val, probably because,
due to the periodic nature of qubit rotations, there are many local minima
init_var_num_to_rads = {1: 1., 2: 3.}
all_var_nums = init_var_num_to_rads.keys()
For convenience, we define a function case()
that creates an
object of MeanHamilMinimizer
, asks that object to minimize the cost function,
and then case()
returns the final result of that minimization process.
As advertised at the beginning of this notebook, we
use various sub-methods of scipy.optimize.minimize
num_samples = 0
print_hiatus = 25
verbose = False
np.random.seed(1234)
emp_mhamil = MeanHamil_native(file_prefix, num_qbits, hamil,
all_var_nums, fun_name_to_fun, simulator_name='SEO_simulator', num_samples=num_samples)
targ_mhamil = MeanHamil_native(file_prefix, num_qbits, hamil,
all_var_nums, fun_name_to_fun, simulator_name='SEO_simulator') # zero samples
def case(**kwargs):
return MeanHamilMinimizer(emp_mhamil, targ_mhamil,
all_var_nums, init_var_num_to_rads,
print_hiatus, verbose).find_min(minlib='scipy', **kwargs)
We focus on two cases, num_samples=0 and num_samples > 0.
When MeanHamil
is told that
num_samples=0, it does no sampling. It just calculates
the mean value of hamil
"exactly', using the exact
final state vector calculated by Qubiter's SEO_simulator
class
min_method = 'Powell'
emp_mhamil.num_samples = 0
case(method=min_method)
x_val~ (#1, #2) iter=0, cost=-0.129256, x_val=1.000000, 3.000000 iter=25, cost=-0.128526, x_val=2.275041, 3.209272 iter=50, cost=-0.247852, x_val=1.105711, 3.540960 iter=75, cost=-0.241892, x_val=0.978839, 3.768145
direc: array([[ 0. , 1. ], [-0.11503823, 0.22534013]]) fun: -0.2479366921073845 message: 'Optimization terminated successfully.' nfev: 98 nit: 4 status: 0 success: True x: array([1.09458962, 3.55198269])
When num_samples > 0, Qubiter calculates the final state vector
of the circuit, then it samples that num_samples times,
obtains an empirical probability distribution from that,
and then it calculates the mean value of hamil
using that empirical distribution
min_method = 'Powell'
emp_mhamil.num_samples = 1000
case(method=min_method)
x_val~ (#1, #2) iter=0, cost=-0.099200, x_val=1.000000, 3.000000 iter=25, cost=-0.196800, x_val=1.248064, 3.472136 iter=50, cost=-0.219200, x_val=1.223266, 3.798374
direc: array([[1., 0.], [0., 1.]]) fun: -0.22179999999999994 message: 'Optimization terminated successfully.' nfev: 61 nit: 2 status: 0 success: True x: array([1.22326557, 3.77708762])
# very sensitive to eps, seems to be hitting discontinuity
min_method='CG'
num_sample= 0
case(options={'disp': True, 'eps':1e-2})
x_val~ (#1, #2) iter=0, cost=-0.139800, x_val=1.000000, 3.000000 iter=25, cost=-0.128000, x_val=0.999909, 2.999462 iter=50, cost=-0.158200, x_val=1.009999, 2.999992 Warning: Desired error not necessarily achieved due to precision loss. Current function value: -0.163000 Iterations: 1 Function evaluations: 75 Gradient evaluations: 16
fun: -0.16300000000000003 hess_inv: array([[1.06504473, 0.10933946], [0.10933946, 0.01122658]]) jac: array([ 1.36, -0.4 ]) message: 'Desired error not necessarily achieved due to precision loss.' nfev: 75 nit: 1 njev: 16 status: 2 success: False x: array([0.99999857, 2.99999152])