XEB fidelity estimator for large circuits

cirq/__init__.py

@@ -85,6 +85,7 @@
 )
 
 from cirq.experiments import (
+    linear_xeb_fidelity,
     generate_supremacy_circuit_google_v2,
     generate_supremacy_circuit_google_v2_bristlecone,
     generate_supremacy_circuit_google_v2_grid,

cirq/experiments/__init__.py

@@ -16,3 +16,6 @@
     build_entangling_layers,
     cross_entropy_benchmarking,
 )
+
+from cirq.experiments.fidelity_estimation import (
+    linear_xeb_fidelity,)

cirq/experiments/fidelity_estimation.py

@@ -0,0 +1,70 @@
+"""Estimation of fidelity associated with experimental circuit executions."""
+
+from typing import Sequence
+
+import numpy as np
+
+from cirq.circuits import Circuit
+from cirq.ops import QubitOrder, QubitOrderOrList
+from cirq.sim import final_wavefunction
+
+
+def linear_xeb_fidelity(
+        circuit: Circuit,
+        bitstrings: Sequence[int],
+        qubit_order: QubitOrderOrList = QubitOrder.DEFAULT,
+) -> float:
+    """Computes fidelity estimate from one circuit using linear XEB estimator.
+
+    Fidelity quantifies the similarity of two quantum states. Here, we estimate
+    the fidelity between the theoretically predicted output state of circuit and
+    the state producted in its experimental realization. Note that we don't know
+    the latter state. Nevertheless, we can estimate the fidelity between the two
+    states from the knowledge of the bitstrings observed in the experiment.
+
+    This estimation procedure makes two assumptions. First, it assumes that the
+    circuit is sufficiently scrambling that its output probabilities follow the
+    Porter-Thomas distribution. This assumption holds for typical instances of
+    random quantum circuits of sufficient depth. Second, it assumes that the
+    circuit uses enough qubits so that the Porter-Thomas distribution can be
+    approximated with the exponential distribution.
+
+    In practice the validity of these assumptions can be confirmed by plotting
+    a histogram of output probabilities and comparing it to the exponential
+    distribution.
+
+    In order to make the estimate more robust one should average the estimates
+    over many random circuits. The API supports per-circuit fidelity estimation
+    to enable users to examine the properties of estimate distribution over
+    many circuits.
+
+    See https://arxiv.org/abs/1608.00263 for more details.
+
+    Args:
+        circuit: Random quantum circuit which has been executed on quantum
+            processor under test
+        bitstrings: Results of terminal all-qubit measurements performed after
+            each circuit execution as integer array where each integer is
+            formed from measured qubit values according to `qubit_order` from
+            most to least significant qubit, i.e. in the order consistent with
+            `cirq.final_wavefunction`.
+        qubit_order: Qubit order used to construct bitstrings enumerating
+            qubits starting with the most sigificant qubit
+    Returns:
+        Estimate of fidelity associated with an experimental realization of
+        circuit which yielded measurements in bitstrings.
+    Raises:
+        ValueError: Circuit is inconsistent with qubit order or one of the
+            bitstrings is inconsistent with the number of qubits.
+    """
+    dim = np.product(circuit.qid_shape())
+
+    for bitstring in bitstrings:
+        if not 0 <= bitstring < dim:
+            raise ValueError(
+                f'Bitstring {bitstring} could not have been observed '
+                f'on {len(circuit.qid_shape())} qubits.')
+
+    output_state = final_wavefunction(circuit, qubit_order=qubit_order)
+    output_probabilities = np.abs(output_state)**2
+    return dim * np.mean(output_probabilities[bitstrings]) - 1

cirq/experiments/fidelity_estimation_test.py

@@ -0,0 +1,85 @@
+import itertools
+
+from typing import Sequence
+
+import numpy as np
+import pytest
+
+import cirq
+
+
+def sample_noisy_bitstrings(circuit: cirq.Circuit,
+                            qubit_order: Sequence[cirq.Qid],
+                            depolarization: float,
+                            repetitions: int) -> np.ndarray:
+    assert 0 <= depolarization <= 1
+    dim = np.product(circuit.qid_shape())
+    n_incoherent = int(depolarization * repetitions)
+    n_coherent = repetitions - n_incoherent
+    incoherent_samples = np.random.randint(dim, size=n_incoherent)
+    circuit_with_measurements = cirq.Circuit.from_ops(
+        circuit, cirq.measure(*qubit_order, key='m'))
+    # TODO(viathor): Remove conditional after #2114.
+    if n_coherent > 0:
+        r = cirq.sample(circuit_with_measurements, repetitions=n_coherent)
+        coherent_samples = r.data['m'].to_numpy()
+        return np.concatenate((coherent_samples, incoherent_samples))
+    return incoherent_samples
+
+
+def make_random_quantum_circuit(qubits: Sequence[cirq.Qid],
+                                depth: int) -> cirq.Circuit:
+    SQ_GATES = [cirq.X**0.5, cirq.Y**0.5, cirq.T]
+    circuit = cirq.Circuit()
+    cz_start = 0
+    for q in qubits:
+        circuit.append(cirq.H(q))
+    for _ in range(depth):
+        for q in qubits:
+            random_gate = SQ_GATES[np.random.randint(len(SQ_GATES))]
+            circuit.append(random_gate(q))
+        for q0, q1 in zip(itertools.islice(qubits, cz_start, None, 2),
+                          itertools.islice(qubits, cz_start + 1, None, 2)):
+            circuit.append(cirq.CNOT(q0, q1))
+        cz_start = 1 - cz_start
+    for q in qubits:
+        circuit.append(cirq.H(q))
+    return circuit
+
+
+@pytest.mark.parametrize('depolarization', (0.0, 0.2, 0.5, 0.7, 1.0))
+def test_linear_xeb_fidelity(depolarization):
+    prng_state = np.random.get_state()
+    np.random.seed(0)
+
+    fs = []
+    for _ in range(10):
+        qubits = cirq.LineQubit.range(5)
+        circuit = make_random_quantum_circuit(qubits, depth=12)
+        bitstrings = sample_noisy_bitstrings(circuit,
+                                             qubits,
+                                             depolarization,
+                                             repetitions=5000)
+        f = cirq.linear_xeb_fidelity(circuit, bitstrings, qubits)
+        fs.append(f)
+    estimated_fidelity = np.mean(fs)
+    expected_fidelity = 1 - depolarization
+    assert np.isclose(estimated_fidelity, expected_fidelity, atol=0.09)
+
+    np.random.set_state(prng_state)
+
+
+def test_linear_xeb_fidelity_invalid_qubits():
+    q0, q1, q2 = cirq.LineQubit.range(3)
+    circuit = cirq.Circuit.from_ops(cirq.H(q0), cirq.CNOT(q0, q1))
+    bitstrings = sample_noisy_bitstrings(circuit, (q0, q1, q2), 0.9, 10)
+    with pytest.raises(ValueError):
+        cirq.linear_xeb_fidelity(circuit, bitstrings, (q0, q2))
+
+
+def test_linear_xeb_fidelity_invalid_bitstrings():
+    q0, q1 = cirq.LineQubit.range(2)
+    circuit = cirq.Circuit.from_ops(cirq.H(q0), cirq.CNOT(q0, q1))
+    bitstrings = [0, 1, 2, 3, 4]
+    with pytest.raises(ValueError):
+        cirq.linear_xeb_fidelity(circuit, bitstrings, (q0, q1))