Add Rust quantum volume function (Qiskit#13238)

* Add Rust quantum volume function

This commit adds a new function quantum_volume used for generating a
quantum volume model circuit. This new function is defined in Rust and
multithreaded to improve the throughput of the circuit generation. This
new function will eventually replace the existing QuantumVolume class as
part of Qiskit#13046. Since quantum volume is a circuit defined by it's
structure using a generator function is inline with the goals of Qiskit#13046.

Right now the performance is bottlenecked by the creation of the
UnitaryGate objects as these are still defined solely in Python. We'll
likely need to port the class to have a rust native representation to
further speed up the construction of the circuit.

* Adjust type hints on python function

Co-authored-by: Julien Gacon <[email protected]>

* Add missing __future__ import

The previous commit was relying on the behavior of the annotations
future import but neglected to add it. This commit corrects the
oversight.

* Add comment on random unitary algorithm

* Reduce allocations random_unitaries

The previous implementation had 4 heap allocations for each random
unitary constructed, this commit uses some fixed sized stack allocated
arrays and reduces that to two allocations one for q and r from the
factorization. We'll always need at least one for the `Array2` that gets
stored in each `UnitaryGate` as a numpy array. But to reduce to just
this we'll need a method of computing the QR factorization without an
allocation for the result space, nalgebtra might be a path for doing
that. While this currently isn't a bottleneck as the `UnitaryGate`
python object creation is the largest source of runtime, but assuming
that's fixed in the future this might have a larger impact.

* Preallocate unitaries for serial path

When executing in the serial path we previously were working directly
with an iterator where the 2q unitaries we're created on the iterator
that we passed directly to circuit constructor. However testing shows
that precomputing all the unitaries into a Vec and passing the iterator
off of that to the circuit constructor is marginally but consistently
faster. So this commit pivots to using that instead.

* Fix determinism and error handling of of qv function

This commit fixes two issues in the reproducibility of the quantum
volume circuit. The first was the output unitary matrices for a fixed
seed would differ between the parallel and serial execution path. This
was because how the RNGs were used was different in the different code
paths. This change results in the serial path being marginally less
efficient, but it shouldn't be a big deal when compared to getting
different results in different contexts. The second was the seed usage
in parallel mode was dependent on the number of threads on the local
system. This was problematic because the exact circuit generated between
two systems would be different even with a fixed seed. This was fixed to
avoid depending on the number of threads to determine how the seeds were
used across multiple threads.

The last fix here was a change to the error handling so that the
CircuitData constructor used to create the circuit object can handle a
fallible iterator. Previously we were throwing away the python error and
panicking if the Python call to generate the UnitaryGate object raised
an error for any reason.

* Mention the new function is multithreaded in docstring

* Update qiskit/circuit/library/quantum_volume.py

Co-authored-by: Julien Gacon <[email protected]>

---------

Co-authored-by: Julien Gacon <[email protected]>
Co-authored-by: Julien Gacon <[email protected]>

crates/accelerate/src/circuit_library/mod.rs

@@ -14,9 +14,11 @@ use pyo3::prelude::*;
 
 mod entanglement;
 mod pauli_feature_map;
+mod quantum_volume;
 
 pub fn circuit_library(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(pauli_feature_map::pauli_feature_map))?;
     m.add_wrapped(wrap_pyfunction!(entanglement::get_entangler_map))?;
+    m.add_wrapped(wrap_pyfunction!(quantum_volume::quantum_volume))?;
     Ok(())
 }

crates/accelerate/src/circuit_library/pauli_feature_map.rs

@@ -138,7 +138,7 @@ fn pauli_evolution(
 ///     insert_barriers: Whether to insert barriers in between the Hadamard and evolution layers.
 ///     data_map_func: An accumulation function that takes as input a vector of parameters the
 ///         current gate acts on and returns a scalar.
-///     
+///
 /// Returns:
 ///     The ``CircuitData`` to construct the Pauli feature map.
 #[pyfunction]
@@ -207,7 +207,13 @@ pub fn pauli_feature_map(
         }
     }
 
-    CircuitData::from_packed_operations(py, feature_dimension, 0, packed_insts, Param::Float(0.0))
+    CircuitData::from_packed_operations(
+        py,
+        feature_dimension,
+        0,
+        packed_insts.into_iter().map(Ok),
+        Param::Float(0.0),
+    )
 }
 
 fn _get_h_layer(feature_dimension: u32) -> impl Iterator<Item = Instruction> {

crates/accelerate/src/circuit_library/quantum_volume.rs

@@ -0,0 +1,176 @@
+// This code is part of Qiskit.
+//
+// (C) Copyright IBM 2024
+//
+// 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.
+
+use pyo3::prelude::*;
+
+use crate::getenv_use_multiple_threads;
+use faer_ext::{IntoFaerComplex, IntoNdarrayComplex};
+use ndarray::prelude::*;
+use num_complex::Complex64;
+use numpy::IntoPyArray;
+use rand::prelude::*;
+use rand_distr::StandardNormal;
+use rand_pcg::Pcg64Mcg;
+use rayon::prelude::*;
+
+use qiskit_circuit::circuit_data::CircuitData;
+use qiskit_circuit::imports::UNITARY_GATE;
+use qiskit_circuit::operations::Param;
+use qiskit_circuit::operations::PyInstruction;
+use qiskit_circuit::packed_instruction::PackedOperation;
+use qiskit_circuit::{Clbit, Qubit};
+use smallvec::{smallvec, SmallVec};
+
+type Instruction = (
+    PackedOperation,
+    SmallVec<[Param; 3]>,
+    Vec<Qubit>,
+    Vec<Clbit>,
+);
+
+#[inline(always)]
+fn random_complex(rng: &mut Pcg64Mcg) -> Complex64 {
+    Complex64::new(rng.sample(StandardNormal), rng.sample(StandardNormal))
+        * std::f64::consts::FRAC_1_SQRT_2
+}
+
+// This function's implementation was modeled off of the algorithm used in the
+// `scipy.stats.unitary_group.rvs()` function defined here:
+//
+// https://github.com/scipy/scipy/blob/v1.14.1/scipy/stats/_multivariate.py#L4224-L4256
+#[inline]
+fn random_unitaries(seed: u64, size: usize) -> impl Iterator<Item = Array2<Complex64>> {
+    let mut rng = Pcg64Mcg::seed_from_u64(seed);
+
+    (0..size).map(move |_| {
+        let raw_numbers: [[Complex64; 4]; 4] = [
+            [
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+            ],
+            [
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+            ],
+            [
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+            ],
+            [
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+                random_complex(&mut rng),
+            ],
+        ];
+
+        let qr = aview2(&raw_numbers).into_faer_complex().qr();
+        let r = qr.compute_r();
+        let diag: [Complex64; 4] = [
+            r[(0, 0)].to_num_complex() / r[(0, 0)].abs(),
+            r[(1, 1)].to_num_complex() / r[(1, 1)].abs(),
+            r[(2, 2)].to_num_complex() / r[(2, 2)].abs(),
+            r[(3, 3)].to_num_complex() / r[(3, 3)].abs(),
+        ];
+        let mut q = qr.compute_q().as_ref().into_ndarray_complex().to_owned();
+        q.axis_iter_mut(Axis(0)).for_each(|mut row| {
+            row.iter_mut()
+                .enumerate()
+                .for_each(|(index, val)| *val *= diag[index])
+        });
+        q
+    })
+}
+
+const UNITARY_PER_SEED: usize = 50;
+
+#[pyfunction]
+pub fn quantum_volume(
+    py: Python,
+    num_qubits: u32,
+    depth: usize,
+    seed: Option<u64>,
+) -> PyResult<CircuitData> {
+    let width = num_qubits as usize / 2;
+    let num_unitaries = width * depth;
+    let mut permutation: Vec<Qubit> = (0..num_qubits).map(Qubit).collect();
+
+    let mut build_instruction = |(unitary_index, unitary_array): (usize, Array2<Complex64>),
+                                 rng: &mut Pcg64Mcg|
+     -> PyResult<Instruction> {
+        let layer_index = unitary_index % width;
+        if layer_index == 0 {
+            permutation.shuffle(rng);
+        }
+        let unitary = unitary_array.into_pyarray_bound(py);
+        let unitary_gate = UNITARY_GATE
+            .get_bound(py)
+            .call1((unitary.clone(), py.None(), false))?;
+        let instruction = PyInstruction {
+            qubits: 2,
+            clbits: 0,
+            params: 1,
+            op_name: "unitary".to_string(),
+            control_flow: false,
+            instruction: unitary_gate.unbind(),
+        };
+        let qubit = layer_index * 2;
+        Ok((
+            PackedOperation::from_instruction(Box::new(instruction)),
+            smallvec![Param::Obj(unitary.unbind().into())],
+            vec![permutation[qubit], permutation[qubit + 1]],
+            vec![],
+        ))
+    };
+
+    let mut per_thread = num_unitaries / UNITARY_PER_SEED;
+    if per_thread == 0 {
+        per_thread = 10;
+    }
+    let mut outer_rng = match seed {
+        Some(seed) => Pcg64Mcg::seed_from_u64(seed),
+        None => Pcg64Mcg::from_entropy(),
+    };
+    let seed_vec: Vec<u64> = rand::distributions::Standard
+        .sample_iter(&mut outer_rng)
+        .take(num_unitaries)
+        .collect();
+
+    let unitaries: Vec<Array2<Complex64>> = if getenv_use_multiple_threads() && num_unitaries > 200
+    {
+        seed_vec
+            .par_chunks(per_thread)
+            .flat_map_iter(|seeds| random_unitaries(seeds[0], seeds.len()))
+            .collect()
+    } else {
+        seed_vec
+            .chunks(per_thread)
+            .flat_map(|seeds| random_unitaries(seeds[0], seeds.len()))
+            .collect()
+    };
+    CircuitData::from_packed_operations(
+        py,
+        num_qubits,
+        0,
+        unitaries
+            .into_iter()
+            .enumerate()
+            .map(|x| build_instruction(x, &mut outer_rng)),
+        Param::Float(0.),
+    )
+}

crates/accelerate/src/two_qubit_decompose.rs

@@ -2176,18 +2176,18 @@ impl TwoQubitBasisDecomposer {
                     .gates
                     .into_iter()
                     .map(|(gate, params, qubits)| match gate {
-                        Some(gate) => (
+                        Some(gate) => Ok((
                             PackedOperation::from_standard(gate),
                             params.into_iter().map(Param::Float).collect(),
                             qubits.into_iter().map(|x| Qubit(x.into())).collect(),
                             Vec::new(),
-                        ),
-                        None => (
+                        )),
+                        None => Ok((
                             kak_gate.operation.clone(),
                             kak_gate.params.clone(),
                             qubits.into_iter().map(|x| Qubit(x.into())).collect(),
                             Vec::new(),
-                        ),
+                        )),
                     }),
                 Param::Float(sequence.global_phase),
             ),

crates/circuit/src/circuit_data.rs

@@ -134,12 +134,12 @@ impl CircuitData {
     ) -> PyResult<Self>
     where
         I: IntoIterator<
-            Item = (
+            Item = PyResult<(
                 PackedOperation,
                 SmallVec<[Param; 3]>,
                 Vec<Qubit>,
                 Vec<Clbit>,
-            ),
+            )>,
         >,
     {
         let instruction_iter = instructions.into_iter();
@@ -150,7 +150,8 @@ impl CircuitData {
             instruction_iter.size_hint().0,
             global_phase,
         )?;
-        for (operation, params, qargs, cargs) in instruction_iter {
+        for item in instruction_iter {
+            let (operation, params, qargs, cargs) = item?;
             let qubits = res.qargs_interner.insert_owned(qargs);
             let clbits = res.cargs_interner.insert_owned(cargs);
             let params = (!params.is_empty()).then(|| Box::new(params));

qiskit/circuit/library/__init__.py

@@ -318,6 +318,7 @@
    HiddenLinearFunction
    IQP
    QuantumVolume
+   quantum_volume
    PhaseEstimation
    GroverOperator
    PhaseOracle
@@ -564,7 +565,7 @@
     StatePreparation,
     Initialize,
 )
-from .quantum_volume import QuantumVolume
+from .quantum_volume import QuantumVolume, quantum_volume
 from .fourier_checking import FourierChecking
 from .graph_state import GraphState
 from .hidden_linear_function import HiddenLinearFunction

qiskit/circuit/library/quantum_volume.py

@@ -12,11 +12,14 @@
 
 """Quantum Volume model circuit."""
 
+from __future__ import annotations
+
 from typing import Optional, Union
 
 import numpy as np
 from qiskit.circuit import QuantumCircuit, CircuitInstruction
 from qiskit.circuit.library.generalized_gates import PermutationGate, UnitaryGate
+from qiskit._accelerate.circuit_library import quantum_volume as qv_rs
 
 
 class QuantumVolume(QuantumCircuit):
@@ -113,3 +116,42 @@ def __init__(
                     base._append(CircuitInstruction(gate, qubits[qubit : qubit + 2]))
         if not flatten:
             self._append(CircuitInstruction(base.to_instruction(), tuple(self.qubits)))
+
+
+def quantum_volume(
+    num_qubits: int,
+    depth: int | None = None,
+    seed: int | np.random.Generator | None = None,
+) -> QuantumCircuit:
+    """A quantum volume model circuit.
+
+    The model circuits are random instances of circuits used to measure
+    the Quantum Volume metric, as introduced in [1].
+
+    The model circuits consist of layers of Haar random
+    elements of SU(4) applied between corresponding pairs
+    of qubits in a random bipartition.
+
+    This function is multithreaded and will launch a thread pool with threads equal to the number
+    of CPUs by default. You can tune the number of threads with the ``RAYON_NUM_THREADS``
+    environment variable. For example, setting ``RAYON_NUM_THREADS=4`` would limit the thread pool
+    to 4 threads.
+
+    **Reference Circuit:**
+
+    .. plot::
+
+       from qiskit.circuit.library import quantum_volume
+       circuit = quantum_volume(5, 6, seed=10)
+       circuit.draw('mpl')
+
+    **References:**
+
+    [1] A. Cross et al. Validating quantum computers using
+    randomized model circuits, Phys. Rev. A 100, 032328 (2019).
+    `arXiv:1811.12926 <https://arxiv.org/abs/1811.12926>`__
+    """
+    if isinstance(seed, np.random.Generator):
+        seed = seed.integers(0, dtype=np.uint64)
+    depth = depth or num_qubits
+    return QuantumCircuit._from_circuit_data(qv_rs(num_qubits, depth, seed))

releasenotes/notes/add-qv-function-a8990e248d5e7e1a.yaml

@@ -0,0 +1,11 @@
+---
+features_circuits:
+  - |
+    Added a new function :func:`.quantum_volume` for generating a quantum volume
+    :class:`.QuantumCircuit` object as defined in A. Cross et al. Validating quantum computers
+    using randomized model circuits, Phys. Rev. A 100, 032328 (2019)
+    `https://link.aps.org/doi/10.1103/PhysRevA.100.032328 <https://link.aps.org/doi/10.1103/PhysRevA.100.032328>`__.
+    This new function differs from the existing :class:`.QuantumVolume` class in that it returns
+    a :class:`.QuantumCircuit` object instead of building a subclass object. The second is
+    that this new function is multithreaded and implemented in rust so it generates the output
+    circuit ~10x faster than the :class:`.QuantumVolume` class.

test/python/circuit/library/test_quantum_volume.py

@@ -16,6 +16,7 @@
 
 from test.utils.base import QiskitTestCase
 from qiskit.circuit.library import QuantumVolume
+from qiskit.circuit.library.quantum_volume import quantum_volume
 
 
 class TestQuantumVolumeLibrary(QiskitTestCase):
@@ -35,6 +36,20 @@ def test_qv_seed_reproducibility(self):
         right = QuantumVolume(4, 4, seed=2024, flatten=True)
         self.assertEqual(left, right)
 
+    def test_qv_function_seed_reproducibility(self):
+        """Test qv circuit."""
+        left = quantum_volume(10, 10, seed=128)
+        right = quantum_volume(10, 10, seed=128)
+        self.assertEqual(left, right)
+
+        left = quantum_volume(10, 10, seed=256)
+        right = quantum_volume(10, 10, seed=256)
+        self.assertEqual(left, right)
+
+        left = quantum_volume(10, 10, seed=4196)
+        right = quantum_volume(10, 10, seed=4196)
+        self.assertEqual(left, right)
+
 
 if __name__ == "__main__":
     unittest.main()