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,
                                           return_mitigation_overhead=return_stddev)
        
        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.

Program metadata

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.

Next we add the needed program data to a dictionary for uploading with our program.

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"}
  ]

Upload the program

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
IBMQ.load_account();
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 *
%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.