Pluggable Clifford synthesis for RB circuits (#1288)

### Summary
Make the Clifford synthesis algorithm for RB circuits pluggable
(implementing it as a `HighLevelSynthesisPlugin`).
Fixes #1279 and #1023.
Change to accept Clifford elements consisting only of instructions
supported by the backend for `interleaved_element` option in
`InterleavedRB`.
Speed up 2Q RB/IRB for backends with unidirectional 2q gates, e.g. IBM's
127Q Eagle processors.

### Details and comments
Previously, for 3Q+ RB circuits, entire circuit is transpiled at once
and hence for each of the resulting Cliffords, the initial and the final
layout may differ, that means sampled Cliffords are changed during the
transpilation. Also in the worst case, the resulting circuit may use
physical qubits not in the supplied `physical_qubits`. To avoid that,
this commit changes to transpile an RB sequence Clifford by Clifford.
The Clifford synthesis algorithm (`rb_default`) is implemented as a
`HighLevelSynthesisPlugin` (see `RBDefaultCliffordSynthesis` in
`clifford_synthesis.py`), which forces the initial layout (i.e.
guarantees the initial layout = the final layout) and physical qubits to
use. As a byproduct, the performance of 2Q RB/IRB for backends with
directed 2q gates (e.g. IBM's 127Q Eagle processors) is drastically
improved. For those cases, previously we had to rely on `transpile`
function to make generated circuits comply with the coupling map,
however, after this commit, we can synthesize Cliffords with considering
the 2q-gate direction and go through the fast path introduced in #982.

Depends on Qiskit/qiskit#10477 (qiskit 0.45)

---------

Co-authored-by: Helena Zhang <[email protected]>

Author: itoko

qiskit_experiments/library/randomized_benchmarking/clifford_synthesis.py

@@ -0,0 +1,101 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""
+Clifford synthesis plugins for randomized benchmarking
+"""
+from __future__ import annotations
+
+from typing import Sequence
+
+from qiskit.circuit import QuantumCircuit, Operation
+from qiskit.circuit.equivalence_library import SessionEquivalenceLibrary as sel
+from qiskit.exceptions import QiskitError
+from qiskit.synthesis.clifford import synth_clifford_full
+from qiskit.transpiler import PassManager, CouplingMap, Layout, Target
+from qiskit.transpiler.passes import (
+    SabreSwap,
+    LayoutTransformation,
+    BasisTranslator,
+    CheckGateDirection,
+    GateDirection,
+    Optimize1qGatesDecomposition,
+)
+from qiskit.transpiler.passes.synthesis.plugin import HighLevelSynthesisPlugin
+
+
+class RBDefaultCliffordSynthesis(HighLevelSynthesisPlugin):
+    """Default Clifford synthesis plugin for randomized benchmarking."""
+
+    def run(
+        self,
+        high_level_object: Operation,
+        coupling_map: CouplingMap | None = None,
+        target: Target | None = None,
+        qubits: Sequence | None = None,
+        **options,
+    ) -> QuantumCircuit:
+        """Run synthesis for the given Clifford.
+
+        Args:
+            high_level_object: The operation to synthesize to a
+                :class:`~qiskit.circuit.QuantumCircuit` object.
+            coupling_map: The reduced coupling map of the backend. For example,
+                if physical qubits [5, 6, 7] to be benchmarked is connected
+                as 5 - 7 - 6 linearly, the reduced coupling map is 0 - 2 - 1.
+            target: A target representing the target backend, which will be ignored in this plugin.
+            qubits: List of physical qubits over which the operation is defined,
+                which will be ignored in this plugin.
+            options: Additional method-specific optional kwargs,
+                which must include ``basis_gates``, basis gates to be used for the synthesis.
+
+        Returns:
+            The quantum circuit representation of the Operation
+                when successful, and ``None`` otherwise.
+
+        Raises:
+            QiskitError: If basis_gates is not supplied.
+        """
+        # synthesize cliffords
+        circ = synth_clifford_full(high_level_object)
+
+        # post processing to comply with basis gates and coupling map
+        if coupling_map is None:  # Sabre does not work with coupling_map=None
+            return circ
+
+        basis_gates = options.get("basis_gates", None)
+        if basis_gates is None:
+            raise QiskitError("basis_gates are required to run this synthesis plugin")
+
+        basis_gates = list(basis_gates)
+
+        # Run Sabre routing and undo the layout change
+        # assuming Sabre routing does not change the initial layout.
+        # And then decompose swap gates, fix 2q-gate direction and optimize 1q gates
+        initial_layout = Layout.generate_trivial_layout(*circ.qubits)
+        undo_layout_change = LayoutTransformation(
+            coupling_map=coupling_map, from_layout="final_layout", to_layout=initial_layout
+        )
+
+        def _direction_condition(property_set):
+            return not property_set["is_direction_mapped"]
+
+        pm = PassManager(
+            [
+                SabreSwap(coupling_map),
+                undo_layout_change,
+                BasisTranslator(sel, basis_gates),
+                CheckGateDirection(coupling_map),
+            ]
+        )
+        pm.append([GateDirection(coupling_map)], condition=_direction_condition)
+        pm.append([Optimize1qGatesDecomposition(basis=basis_gates)])
+        return pm.run(circ)

qiskit_experiments/library/randomized_benchmarking/clifford_utils.py

@@ -29,8 +29,12 @@
 from qiskit.compiler import transpile
 from qiskit.exceptions import QiskitError
 from qiskit.quantum_info import Clifford, random_clifford
+from qiskit.transpiler import CouplingMap, PassManager
+from qiskit.transpiler.passes.synthesis.high_level_synthesis import HLSConfig, HighLevelSynthesis
 from qiskit.utils.deprecation import deprecate_func
 
+DEFAULT_SYNTHESIS_METHOD = "rb_default"
+
 _DATA_FOLDER = os.path.join(os.path.dirname(__file__), "data")
 
 _CLIFFORD_COMPOSE_1Q = np.load(f"{_DATA_FOLDER}/clifford_compose_1q.npz")["table"]
@@ -110,37 +114,141 @@ def _circuit_compose(
     return self
 
 
-def _truncate_inactive_qubits(
-    circ: QuantumCircuit, active_qubits: Sequence[Qubit]
+def _synthesize_clifford(
+    clifford: Clifford,
+    basis_gates: Optional[Tuple[str]],
+    coupling_tuple: Optional[Tuple[Tuple[int, int]]] = None,
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
 ) -> QuantumCircuit:
-    res = QuantumCircuit(active_qubits, name=circ.name, metadata=circ.metadata)
-    for inst in circ:
-        if all(q in active_qubits for q in inst.qubits):
-            res.append(inst)
-    res.calibrations = circ.calibrations
-    return res
+    """Synthesize a circuit of a Clifford element. The resulting circuit contains only
+    ``basis_gates`` and it complies with ``coupling_tuple``.
+
+    Args:
+        clifford: Clifford element to be converted
+        basis_gates: basis gates to use in the conversion
+        coupling_tuple: coupling map to use in the conversion in the form of tuple of edges
+        synthesis_method: conversion algorithm name
+
+    Returns:
+        Synthesized circuit
+    """
+    qc = QuantumCircuit(clifford.num_qubits, name=str(clifford))
+    qc.append(clifford, qc.qubits)
+    return _synthesize_clifford_circuit(
+        qc,
+        basis_gates=basis_gates,
+        coupling_tuple=coupling_tuple,
+        synthesis_method=synthesis_method,
+    )
 
 
 def _synthesize_clifford_circuit(
-    circuit: QuantumCircuit, basis_gates: Tuple[str]
+    circuit: QuantumCircuit,
+    basis_gates: Optional[Tuple[str]],
+    coupling_tuple: Optional[Tuple[Tuple[int, int]]] = None,
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
 ) -> QuantumCircuit:
-    # synthesizes clifford circuits using given basis gates, for use during
-    # custom transpilation during RB circuit generation.
-    return transpile(circuit, basis_gates=list(basis_gates), optimization_level=1)
+    """Convert a Clifford circuit into one composed of ``basis_gates`` with
+    satisfying ``coupling_tuple`` using the specified synthesis method.
+
+    Args:
+        circuit: Clifford circuit to be converted
+        basis_gates: basis gates to use in the conversion
+        coupling_tuple: coupling map to use in the conversion in the form of tuple of edges
+        synthesis_method: name of Clifford synthesis algorithm to use
+
+    Returns:
+        Synthesized circuit
+    """
+    if basis_gates:
+        basis_gates = list(basis_gates)
+    coupling_map = CouplingMap(coupling_tuple) if coupling_tuple else None
+
+    # special handling for 1q or 2q case for speed
+    if circuit.num_qubits <= 2:
+        if synthesis_method == DEFAULT_SYNTHESIS_METHOD:
+            return transpile(
+                circuit,
+                basis_gates=basis_gates,
+                coupling_map=coupling_map,
+                optimization_level=1,
+            )
+        else:
+            # Provided custom synthesis method, re-synthesize Clifford circuit
+            # convert the circuit back to a Clifford object and then call the synthesis plugin
+            new_circuit = QuantumCircuit(circuit.num_qubits, name=circuit.name)
+            new_circuit.append(Clifford(circuit), new_circuit.qubits)
+            circuit = new_circuit
+
+    # for 3q+ or custom synthesis method, synthesizes clifford circuit
+    hls_config = HLSConfig(clifford=[(synthesis_method, {"basis_gates": basis_gates})])
+    pm = PassManager([HighLevelSynthesis(hls_config=hls_config, coupling_map=coupling_map)])
+    circuit = pm.run(circuit)
+    return circuit
 
 
-@lru_cache(maxsize=None)
+@lru_cache(maxsize=256)
 def _clifford_1q_int_to_instruction(
-    num: Integral, basis_gates: Optional[Tuple[str]]
+    num: Integral,
+    basis_gates: Optional[Tuple[str]],
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
 ) -> Instruction:
-    return CliffordUtils.clifford_1_qubit_circuit(num, basis_gates).to_instruction()
+    return CliffordUtils.clifford_1_qubit_circuit(
+        num, basis_gates=basis_gates, synthesis_method=synthesis_method
+    ).to_instruction()
 
 
 @lru_cache(maxsize=11520)
 def _clifford_2q_int_to_instruction(
-    num: Integral, basis_gates: Optional[Tuple[str]]
+    num: Integral,
+    basis_gates: Optional[Tuple[str]],
+    coupling_tuple: Optional[Tuple[Tuple[int, int]]],
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
 ) -> Instruction:
-    return CliffordUtils.clifford_2_qubit_circuit(num, basis_gates).to_instruction()
+    return CliffordUtils.clifford_2_qubit_circuit(
+        num,
+        basis_gates=basis_gates,
+        coupling_tuple=coupling_tuple,
+        synthesis_method=synthesis_method,
+    ).to_instruction()
+
+
+def _hash_cliff(cliff):
+    return cliff.tableau.tobytes(), cliff.tableau.shape
+
+
+def _dehash_cliff(cliff_hash):
+    tableau = np.frombuffer(cliff_hash[0], dtype=bool).reshape(cliff_hash[1])
+    return Clifford(tableau)
+
+
+def _clifford_to_instruction(
+    clifford: Clifford,
+    basis_gates: Optional[Tuple[str]],
+    coupling_tuple: Optional[Tuple[Tuple[int, int]]],
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
+) -> Instruction:
+    return _cached_clifford_to_instruction(
+        _hash_cliff(clifford),
+        basis_gates=basis_gates,
+        coupling_tuple=coupling_tuple,
+        synthesis_method=synthesis_method,
+    )
+
+
+@lru_cache(maxsize=256)
+def _cached_clifford_to_instruction(
+    cliff_hash: Tuple[str, Tuple[int, int]],
+    basis_gates: Optional[Tuple[str]],
+    coupling_tuple: Optional[Tuple[Tuple[int, int]]],
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
+) -> Instruction:
+    return _synthesize_clifford(
+        _dehash_cliff(cliff_hash),
+        basis_gates=basis_gates,
+        coupling_tuple=coupling_tuple,
+        synthesis_method=synthesis_method,
+    ).to_instruction()
 
 
 # The classes VGate and WGate are not actually used in the code - we leave them here to give
@@ -254,7 +362,12 @@ def random_clifford_circuits(
 
     @classmethod
     @lru_cache(maxsize=24)
-    def clifford_1_qubit_circuit(cls, num, basis_gates: Optional[Tuple[str, ...]] = None):
+    def clifford_1_qubit_circuit(
+        cls,
+        num,
+        basis_gates: Optional[Tuple[str, ...]] = None,
+        synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
+    ):
         """Return the 1-qubit clifford circuit corresponding to ``num``,
         where ``num`` is between 0 and 23.
         """
@@ -275,20 +388,28 @@ def clifford_1_qubit_circuit(cls, num, basis_gates: Optional[Tuple[str, ...]] =
             qc.z(0)
 
         if basis_gates:
-            qc = _synthesize_clifford_circuit(qc, basis_gates)
+            qc = _synthesize_clifford_circuit(qc, basis_gates, synthesis_method=synthesis_method)
 
         return qc
 
     @classmethod
     @lru_cache(maxsize=11520)
-    def clifford_2_qubit_circuit(cls, num, basis_gates: Optional[Tuple[str, ...]] = None):
+    def clifford_2_qubit_circuit(
+        cls,
+        num,
+        basis_gates: Optional[Tuple[str, ...]] = None,
+        coupling_tuple: Optional[Tuple[Tuple[int, int]]] = None,
+        synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
+    ):
         """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)):
             if basis_gates:
-                layer_circ = _transformed_clifford_layer(layer, idx, basis_gates)
+                layer_circ = _transformed_clifford_layer(
+                    layer, idx, basis_gates, coupling_tuple, synthesis_method=synthesis_method
+                )
             else:
                 layer_circ = _CLIFFORD_LAYER[layer][idx]
             _circuit_compose(qc, layer_circ, qubits=(0, 1))
@@ -578,13 +699,22 @@ def _clifford_2q_nums_from_2q_circuit(qc: QuantumCircuit) -> Iterable[Integral]:
 ]
 
 
-@lru_cache(maxsize=None)
+@lru_cache(maxsize=256)
 def _transformed_clifford_layer(
-    layer: int, index: Integral, basis_gates: Tuple[str, ...]
+    layer: int,
+    index: Integral,
+    basis_gates: Tuple[str, ...],
+    coupling_tuple: Optional[Tuple[Tuple[int, int]]],
+    synthesis_method: str = DEFAULT_SYNTHESIS_METHOD,
 ) -> QuantumCircuit:
     # Return the index-th quantum circuit of the layer translated with the basis_gates.
     # The result is cached for speed.
-    return _synthesize_clifford_circuit(_CLIFFORD_LAYER[layer][index], basis_gates)
+    return _synthesize_clifford_circuit(
+        _CLIFFORD_LAYER[layer][index],
+        basis_gates=basis_gates,
+        coupling_tuple=coupling_tuple,
+        synthesis_method=synthesis_method,
+    )
 
 
 def _num_from_layer_indices(triplet: Tuple[Integral, Integral, Integral]) -> Integral: