Cross Entropy Benchmarking Theory

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
    print("installed cirq.")

Cross entropy benchmarking uses the properties of random quantum programs to determine the fidelity of a wide variety of circuits. When applied to circuits with many qubits, XEB can characterize the performance of a large device. When applied to deep, two-qubit circuits it can be used to accurately characterize a two-qubit interaction potentially leading to better calibration.

# Standard imports
import numpy as np
import cirq

from cirq.contrib.svg import SVGCircuit

The action of random circuits with noise

An XEB experiment collects data from the execution of random circuits subject to noise. The effect of applying a random circuit with unitary $U$ is modeled as $U$ followed by a depolarizing channel. The result is that the initial state $|πœ“βŸ©$ is mapped to a density matrix $ρ_U$ as follows:

$$ |πœ“βŸ© β†’ ρ_U = f |πœ“_UβŸ©βŸ¨πœ“_U| + (1 - f) I / D $$

where $|πœ“_U⟩ = U|πœ“βŸ©$, $D$ is the dimension of the Hilbert space, $I / D$ is the maximally mixed state, and $f$ is the fidelity with which the circuit is applied.

For this model to be accurate, we require $U$ to be a random circuit that scrambles errors. In practice, we use a particular circuit ansatz consisting of random single-qubit rotations interleaved with entangling gates.

Possible single-qubit rotations

These 8*8 possible rotations are chosen randomly when constructing the circuit.

Geometrically, we choose 8 axes in the XY plane to perform a quarter-turn (pi/2 rotation) around. This is followed by a rotation around the Z axis of 8 different magnitudes.

exponents = np.linspace(0, 7/4, 8)
exponents
array([0.  , 0.25, 0.5 , 0.75, 1.  , 1.25, 1.5 , 1.75])
import itertools
SINGLE_QUBIT_GATES = [
    cirq.PhasedXZGate(x_exponent=0.5, z_exponent=z, axis_phase_exponent=a)
    for a, z in itertools.product(exponents, repeat=2)
]
SINGLE_QUBIT_GATES[:10], '...'
([cirq.PhasedXZGate(axis_phase_exponent=0.0, x_exponent=0.5, z_exponent=0.0),
  cirq.PhasedXZGate(axis_phase_exponent=0.0, x_exponent=0.5, z_exponent=0.25),
  cirq.PhasedXZGate(axis_phase_exponent=0.0, x_exponent=0.5, z_exponent=0.5),
  cirq.PhasedXZGate(axis_phase_exponent=0.0, x_exponent=0.5, z_exponent=0.75),
  cirq.PhasedXZGate(axis_phase_exponent=0.0, x_exponent=0.5, z_exponent=1.0),
  cirq.PhasedXZGate(axis_phase_exponent=0.0, x_exponent=0.5, z_exponent=1.25),
  cirq.PhasedXZGate(axis_phase_exponent=0.0, x_exponent=0.5, z_exponent=1.5),
  cirq.PhasedXZGate(axis_phase_exponent=0.0, x_exponent=0.5, z_exponent=1.75),
  cirq.PhasedXZGate(axis_phase_exponent=0.25, x_exponent=0.5, z_exponent=0.0),
  cirq.PhasedXZGate(axis_phase_exponent=0.25, x_exponent=0.5, z_exponent=0.25)],
 '...')

Random circuit

We use random_rotations_between_two_qubit_circuit to generate a random two-qubit circuit. Note that we provide the possible single-qubit rotations from above and declare that our two-qubit operation is the $\sqrt{i\mathrm{SWAP} }$ gate.

import cirq.google as cg
from cirq.experiments import random_quantum_circuit_generation as rqcg

SQRT_ISWAP = cirq.ISWAP**0.5

q0, q1 = cirq.LineQubit.range(2)
circuit = rqcg.random_rotations_between_two_qubit_circuit(
    q0, q1, 
    depth=4, 
    two_qubit_op_factory=lambda a, b, _: SQRT_ISWAP(a, b), 
    single_qubit_gates=SINGLE_QUBIT_GATES
)
SVGCircuit(circuit)
findfont: Font family ['Arial'] not found. Falling back to DejaVu Sans.

svg

Estimating fidelity

Let $O_U$ be an observable that is diagonal in the computational basis. Then the expectation value of $O_U$ on $ρ_U$ is given by

$$ Tr(ρ_U O_U) = f βŸ¨πœ“_U|O_U|πœ“_U⟩ + (1 - f) Tr(O_U / D). $$

This equation shows how $f$ can be estimated, since $Tr(ρ_U O_U)$ can be estimated from experimental data, and $βŸ¨πœ“_U|O_U|πœ“_U⟩$ and $Tr(O_U / D)$ can be computed.

Let $e_U = βŸ¨πœ“_U|O_U|πœ“_U⟩$, $u_U = Tr(O_U / D)$, and $m_U$ denote the experimental estimate of $Tr(ρ_U O_U)$. We can write the following linear equation (equivalent to the expression above):

$$ m_U = f e_U + (1-f) u_U \\ m_U - u_U = f (e_U - u_U) $$
# Make long circuits (which we will truncate)
MAX_DEPTH = 100
circuits = [
    rqcg.random_rotations_between_two_qubit_circuit(
        q0, q1, 
        depth=MAX_DEPTH, 
        two_qubit_op_factory=lambda a, b, _: SQRT_ISWAP(a, b), 
        single_qubit_gates=SINGLE_QUBIT_GATES)
    for _ in range(10)
]
# We will truncate to these lengths
cycle_depths = np.arange(3, MAX_DEPTH, 9)
cycle_depths
array([ 3, 12, 21, 30, 39, 48, 57, 66, 75, 84, 93])

Execute circuits

Cross entropy benchmarking requires sampled bitstrings from the device being benchmarked as well as the true probabilities from a noiseless simulation. We find these quantities for all (cycle_depth, circuit) permutations.

pure_sim = cirq.Simulator()

P_DEPOL = 5e-3
noisy_sim = cirq.DensityMatrixSimulator(noise=cirq.depolarize(P_DEPOL))

# These two qubit circuits have 2^2 = 4 probabilities
DIM = 4

records = []
for cycle_depth in cycle_depths:
    for circuit_i, circuit in enumerate(circuits):

        # Truncate the long circuit to the requested cycle_depth
        circuit_depth = cycle_depth * 2 + 1
        assert circuit_depth <= len(circuit)
        trunc_circuit = circuit[:circuit_depth]

        # Pure-state simulation
        psi = pure_sim.simulate(trunc_circuit)
        psi = psi.final_state_vector
        pure_probs = np.abs(psi)**2

        # Noisy execution
        meas_circuit = trunc_circuit + cirq.measure(q0, q1)
        sampled_inds = noisy_sim.sample(meas_circuit, repetitions=10_000).values[:,0]
        sampled_probs = np.bincount(sampled_inds, minlength=DIM) / len(sampled_inds)

        # Save the results
        records += [{
            'circuit_i': circuit_i,
            'cycle_depth': cycle_depth,
            'circuit_depth': circuit_depth,
            'pure_probs': pure_probs,
            'sampled_probs': sampled_probs,
        }]
        print('.', end='', flush=True)
..............................................................................................................

What's the observable

What is $O_U$? Let's define it to be the observable that gives the sum of all probabilities, i.e.

$$ O_U |x \rangle = p(x) |x \rangle $$

for any bitstring $x$. We can use this to derive expressions for our quantities of interest.

$$ e_U = \langle \psi_U | O_U | \psi_U \rangle \\ = \sum_x a_x^* \langle x | O_U | x \rangle a_x \\ = \sum_x p(x) \langle x | O_U | x \rangle \\ = \sum_x p(x) p(x) $$

$e_U$ is simply the sum of squared ideal probabilities. $u_U$ is a normalizing factor that only depends on the operator. Since this operator has the true probabilities in the definition, they show up here anyways.

$$ u_U = \mathrm{Tr}[O_U / D] \\ = 1/D \sum_x \langle x | O_U | x \rangle \\ = 1/D \sum_x p(x) $$

For the measured values, we use the definition of an expectation value

$$ \langle f(x) \rangle_\rho = \sum_x p(x) f(x) $$

It becomes notationally confusing because remember: our operator on basis states returns the ideal probability of that basis state $p(x)$. The probability of observing a measured basis state is estimated from samples and denoted $p_\mathrm{est}(x)$ here.

$$ m_U = \mathrm{Tr}[\rho_U O_U] \\ = \langle O_U \rangle_{\rho_U} = \sum_{x} p_\mathrm{est}(x) p(x) $$
for record in records:
    e_u = np.sum(record['pure_probs']**2)
    u_u = np.sum(record['pure_probs']) / DIM
    m_u = np.sum(record['pure_probs'] * record['sampled_probs'])
    record.update(
        e_u=e_u,
        u_u=u_u,
        m_u=m_u,    
    )

Remember:

$$ m_U - u_U = f (e_U - u_U) $$

We estimate f by performing least squares minimization of the quantity

$$ f (e_U - u_U) - (m_U - u_U) $$

over different random circuits. The solution to the least squares problem is given by

$$ f = (βˆ‘_U (m_U - u_U) * (e_U - u_U)) / (βˆ‘_U (e_U - u_U)^2) $$
import pandas as pd
df = pd.DataFrame(records)
df['y'] = df['m_u'] - df['u_u']
df['x'] = df['e_u'] - df['u_u']

df['numerator'] = df['x'] * df['y']
df['denominator'] = df['x'] ** 2
df.head()

Fit

We'll plot the linear relationship and least-squares fit while we transform the raw DataFrame into one containing fidelities.

%matplotlib inline
from matplotlib import pyplot as plt

# Color by cycle depth
import seaborn as sns
colors = sns.cubehelix_palette(n_colors=len(cycle_depths))        
colors = {k: colors[i] for i, k in enumerate(cycle_depths)}

_lines = []
def per_cycle_depth(df):
    fid_lsq = df['numerator'].sum() / df['denominator'].sum()

    cycle_depth = df.name    
    xx = np.linspace(0, df['x'].max())
    l, = plt.plot(xx, fid_lsq*xx, color=colors[cycle_depth])
    plt.scatter(df['x'], df['y'], color=colors[cycle_depth])

    global _lines
    _lines += [l] # for legend
    return pd.Series({'fid_lsq': fid_lsq})

fids = df.groupby('cycle_depth').apply(per_cycle_depth).reset_index()
plt.xlabel(r'$e_U - u_U$', fontsize=18)
plt.ylabel(r'$m_U - u_U$', fontsize=18)
_lines = np.asarray(_lines)
plt.legend(_lines[[0,-1]], cycle_depths[[0,-1]], loc='best', title='Cycle depth')
plt.tight_layout()

png

Fidelities

plt.plot(fids['cycle_depth'], fids['fid_lsq'], label='LSq')

xx = np.linspace(0, fids['cycle_depth'].max())
plt.plot(xx, (1-P_DEPOL)**(4*xx), label=r'$(1-\mathrm{depol})^{4d}$')

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

png