Parallel XEB

View on QuantumAI Run in Google Colab View source on GitHub Download notebook 
try:
    import cirq
except ImportError:
    print("installing cirq...")
    !pip install --quiet cirq~=1.0.dev
    print("installed cirq.")

This notebook demonstrates how to use the functionality in cirq.experiments to run parallel XEB end-to-end. "Parallel" means we characterize multiple pairs simultaneously.

import cirq
import numpy as np

%matplotlib inline
from matplotlib import pyplot as plt

Parallel XEB with library functions

The entire XEB workflow can be run by calling cirq.experiments.parallel_two_qubit_xeb and the combined single-qubit randomized benchmarking (RB) and XEB workflows can be run by calling cirq.experiments.run_rb_and_xeb.

# Simulation
qubits = cirq.GridQubit.rect(3, 2, 4, 3)
result = cirq.experiments.parallel_two_qubit_xeb(
    sampler=cirq.DensityMatrixSimulator(noise=cirq.depolarize(5e-3), dtype=np.complex128),  # Any simulator or a ProcessorSampler.
    qubits=qubits
)

100%|██████████| 243/243 [02:59<00:00,  1.36it/s]

# The returned result is an instance of the `TwoQubitXEBResult` class which provides visualization methods like 
result.plot_heatmap();  # plot the heatmap of XEB errors
result.plot_fitted_exponential(*qubits[:2]);  # plot the fitted model of xeb error of a qubit pair.
result.plot_histogram();  # plot a histogram of all xeb errors.

png

png

png

# `TwoQubitXEBResult` also has methods to retrieve errors.
print('pauli errors:', result.pauli_error())
print('xeb errors:', result.xeb_error(*qubits[:2]))
print('xeb fidelity:', result.xeb_fidelity(*qubits[:2]))

pauli errors: {(cirq.GridQubit(4, 3), cirq.GridQubit(4, 4)): 0.020117979533971764, (cirq.GridQubit(4, 3), cirq.GridQubit(5, 3)): 0.019975326456336887, (cirq.GridQubit(4, 4), cirq.GridQubit(5, 4)): 0.019313445323787537, (cirq.GridQubit(5, 3), cirq.GridQubit(5, 4)): 0.01938205020704624, (cirq.GridQubit(5, 3), cirq.GridQubit(6, 3)): 0.019807522577596597, (cirq.GridQubit(5, 4), cirq.GridQubit(6, 4)): 0.02039249169943736, (cirq.GridQubit(6, 3), cirq.GridQubit(6, 4)): 0.01975950793158498}
xeb errors: 0.01609438362717741
xeb fidelity: 0.9839056163728226

The run_rb_and_xeb method returns an object of type InferredXEBResult which is like TwoQubitXEBResult except that it removes the single-qubit errors obtained from the single-qubit randomized benchmarking (RB) experiment to isolate the error from the two qubit gate.

Step by step XEB

The rest of this notebook explains how the parallel_two_qubit_xeb works internally. Note that the notebook uses SQRT_ISWAP as the entangling gate while parallel_two_qubit_xeb and run_rb_and_xeb default to CZ.

Set up Random Circuits

We create a library of 10 random, two-qubit circuits using the sqrt(ISWAP) gate. These library circuits will be mixed-and-matched among all the pairs on the device we aim to characterize.

from cirq.experiments import random_quantum_circuit_generation as rqcg

circuit_library = rqcg.generate_library_of_2q_circuits(
    n_library_circuits=20, 
    two_qubit_gate=cirq.ISWAP**0.5,
    random_state=52,
)
print(len(circuit_library))

20

# We will truncate to these lengths
max_depth = 100
cycle_depths = np.arange(3, max_depth, 20)
cycle_depths

array([ 3, 23, 43, 63, 83])

Determine the device topology

We will run on all pairs from a given device topology. Below, you can supply a device_name if you're authenticated to run on Google QCS. In that case, we will get the device object from the cloud endpoint and turn it into a graph of qubits. Otherwise, we mock a device graph by allocating arbitrary cirq.GridQubits to turn into a graph.

device_name = None  # change me!

import cirq.contrib.routing as ccr
import networkx as nx

if device_name is None:
    qubits = cirq.GridQubit.rect(3, 2, 4, 3)
    # Delete one qubit from the rectangular arangement to
    # 1) make it irregular 2) simplify simulation.
    qubits = qubits[:-1]
    sampler = cirq.DensityMatrixSimulator(noise=cirq.depolarize(5e-3))
    graph = ccr.gridqubits_to_graph_device(qubits)
else:
    import cirq_google as cg
    sampler = cg.get_engine_sampler(device_name, gate_set_name='sqrt_iswap')
    device = cg.get_engine_device(device_name)
    qubits = sorted(device.qubits)
    graph = ccr.gridqubits_to_graph_device(device.qubits)


pos = {q: (q.row, q.col) for q in qubits}
nx.draw_networkx(graph, pos=pos)

png

Set up our combinations

We take the library of two-qubit circuits in circuit_library and mix-and-match to sampled in parallel.

We will pass combs_by_layer and circuit_library to the sampling function which will "zip" the circuits according to these combinations. The outer list corresponds to the four cirq.GridInteractionLayers (one of four for the degree-four GridQubit-implied graph). The inner combinations matrix is a (n_combinations, n_pairs) ndarray of integers which index into the circuit library.

combs_by_layer = rqcg.get_random_combinations_for_device(
    n_library_circuits=len(circuit_library),
    n_combinations=10,
    device_graph=graph,
    random_state=53,
)
combs_by_layer

[CircuitLibraryCombination(layer=cirq.experiments.GridInteractionLayer(col_offset=0, vertical=True, stagger=True), combinations=array([[ 5, 16],
        [12,  9],
        [ 5, 18],
        [11,  3],
        [ 6,  9],
        [13,  3],
        [11,  6],
        [14, 12],
        [16, 10],
        [18, 15]]), pairs=[(cirq.GridQubit(4, 4), cirq.GridQubit(5, 4)), (cirq.GridQubit(5, 3), cirq.GridQubit(6, 3))]),
 CircuitLibraryCombination(layer=cirq.experiments.GridInteractionLayer(col_offset=1, vertical=True, stagger=True), combinations=array([[16],
        [ 3],
        [12],
        [ 0],
        [ 5],
        [ 5],
        [ 0],
        [ 7],
        [ 9],
        [12]]), pairs=[(cirq.GridQubit(4, 3), cirq.GridQubit(5, 3))]),
 CircuitLibraryCombination(layer=cirq.experiments.GridInteractionLayer(col_offset=1, vertical=False, stagger=True), combinations=array([[13],
        [ 8],
        [ 8],
        [13],
        [ 1],
        [11],
        [11],
        [ 8],
        [14],
        [14]]), pairs=[(cirq.GridQubit(4, 3), cirq.GridQubit(4, 4))]),
 CircuitLibraryCombination(layer=cirq.experiments.GridInteractionLayer(col_offset=0, vertical=False, stagger=True), combinations=array([[18],
        [16],
        [ 2],
        [ 6],
        [ 1],
        [ 9],
        [11],
        [ 1],
        [ 9],
        [ 6]]), pairs=[(cirq.GridQubit(5, 3), cirq.GridQubit(5, 4))])]

Visualize

Here, we draw the four layers' active pairs.

fig, axes = plt.subplots(2,2, figsize=(9,6))
for comb_layer, ax in zip(combs_by_layer, axes.reshape(-1)):
    active_qubits = np.array(comb_layer.pairs).reshape(-1)
    colors = ['red' if q in active_qubits else 'blue' for q in graph.nodes]
    nx.draw_networkx(graph, pos=pos, node_color=colors, ax=ax)
    nx.draw_networkx_edges(graph, pos=pos, edgelist=comb_layer.pairs, width=3, edge_color='red', ax=ax)

plt.tight_layout()

png

Take Data

The following call will execute the zipped circuits and sample bitstrings.

from cirq.experiments.xeb_sampling import sample_2q_xeb_circuits
sampled_df = sample_2q_xeb_circuits(
    sampler=sampler,
    circuits=circuit_library,
    cycle_depths=cycle_depths,
    combinations_by_layer=combs_by_layer,
    shuffle=np.random.RandomState(52),
    repetitions=10_000,
)
sampled_df

100%|██████████| 207/207 [00:49<00:00,  4.22it/s]

Benchmark Fidelities

from cirq.experiments.xeb_fitting import benchmark_2q_xeb_fidelities
fids = benchmark_2q_xeb_fidelities(
    sampled_df=sampled_df,
    circuits=circuit_library,
    cycle_depths=cycle_depths,
)
fids

from cirq.experiments.xeb_fitting import fit_exponential_decays, exponential_decay
fidelities = fit_exponential_decays(fids)

heatmap_data = {}

for (_, _, pair), fidelity in fidelities.layer_fid.items():
    heatmap_data[pair] = 1.0 - fidelity

cirq.TwoQubitInteractionHeatmap(heatmap_data).plot();

png

for i, record in fidelities.iterrows():
    plt.axhline(1, color='grey', ls='--')
    plt.plot(record['cycle_depths'], record['fidelities'], 'o')
    xx = np.linspace(0, np.max(record['cycle_depths']))
    plt.plot(xx, exponential_decay(xx, a=record['a'], layer_fid=record['layer_fid']))
    plt.show()

png

png

png

png

png

import seaborn as sns

# Give each pair its own color
colors = sns.cubehelix_palette(n_colors=graph.number_of_edges())
colors = dict(zip(graph.edges, colors))

# Exponential reference
xx = np.linspace(0, fids['cycle_depth'].max())
plt.plot(xx, (1-5e-3)**(4*xx), label=r'Exponential Reference', color='black')

# Plot each pair
def _p(fids):
    q0, q1 = fids.name
    plt.plot(fids['cycle_depth'], fids['fidelity'], 
             'o-', label=f'{q0}-{q1}', color=colors[fids.name], alpha=0.5)
fids.groupby('pair').apply(_p)

plt.ylabel('Circuit fidelity')
plt.xlabel('Cycle Depth $d$')
plt.legend(loc='best')
plt.tight_layout()

/tmpfs/tmp/ipykernel_49834/2578408089.py:16: DeprecationWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.
  fids.groupby('pair').apply(_p)

png

Optimize PhasedFSimGate parameters

We know what circuits we requested, and in this simulated example, we know what coherent error has happened. But in a real experiment, there is likely unknown coherent error that you would like to characterize. Therefore, we make the five angles in PhasedFSimGate free parameters and use a classical optimizer to find which set of parameters best describes the data we collected from the noisy simulator (or device, if this was a real experiment).

import multiprocessing
pool = multiprocessing.get_context('spawn').Pool()

from cirq.experiments.xeb_fitting import (
    parameterize_circuit, 
    characterize_phased_fsim_parameters_with_xeb_by_pair, 
    SqrtISwapXEBOptions,
)

# Set which angles we want to characterize (all)
options = SqrtISwapXEBOptions(
    characterize_theta = True,
    characterize_zeta = True,
    characterize_chi = True,
    characterize_gamma = True,
    characterize_phi = True
)
# Parameterize the sqrt(iswap)s in our circuit library
pcircuits = [parameterize_circuit(circuit, options) for circuit in circuit_library]

# Run the characterization loop
characterization_result = characterize_phased_fsim_parameters_with_xeb_by_pair(
    sampled_df,
    pcircuits,
    cycle_depths,
    options,
    pool=pool,
    # ease tolerance so it converges faster:
    fatol=1e-2, 
    xatol=1e-2
)

characterization_result.final_params

{(cirq.GridQubit(5, 3), cirq.GridQubit(5, 4)): {'theta': -0.8321148131969989,
  'zeta': 0.10511205256930206,
  'chi': -0.0772924518567982,
  'gamma': -0.0077367253693801895,
  'phi': -0.002721022112115848},
 (cirq.GridQubit(4, 3), cirq.GridQubit(4, 4)): {'theta': -0.785265785514893,
  'zeta': 0.025935985318611354,
  'chi': -0.0304168347183542,
  'gamma': 0.004015471594176709,
  'phi': 0.0681226613898939},
 (cirq.GridQubit(4, 3), cirq.GridQubit(5, 3)): {'theta': -0.7549372852935057,
  'zeta': 0.029798183701181105,
  'chi': -0.03689390288774589,
  'gamma': 0.029746257267790458,
  'phi': -0.0642942609274198},
 (cirq.GridQubit(4, 4), cirq.GridQubit(5, 4)): {'theta': -0.779049485025755,
  'zeta': 0.032414929019303595,
  'chi': -0.019544238644370535,
  'gamma': 0.009969679046614537,
  'phi': 0.032517943016781486},
 (cirq.GridQubit(5, 3), cirq.GridQubit(6, 3)): {'theta': -0.7599534579368057,
  'zeta': 0.0016406747093006828,
  'chi': -0.045759826851075325,
  'gamma': -0.0026306173794372906,
  'phi': 0.0432584725364759} }

characterization_result.fidelities_df

from cirq.experiments.xeb_fitting import before_and_after_characterization
before_after_df = before_and_after_characterization(fids, characterization_result)
before_after_df

for i, row in before_after_df.iterrows():
    plt.axhline(1, color='grey', ls='--')
    plt.plot(row['cycle_depths_0'], row['fidelities_0'], '*', color='red')
    plt.plot(row['cycle_depths_c'], row['fidelities_c'], 'o', color='blue')

    xx = np.linspace(0, np.max(row['cycle_depths_0']))
    plt.plot(xx, exponential_decay(xx, a=row['a_0'], layer_fid=row['layer_fid_0']), 
             color='red', label=f'f_0 = {row["layer_fid_0"]:.3f}')
    plt.plot(xx, exponential_decay(xx, a=row['a_c'], layer_fid=row['layer_fid_c']),
             color='blue', label=f'f_c = {row["layer_fid_c"]:.3f}')

    plt.xlabel('Cycle Depth')
    plt.ylabel('Fidelity')
    plt.legend(loc='best')
    plt.tight_layout()
    plt.show()

png

png

png

png

png