# Custom Expectation Value Program for the Qiskit Runtime¶

Paul Nation

IBM Quantum Partners Technical Enablement Team

Here we will show how to make a program that takes a circuit, or list of circuits, and computes the expectation values of one or more diagonal operators.

## Prerequisites¶

• You must have Qiskit 0.30+ installed.

• You must have an IBM Quantum Experience account with the ability to upload a Runtime program. Currently there is no way to know if you have Runtime upload ability outside of an email from IBM.

## Background¶

The primary method by which information is obtained from quantum computers is via expectation values. Indeed, the samples that come from executing a quantum circuit multiple times, once converted to probabilities, can be viewed as just a finite sample approximation to the expectation value for the projection operators corresponding to each bitstring. More practically, many quantum algorithms require computing expectation values over Pauli operators, e.g. Variational Quantum Eigensolvers, and thus having a runtime program that computes these quantities is of fundamental importance. Here we look at one such example, where an user passes one or more circuits and expectation operators and gets back the computed expectation values, and possibly error bounds.

### Expectation value of a diagonal operator¶

Consider a generic observable given by the tensor product of diagonal operators over $$N$$ qubits $$O = O_{N-1}\dots O_{0}$$ where the subscript indicates the qubit on which the operator acts. Then for a set of observed $$M$$ bitstrings $$\{b_{0}, \dots b_{M-1}\}$$, where $$M \leq 2^N$$, with corresponding approximate probabilites $$p_{m}$$ the expectation value is given by

$\langle O\rangle \simeq \sum_{m=0}^{M-1} p_{m}\prod_{n=0}^{N-1}O_{n}[b_{m}[N-n-1], b_{m}[N-n-1]],$

where $$O_{n}[b_{m}[N-n-1], b_{m}[N-n-1]]$$ is the diagonal element of $$O_{n}$$ specified by the $$N-n-1$$th bit in bitstring $$b_{m}$$. The reason for the complicated indexing in b_{m} is because Qiskit uses least-sginificant bit indexing where the zeroth element of the bit-strings is given by the right-most bit.

Here we will use built-in routines to compute these expectation values. However, it is not hard to do yourself, with plenty of examples to be found.

## Main program¶

Here we define our main function for the expectation value runtime program. As always, our program must start with the backend, and user_messenger arguements, followed by the actual inputs we pass to the program. Here our options are quite simple:

• circuits: A single QuantumCircuit or list of QuantumCircuits to be executed on the target backend.

• expectation_operators: The operators we want to evaluate. These can be strings of diagonal Pauli’s, eg, ZIZZ, or custom operators defined by dictionarys. For example, the projection operator on the all ones state of 4 qubits is {'1111': 1}.

• shots: Howe many times to sample each circuit.

• transpiler_config: A dictionary that passes additional arguments on to the transpile function, eg. optimization_level.

• run_config: A dictionary that passes additional arguments on to backend.run().

• skip_transpilation: A flag to skip transpilation altogether and just run the circuits. This is useful for situations where you need to transpile parameterized circuits once, but must bind parameters multiple times and evaluate.

• return_stddev: Flag to return bound on standard deviation. If using measurement mitigation this adds some overhead to the computation.

• use_measurement_mitigation: Use M3 measurement mitigation and compute expecation value and standard deviation bound from quasi-probabilities.

At the top of the cell below you will see a commented out %%writefile sample_expval.py. We will use this to convert the cell to a Python module named sample_expval.py to upload.

#%%writefile sample_expval.py
import mthree
from qiskit import transpile

# The entrypoint for our Runtime Program
def main(backend, user_messenger,
circuits,
expectation_operators='',
shots = 8192,
transpiler_config={},
run_config={},
skip_transpilation=False,
return_stddev=False,
use_measurement_mitigation=False,
):

"""Compute expectation
values for a list of operators after
executing a list of circuits on the target backend.

Parameters:
backend (ProgramBackend): Qiskit backend instance.
user_messenger (UserMessenger): Used to communicate with the program user.
circuits: (QuantumCircuit or list): A single list of QuantumCircuits.
expectation_operators (str or dict or list): Expectation values to evaluate.
shots (int): Number of shots to take per circuit.
transpiler_config (dict): A collection of kwargs passed to transpile().
run_config (dict): A collection of kwargs passed to backend.run().
skip_transpilation (bool): Skip transpiling of circuits, default=False.
return_stddev (bool): Return upper bound on standard devitation,
default=False.
use_measurement_mitigation (bool): Improve resulting using measurement
error mitigation, default=False.

Returns:
array_like: Returns array of expectation values or a list of (expval, stddev)
tuples if return_stddev=True.
"""

# transpiling the circuits using given transpile options
if not skip_transpilation:
trans_circuits = transpile(circuits, backend=backend,
**transpiler_config)
# Make sure everything is a list
if not isinstance(trans_circuits, list):
trans_circuits = [trans_circuits]
# If skipping set circuits -> trans_circuits
else:
if not isinstance(circuits, list):
trans_circuits = [circuits]
else:
trans_circuits = circuits

# If we are given a single circuit but requesting multiple expectation
# values, then set flag to make multiple pointers to same result.
duplicate_results = False
if isinstance(expectation_operators, list):
if len(expectation_operators) and len(trans_circuits) == 1:
duplicate_results = True

# If doing measurement mitigation we must build and calibrate a
# mitigator object.  Will also determine which qubits need to be
# calibrated.
if use_measurement_mitigation:
# Get an the measurement mappings at end of circuits
meas_maps = mthree.utils.final_measurement_mapping(trans_circuits)
# Get an M3 mitigator
mit = mthree.M3Mitigation(backend)
# Calibrate over the set of qubits measured in the transpiled circuits.
mit.cals_from_system(meas_maps)

# Compute raw results
result = backend.run(trans_circuits, shots=shots, **run_config).result()
raw_counts = result.get_counts()

# When using measurement mitigation we need to apply the correction and then
# compute the expectation values from the computed quasi-probabilities.
if use_measurement_mitigation:
quasi_dists = mit.apply_correction(raw_counts, meas_maps,

if duplicate_results:
quasi_dists = mthree.classes.QuasiCollection(
[quasi_dists]*len(expectation_operators))
# There are two different calls depending on what we want returned.
if return_stddev:
return quasi_dists.expval_and_stddev(expectation_operators)
return quasi_dists.expval(expectation_operators)

# If the program didn't return in the mitigation loop above it means
# we are processing the raw_counts data.  We do so here using the
# mthree utilities
if duplicate_results:
raw_counts = [raw_counts]*len(expectation_operators)
if return_stddev:
return mthree.utils.expval_and_stddev(raw_counts, expectation_operators)
return mthree.utils.expval(raw_counts, expectation_operators)


## Local testing¶

Here we test with a local “Fake” backend that mimics the noise properties of a real system and a 4-qubit GHZ state.

from qiskit import QuantumCircuit
from qiskit.test.mock import FakeSantiago
from qiskit.providers.ibmq.runtime import UserMessenger
msg = UserMessenger()
backend = FakeSantiago()

qc = QuantumCircuit(4)
qc.h(2)
qc.cx(2, 1)
qc.cx(1, 0)
qc.cx(2, 3)
qc.measure_all()

main(backend, msg,
qc,
expectation_operators=['ZZZZ', 'IIII', 'IZZZ'],
transpiler_config={'optimization_level':3, 'layout_method': 'sabre',
'routing_method': 'sabre'},
run_config={},
skip_transpilation=False,
return_stddev=False,
use_measurement_mitigation=True
)

array([0.974614  , 1.        , 0.02428596])


If we have done our job correctly, the above should print out two expectation values close to one and a final expectation value close to zero.

Important: At this time all program names must be unique amongst all users. Thus you need to change the metadata name value to soemthing other than sample-expval if you want to upload.

meta = {
"name": "sample-expval",
"description": "A sample expectation value program.",
"max_execution_time": 1000,
"version": "1.0",
}

meta["parameters"] = [
{"name": "circuits", "description": "A single or list of QuantumCircuits.",
"type": "QuantumCircuit or list", "required": True},
{"name": "expectation_operators", "description": "One or more expectation values to evaluate",
"type": "str or dict or list", "required": False},
{"name": "shots", "description": "Number of shots per circuit.", "type": "int", "required": False},
{"name": "transpiler_config", "description": "A collection of kwargs passed to transpile.",
"type": "dict", "required": False},
{"name": "run_config", "description": "A collection of kwargs passed to backend.run.",
"type": "dict", "required": False},
{"name": "return_stddev", "description": "Return upper-bound on standard deviation.",
"type": "bool", "required": False},
{"name": "skip_transpilation", "description": "Skip circuit transpilation.",
"type": "bool", "required": False},
{"name": "use_measurement_mitigation", "description": "Use measurement mitigation to improve results.",
"type": "bool", "required": False}
]

meta["return_values"]: [
{"name": "result", "description": "A list of expectation values and optionally standard deviations.", "type": "list"}
]


We are now in a position to upload the program. To do so we first uncomment and excute the line %%writefile sample_expval.py giving use the sample_expval.py file we need to upload.

from qiskit import IBMQ

provider = IBMQ.get_provider(group='deployed')

program_id = provider.runtime.upload_program(data='sample_expval.py', metadata=meta)
program_id

'sample-expval'


### Delete program if needed¶

#provider.runtime.delete_program(program_id)


## Making the program public¶

If you want to make your program(s) public for the whole world to see you can do so by calling

#provider.runtime.set_program_visibility(program_id=program_id, public=True)


## Wrapping the runtime program¶

As always, it is best to wrap the call to the runtime program with a function (or possibly a class) that makes input easier and does some validation.

def expectation_value_runner(backend,
circuits,
expectation_operators='',
shots = 8192,
transpiler_config={},
run_config={},
skip_transpilation=False,
return_stddev=False,
use_measurement_mitigation=False):

"""Compute expectation values for a list of operators after
executing a list of circuits on the target backend.

Parameters:
backend (Backend or str): Qiskit backend instance or name.
circuits: (QuantumCircuit or list): A single or list of QuantumCircuits.
expectation_operators (str or dict or list): Expectation values to evaluate.
shots (int): Number of shots to take per circuit.
transpiler_config (dict): A collection of kwargs passed to transpile().
run_config (dict): A collection of kwargs passed to backend.run().
return_stddev (bool): Return upper bound on standard devitation,
default=False.
skip_transpilation (bool): Skip transpiling of circuits, default=False.
use_measurement_mitigation (bool): Improve resulting using measurement
error mitigation, default=False.

Returns:
array_like: Returns array of expectation values or a list of (expval, stddev)
pairs if return_stddev=True.
"""
if not isinstance(backend, str):
backend = backend.name()
options = {'backend_name': backend}

if isinstance(circuits, list) and len(circuits) != 1:
if isinstance(expectation_operators, list):
if len(circuits) != 1 and  len(expectation_operators) == 1:
expectation_operators = expectation_operators*len(circuits)
elif len(circuits) != len(expectation_operators):
raise ValueError('Number of circuits must match number of expectation \
values if more than one of each')
inputs = {}
inputs['circuits'] = circuits
inputs['expectation_operators'] = expectation_operators
inputs['shots'] = shots
inputs['transpiler_config'] = transpiler_config
inputs['run_config'] = run_config
inputs['return_stddev'] = return_stddev
inputs['skip_transpilation'] = skip_transpilation
inputs['use_measurement_mitigation'] = use_measurement_mitigation

return provider.runtime.run('sample-expval', options=options, inputs=inputs)


### Trying it out¶

Because we made our program public anyone can try it out. Lets do so here with our previously made GHZ state and running on the simulator.

backend = provider.backend.ibmq_qasm_simulator

all_zeros_proj = {'0000': 1}
all_ones_proj = {'1111': 1}
job = expectation_value_runner(backend, qc, [all_zeros_proj, all_ones_proj, 'ZZZZ'])

job.result()

array([0.50012207, 0.49987793, 1.        ])


The first two projectors should be nearly $$0.50$$ as they tell use the probability of being in the all zeros and ones states, respectively, which should be 50/50 for our GHZ state. The final expectation value of ZZZZ should be one since this is a GHZ over an even number of qubits. It should go close to zero for an odd number.

qc2 = QuantumCircuit(3)
qc2.h(2)
qc2.cx(2, 1)
qc2.cx(1, 0)
qc2.measure_all()

all_zeros_proj = {'000': 1}
all_ones_proj = {'111': 1}
job2 = expectation_value_runner(backend, qc2, [all_zeros_proj, all_ones_proj, 'ZZZ'])

job2.result()

array([0.50415039, 0.49584961, 0.00830078])


## Quantum Volume as an expectation value¶

Here we formulate QV as an expectation value of a projector onto the heavy-output elements on a distribution. We can then use our expectation value routine to compute whether a given circuit has passed the QV metric.

QV is defined in terms of heavy-ouputs of a distribution. Heavy-outputs are those bit-strings that are those that have probabilities above the median value of the distribution. Below we define the projection operator onto the set of bit-strings that are heavy-outputs for a given distribution.

def heavy_projector(qv_probs):
"""Forms the projection operator onto the heavy-outputs of a given probability distribution.

Parameters:
qv_probs (dict): A dictionary of bitstrings and associated probabilities.

Returns:
dict: Projector onto the heavy-set.
"""
median_prob = np.median(list(qv_probs.values()))
heavy_strs = {}
for key, val in qv_probs.items():
if val > median_prob:
heavy_strs[key] = 1
return heavy_strs


Now we generate 10 QV circuits as our dataset.

import numpy as np
from qiskit.quantum_info import Statevector
from qiskit.circuit.library import QuantumVolume
# Generate QV circuits
N = 10
qv_circs = [QuantumVolume(5) for _ in range(N)]


Next, we have to determine the heavy-set of each circuit from the ideal answer, and then pass this along to our heavy-set projector function that we defined above.

ideal_probs = [Statevector.from_instruction(circ).probabilities_dict() for circ in qv_circs]
heavy_projectors = [heavy_projector(probs) for probs in ideal_probs]


QV circuits have no meaasurements on them so need to add them:

circs = [circ.measure_all(inplace=False) for circ in qv_circs]


With a list of circuits and projection operators we now need only to pass both sets to our above expection value runner targeting the desired backend. We will also set the best transpiler arguments to give us a sporting chance of getting some passing scores.

backend = provider.backend.ibmq_manila

job3 = expectation_value_runner(backend, circs, heavy_projectors,
transpiler_config={'optimization_level':3, 'layout_method': 'sabre',
'routing_method': 'sabre'})

qv_scores = job3.result()
qv_scores

array([0.52392578, 0.60400391, 0.57189941, 0.57897949, 0.7734375 ,
0.65844727, 0.56225586, 0.73706055, 0.69030762, 0.61193848])


A passing QV score is one where the value of the heavy-set projector is above $$2/3$$. So let us see who passed:

qv_scores > 2/3

array([False, False, False, False,  True, False, False,  True,  True,
False])

from qiskit.tools.jupyter import *