Improve performance and randomness of QuantumVolume (Qiskit#12097)

* Improve performance and randomness of `QuantumVolume`

This hugely improves the construction time of `QuantumVolume` circuits,
in part by removing the previous behaviour of using low-entropy
indiviual RNG instances for each SU4 matrix.  Now that we need larger
circuits, this would already have been a non-trivial biasing of the
random outputs, but also, calling Scipy random variates in a loop is
_much_ less efficient than vectorising the entire process in one go.

Along with changes to the RNG, this commit also adds a faster path to
`UnitaryGate` to shave off some useless repeated calculations (this can
likely be used elsewhere in Qiskit) and reworks the `QuantumVolume`
builder to use more efficient circuit construction methods.

* Make test explicitly test reproducibility

* Protect best-effort seed retrieval against old Numpy

qiskit/circuit/library/generalized_gates/unitary.py

@@ -70,6 +70,8 @@ def __init__(
         data: numpy.ndarray | Gate | BaseOperator,
         label: str | None = None,
         check_input: bool = True,
+        *,
+        num_qubits: int | None = None,
     ) -> None:
         """Create a gate from a numeric unitary matrix.
 
@@ -81,6 +83,7 @@ def __init__(
                 be skipped. This should only ever be used if you know the
                 input is unitary, setting this to ``False`` and passing in
                 a non-unitary matrix will result unexpected behavior and errors.
+            num_qubits: If given, the number of qubits in the matrix.  If not given, it is inferred.
 
         Raises:
             ValueError: If input data is not an N-qubit unitary operator.
@@ -97,7 +100,7 @@ def __init__(
         # Convert to numpy array in case not already an array
         data = numpy.asarray(data, dtype=complex)
         input_dim, output_dim = data.shape
-        num_qubits = int(math.log2(input_dim))
+        num_qubits = num_qubits if num_qubits is not None else int(math.log2(input_dim))
         if check_input:
             # Check input is unitary
             if not is_unitary_matrix(data):

qiskit/circuit/library/quantum_volume.py

@@ -15,9 +15,8 @@
 from typing import Optional, Union
 
 import numpy as np
-from qiskit.quantum_info.random import random_unitary
-from qiskit.circuit import QuantumCircuit
-from qiskit.circuit.library.generalized_gates.permutation import Permutation
+from qiskit.circuit import QuantumCircuit, CircuitInstruction
+from qiskit.circuit.library.generalized_gates import PermutationGate, UnitaryGate
 
 
 class QuantumVolume(QuantumCircuit):
@@ -60,6 +59,8 @@ def __init__(
         depth: Optional[int] = None,
         seed: Optional[Union[int, np.random.Generator]] = None,
         classical_permutation: bool = True,
+        *,
+        flatten: bool = False,
     ) -> None:
         """Create quantum volume model circuit of size num_qubits x depth.
 
@@ -69,46 +70,46 @@ def __init__(
             seed: Random number generator or generator seed.
             classical_permutation: use classical permutations at every layer,
                 rather than quantum.
+            flatten: If ``False`` (the default), construct a circuit that contains a single
+                instruction, which in turn has the actual volume structure.  If ``True``, construct
+                the volume structure directly.
         """
-        # Initialize RNG
-        if seed is None:
-            rng_set = np.random.default_rng()
-            seed = rng_set.integers(low=1, high=1000)
-        if isinstance(seed, np.random.Generator):
-            rng = seed
-        else:
-            rng = np.random.default_rng(seed)
+        import scipy.stats
 
         # Parameters
         depth = depth or num_qubits  # how many layers of SU(4)
-        width = int(np.floor(num_qubits / 2))  # how many SU(4)s fit in each layer
-        name = "quantum_volume_" + str([num_qubits, depth, seed]).replace(" ", "")
+        width = num_qubits // 2  # how many SU(4)s fit in each layer
+        rng = seed if isinstance(seed, np.random.Generator) else np.random.default_rng(seed)
+        if seed is None:
+            # Get the internal entropy used to seed the default RNG, if no seed was given.  This
+            # stays in the output name, so effectively stores a way of regenerating the circuit.
+            # This is just best-effort only, for backwards compatibility, and isn't critical (if
+            # someone needs full reproducibility, they should be manually controlling the seeding).
+            seed = getattr(getattr(rng.bit_generator, "seed_seq", None), "entropy", None)
 
-        # Generator random unitary seeds in advance.
-        # Note that this means we are constructing multiple new generator
-        # objects from low-entropy integer seeds rather than pass the shared
-        # generator object to the random_unitary function. This is done so
-        # that we can use the integer seed as a label for the generated gates.
-        unitary_seeds = rng.integers(low=1, high=1000, size=[depth, width])
+        super().__init__(
+            num_qubits, name="quantum_volume_" + str([num_qubits, depth, seed]).replace(" ", "")
+        )
+        base = self if flatten else QuantumCircuit(num_qubits, name=self.name)
 
         # For each layer, generate a permutation of qubits
         # Then generate and apply a Haar-random SU(4) to each pair
-        circuit = QuantumCircuit(num_qubits, name=name)
-        perm_0 = list(range(num_qubits))
-        for d in range(depth):
-            perm = rng.permutation(perm_0)
-            if not classical_permutation:
-                layer_perm = Permutation(num_qubits, perm)
-                circuit.compose(layer_perm, inplace=True)
-            for w in range(width):
-                seed_u = unitary_seeds[d][w]
-                su4 = random_unitary(4, seed=seed_u).to_instruction()
-                su4.label = "su4_" + str(seed_u)
-                if classical_permutation:
-                    physical_qubits = int(perm[2 * w]), int(perm[2 * w + 1])
-                    circuit.compose(su4, [physical_qubits[0], physical_qubits[1]], inplace=True)
-                else:
-                    circuit.compose(su4, [2 * w, 2 * w + 1], inplace=True)
-
-        super().__init__(*circuit.qregs, name=circuit.name)
-        self.compose(circuit.to_instruction(), qubits=self.qubits, inplace=True)
+        unitaries = scipy.stats.unitary_group.rvs(4, depth * width, rng).reshape(depth, width, 4, 4)
+        qubits = tuple(base.qubits)
+        for row in unitaries:
+            perm = rng.permutation(num_qubits)
+            if classical_permutation:
+                for w, unitary in enumerate(row):
+                    gate = UnitaryGate(unitary, check_input=False, num_qubits=2)
+                    qubit = 2 * w
+                    base._append(
+                        CircuitInstruction(gate, (qubits[perm[qubit]], qubits[perm[qubit + 1]]))
+                    )
+            else:
+                base._append(CircuitInstruction(PermutationGate(perm), qubits))
+                for w, unitary in enumerate(row):
+                    gate = UnitaryGate(unitary, check_input=False, num_qubits=2)
+                    qubit = 2 * w
+                    base._append(CircuitInstruction(gate, qubits[qubit : qubit + 2]))
+        if not flatten:
+            self._append(CircuitInstruction(base.to_instruction(), tuple(self.qubits)))

releasenotes/notes/qv-perf-be76290f472e4777.yaml

@@ -0,0 +1,27 @@
+---
+features_circuits:
+  - |
+    Construction time for :class:`.QuantumVolume` circuits has been significantly improved, on the
+    order of 10x or a bit more.  The internal SU4 gates will now also use more bits of randomness
+    during their generation, leading to more representative volume circuits, especially at large
+    widths and depths.
+  - |
+    :class:`.QuantumVolume` now has a ``flatten`` keyword argument.  This defaults to ``False``,
+    where the constructed circuit contains a single instruction that in turn contains the actual
+    volume structure.  If set ``True``, the circuit will directly have the volumetric SU4 matrices.
+  - |
+    :class:`.UnitaryGate` now accepts an optional ``num_qubits`` argument.  The only effect of this
+    is to skip the inference of the qubit count, which can be helpful for performance when many
+    gates are being constructed.
+upgrade_circuits:
+  - |
+    The random-number usage of :class:`.QuantumVolume` has changed, so you will get a different
+    circuit for a fixed seed between older versions of Qiskit and this version.  The random-unitary
+    generation now uses more bits of entropy, so large circuits will be less biased.
+  - |
+    The internal :class:`.UnitaryGate` instances in the definition of a :class:`.QuantumVolume`
+    circuit will no longer have a :attr:`~.Instruction.label` field set.  Previously this was set
+    to the string ``su4_<seed>`` where ``<seed>`` was a three-digit number denoting the seed of an
+    internal Numpy pRNG instance for that gate.  Doing this was a serious performance problem, and
+    the seed ought not to have been useful; if you need to retrieve the matrix from the gate, simply
+    use the :meth:`.Gate.to_matrix` method.

test/python/circuit/library/test_quantum_volume.py

@@ -15,25 +15,25 @@
 import unittest
 
 from test.utils.base import QiskitTestCase
-from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.library import QuantumVolume
-from qiskit.quantum_info import Operator
-from qiskit.quantum_info.random import random_unitary
 
 
 class TestQuantumVolumeLibrary(QiskitTestCase):
     """Test library of quantum volume quantum circuits."""
 
-    def test_qv(self):
+    def test_qv_seed_reproducibility(self):
         """Test qv circuit."""
-        circuit = QuantumVolume(2, 2, seed=2, classical_permutation=False)
-        expected = QuantumCircuit(2)
-        expected.swap(0, 1)
-        expected.append(random_unitary(4, seed=837), [0, 1])
-        expected.append(random_unitary(4, seed=262), [0, 1])
-        expected = Operator(expected)
-        simulated = Operator(circuit)
-        self.assertTrue(expected.equiv(simulated))
+        left = QuantumVolume(4, 4, seed=28, classical_permutation=False)
+        right = QuantumVolume(4, 4, seed=28, classical_permutation=False)
+        self.assertEqual(left, right)
+
+        left = QuantumVolume(4, 4, seed=3, classical_permutation=True)
+        right = QuantumVolume(4, 4, seed=3, classical_permutation=True)
+        self.assertEqual(left, right)
+
+        left = QuantumVolume(4, 4, seed=2024, flatten=True)
+        right = QuantumVolume(4, 4, seed=2024, flatten=True)
+        self.assertEqual(left, right)
 
 
 if __name__ == "__main__":