Compressive Gate Set Tomography (GST)#

This notebook gives an introduction on how to run compressive gate set tomography, what the input parameters mean and how to display different observations and plots.

This notebook and the compressive GST functionality requires the optional dependency “mgst” to be installed.
For background, see the journal reference at https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.4.010325

%load_ext autoreload
from iqm.benchmarks.compressive_gst.compressive_gst import GSTConfiguration, CompressiveGST

import json
import numpy as np
import matplotlib.pyplot as plt

Choose (or define) a backend#

#backend = "IQMFakeAdonis"
backend = "IQMFakeApollo"
#backend = "garnet"

Minimal GST Experiment configurations#

The most important parameters are the following:

qubits (List[int]): The qubits on the backend where the experiment is performed.
gate_set (Union[str, List[Type[QuantumCircuit]]]): Either one of the currently predefined gate sets "1QXYI", "2QXYCZ", "2QXYCZ_extended", "3QXYCZ", or a list of quantum circuits.
num_circuits (int): The number of circuits for the experiment. Recommended are at least 50 for single qubit GST, at least 300 for two-qubit GST and at least 2000 for 3-qubit GST.
shots (int): The number of shots per circuit.
rank (int): The Kraus rank of the reconstruction, i.e. the number of Kraus operators for each gate. $ $G (ρ) = \sum_{i = 1}^{rank} K_{i} ρ K_{i}^{†}$ $ Setting rank=1 will trigger a unitary gate fit, leading to a gate parametrization output in terms Hamiltonian parameters. The maximum rank is given by the physical dimension squared. For fully rigorous performance computation of metrics such as the average gate fidelity or the diamond distance it is recommended to choose full rank. A low rank has the benefit of needing less circuits and less computing time, while still capturing the dominant error sources.
bootstrap_samples (int): If bootstrapping error bars are to be generated, this variable sets the number of bootstrap runs. Recommended for trustworthy error bars are 50 samples. The default is 0, since computing bootstrap error bars is very time consuming.

Minimal_1Q_GST = GSTConfiguration(
    qubit_layouts=[[0]],
    gate_set="1QXYI",
    num_circuits=50,
    shots=1000,
    rank=4,
    bootstrap_samples=0,
)

Minimal_2Q_GST = GSTConfiguration(
    qubit_layouts=[[0,1]],
    gate_set="2QXYCZ",
    num_circuits=1000,
    shots=1000,
    rank=16,
    bootstrap_samples=0,
)

Execute GST Experiment(s)#

Be prepared that the first execution on a new system will take an extra 1-2 minutes to compile the lower level optimization code.

benchmark = CompressiveGST(backend, Minimal_1Q_GST)
benchmark.run()
result = benchmark.analyze()

2025-04-01 16:06:42,999 - iqm.benchmarks.logging_config - INFO - Now generating 50 random GST circuits...
2025-04-01 16:06:43,430 - iqm.benchmarks.logging_config - INFO - Will transpile all 50 circuits according to fixed physical layout
2025-04-01 16:06:43,430 - iqm.benchmarks.logging_config - INFO - Transpiling for backend IQMFakeApolloBackend with optimization level 0, sabre routing method all circuits
2025-04-01 16:06:44,259 - iqm.benchmarks.logging_config - INFO - Submitting batch with 50 circuits corresponding to qubits [0]
2025-04-01 16:06:44,338 - iqm.benchmarks.logging_config - INFO - Now executing the corresponding circuit batch
2025-04-01 16:06:44,434 - iqm.benchmarks.logging_config - INFO - Retrieving all counts
2025-04-01 16:06:46,644 - iqm.benchmarks.logging_config - INFO - Starting mGST optimization...
 26%|██▌       | 64/250 [00:01<00:05, 36.20it/s]
2025-04-01 16:06:48,462 - iqm.benchmarks.logging_config - INFO - Convergence criterion satisfied
2025-04-01 16:06:48,462 - iqm.benchmarks.logging_config - INFO - Final objective 1.56e-4 in time 1.82s

To get a unitary model of the gate set from the same data, the rank parameter of the benchmark object can be set to 1. The analysis will give a Hamiltonian parametrization of the gate set and produce different plots.

from copy import deepcopy
benchmark_rK1 = deepcopy(benchmark)
benchmark_rK1.runs[0].dataset.attrs["rank"] = 1
result_rK1 = benchmark_rK1.analyze()

2025-04-01 16:07:55,034 - iqm.benchmarks.logging_config - INFO - Starting mGST optimization...
 11%|█         | 27/250 [00:01<00:08, 24.83it/s]
2025-04-01 16:07:56,147 - iqm.benchmarks.logging_config - INFO - Convergence criterion satisfied
2025-04-01 16:07:56,148 - iqm.benchmarks.logging_config - INFO - Final objective 1.99e-4 in time 1.11s

Examine the results#

High level results stored at result.observations#

For uncertaintites on the observations set bootstrap_samples $\geq 10$ . The high level results are stored in a list of Observations under restult.observations. To access only those observations corresponding to a specific qubit layout, one can use the identifier attribute:

qubit_layout = [0]
for observation in result.observations:
    if observation.identifier.string_identifier == str(qubit_layout):
        print(f"{observation.name}: {observation.value} +/- {observation.uncertainty}")

Avg. gate fidelity - Idle:QB1: 0.99693 +/- nan
Diamond distance - Idle:QB1: 0.00703 +/- nan
Unitarity - Idle:QB1: 0.98781 +/- nan
Avg. gate fidelity - Rx(pi/2):QB1: 0.99918 +/- nan
Diamond distance - Rx(pi/2):QB1: 0.0013 +/- nan
Unitarity - Rx(pi/2):QB1: 0.99674 +/- nan
Avg. gate fidelity - Ry(pi/2):QB1: 0.99509 +/- nan
Diamond distance - Ry(pi/2):QB1: 0.01332 +/- nan
Unitarity - Ry(pi/2):QB1: 0.98074 +/- nan
Mean TVD: estimate - data: 0.00698 +/- nan
Mean TVD: target - data: 0.03296 +/- nan
POVM - diamond dist.: 0.00601 +/- nan
State - trace dist.: 0.04366 +/- nan

Accessing the final gate set estimates for further processing#

In addition to the high level observations above, the full process matrices for each gate and the full parametrizations for initial state and measurement are stored. They can be accessed under result.dataset.attrs[f"results_layout_{qubit_layout}"] as follows.

qubit_layout = [0]
print(result.dataset.attrs[f"results_layout_{qubit_layout}"].keys())

dict_keys(['raw_Kraus_operators', 'raw_gates', 'raw_POVM', 'raw_state', 'gauge_opt_gates', 'gauge_opt_gates_Pauli_basis', 'gauge_opt_POVM', 'gauge_opt_state', 'main_mGST_time', 'gauge_optimization_time', 'choi_evals', 'full_metrics'])

Two gate sets are saved, the raw gate set and the gauge-optimized gate set. In most instances the gauge-optimized gate set should be used for further processing, since it gives the gate set in the reference frame in which the target gates are defined.

raw_gates = result.dataset.attrs[f"results_layout_{qubit_layout}"]["raw_gates"]
gauge_opt_gates = result.dataset.attrs[f"results_layout_{qubit_layout}"]["gauge_opt_gates"]

The "raw_gates" and "gauge_opt_gates" keys in the outcome dictionary contain a 3D numpy array, where i.e. gate #1 is accessed with raw_results["gauge_opt_state"][0], and so on.

print(np.array_str(gauge_opt_gates[0], precision=3, suppress_small=True))

[[ 0.995+0.j    -0.002-0.002j -0.002+0.002j  0.002+0.j   ]
 [ 0.005-0.002j  0.994-0.003j  0.002-0.002j -0.005+0.004j]
 [ 0.005+0.002j  0.002+0.002j  0.994+0.003j -0.005-0.004j]
 [ 0.005+0.j     0.002+0.002j  0.002-0.002j  0.998+0.j   ]]

Plots#

The plots can be accessed in the notebook via result.plots, a dictionary containing all figure objects. These can then be diplayed as shown below, or saved to disc from the notebook. Currently all gauge optimized gate superoperators are shown as matrix plots in their Pauli basis representation, while the state preparation and measurement outcomes are shown as matrix plots in standard basis. For reference, a sinlge qubit superoperator for gate $G$ in the Pauli-basis has entries $G_{i j}$ defined via $ $G (ρ) = \frac{1}{2} \sum_{i, j = 1}^{4} G_{i j} P_{i} T r (P_{j} ρ),$ $w h e r e$ P_i, P_j$ are Pauli matrices.

In addition to gate plots, selections of gate error measures and gate parameters are also stored in figure objects.

Stored plot names can be displayed and individually plotted as follows.

print(list(result.plots.keys()))
result.plot('layout_[0]_gate_metrics')

['layout_[0]_choi_eigenvalues', 'layout_[0]_gate_metrics', 'layout_[0]_other_metrics', 'layout_[0]_process_matrix_0', 'layout_[0]_process_matrix_1', 'layout_[0]_process_matrix_2', 'layout_[0]_SPAM_matrices_real', 'layout_[0]_SPAM_matrices_imag']

../_images/85973b8a8da2e8539b4c570bd26a3cf94b41010b39267c1e8e3108bc9f05b932.png

Rank 1 results (for example the gate Hamiltonians in the Pauli basis)

print(list(result_rK1.plots.keys()))
result_rK1.plot('layout_[0]_hamiltonian_parameters')

['layout_[0]_hamiltonian_parameters', 'layout_[0]_gate_metrics', 'layout_[0]_other_metrics', 'layout_[0]_process_matrix_0', 'layout_[0]_process_matrix_1', 'layout_[0]_process_matrix_2', 'layout_[0]_SPAM_matrices_real', 'layout_[0]_SPAM_matrices_imag']

../_images/fb8fc196a1797b6e27812197746f860defc9cc60b0468dcac5a05751c2608120.png

Alternatively, all plots (of all layouts) can be displayed via result.plot_all().

Creating a custom gate set#

In this example we define a gate set ourselves as a list of quantum circuits. The example gate set chosen here is the “XYI” gate set with the addition of $π$ -rotations around the X- and Y-axis. Note that we also increased the number of random GST sequences in the GST configuration from 50 to 100 to account for the larger gate set.

from qiskit import QuantumCircuit
from qiskit.circuit.library import CZGate, RGate

# Define a list of gate instructions
gate_list = [RGate(1e-10, 0), RGate(0.5 * np.pi, 0), RGate(0.5 * np.pi, np.pi / 2), RGate(np.pi, 0), RGate(np.pi, np.pi / 2)]
# Define the gate set as a list of circuits, each with one gate instruction
gate_set = [QuantumCircuit(1, 0) for _ in range(len(gate_list))]
for i, gate in enumerate(gate_list):
    gate_set[i].append(gate, [0])
# Define the gate names (to be used for indentification and in plots/tables)
gate_labels = ["Idle", "Rx(pi/2)", "Ry(pi/2)", "Rx(pi)", "Ry(pi)"]

# Checking the gate set for correctness
for gate in gate_set:
    print(gate)

   ┌────────┐
q: ┤ R(0,0) ├
   └────────┘
   ┌──────────┐
q: ┤ R(π/2,0) ├
   └──────────┘
   ┌────────────┐
q: ┤ R(π/2,π/2) ├
   └────────────┘
   ┌────────┐
q: ┤ R(π,0) ├
   └────────┘
   ┌──────────┐
q: ┤ R(π,π/2) ├
   └──────────┘

CUSTOM_1Q_GST = GSTConfiguration(
    qubit_layouts=[[0]],
    gate_set=gate_set,
    gate_labels = gate_labels,
    num_circuits=100,
    shots=1000,
    rank=4,
    verbose=True,
    bootstrap_samples=0,
)

benchmark = CompressiveGST(backend, CUSTOM_1Q_GST)
benchmark.run()
result = benchmark.analyze()

2025-01-08 12:51:46,758 - iqm.benchmarks.logging_config - INFO - Now generating 100 random GST circuits...
2025-01-08 12:51:46,806 - iqm.benchmarks.logging_config - INFO - Will transpile all 100 circuits according to fixed physical layout
2025-01-08 12:51:46,806 - iqm.benchmarks.logging_config - INFO - Transpiling for backend IQMFakeApolloBackend with optimization level 0, sabre routing method all circuits
2025-01-08 12:51:47,075 - iqm.benchmarks.logging_config - INFO - Submitting batch with 100 circuits corresponding to qubits [0]
2025-01-08 12:51:47,075 - iqm.benchmarks.logging_config - INFO - Now executing the corresponding circuit batch
2025-01-08 12:51:47,391 - iqm.benchmarks.logging_config - INFO - Retrieving all counts
2025-01-08 12:51:49,408 - iqm.benchmarks.logging_config - INFO - Starting mGST optimization...
 12%|█████████▋                                                                       | 12/100 [00:19<02:22,  1.61s/it]
2025-01-08 12:52:08,800 - iqm.benchmarks.logging_config - INFO - Convergence criterion satisfied
2025-01-08 12:52:08,803 - iqm.benchmarks.logging_config - INFO - Final objective 2.42e-4 in time 19.39s

Checking out the process matrix corresponding to the pi rotation around the y-axis#

gauge_opt_gates = result.dataset.attrs[f"results_layout_{str([0])}"]["gauge_opt_gates"]
print(np.array_str(gauge_opt_gates[4], precision=3, suppress_small=True))

[[ 0.011+0.j    -0.001-0.002j -0.001+0.002j  0.998+0.j   ]
 [-0.003-0.008j -0.   -0.004j -0.989-0.001j  0.003-0.001j]
 [-0.003+0.008j -0.989+0.001j -0.   +0.004j  0.003+0.001j]
 [ 0.989+0.j     0.001+0.002j  0.001-0.002j  0.002+0.j   ]]