Randomized benchmarking Clifford synthesis optimization

[wshanks]

The randomized benchmarking circuit generation takes sequences of Cliffords and converts them to basis gates by using a synthesis routine (that can optionally be changed) and then transpiling to a given basis gates list with optimization level 1. Transpilation is done on individual Cliffords so that the results can be cached rather than transpiling entire circuits. So the standard transpilation step is not performed. Recently in #1567, the two qubit circuit generation was improved. In that case, the default synthesis routine produced two qubit Cliffords in three steps and the transpilation was done on each step but not on the combination so there was room left to combine single qubit gates across the three layers.

While looking into #1567, it was noticed that the transpiled randomized benchmarking sequences were still not optimized. In particular, the transpilation uses only a list of gates. For IBM backends, rz is a virtual gate with zero error. Ideally, the sx and x gate counts would be minimized while allowing as many rz gates as necessary. Doing this would allow to keep the single qubit gate count as small as possible compared to the two qubit gate count. Usually, one wants to infer the two qubit gate error which is larger than the single qubit gate error. The smaller the single qubit gate count the more accurate it is just to ignore the single qubit gate error contribution.

The following code transpiles every two qubit Clifford for the (sx, x, rz, cz) basis set and looks at the averge gate counts for different optimization levels and with and without setting errors on the gates (with 0 error on rz):

import itertools
import time
from functools import lru_cache
from numbers import Integral
from typing import Optional, Sequence, Tuple, Union

import numpy

from qiskit import QuantumCircuit, transpile
from qiskit.circuit import CircuitInstruction, Qubit
from qiskit.circuit.library import XGate, YGate, ZGate


samples = 10_000
seed = 1234

rng = numpy.random.default_rng(seed=seed)

NUM_CLIFFORD_2_QUBIT = 11520

clifford_indices = rng.integers(NUM_CLIFFORD_2_QUBIT, size=samples)


def _append_v_w(qc, vw0, vw1):
    if vw0 == "v":
        qc.sdg(0)
        qc.h(0)
    elif vw0 == "w":
        qc.h(0)
        qc.s(0)
    if vw1 == "v":
        qc.sdg(1)
        qc.h(1)
    elif vw1 == "w":
        qc.h(1)
        qc.s(1)


def _create_cliff_2q_layer_0():
    """Layer 0 consists of 0 or 1 H gates on each qubit, followed by 0/1/2 V gates on each qubit.
    Number of Cliffords == 36."""
    circuits = []
    num_h = [0, 1]
    v_w_gates = ["i", "v", "w"]
    for h0, h1, v0, v1 in itertools.product(num_h, num_h, v_w_gates, v_w_gates):
        qc = QuantumCircuit(2)
        for _ in range(h0):
            qc.h(0)
        for _ in range(h1):
            qc.h(1)
        _append_v_w(qc, v0, v1)
        circuits.append(qc)
    return circuits


def _create_cliff_2q_layer_1():
    """Layer 1 consists of one of the following:
    - nothing
    - cx(0,1) followed by 0/1/2 V gates on each qubit
    - cx(0,1), cx(1,0) followed by 0/1/2 V gates on each qubit
    - cx(0,1), cx(1,0), cx(0,1)
    Number of Cliffords == 20."""
    circuits = [QuantumCircuit(2)]  # identity at the beginning

    v_w_gates = ["i", "v", "w"]
    for v0, v1 in itertools.product(v_w_gates, v_w_gates):
        qc = QuantumCircuit(2)
        qc.cx(0, 1)
        _append_v_w(qc, v0, v1)
        circuits.append(qc)

    for v0, v1 in itertools.product(v_w_gates, v_w_gates):
        qc = QuantumCircuit(2)
        qc.cx(0, 1)
        qc.cx(1, 0)
        _append_v_w(qc, v0, v1)
        circuits.append(qc)

    qc = QuantumCircuit(2)  # swap at the end
    qc.cx(0, 1)
    qc.cx(1, 0)
    qc.cx(0, 1)
    circuits.append(qc)
    return circuits


def _create_cliff_2q_layer_2():
    """Layer 2 consists of a Pauli gate on each qubit {Id, X, Y, Z}.
    Number of Cliffords == 16."""
    circuits = []
    pauli = ("i", XGate(), YGate(), ZGate())
    for p0, p1 in itertools.product(pauli, pauli):
        qc = QuantumCircuit(2)
        if p0 != "i":
            qc.append(p0, [0])
        if p1 != "i":
            qc.append(p1, [1])
        circuits.append(qc)
    return circuits


_CLIFFORD_LAYER = (
    _create_cliff_2q_layer_0(),
    _create_cliff_2q_layer_1(),
    _create_cliff_2q_layer_2(),
)
_NUM_LAYER_1 = 20
_NUM_LAYER_2 = 16

basis_gates = ["sx", "x", "rz", "cz"]
coupling_tuple = ((0, 1),)
synthesis_method = "rb_default"


def _layer_indices_from_num(num: Integral) -> Tuple[Integral, Integral, Integral]:
    """Return the triplet of layer indices corresponding to the input number."""
    idx2 = num % _NUM_LAYER_2
    num = num // _NUM_LAYER_2
    idx1 = num % _NUM_LAYER_1
    idx0 = num // _NUM_LAYER_1
    return idx0, idx1, idx2


def _circuit_compose(
    self: QuantumCircuit, other: QuantumCircuit, qubits: Sequence[Union[Qubit, int]]
) -> QuantumCircuit:
    # Simplified QuantumCircuit.compose with clbits=None, front=False, inplace=True, wrap=False
    # without any validation, parameter_table/calibrations updates and copy of operations
    # The input circuit `self` is changed inplace.
    qubit_map = {
        other.qubits[i]: (self.qubits[q] if isinstance(q, int) else q) for i, q in enumerate(qubits)
    }
    for instr in other:
        self._data.append(
            CircuitInstruction(
                operation=instr.operation,
                qubits=[qubit_map[q] for q in instr.qubits],
                clbits=instr.clbits,
            ),
        )
    self.global_phase += other.global_phase
    return self


@lru_cache(maxsize=256)
def _transformed_clifford_layer(
    layer: int,
    index: Integral,
    basis_gates: Tuple[str, ...],
) -> QuantumCircuit:
    # Return the index-th quantum circuit of the layer translated with the basis_gates.
    # The result is cached for speed.
    return transpile(
        _CLIFFORD_LAYER[layer][index],
        basis_gates=list(basis_gates),
        coupling_map=None,
        optimization_level=1,
    )


def clifford_2_qubit_circuit_no_error(
    num,
    basis_gates: list[str, ...],
    final_transpile_opt_level: int,
):
    """Return the 2-qubit clifford circuit corresponding to `num`
    where `num` is between 0 and 11519.
    """
    qc = QuantumCircuit(2, name=f"Clifford-2Q({num})")
    for layer, idx in enumerate(_layer_indices_from_num(num)):
        layer_circ = _transformed_clifford_layer(
            layer, idx, tuple(basis_gates)
        )
        _circuit_compose(qc, layer_circ, qubits=(0, 1))

    # This transpile step is not included in qiskit-experiments
    if final_transpile_opt_level > 0:
        qc = transpile(
            qc,
            basis_gates=tuple(basis_gates),
            coupling_map=None,
            optimization_level=final_transpile_opt_level,
        )
    return qc


def clifford_2_qubit_circuit(
    num,
    basis_gates: list[str, ...],
    final_transpile_opt_level: int,
    final_transpile_errors: bool,
):
    if final_transpile_errors:
        return clifford_2_qubit_circuit_with_error(num, basis_gates, final_transpile_opt_level)
    else:
        return clifford_2_qubit_circuit_no_error(num, basis_gates, final_transpile_opt_level)


def clifford_2_qubit_count_gates(
    num,
    basis_gates: list[str, ...],
    final_transpile_opt_level: int,
    final_transpile_errors: bool,
):
    """Return the gate count for 2-qubit clifford circuit corresponding to `num`
    where `num` is between 0 and 11519.
    """
    qc = clifford_2_qubit_circuit(num, basis_gates, final_transpile_opt_level, final_transpile_errors)
    counts1q = 0
    countsrz = 0
    counts2q = 0
    for inst in qc.data:
        if inst.operation.num_qubits == 1 and inst.operation.name != "rz":
            counts1q += 1
        elif inst.operation.name == "rz":
            countsrz += 1
        elif inst.operation.num_qubits == 2:
            counts2q += 1
    return counts1q, countsrz, counts2q


def clifford_2_qubit_circuit_with_error(
    num,
    basis_gates: list[str, ...],
    final_transpile_opt_level: int,
):
    """Return the 2-qubit clifford circuit corresponding to `num`
    where `num` is between 0 and 11519.
    """
    qc = QuantumCircuit(2, name=f"Clifford-2Q({num})")
    for layer, idx in enumerate(_layer_indices_from_num(num)):
        layer_circ = _transformed_clifford_layer(
            layer, idx, tuple(basis_gates)
        )
        _circuit_compose(qc, layer_circ, qubits=(0, 1))

    from qiskit.transpiler import Target, InstructionProperties
    from qiskit.circuit.library import CZGate, SXGate, XGate, RZGate
    from qiskit.circuit import Parameter
    target = Target(num_qubits=2)
    target.add_instruction(
        CZGate(),
        {(0, 1): InstructionProperties(duration=68e-9, error=3e-3)},
    )
    for gate in (SXGate(), XGate()):
        target.add_instruction(
            gate,
            {
                (q,): InstructionProperties(duration=68e-9, error=4e-4)
                for q in range(2)
            },
        )
    target.add_instruction(
        RZGate(Parameter("t")),
        {
            (0,): InstructionProperties(duration=0, error=0),
            (1,): InstructionProperties(duration=0, error=0),
        },
    )

    # This transpile step is not included in qiskit-experiments
    if final_transpile_opt_level > 0:
        qc = transpile(
            qc,
            target=target,
            optimization_level=final_transpile_opt_level,
        )
    return qc


for opt_level in range(3):
    start = time.perf_counter()
    for errors in (True, False):
        counts = [
            clifford_2_qubit_count_gates(i, ["cz", "rz", "sx", "x"], opt_level, errors)
            for i in range(NUM_CLIFFORD_2_QUBIT)
        ]
    gen_time = time.perf_counter() - start
    counts1q, countsrz, counts2q = numpy.average(counts, axis=0)
    print(f"With{'out' if not errors else ''} errors:")
    print(f"Found {counts1q:.2f} 1Q, {countsrz:.2f} rz, {counts2q} 2Q in {gen_time:0.1f} sec for opt_level={opt_level}")

The output of the script (with Qiskit 2.0.3) is:

With errors:
Found 5.63 1Q, 10.27 rz, 1.5 2Q in 7.0 sec for opt_level=0
Without errors:
Found 5.63 1Q, 10.27 rz, 1.5 2Q in 9.7 sec for opt_level=0
With errors:
Found 4.03 1Q, 6.92 rz, 1.5 2Q in 104.3 sec for opt_level=1
Without errors:
Found 4.10 1Q, 6.84 rz, 1.5 2Q in 173.4 sec for opt_level=1
With errors:
Found 3.91 1Q, 4.36 rz, 1.4302083333333333 2Q in 164.4 sec for opt_level=2
Without errors:
Found 4.38 1Q, 4.39 rz, 1.4302083333333333 2Q in 304.7 sec for opt_level=2

With #1567, two qubit StandardRB now corresponds to the optimization level 1 without error case, which has 4.1 sx or x gates (down from 5.6 before #1567). We see that without considering errors optimization level 2 actually produces more sx and x gates in order to reduce the rz (and total) gate count. If errors were considered the sx+x gate count could be reduced further to 4.03 (optimization level 2 finds a way to reduce gate counts further but keeps the ratio of single to two qubit gates the same, so it does not provide any benefit for benchmarking). Further optimizing the gate count by doing a brute force search through the possible gate combinations to produce the Cliffords can produce a set with an average count of sx+x gates of 3.3 (code for this search not shown here).

We can choose not to take any action on this information. I wanted to document it for now in case we want to do something to change the gate counts in the future.

[wshanks]

Here is a variant where only the optimization stage of the preset pass manager is used instead of calling transpile() and using all the stages. Additionally, this version drops the first transpile and only transpiles the Clifford circuit with all the layers appended.

import itertools
import time
from functools import lru_cache
from numbers import Integral
from typing import Sequence, Tuple, Union

import numpy

from qiskit import QuantumCircuit, transpile, generate_preset_pass_manager
from qiskit.circuit import CircuitInstruction, Qubit
from qiskit.circuit.library import XGate, YGate, ZGate


samples = 10_000
seed = 1234

rng = numpy.random.default_rng(seed=seed)

NUM_CLIFFORD_2_QUBIT = 11520

clifford_indices = rng.integers(NUM_CLIFFORD_2_QUBIT, size=samples)


def _append_v_w(qc, vw0, vw1):
    if vw0 == "v":
        qc.sdg(0)
        qc.h(0)
    elif vw0 == "w":
        qc.h(0)
        qc.s(0)
    if vw1 == "v":
        qc.sdg(1)
        qc.h(1)
    elif vw1 == "w":
        qc.h(1)
        qc.s(1)


def _create_cliff_2q_layer_0():
    """Layer 0 consists of 0 or 1 H gates on each qubit, followed by 0/1/2 V gates on each qubit.
    Number of Cliffords == 36."""
    circuits = []
    num_h = [0, 1]
    v_w_gates = ["i", "v", "w"]
    for h0, h1, v0, v1 in itertools.product(num_h, num_h, v_w_gates, v_w_gates):
        qc = QuantumCircuit(2)
        for _ in range(h0):
            qc.h(0)
        for _ in range(h1):
            qc.h(1)
        _append_v_w(qc, v0, v1)
        circuits.append(qc)
    return circuits


def _create_cliff_2q_layer_1():
    """Layer 1 consists of one of the following:
    - nothing
    - cx(0,1) followed by 0/1/2 V gates on each qubit
    - cx(0,1), cx(1,0) followed by 0/1/2 V gates on each qubit
    - cx(0,1), cx(1,0), cx(0,1)
    Number of Cliffords == 20."""
    circuits = [QuantumCircuit(2)]  # identity at the beginning

    v_w_gates = ["i", "v", "w"]
    for v0, v1 in itertools.product(v_w_gates, v_w_gates):
        qc = QuantumCircuit(2)
        qc.cx(0, 1)
        _append_v_w(qc, v0, v1)
        circuits.append(qc)

    for v0, v1 in itertools.product(v_w_gates, v_w_gates):
        qc = QuantumCircuit(2)
        qc.cx(0, 1)
        qc.cx(1, 0)
        _append_v_w(qc, v0, v1)
        circuits.append(qc)

    qc = QuantumCircuit(2)  # swap at the end
    qc.cx(0, 1)
    qc.cx(1, 0)
    qc.cx(0, 1)
    circuits.append(qc)
    return circuits


def _create_cliff_2q_layer_2():
    """Layer 2 consists of a Pauli gate on each qubit {Id, X, Y, Z}.
    Number of Cliffords == 16."""
    circuits = []
    pauli = ("i", XGate(), YGate(), ZGate())
    for p0, p1 in itertools.product(pauli, pauli):
        qc = QuantumCircuit(2)
        if p0 != "i":
            qc.append(p0, [0])
        if p1 != "i":
            qc.append(p1, [1])
        circuits.append(qc)
    return circuits


_CLIFFORD_LAYER = (
    _create_cliff_2q_layer_0(),
    _create_cliff_2q_layer_1(),
    _create_cliff_2q_layer_2(),
)
_NUM_LAYER_1 = 20
_NUM_LAYER_2 = 16

basis_gates = ["sx", "x", "rz", "cz"]
coupling_tuple = ((0, 1),)
synthesis_method = "rb_default"


def _layer_indices_from_num(num: Integral) -> Tuple[Integral, Integral, Integral]:
    """Return the triplet of layer indices corresponding to the input number."""
    idx2 = num % _NUM_LAYER_2
    num = num // _NUM_LAYER_2
    idx1 = num % _NUM_LAYER_1
    idx0 = num // _NUM_LAYER_1
    return idx0, idx1, idx2


def _circuit_compose(
    self: QuantumCircuit, other: QuantumCircuit, qubits: Sequence[Union[Qubit, int]]
) -> QuantumCircuit:
    # Simplified QuantumCircuit.compose with clbits=None, front=False, inplace=True, wrap=False
    # without any validation, parameter_table/calibrations updates and copy of operations
    # The input circuit `self` is changed inplace.
    qubit_map = {
        other.qubits[i]: (self.qubits[q] if isinstance(q, int) else q) for i, q in enumerate(qubits)
    }
    for instr in other:
        self._data.append(
            CircuitInstruction(
                operation=instr.operation,
                qubits=[qubit_map[q] for q in instr.qubits],
                clbits=instr.clbits,
            ),
        )
    self.global_phase += other.global_phase
    return self


@lru_cache(maxsize=256)
def _transformed_clifford_layer(
    layer: int,
    index: Integral,
    basis_gates: Tuple[str, ...],
) -> QuantumCircuit:
    # Return the index-th quantum circuit of the layer translated with the basis_gates.
    # The result is cached for speed.
    return _CLIFFORD_LAYER[layer][index]
    # return transpile(
    #     _CLIFFORD_LAYER[layer][index],
    #     basis_gates=list(basis_gates),
    #     coupling_map=None,
    #     optimization_level=1,
    # )


def clifford_2_qubit_circuit_no_error(
    num,
    basis_gates: list[str, ...],
    final_transpile_opt_level: int,
):
    """Return the 2-qubit clifford circuit corresponding to `num`
    where `num` is between 0 and 11519.
    """
    qc = QuantumCircuit(2, name=f"Clifford-2Q({num})")
    for layer, idx in enumerate(_layer_indices_from_num(num)):
        layer_circ = _transformed_clifford_layer(
            layer, idx, tuple(basis_gates)
        )
        _circuit_compose(qc, layer_circ, qubits=(0, 1))

    # This transpile step is not included in qiskit-experiments
    if final_transpile_opt_level > 0:
        pm = generate_preset_pass_manager(basis_gates=basis_gates, optimization_level=final_transpile_opt_level)
        qc = pm.optimization.run(qc)
    return qc


def clifford_2_qubit_circuit(
    num,
    basis_gates: list[str, ...],
    final_transpile_opt_level: int,
    final_transpile_errors: bool,
):
    if final_transpile_errors:
        return clifford_2_qubit_circuit_with_error(num, basis_gates, final_transpile_opt_level)
    else:
        return clifford_2_qubit_circuit_no_error(num, basis_gates, final_transpile_opt_level)


def clifford_2_qubit_count_gates(
    num,
    basis_gates: list[str, ...],
    final_transpile_opt_level: int,
    final_transpile_errors: bool,
):
    """Return the gate count for 2-qubit clifford circuit corresponding to `num`
    where `num` is between 0 and 11519.
    """
    qc = clifford_2_qubit_circuit(num, basis_gates, final_transpile_opt_level, final_transpile_errors)
    counts1q = 0
    countsrz = 0
    counts2q = 0
    for inst in qc.data:
        if inst.operation.num_qubits == 1 and inst.operation.name != "rz":
            counts1q += 1
        elif inst.operation.name == "rz":
            countsrz += 1
        elif inst.operation.num_qubits == 2:
            counts2q += 1
    return counts1q, countsrz, counts2q


def clifford_2_qubit_circuit_with_error(
    num,
    basis_gates: list[str, ...],
    final_transpile_opt_level: int,
):
    """Return the 2-qubit clifford circuit corresponding to `num`
    where `num` is between 0 and 11519.
    """
    qc = QuantumCircuit(2, name=f"Clifford-2Q({num})")
    for layer, idx in enumerate(_layer_indices_from_num(num)):
        layer_circ = _transformed_clifford_layer(
            layer, idx, tuple(basis_gates)
        )
        _circuit_compose(qc, layer_circ, qubits=(0, 1))

    from qiskit.transpiler import Target, InstructionProperties
    from qiskit.circuit.library import CZGate, SXGate, XGate, RZGate
    from qiskit.circuit import Parameter
    target = Target(num_qubits=2)
    target.add_instruction(
        CZGate(),
        {(0, 1): InstructionProperties(duration=68e-9, error=3e-3)},
    )
    for gate in (SXGate(), XGate()):
        target.add_instruction(
            gate,
            {
                (q,): InstructionProperties(duration=68e-9, error=4e-4)
                for q in range(2)
            },
        )
    target.add_instruction(
        RZGate(Parameter("t")),
        {
            (0,): InstructionProperties(duration=0, error=0),
            (1,): InstructionProperties(duration=0, error=0),
        },
    )

    # This transpile step is not included in qiskit-experiments
    if final_transpile_opt_level > 0:
        pm = generate_preset_pass_manager(target=target, optimization_level=final_transpile_opt_level)
        qc = pm.optimization.run(qc)
    return qc


for opt_level in range(3):
    start = time.perf_counter()
    for errors in (True, False):
        counts = [
            clifford_2_qubit_count_gates(i, ["cz", "rz", "sx", "x"], opt_level, errors)
            for i in range(NUM_CLIFFORD_2_QUBIT)
        ]
        gen_time = time.perf_counter() - start
        counts1q, countsrz, counts2q = numpy.average(counts, axis=0)
        print(f"With{'out' if not errors else ''} errors:")
        print(f"Found {counts1q:.2f} 1Q, {countsrz:.2f} rz, {counts2q} 2Q in {gen_time:0.1f} sec for opt_level={opt_level}")

The output is:

With errors:
Found 7.57 1Q, 0.00 rz, 1.5 2Q in 2.5 sec for opt_level=0
Without errors:
Found 7.57 1Q, 0.00 rz, 1.5 2Q in 3.9 sec for opt_level=0
With errors:
Found 4.03 1Q, 6.92 rz, 1.5 2Q in 101.5 sec for opt_level=1
Without errors:
Found 4.03 1Q, 6.92 rz, 1.5 2Q in 180.3 sec for opt_level=1
With errors:
Found 4.31 1Q, 4.52 rz, 1.5 2Q in 140.7 sec for opt_level=2
Without errors:
Found 4.31 1Q, 4.52 rz, 1.5 2Q in 270.8 sec for opt_level=2

We are not too interested in optimization level 0. It is not getting transpiled at all with this version of the code.

Oddly, for levels 1 and 2, we get the same results whether or not we consider gate error (but with vastly different run times). The sx/x count gets better for optimization level 1 without errors this way. It seems like something might be off though. Somehow the sx/x count goes up for optimization level 2 even when gate error is considered compared to level 1. This version skipping the layout and routing stages produces 1.5 cz per Clifford which is what should be expected. Maybe steps in other stages (like init?) are also important for overall gate count optimization.