Custom VQE Program for Qiskit Runtime¶
Here we will demonstrate how to create, upload, and use a custom Program for Qiskit Runtime. As the utility of the Runtime execution engine lies in its ability to execute many quantum circuits with low latencies, this tutorial will show how to create your own Variational Quantum Eigensolver (VQE) program from scratch.
Prerequisites¶
You must have Qiskit 0.30+ installed.
You must have an IBM Quantum account with the ability to upload a runtime program. You have this ability if you belong to more than one provider.
Current limitations¶
The runtime execution engine currently has the following limitations that must be kept in mind:
The Docker images used by the runtime include only Qiskit and its dependencies, with few exceptions. One exception is the inclusion of the
mthreemeasurement mitigation package.For security reasons, the runtime cannot make internet calls outside of the environment.
Your runtime program name must not contain an underscore
_, otherwise it will cause an error when you try to execute it.
As Qiskit Runtime matures, these limitations will be removed.
Simple VQE¶
VQE is an hybrid quantum-classical optimization procedure that finds the lowest eigenstate and eigenenergy of a linear system defined by a given Hamiltonian of Pauli operators. For example, consider the following two-qubit Hamiltonian:
where \(A\) is numerical coefficient and the subscripts label the qubits on which the operators act. The zero index being farthest right is the ordering used in Qiskit. The Pauli operators tell us which measurement basis to to use when measuring each of the qubits.
We want to find the ground state (lowest energy state) of this Hamiltonian, and the associated eigenvector. To do this we must start at a given initial state and iteratively vary the parameters that define this state using a classical optimizer, such that the computed energies of subsequent steps are nominally lower than those previously. The parameterized state of the system is defined by an ansatz quantum circuit that should have non-zero support in the direction of the ground state. Because in general we do not know the solution, the choice of ansatz circuit can be highly problem-specific with a form dictated by additional information. For further information about variational algorithms, we point the reader to Nature Reviews Physics volume 3, 625 (2021).
Thus we need at least the following inputs to create our VQE quantum program:
A representation of the Hamiltonian that specifies the problem.
A choice of parameterized ansatz circuit, and the ability to pass configuration options, if any.
However, the following are also beneficial inputs that users might want to have:
Add the ability to pass an initial state.
Vary the number of shots that are taken.
Ability to select which classical optimizer is used, and set configuraton values, if any.
Ability to turn on and off measurement mitigation.
Specifying the form of the input values¶
All inputs to runtime programs must be serializable objects. That is to say, whatever you pass into a runtime program must be able to be converted to JSON format. It is thus beneficial to keep inputs limited to basic data types and structures unless you have experience with custom object serialization, or they are common Qiskit types such as QuantumCircuit etc that the built-in RuntimeEncoder can handle. Fortunately, the VQE program described above can be made out of simple Python components.
First, it is possible to represent any Hamiltonian using a list of values with each containing the numerical coefficeint for each term and the string representation for the Pauli operators. For the above example, the ground state energy with \(A=1\) is \(-3\) and we can write it as:
H = [(1, 'XX'), (1, 'YY'), (1, 'ZZ')]
Next we have to provide the ability to specify the parameterized Ansatz circuit. Here we will take advange of the fact that many ansatz circuits are pre-defined in the Qiskit Circuit Library. Examples can be found in the N-local circuits section.
We would like the user to be able to select between ansatz options such as: NLocal, TwoLocal, and EfficientSU2. We could have the user pass the whole ansatz circuit to the program; however, in order to reduce the size of the upload we will pass the ansatz by name. In the runtime program, we can take this name and get the class that it corresponds to from the library using, for example,
import qiskit.circuit.library.n_local as lib_local
ansatz = getattr(lib_local, 'EfficientSU2')
For the ansatz configuration, we will pass a simple dict of values.
Optionals¶
If we want to add the ability to pass an initial state, then we will need to add the ability to pass a 1D list/ NumPy array. Because the number of parameters depends on the ansatz and its configuration, the user would have to know what ansatz they are doing ahead of time.
Selecting a number of shots requires simply passing an integer value.
Here we will allow selecting a classical optimizer by name from those in SciPy, and a
dictof configuration parameters. Note that for execution on an actual system, the noise inherent in today’s quantum systems makes having a stochastic optimizer crucial to success. SciPy does not have such a choice, and the one built into Qiskit is wrapped in such a manner as to make it difficult to use elsewhere. As such, here we will use an SPSA optimizer written to match the style of those in SciPy. This function is given in Appendix A.
Finally, for measurement error mitigation we can simply pass a boolean (True/False) value.
Main program¶
We are now in a position to start building our main program. However, before doing so we point out that it makes the code cleaner to make a separate fuction that takes strings of Pauli operators that define our Hamiltonian and convert them to a list of circuits with single-qubit gates that change the measurement basis for each qubit, if needed. This function is given in Appendix B.
Required signature¶
Every runtime program is defined via the main function, and must have the following input signature:
main(backend, user_message, *args, **kwargs)
where backend is the backend that the program is to be executed on, and user_message is the class by which interim (and possibly final) results are communicated back to the user. After these two items, we add our program-specific arguments and keyword arguments.
The main VQE program¶
Here is the main program for our sample VQE. What each element of the function does is written in the comments before the element appears.
# Grab functions and modules from dependencies
import numpy as np
import scipy.optimize as opt
from scipy.optimize import OptimizeResult
import mthree
# Grab functions and modules from Qiskit needed
from qiskit import QuantumCircuit, transpile
import qiskit.circuit.library.n_local as lib_local
# The entrypoint for our Runtime Program
def main(backend, user_messenger,
hamiltonian,
ansatz='EfficientSU2',
ansatz_config={},
x0=None,
optimizer='SPSA',
optimizer_config={'maxiter': 100},
shots = 8192,
use_measurement_mitigation=False
):
"""
The main sample VQE program.
Parameters:
backend (ProgramBackend): Qiskit backend instance.
user_messenger (UserMessenger): Used to communicate with the
program user.
hamiltonian (list): Hamiltonian whose ground state we want to find.
ansatz (str): Optional, name of ansatz quantum circuit to use,
default='EfficientSU2'
ansatz_config (dict): Optional, configuration parameters for the
ansatz circuit.
x0 (array_like): Optional, initial vector of parameters.
optimizer (str): Optional, string specifying classical optimizer,
default='SPSA'.
optimizer_config (dict): Optional, configuration parameters for the
optimizer.
shots (int): Optional, number of shots to take per circuit.
use_measurement_mitigation (bool): Optional, use measurement mitigation,
default=False.
Returns:
OptimizeResult: The result in SciPy optimization format.
"""
# Split the Hamiltonian into two arrays, one for coefficients, the other for
# operator strings
coeffs = np.array([item[0] for item in hamiltonian], dtype=complex)
op_strings = [item[1] for item in hamiltonian]
# The number of qubits needed is given by the number of elements in the strings
# the defiune the Hamiltonian. Here we grab this data from the first element.
num_qubits = len(op_strings[0])
# We grab the requested ansatz circuit class from the Qiskit circuit library
# n_local module and configure it using the number of qubits and options
# passed in the ansatz_config.
ansatz_instance = getattr(lib_local, ansatz)
ansatz_circuit = ansatz_instance(num_qubits, **ansatz_config)
# Here we use our convenence function from Appendix B to get measurement circuits
# with the correct single-qubit rotation gates.
meas_circs = opstr_to_meas_circ(op_strings)
# When computing the expectation value for the energy, we need to know if we
# evaluate a Z measurement or and identity measurement. Here we take and X and Y
# operator in the strings and convert it to a Z since we added the rotations
# with the meas_circs.
meas_strings = [string.replace('X', 'Z').replace('Y', 'Z') for string in op_strings]
# Take the ansatz circuits, add the single-qubit measurement basis rotations from
# meas_circs, and finally append the measurements themselves.
full_circs = [ansatz_circuit.compose(mcirc).measure_all(inplace=False) for mcirc in meas_circs]
# Get the number of parameters in the ansatz circuit.
num_params = ansatz_circuit.num_parameters
# Use a given initial state, if any, or do random initial state.
if x0:
x0 = np.asarray(x0, dtype=float)
if x0.shape[0] != num_params:
raise ValueError('Number of params in x0 ({}) does not match number \
of ansatz parameters ({})'. format(x0.shape[0],
num_params))
else:
x0 = 2*np.pi*np.random.rand(num_params)
# Because we are in general targeting a real quantum system, our circuits must be transpiled
# to match the system topology and, hopefully, optimize them.
# Here we will set the transpiler to the most optimal settings where 'sabre' layout and
# routing are used, along with full O3 optimization.
# This works around a bug in Qiskit where Sabre routing fails for simulators (Issue #7098)
trans_dict = {}
if not backend.configuration().simulator:
trans_dict = {'layout_method': 'sabre', 'routing_method': 'sabre'}
trans_circs = transpile(full_circs, backend, optimization_level=3, **trans_dict)
# If using measurement mitigation we need to find out which physical qubits our transpiled
# circuits actually measure, construct a mitigation object targeting our backend, and
# finally calibrate our mitgation by running calibration circuits on the backend.
if use_measurement_mitigation:
maps = mthree.utils.final_measurement_mapping(trans_circs)
mit = mthree.M3Mitigation(backend)
mit.cals_from_system(maps)
# Here we define a callback function that will stream the optimizer parameter vector
# back to the user after each iteration. This uses the `user_messenger` object.
# Here we convert to a list so that the return is user readable locally, but
# this is not required.
def callback(xk):
user_messenger.publish(list(xk))
# This is the primary VQE function executed by the optimizer. This function takes the
# parameter vector as input and returns the energy evaluated using an ansatz circuit
# bound with those parameters.
def vqe_func(params):
# Attach (bind) parameters in params vector to the transpiled circuits.
bound_circs = [circ.bind_parameters(params) for circ in trans_circs]
# Submit the job and get the resultant counts back
counts = backend.run(bound_circs, shots=shots).result().get_counts()
# If using measurement mitigation apply the correction and
# compute expectation values from the resultant quasiprobabilities
# using the measurement strings.
if use_measurement_mitigation:
quasi_collection = mit.apply_correction(counts, maps)
expvals = quasi_collection.expval(meas_strings)
# If not doing any mitigation just compute expectation values
# from the raw counts using the measurement strings.
# Since Qiskit does not have such functionality we use the convenence
# function from the mthree mitigation module.
else:
expvals = mthree.utils.expval(counts, meas_strings)
# The energy is computed by simply taking the product of the coefficients
# and the computed expectation values and summing them. Here we also
# take just the real part as the coefficients can possibly be complex,
# but the energy (eigenvalue) of a Hamiltonian is always real.
energy = np.sum(coeffs*expvals).real
return energy
# Here is where we actually perform the computation. We begin by seeing what
# optimization routine the user has requested, eg. SPSA verses SciPy ones,
# and dispatch to the correct optimizer. The selected optimizer starts at
# x0 and calls 'vqe_func' everytime the optimizer needs to evaluate the cost
# function. The result is returned as a SciPy OptimizerResult object.
# Additionally, after every iteration, we use the 'callback' function to
# publish the interm results back to the user. This is important to do
# so that if the Program terminates unexpectedly, the user can start where they
# left off.
# Since SPSA is not in SciPy need if statement
if optimizer == 'SPSA':
res = fmin_spsa(vqe_func, x0, args=(), **optimizer_config,
callback=callback)
# All other SciPy optimizers here
else:
res = opt.minimize(vqe_func, x0, method=optimizer,
options=optimizer_config, callback=callback)
# Return result. OptimizeResult is a subclass of dict.
return res
Local testing¶
We can test whether our routine works by simply calling the main function with a backend instance, a UserMessenger, and sample arguments.
from qiskit.providers.ibmq.runtime import UserMessenger
msg = UserMessenger()
# Use the local Aer simulator
from qiskit import Aer
backend = Aer.get_backend('qasm_simulator')
# Execute the main routine for our simple two-qubit Hamiltonian H, and perform 5 iterations of the SPSA solver.
main(backend, msg, H, optimizer_config={'maxiter': 5})
[1.419780432710152, 2.3984284215892018, 1.1306533554149105, 1.8357672762510684, 5.414120644000338, 6.107301966755861, -0.013391355872252708, 5.615586607539193, 4.211781149943555, 1.792388243059789, 4.203949657158362, 0.1038271369149637, 2.4220098073658884, 4.617958787629208, 2.9969591661895865, 1.5490655190231735]
[2.1084925021737537, 3.0871404910528035, 0.4419412859513089, 2.52447934571467, 4.725408574536736, 5.418589897292259, -0.7021034253358543, 6.3042986770027944, 3.523069080479953, 1.1036761735961873, 3.5152375876947604, 0.7925392063785653, 3.11072187682949, 5.30667085709281, 3.685671235653188, 0.8603534495595718]
[1.7365578685005831, 3.459075124725974, 0.8138759196244794, 2.8964139793878405, 4.353473940863566, 5.046655263619089, -1.0740380590090248, 5.932364043329624, 3.1511344468067826, 1.475610807269358, 3.8871722213679307, 1.1644738400517358, 2.73878724315632, 4.934736223419639, 4.057605869326359, 1.2322880832327423]
[1.7839871181735734, 3.4116458750529834, 0.766446669951489, 2.84898472971485, 4.306044691190576, 5.094084513292079, -1.0266088093360346, 5.884934793656634, 3.198563696479773, 1.5230400569423481, 3.8397429716949403, 1.1170445903787456, 2.6913579934833294, 4.887306973746649, 4.105035118999349, 1.2797173329057325]
[1.122687940285629, 4.072945052940928, 1.4277458478394336, 2.1876855518269056, 3.6447455133026314, 5.755383691180024, -1.687907987223979, 6.546233971544579, 2.5372645185918286, 2.1843392348302926, 4.501042149582885, 1.7783437682666903, 3.352657171371274, 4.226007795858704, 4.766334296887294, 0.618418155017788]
fun: -1.72705078125
message: 'Optimization terminated successfully.'
nfev: 10
nit: 5
success: True
x: array([ 1.12268794, 4.07294505, 1.42774585, 2.18768555, 3.64474551,
5.75538369, -1.68790799, 6.54623397, 2.53726452, 2.18433923,
4.50104215, 1.77834377, 3.35265717, 4.2260078 , 4.7663343 ,
0.61841816])
Having executed the above, we see that there are 5 parameter arrays returned, one for each callback, along with the final optimization result. The parameter arrays are the interim results, and the UserMessenger object prints these values to the cell output. The output itself is the answer we obtained, expressed as a SciPy OptimizerResult object.
Program metadata¶
Program metadata is essentially the docstring for a runtime program. It describes overall program information such as the program name, description, version, and the max_execution_time the program is allowed to run, as well as details the inputs and the outputs the program expects. At a bare minimum the values described above are required
import uuid
meta = {
"name": "sample-vqe-{}".format(uuid.uuid4()),
"description": "A sample VQE program.",
"max_execution_time": 100000,
"version": "1.0",
}
It is important to set the max_execution_time high enough so that your program does not get terminated unexpectedly. Additionally, one should make sure that interim results are sent back to the user so that, if something does happen, the user can start where they left off.
It is, however, good form to detail the parameters and return types, as well as interim results. That being said, if making a runtime intended to be used by others, this information would also likely be mirrored in the signature of a function or class that the user would interact with directly; end users should not directly call runtime programs. We will see why below. Nevertheless, let us add to our metadata. First, the parameters section details the inputs the user is able to pass:
meta["parameters"] = [
{"name": "hamiltonian", "description": "Hamiltonian whose ground state we want to find.", "type": "list", "required": True},
{"name": "ansatz", "description": "Name of ansatz quantum circuit to use, default='EfficientSU2'", "type": "str", "required": False},
{"name": "ansatz_config", "description": "Configuration parameters for the ansatz circuit.", "type": "dict", "required": False},
{"name": "x0", "description": "Initial vector of parameters.", "type": "ndarray", "required": False},
{"name": "optimizer", "description": "Classical optimizer to use, default='SPSA'.", "type": "str", "required": False},
{"name": "optimizer_config", "description": "Configuration parameters for the optimizer.", "type": "dict", "required": False},
{"name": "shots", "description": "Number of shots to take per circuit.", "type": "int", "required": False},
{"name": "use_measurement_mitigation", "description": "Use measurement mitigation, default=False.", "type": "bool", "required": False}
]
Next, the return_values section tells about the return types:
meta['return_values'] = [
{"name": "result", "description": "Final result in SciPy optimizer format.", "type": "OptimizeResult"}
]
and finally let us specify what comes back when an interim result is returned:
meta["interim_results"] = [
{"name": "params", "description": "Parameter vector at current optimization step", "type": "ndarray"},
]
Uploading the program¶
We now have all the ingredients needed to upload our program. To do so we need to collect all of our code in one file, here called sample_vqe.py for uploading. This limitation will be removed in later versions of Qiskit Runtime. Alternatively, if the entire code is contained within a single Jupyter notebook cell, this can be done using the magic function
%%writefile my_program.py
To actually upload the program we need to get a Provider from our IBM Quantum account:
from qiskit import IBMQ
IBMQ.load_account()
<AccountProvider for IBMQ(hub='ibm-q', group='open', project='main')>
provider = IBMQ.get_provider(group='deployed')
Program upload¶
The call to program_upload takes the target Python file as data and the metadata as inputs. If you have already uploaded the program this will raise an error and you must delete it first to continue.
program_id = provider.runtime.upload_program(data='sample_vqe.py', metadata=meta)
program_id
'sample-vqe-1c65cfd4-9551-42fc-bfbe-099e7fd2574f'
Here the returned program_id is the same as the program name given in the metadata. You cannot have more than one program with the same name and program_id. The program_id is how you should reference your program.
Program information¶
We can query the program for information and see that our metadata is correctly being attached:
prog = provider.runtime.program(program_id)
print(prog)
sample-vqe-1c65cfd4-9551-42fc-bfbe-099e7fd2574f:
Name: sample-vqe-1c65cfd4-9551-42fc-bfbe-099e7fd2574f
Description: A sample VQE program.
Version: 1.0
Creation date: 2021-10-06T22:21:05.000000
Max execution time: 100000
Input parameters:
- hamiltonian:
Description: Hamiltonian whose ground state we want to find.
Type: list
Required: True
- ansatz:
Description: Name of ansatz quantum circuit to use, default='EfficientSU2'
Type: str
Required: False
- ansatz_config:
Description: Configuration parameters for the ansatz circuit.
Type: dict
Required: False
- x0:
Description: Initial vector of parameters.
Type: ndarray
Required: False
- optimizer:
Description: Classical optimizer to use, default='SPSA'.
Type: str
Required: False
- optimizer_config:
Description: Configuration parameters for the optimizer.
Type: dict
Required: False
- shots:
Description: Number of shots to take per circuit.
Type: int
Required: False
- use_measurement_mitigation:
Description: Use measurement mitigation, default=False.
Type: bool
Required: False
Interim results:
- params:
Description: Parameter vector at current optimization step
Type: ndarray
Returns:
- result:
Description: Final result in SciPy optimizer format.
Type: OptimizeResult
Deleting a program¶
If you make a mistake and need to delete and/or re-upload the program, you can run the following, passing the program_id:
#provider.runtime.delete_program(program_id)
Running the program¶
Specify parameters¶
To run the program we need to specify the options that are used in the runtime environment (not the program variables). At present, only the backend_name is required.
backend = provider.backend.ibmq_qasm_simulator
options = {'backend_name': backend.name()}
The inputs dictionary is used to pass arguments to the main function itself. For example:
inputs = {}
inputs['hamiltonian'] = H
inputs['optimizer_config']={'maxiter': 10}
Execute the program¶
We now can execute the program and grab the result.
job = provider.runtime.run(program_id, options=options, inputs=inputs)
job.result()
{'fun': -1.93994140625,
'x': array([-1.28984461, 3.73974929, 3.52327612, 1.74979783, 3.13519544,
2.43577395, 1.30425595, 0.04847941, 6.17766827, 1.92879213,
1.95707213, 2.8097762 , 1.95108352, 1.20067124, 7.01868106,
4.36507161]),
'nit': 10,
'nfev': 20,
'message': 'Optimization terminated successfully.',
'success': True}
A few things need to be pointed out. First, we did not get back any interim results, and second, the return object is a plain dictionary. This is because we did not listen for the return results, and we did not tell the job how to format the return result.
Listening for interim results¶
To listen for interm results we need to pass a callback function to provider.runtime.run that stores the results. The callback takes two arguments job_id and the returned data:
interm_results = []
def vqe_callback(job_id, data):
interm_results.append(data)
Executing again we get:
job2 = provider.runtime.run(program_id, options=options, inputs=inputs, callback=vqe_callback)
job2.result()
{'fun': -2.635986328125,
'x': array([ 1.39625003, 3.10967996, 2.46291361, -0.09150619, 1.89013366,
0.48872864, 5.60656903, 1.12770301, 4.04603538, 2.85551118,
0.45677689, 3.46054286, 4.10740117, 4.163728 , 1.53949656,
3.46634995]),
'nit': 10,
'nfev': 20,
'message': 'Optimization terminated successfully.',
'success': True}
print(interm_results)
[[1.1839280526666394, 2.391820224610454, 2.7491281736833244, 0.5771768054969294, 2.349087960882593, 0.20251406828095217, 5.3527505036344865, 1.80726551800796, 2.8686317344166947, 2.4545878612072003, -0.04047464122825306, 4.2780676963333795, 3.27599724292225, 3.5527489679560844, 2.1472927005219273, 3.1637626657075555], [1.1855194978035488, 2.3902287794735444, 2.750719618820234, 0.5755853603600198, 2.3506794060195024, 0.20092262314404263, 5.351159058497577, 1.8088569631448694, 2.870223179553604, 2.452996416070291, -0.04206608636516258, 4.27647625119647, 3.2775886880591596, 3.554340413092994, 2.148884145658837, 3.165354110844465], [1.0411904999135912, 2.534557777363502, 2.8950486167101914, 0.7199143582499773, 2.206350408129545, 0.05659362525408518, 5.206830060607619, 1.664527965254912, 3.0145521774435617, 2.5973254139602484, 0.10226291152479487, 4.420805249086427, 3.133259690169202, 3.6986694109829514, 2.004555147768879, 3.0210251129545074], [1.005580093753927, 2.5701681835231662, 2.9306590228698557, 0.7555247644096416, 2.241960814289209, 0.020983219094420913, 5.242440466767284, 1.7001383714145764, 3.050162583603226, 2.561715007800584, 0.13787331768445915, 4.456415655246091, 3.0976492840095378, 3.663059004823287, 2.0401655539285435, 3.0566355191141716], [1.07047876838977, 2.6350668581590093, 2.8657603482340126, 0.8204234390454845, 2.177062139653366, 0.08588189373026392, 5.307339141403126, 1.6352396967787333, 2.985263908967383, 2.496816333164741, 0.20277199232030216, 4.521314329881934, 3.162547958645381, 3.7279576794591303, 1.9752668792927004, 2.9917368444783285], [1.3994411335364108, 2.96402922330565, 3.1947227133806533, 0.4914610738988439, 2.5060245048000067, -0.2430804714163767, 5.636301506549767, 1.3062773316320926, 3.3142262741140236, 2.8257786983113817, -0.12619037282633846, 4.192351964735293, 3.4915103237920215, 3.3989953143124896, 2.304229244439341, 3.3206992096249692], [1.325020213130704, 3.0384501437113567, 3.1203017929749466, 0.5658819943045507, 2.5804454252057134, -0.16865955101066996, 5.710722426955474, 1.231856411226386, 3.3886471945197303, 2.751357777905675, -0.2006112932320452, 4.117931044329586, 3.417089403386315, 3.4734162347181963, 2.2298083240336344, 3.395120130030676], [1.031941029864989, 2.7453709604456416, 2.8272226097092314, 0.2728028110388356, 2.2873662419399983, 0.12441963225504513, 6.003801610221189, 1.524935594492101, 3.6817263777854454, 2.45827859463996, 0.09246789003366987, 3.8248518610638707, 3.71016858665203, 3.7664954179839114, 1.9367291407679192, 3.102040946764961], [1.4127118235825624, 3.126141754163215, 2.446451815991658, -0.10796798267873797, 1.9065954482224248, 0.5051904259726187, 5.623030816503616, 1.1441648007745275, 4.062497171503019, 2.8390493883575334, 0.47323868375124345, 3.444081067346297, 4.090939380369604, 4.147266211701485, 1.5559583470503457, 3.4828117404825343], [1.3962500340466297, 3.1096799646272824, 2.4629136055275906, -0.09150619314280523, 1.890133658686492, 0.4887286364366859, 5.606569026967683, 1.1277030112385948, 4.046035381967086, 2.855511177893466, 0.4567768942153107, 3.46054285688223, 4.107401169905537, 4.163728001237418, 1.539496557514413, 3.4663499509466016]]
Formatting the returned results¶
In order to format the return results into the desired format, we need to specify a decoder. This decoder must have a decode method that gets called to do the actual conversion. In our case OptimizeResult is a simple sub-class of dict so the formatting is simple.
from qiskit.providers.ibmq.runtime import ResultDecoder
from scipy.optimize import OptimizeResult
class VQEResultDecoder(ResultDecoder):
@classmethod
def decode(cls, data):
data = super().decode(data) # This is required to preformat the data returned.
return OptimizeResult(data)
We can then use this when returning the job result:
job3 = provider.runtime.run(program_id, options=options, inputs=inputs)
job3.result(decoder=VQEResultDecoder)
fun: -0.645751953125
message: 'Optimization terminated successfully.'
nfev: 20
nit: 10
success: True
x: array([ 5.72140052, 2.29687026, 4.13837683, 3.22216958, 4.76184762,
1.20943004, 5.74244574, 2.22665936, 4.34308411, 3.8390838 ,
-0.50949471, 2.15587397, 3.19045035, 5.82751179, 1.95972168,
3.75821819])
Simplifying program execution with wrapping functions¶
While runtime programs are powerful and flexible, they are not the most friendly things to interact with. Therefore, if your program is intended to be used by others, it is best to make wrapper functions and/or classes that simplify the user experience. Moreover, such wrappers allow for validation of user inputs on the client side, which can quickly find errors that would otherwise be raised later during the execution process - something that might have taken hours waiting in queue to get to.
Here we will make two helper routines. First, a job wrapper that allows us to attach and retrieve the interim results directly from the job object itself, as well as decodes for us so that the end user need not worry about formatting the results themselves.
class RuntimeJobWrapper():
"""A simple Job wrapper that attaches interm results directly to the job object itself
in the `interm_results attribute` via the `_callback` function.
"""
def __init__(self):
self._job = None
self._decoder = VQEResultDecoder
self.interm_results = []
def _callback(self, job_id, xk):
"""The callback function that attaches interm results:
Parameters:
job_id (str): The job ID.
xk (array_like): A list or NumPy array to attach.
"""
self.interm_results.append(xk)
def __getattr__(self, attr):
if attr == 'result':
return self.result
else:
if attr in dir(self._job):
return getattr(self._job, attr)
raise AttributeError("Class does not have {}.".format(attr))
def result(self):
"""Get the result of the job as a SciPy OptimizerResult object.
This blocks until job is done, cancelled, or errors.
Returns:
OptimizerResult: A SciPy optimizer result object.
"""
return self._job.result(decoder=self._decoder)
Next, we create the actual function we want users to call to execute our program. To this function we will add a series of simple validation checks (not all checks will be done for simplicity), as well as use the job wrapper defined above to simply the output.
import qiskit.circuit.library.n_local as lib_local
def vqe_runner(backend, hamiltonian,
ansatz='EfficientSU2', ansatz_config={},
x0=None, optimizer='SPSA',
optimizer_config={'maxiter': 100},
shots = 8192,
use_measurement_mitigation=False):
"""Routine that executes a given VQE problem via the sample-vqe program on the target backend.
Parameters:
backend (ProgramBackend): Qiskit backend instance.
hamiltonian (list): Hamiltonian whose ground state we want to find.
ansatz (str): Optional, name of ansatz quantum circuit to use, default='EfficientSU2'
ansatz_config (dict): Optional, configuration parameters for the ansatz circuit.
x0 (array_like): Optional, initial vector of parameters.
optimizer (str): Optional, string specifying classical optimizer, default='SPSA'.
optimizer_config (dict): Optional, configuration parameters for the optimizer.
shots (int): Optional, number of shots to take per circuit.
use_measurement_mitigation (bool): Optional, use measurement mitigation, default=False.
Returns:
OptimizeResult: The result in SciPy optimization format.
"""
options = {'backend_name': backend.name()}
inputs = {}
# Validate Hamiltonian is correct
num_qubits = len(H[0][1])
for idx, ham in enumerate(hamiltonian):
if len(ham[1]) != num_qubits:
raise ValueError('Number of qubits in Hamiltonian term {} does not match {}'.format(idx,
num_qubits))
inputs['hamiltonian'] = hamiltonian
# Validate ansatz is in the module
ansatz_circ = getattr(lib_local, ansatz, None)
if not ansatz_circ:
raise ValueError('Ansatz {} not in n_local circuit library.'.format(ansatz))
inputs['ansatz'] = ansatz
inputs['ansatz_config'] = ansatz_config
# If given x0, validate its length against num_params in ansatz:
if x0:
x0 = np.asarray(x0)
ansatz_circ = ansatz_circ(num_qubits, **ansatz_config)
num_params = ansatz_circ.num_parameters
if x0.shape[0] != num_params:
raise ValueError('Length of x0 {} does not match number of params in ansatz {}'.format(x0.shape[0],
num_params))
inputs['x0'] = x0
# Set the rest of the inputs
inputs['optimizer'] = optimizer
inputs['optimizer_config'] = optimizer_config
inputs['shots'] = shots
inputs['use_measurement_mitigation'] = use_measurement_mitigation
rt_job = RuntimeJobWrapper()
job = provider.runtime.run(program_id, options=options, inputs=inputs, callback=rt_job._callback)
rt_job._job = job
return rt_job
We can now execute our runtime program via this runner function:
job4 = vqe_runner(backend, H, optimizer_config={'maxiter': 15})
job4.result()
fun: -1.349853515625
message: 'Optimization terminated successfully.'
nfev: 30
nit: 15
success: True
x: array([0.09925502, 1.40473727, 1.61291267, 3.45519813, 2.65167136,
4.4163485 , 1.98523376, 5.94459488, 6.46103911, 1.76878845,
1.96124064, 3.31830748, 2.06192779, 4.28293342, 3.2448137 ,
1.63457609])
The interim results are now attached to the job interm_results attribute and, as expected, we see that the length matches the number of iterations performed.
len(job4.interm_results)
15
Conclusion¶
We have demonstrated how to create, upload, and use a custom Qiskit Runtime by creating our own VQE solver from scratch. This tutorial was meant to touch upon every aspect of the process for a real-world example. Within the current limitations of the runtime environment, this example should enable readers to develop their own single-file runtime program. This program is also a good starting point for exploring additional flavors of VQE runtime. For example, it is straightforward to vary the number of shots per iteration, increasing shots as the number of iterations increases. Those looking to go deeper can consider implimenting an adaptive VQE, where the ansatz is not fixed at initialization.
Appendix A¶
Here we code a simple simultaneous perturbation stochastic approximation (SPSA) optimizer for use on noisy quantum systems. Most optimizers do not handle fluctuating cost functions well, so this is a needed addition for executing on real quantum hardware.
import numpy as np
from scipy.optimize import OptimizeResult
def fmin_spsa(func, x0, args=(), maxiter=100,
a=1.0, alpha=0.602, c=1.0, gamma=0.101,
callback=None):
"""
Minimization of scalar function of one or more variables using simultaneous
perturbation stochastic approximation (SPSA).
Parameters:
func (callable): The objective function to be minimized.
``fun(x, *args) -> float``
where x is an 1-D array with shape (n,) and args is a
tuple of the fixed parameters needed to completely
specify the function.
x0 (ndarray): Initial guess. Array of real elements of size (n,),
where ‘n’ is the number of independent variables.
maxiter (int): Maximum number of iterations. The number of function
evaluations is twice as many. Optional.
a (float): SPSA gradient scaling parameter. Optional.
alpha (float): SPSA gradient scaling exponent. Optional.
c (float): SPSA step size scaling parameter. Optional.
gamma (float): SPSA step size scaling exponent. Optional.
callback (callable): Function that accepts the current parameter vector
as input.
Returns:
OptimizeResult: Solution in SciPy Optimization format.
Notes:
See the `SPSA homepage <https://www.jhuapl.edu/SPSA/>`_ for usage and
additional extentions to the basic version implimented here.
"""
A = 0.01 * maxiter
x0 = np.asarray(x0)
x = x0
for kk in range(maxiter):
ak = a*(kk+1.0+A)**-alpha
ck = c*(kk+1.0)**-gamma
# Bernoulli distribution for randoms
deltak = 2*np.random.randint(2, size=x.shape[0])-1
grad = (func(x + ck*deltak, *args) - func(x - ck*deltak, *args))/(2*ck*deltak)
x -= ak*grad
if callback is not None:
callback(x)
return OptimizeResult(fun=func(x, *args), x=x, nit=maxiter, nfev=2*maxiter,
message='Optimization terminated successfully.',
success=True)
Appendix B¶
This is a helper function that converts the Pauli operators in the strings that define the Hamiltonian operators into the appropriate measurements at the end of the circuits. For \(X\) operators this involves adding an \(H\) gate to the qubits to be measured, whereas a \(Y\) operator needs \(S^{+}\) followed by a \(H\). Other choices of Pauli operators require no additional gates prior to measurement.
def opstr_to_meas_circ(op_str):
"""Takes a list of operator strings and makes circuit with the correct post-rotations for measurements.
Parameters:
op_str (list): List of strings representing the operators needed for measurements.
Returns:
list: List of circuits for measurement post-rotations
"""
num_qubits = len(op_str[0])
circs = []
for op in op_str:
qc = QuantumCircuit(num_qubits)
for idx, item in enumerate(op):
if item == 'X':
qc.h(num_qubits-idx-1)
elif item == 'Y':
qc.sdg(num_qubits-idx-1)
qc.h(num_qubits-idx-1)
circs.append(qc)
return circs
from qiskit.tools.jupyter import *
%qiskit_copyright
This code is a part of Qiskit
© Copyright IBM 2017, 2021.
This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.