Improve performance of multiplication of cliffords (#1275)

### Summary

The circuit generation for 2-qubit RB experiments is slow due to the
Clifford multiplication. Internally Clifford numbers are converted to
Qiskit circuits, decomposed and then multiplied. We can avoid the
conversion to circuits and indexing in a sparse table by creating a
dense version of `_CLIFFORD_COMPOSE_2Q`.

### Details and comments

Bechmark:
```
import time
import qiskit
from qiskit_experiments.library import InterleavedRB

gate=qiskit.circuit.library.CXGate()
rb=InterleavedRB(gate, physical_qubits= (0,1),lengths= [1,10,50,100,200, 1000], num_samples= 100)

t0=time.time()
c=rb.circuits()
dt=time.time()-t0
print(dt)
```
improved from 21.2  to 7.9 seconds.

Micro benchmark
```
from qiskit_experiments.library.randomized_benchmarking.clifford_utils import compose_2q
%timeit compose_2q(15, 10412)
#main: 56.2 µs ± 11.4 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
#PR: 3.58 µs ± 362 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
```

### PR checklist (delete when all criteria are met)

- [x] I have read the contributing guide `CONTRIBUTING.md`.
- [x] I have added the tests to cover my changes. (no changes to public
methods)
- [x] I have updated the documentation accordingly. (no changes to
public methods)
- [x] I have added a release note file using `reno` if this change needs
to be documented in the release notes. (no release note required)

qiskit_experiments/library/randomized_benchmarking/clifford_utils.py

@@ -17,10 +17,9 @@
 import os
 from functools import lru_cache
 from numbers import Integral
-from typing import Optional, Union, Tuple, Sequence
+from typing import Optional, Union, Tuple, Sequence, Iterable
 
 import numpy as np
-import scipy.sparse
 from numpy.random import Generator, default_rng
 
 from qiskit.circuit import CircuitInstruction, Qubit
@@ -36,11 +35,13 @@
 
 _CLIFFORD_COMPOSE_1Q = np.load(f"{_DATA_FOLDER}/clifford_compose_1q.npz")["table"]
 _CLIFFORD_INVERSE_1Q = np.load(f"{_DATA_FOLDER}/clifford_inverse_1q.npz")["table"]
-_CLIFFORD_COMPOSE_2Q = scipy.sparse.lil_matrix(
-    scipy.sparse.load_npz(f"{_DATA_FOLDER}/clifford_compose_2q_sparse.npz")
-)
 _CLIFFORD_INVERSE_2Q = np.load(f"{_DATA_FOLDER}/clifford_inverse_2q.npz")["table"]
-
+_clifford_compose_2q_data = np.load(f"{_DATA_FOLDER}/clifford_compose_2q_dense_selected.npz")
+_CLIFFORD_COMPOSE_2Q_DENSE = _clifford_compose_2q_data["table"]
+# valid indices for the columns of the _CLIFFORD_COMPOSE_2Q_DENSE table
+_valid_sparse_indices = _clifford_compose_2q_data["valid_sparse_indices"]
+# map a clifford number to the index of _CLIFFORD_COMPOSE_2Q_DENSE
+_clifford_num_to_dense_index = {idx: ii for ii, idx in enumerate(_valid_sparse_indices)}
 
 # Transpilation utilities
 def _transpile_clifford_circuit(
@@ -422,8 +423,9 @@ def compose_2q(lhs: Integral, rhs: Integral) -> Integral:
     """Return the composition of 2-qubit clifford integers."""
     num = lhs
     for layer, idx in enumerate(_layer_indices_from_num(rhs)):
-        circ = _CLIFFORD_LAYER[layer][idx]
-        num = _compose_num_with_circuit_2q(num, circ)
+        gate_numbers = _CLIFFORD_LAYER_NUMS[layer][idx]
+        for n in gate_numbers:
+            num = _CLIFFORD_COMPOSE_2Q_DENSE[num, _clifford_num_to_dense_index[n]]
     return num
 
 
@@ -434,17 +436,9 @@ def inverse_2q(num: Integral) -> Integral:
 
 def num_from_2q_circuit(qc: QuantumCircuit) -> Integral:
     """Convert a given 2-qubit Clifford circuit to the corresponding integer."""
-    return _compose_num_with_circuit_2q(0, qc)
-
-
-def _compose_num_with_circuit_2q(num: Integral, qc: QuantumCircuit) -> Integral:
-    """Compose a number that represents a Clifford, with a Clifford circuit, and return the
-    number that represents the resulting Clifford."""
-    lhs = num
-    for inst in qc:
-        qubits = tuple(qc.find_bit(q).index for q in inst.qubits)
-        rhs = _num_from_2q_gate(op=inst.operation, qubits=qubits)
-        lhs = _CLIFFORD_COMPOSE_2Q[lhs, rhs]
+    lhs = 0
+    for rhs in _clifford_2q_nums_from_2q_circuit(qc):
+        lhs = _CLIFFORD_COMPOSE_2Q_DENSE[lhs, _clifford_num_to_dense_index[rhs]]
     return lhs
 
 
@@ -570,6 +564,20 @@ def _create_cliff_2q_layer_2():
 _NUM_LAYER_2 = 16
 
 
+def _clifford_2q_nums_from_2q_circuit(qc: QuantumCircuit) -> Iterable[Integral]:
+    """Yield Clifford numbers that represents the 2Q Clifford circuit."""
+    for inst in qc:
+        qubits = tuple(qc.find_bit(q).index for q in inst.qubits)
+        yield _num_from_2q_gate(op=inst.operation, qubits=qubits)
+
+
+# Construct mapping from Clifford layers to series of Clifford numbers
+_CLIFFORD_LAYER_NUMS = [
+    [tuple(_clifford_2q_nums_from_2q_circuit(qc)) for qc in _CLIFFORD_LAYER[layer]]
+    for layer in [0, 1, 2]
+]
+
+
 @lru_cache(maxsize=None)
 def _transformed_clifford_layer(
     layer: int, index: Integral, basis_gates: Tuple[str, ...]

qiskit_experiments/library/randomized_benchmarking/data/generate_clifford_data.py

@@ -109,7 +109,29 @@ def gen_clifford_compose_2q_gate():
         for rhs in _CLIFF_SINGLE_GATE_MAP_2Q.values():
             composed = cliff_lhs.compose(_CLIFF_2Q[rhs])
             products[lhs, rhs] = _TO_INT_2Q[_hash_cliff(composed)]
-    return products.tocsr()
+    return products.tocsc()
+
+
+def gen_clifford_compose_2q_dense() -> tuple[np.typing.NDArray[int], list[int]]:
+    """Generate a dense multiplication table for 2-qubit Clifford numbers
+
+    The multiplication table is generated from the sparse table generated by :meth:`gen_clifford_compose_2q_gate`.
+    Each column contains the full set of Clifford numbers . Each row contains a subset of Clifford multiplications
+    corresponding to the values of entries of `_CLIFF_SINGLE_GATE_MAP_2Q`.
+
+    Returns:
+        Tuple with a dense multiplication table and the valid indices for the columns
+    """
+    _CLIFFORD_COMPOSE_2Q = gen_clifford_compose_2q_gate()
+    number_of_cliffords = _CLIFFORD_COMPOSE_2Q.shape[0]
+    valid_sparse_indices = [
+        num
+        for num in range(number_of_cliffords)
+        if _CLIFFORD_COMPOSE_2Q[:, num].nnz == number_of_cliffords - 1
+    ]
+    _CLIFFORD_COMPOSE_2Q_DENSE = _CLIFFORD_COMPOSE_2Q[:, valid_sparse_indices].toarray()
+
+    return _CLIFFORD_COMPOSE_2Q_DENSE, valid_sparse_indices
 
 
 _GATE_LIST_1Q = [
@@ -184,4 +206,11 @@ def gen_cliff_single_2q_gate_map():
             "_CLIFF_SINGLE_GATE_MAP_2Q must be generated by gen_cliff_single_2q_gate_map()"
         )
     np.savez_compressed("clifford_inverse_2q.npz", table=gen_clifford_inverse_2q())
-    scipy.sparse.save_npz("clifford_compose_2q_sparse.npz", gen_clifford_compose_2q_gate())
+
+    _CLIFFORD_COMPOSE_2Q_DENSE, valid_sparse_indices = gen_clifford_compose_2q_dense()
+
+    np.savez_compressed(
+        "clifford_compose_2q_dense_selected.npz",
+        table=_CLIFFORD_COMPOSE_2Q_DENSE,
+        valid_sparse_indices=valid_sparse_indices,
+    )

test/library/randomized_benchmarking/test_clifford_utils.py

@@ -20,6 +20,7 @@
 from numpy.random import default_rng
 
 from qiskit import QuantumCircuit
+from qiskit.exceptions import QiskitError
 from qiskit.circuit.library import (
     IGate,
     XGate,
@@ -43,6 +44,7 @@
     _num_from_layer_indices,
     _layer_indices_from_num,
     _CLIFFORD_LAYER,
+    _CLIFFORD_INVERSE_2Q,
 )
 
 
@@ -195,3 +197,27 @@ def test_num_from_layer(self):
                 circ.compose(_CLIFFORD_LAYER[layer][idx], inplace=True)
             layered = Clifford(circ)
             self.assertEqual(standard, layered)
+
+    def test_num_from_2q_circuit(self):
+        """Check conversion of circuits to integers with num_from_2q_circuit"""
+        qc = QuantumCircuit(2)
+        qc.h(0)
+        num = num_from_2q_circuit(qc)
+        self.assertEqual(num, 5760)
+        qc = QuantumCircuit(2)
+        qc.u(0, 0, np.pi, 0)
+        with self.assertRaises(QiskitError):
+            # raising an error is ok, num_from_2q_circuit does not support all 2-qubit gates
+            num_from_2q_circuit(qc)
+
+        # regression test for using the dense multiplication table
+        qc = QuantumCircuit(2)
+        qc.cz(1, 0)
+        num = num_from_2q_circuit(qc)
+        self.assertEqual(num, 368)
+
+    def test_clifford_inverse_table(self):
+        """Check correctness of the Clifford inversion table"""
+        for lhs, rhs in enumerate(_CLIFFORD_INVERSE_2Q):
+            c = compose_2q(lhs, rhs)
+            self.assertEqual(c, 0)