Add Rust quantum volume function

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/quantum_volume.rs

@@ -0,0 +1,185 @@
+// 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 std::thread::available_parallelism;
+
+use qiskit_circuit::imports::UNITARY_GATE;
+
+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::operations::Param;
+use qiskit_circuit::operations::PyInstruction;
+use qiskit_circuit::packed_instruction::PackedOperation;
+use qiskit_circuit::Qubit;
+use smallvec::smallvec;
+
+#[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
+    })
+}
+
+#[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| {
+        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))
+            .unwrap();
+        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;
+        (
+            PackedOperation::from_instruction(Box::new(instruction)),
+            smallvec![Param::Obj(unitary.unbind().into())],
+            vec![permutation[qubit], permutation[qubit + 1]],
+            vec![],
+        )
+    };
+
+    if getenv_use_multiple_threads() {
+        let mut per_thread = num_unitaries / available_parallelism().unwrap();
+        if per_thread == 0 {
+            if num_unitaries > 10 {
+                per_thread = 10
+            } else {
+                per_thread = num_unitaries
+            }
+        }
+
+        let mut outer_rng = match seed {
+            Some(seed) => Pcg64Mcg::seed_from_u64(seed),
+            None => Pcg64Mcg::from_entropy(),
+        };
+        let seed_vec: Vec<u64> = outer_rng
+            .clone()
+            .sample_iter(&rand::distributions::Standard)
+            .take(num_unitaries)
+            .collect();
+        let unitaries: Vec<Array2<Complex64>> = seed_vec
+            .into_par_iter()
+            .chunks(per_thread)
+            .flat_map_iter(|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.),
+        )
+    } else {
+        let mut outer_rng = match seed {
+            Some(seed) => Pcg64Mcg::seed_from_u64(seed),
+            None => Pcg64Mcg::from_entropy(),
+        };
+        let seed: u64 = outer_rng.sample(rand::distributions::Standard);
+        CircuitData::from_packed_operations(
+            py,
+            num_qubits,
+            0,
+            random_unitaries(seed, num_unitaries)
+                .enumerate()
+                .map(|x| build_instruction(x, &mut outer_rng)),
+            Param::Float(0.),
+        )
+    }
+}

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,37 @@ 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.
+
+    **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()