Oxidize random clifford (#12695)

* starting to experiment

* porting code

* messy code porting

* printing statements to enable debugging

* fixes

* fixing phase

* removing some of the printing statements

* fixing inaccuracy for cost computation

* Moving some of the functionality to SymplecticMatrix class

* reducing the number of warnings

* formatting

* replacing expensive adjoint and compose operations for symplectic matrices by significantly cheaper in-place prepend operations

* resolving merge conflicts

* cleanup

* code cleanup

* cleanup

* cleanup

* cleanup

* cleanup

* cleanup

* cleanup

* cleanup

* using fast lookup

* cleanup

* clippy

* including params in gate_seq to avoid mapping

* removing unnecessary inner function

* cleanup

* renaming

* changes on the python side

* reno

* adding error handling

* improved error handling

* removing redundant Ok(Some(...))

* using random_clifford in tests

* reorganizing clifford code

* fixes

* formatting

* improved error handling

* do not panic

* formatting

* initial implementation

* cleanup

* changing return to be the full tableau

* parallelized computation for binary matrix multiplication for matrices with >= 10 qubits

* updating python code and adding release notes

* lint, test, reno

* removing obsolete python code

* removing todo (after experimenting with alternative implementation)

* removing another obsolete todo

* unused import

* adding comment

crates/accelerate/src/synthesis/clifford/mod.rs

@@ -12,12 +12,13 @@
 
 mod bm_synthesis;
 mod greedy_synthesis;
+mod random_clifford;
 mod utils;
 
 use crate::synthesis::clifford::bm_synthesis::synth_clifford_bm_inner;
 use crate::synthesis::clifford::greedy_synthesis::GreedyCliffordSynthesis;
 use crate::QiskitError;
-use numpy::PyReadonlyArray2;
+use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2};
 use pyo3::prelude::*;
 use qiskit_circuit::circuit_data::CircuitData;
 use qiskit_circuit::operations::Param;
@@ -43,6 +44,28 @@ fn synth_clifford_greedy(py: Python, clifford: PyReadonlyArray2<bool>) -> PyResu
     CircuitData::from_standard_gates(py, num_qubits as u32, clifford_gates, Param::Float(0.0))
 }
 
+/// Generate a random Clifford tableau.
+///
+/// The Clifford is sampled using the method of the paper "Hadamard-free circuits
+/// expose the structure of the Clifford group" by S. Bravyi and D. Maslov (2020),
+/// `https://arxiv.org/abs/2003.09412`__.
+///
+/// Args:
+///     num_qubits: the number of qubits.
+///     seed: an optional random seed.
+/// Returns:
+///     result: a random clifford tableau.
+#[pyfunction]
+#[pyo3(signature = (num_qubits, seed=None))]
+fn random_clifford_tableau(
+    py: Python,
+    num_qubits: usize,
+    seed: Option<u64>,
+) -> PyResult<Py<PyArray2<bool>>> {
+    let tableau = random_clifford::random_clifford_tableau_inner(num_qubits, seed);
+    Ok(tableau.into_pyarray_bound(py).unbind())
+}
+
 /// Create a circuit that optimally synthesizes a given Clifford operator represented as
 /// a tableau for Cliffords up to 3 qubits.
 ///
@@ -60,5 +83,6 @@ fn synth_clifford_bm(py: Python, clifford: PyReadonlyArray2<bool>) -> PyResult<C
 pub fn clifford(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(synth_clifford_greedy, m)?)?;
     m.add_function(wrap_pyfunction!(synth_clifford_bm, m)?)?;
+    m.add_function(wrap_pyfunction!(random_clifford_tableau, m)?)?;
     Ok(())
 }

crates/accelerate/src/synthesis/clifford/random_clifford.rs

@@ -0,0 +1,162 @@
+// 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 crate::synthesis::linear::utils::{
+    binary_matmul_inner, calc_inverse_matrix_inner, replace_row_inner, swap_rows_inner,
+};
+use ndarray::{concatenate, s, Array1, Array2, ArrayView2, ArrayViewMut2, Axis};
+use rand::{Rng, SeedableRng};
+use rand_pcg::Pcg64Mcg;
+
+/// Sample from the quantum Mallows distribution.
+fn sample_qmallows(n: usize, rng: &mut Pcg64Mcg) -> (Array1<bool>, Array1<usize>) {
+    // Hadamard layer
+    let mut had = Array1::from_elem(n, false);
+
+    // Permutation layer
+    let mut perm = Array1::from_elem(n, 0);
+    let mut inds: Vec<usize> = (0..n).collect();
+
+    for i in 0..n {
+        let m = n - i;
+        let eps: f64 = 4f64.powi(-(m as i32));
+        let r: f64 = rng.gen();
+        let index: usize = -((r + (1f64 - r) * eps).log2().ceil() as isize) as usize;
+        had[i] = index < m;
+        let k = if index < m { index } else { 2 * m - index - 1 };
+        perm[i] = inds[k];
+        inds.remove(k);
+    }
+    (had, perm)
+}
+
+/// Add symmetric random boolean value to off diagonal entries.
+fn fill_tril(mut mat: ArrayViewMut2<bool>, rng: &mut Pcg64Mcg, symmetric: bool) {
+    let n = mat.shape()[0];
+    for i in 0..n {
+        for j in 0..i {
+            mat[[i, j]] = rng.gen();
+            if symmetric {
+                mat[[j, i]] = mat[[i, j]];
+            }
+        }
+    }
+}
+
+/// Invert a lower-triangular matrix with unit diagonal.
+fn inverse_tril(mat: ArrayView2<bool>) -> Array2<bool> {
+    calc_inverse_matrix_inner(mat, false).unwrap()
+}
+
+/// Generate a random Clifford tableau.
+///
+/// The Clifford is sampled using the method of the paper "Hadamard-free circuits
+/// expose the structure of the Clifford group" by S. Bravyi and D. Maslov (2020),
+/// `https://arxiv.org/abs/2003.09412`__.
+///
+/// The function returns a random clifford tableau.
+pub fn random_clifford_tableau_inner(num_qubits: usize, seed: Option<u64>) -> Array2<bool> {
+    let mut rng = match seed {
+        Some(seed) => Pcg64Mcg::seed_from_u64(seed),
+        None => Pcg64Mcg::from_entropy(),
+    };
+
+    let (had, perm) = sample_qmallows(num_qubits, &mut rng);
+
+    let mut gamma1: Array2<bool> = Array2::from_elem((num_qubits, num_qubits), false);
+    for i in 0..num_qubits {
+        gamma1[[i, i]] = rng.gen();
+    }
+    fill_tril(gamma1.view_mut(), &mut rng, true);
+
+    let mut gamma2: Array2<bool> = Array2::from_elem((num_qubits, num_qubits), false);
+    for i in 0..num_qubits {
+        gamma2[[i, i]] = rng.gen();
+    }
+    fill_tril(gamma2.view_mut(), &mut rng, true);
+
+    let mut delta1: Array2<bool> = Array2::from_shape_fn((num_qubits, num_qubits), |(i, j)| i == j);
+    fill_tril(delta1.view_mut(), &mut rng, false);
+
+    let mut delta2: Array2<bool> = Array2::from_shape_fn((num_qubits, num_qubits), |(i, j)| i == j);
+    fill_tril(delta2.view_mut(), &mut rng, false);
+
+    // Compute stabilizer table
+    let zero = Array2::from_elem((num_qubits, num_qubits), false);
+    let prod1 = binary_matmul_inner(gamma1.view(), delta1.view()).unwrap();
+    let prod2 = binary_matmul_inner(gamma2.view(), delta2.view()).unwrap();
+    let inv1 = inverse_tril(delta1.view()).t().to_owned();
+    let inv2 = inverse_tril(delta2.view()).t().to_owned();
+
+    let table1 = concatenate(
+        Axis(0),
+        &[
+            concatenate(Axis(1), &[delta1.view(), zero.view()])
+                .unwrap()
+                .view(),
+            concatenate(Axis(1), &[prod1.view(), inv1.view()])
+                .unwrap()
+                .view(),
+        ],
+    )
+    .unwrap();
+
+    let table2 = concatenate(
+        Axis(0),
+        &[
+            concatenate(Axis(1), &[delta2.view(), zero.view()])
+                .unwrap()
+                .view(),
+            concatenate(Axis(1), &[prod2.view(), inv2.view()])
+                .unwrap()
+                .view(),
+        ],
+    )
+    .unwrap();
+
+    // Compute the full stabilizer tableau
+
+    // The code below is identical to the Python implementation, but is based on the original
+    // code in the paper.
+
+    let mut table = Array2::from_elem((2 * num_qubits, 2 * num_qubits), false);
+
+    // Apply qubit permutation
+    for i in 0..num_qubits {
+        replace_row_inner(table.view_mut(), i, table2.slice(s![i, ..]));
+        replace_row_inner(
+            table.view_mut(),
+            perm[i] + num_qubits,
+            table2.slice(s![perm[i] + num_qubits, ..]),
+        );
+    }
+
+    // Apply layer of Hadamards
+    for i in 0..num_qubits {
+        if had[i] {
+            swap_rows_inner(table.view_mut(), i, i + num_qubits);
+        }
+    }
+
+    // Apply table
+    let random_symplectic_mat = binary_matmul_inner(table1.view(), table.view()).unwrap();
+
+    // Generate random phases
+    let random_phases: Array2<bool> = Array2::from_shape_fn((2 * num_qubits, 1), |_| rng.gen());
+
+    let random_tableau: Array2<bool> = concatenate(
+        Axis(1),
+        &[random_symplectic_mat.view(), random_phases.view()],
+    )
+    .unwrap();
+    random_tableau
+}

crates/accelerate/src/synthesis/linear/utils.rs

@@ -10,9 +10,17 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
-use ndarray::{concatenate, s, Array2, ArrayView2, ArrayViewMut2, Axis};
+use ndarray::{azip, concatenate, s, Array2, ArrayView1, ArrayView2, ArrayViewMut2, Axis, Zip};
 use rand::{Rng, SeedableRng};
 use rand_pcg::Pcg64Mcg;
+use rayon::iter::{IndexedParallelIterator, ParallelIterator};
+use rayon::prelude::IntoParallelIterator;
+
+use crate::getenv_use_multiple_threads;
+
+/// Specifies the minimum number of qubits in order to parallelize computations
+/// (this number is chosen based on several local experiments).
+const PARALLEL_THRESHOLD: usize = 10;
 
 /// Binary matrix multiplication
 pub fn binary_matmul_inner(
@@ -30,11 +38,29 @@ pub fn binary_matmul_inner(
         ));
     }
 
-    Ok(Array2::from_shape_fn((n1_rows, n2_cols), |(i, j)| {
-        (0..n2_rows)
-            .map(|k| mat1[[i, k]] & mat2[[k, j]])
-            .fold(false, |acc, v| acc ^ v)
-    }))
+    let run_in_parallel = getenv_use_multiple_threads();
+
+    if n1_rows < PARALLEL_THRESHOLD || !run_in_parallel {
+        Ok(Array2::from_shape_fn((n1_rows, n2_cols), |(i, j)| {
+            (0..n2_rows)
+                .map(|k| mat1[[i, k]] & mat2[[k, j]])
+                .fold(false, |acc, v| acc ^ v)
+        }))
+    } else {
+        let mut result = Array2::from_elem((n1_rows, n2_cols), false);
+        result
+            .axis_iter_mut(Axis(0))
+            .into_par_iter()
+            .enumerate()
+            .for_each(|(i, mut row)| {
+                for j in 0..n2_cols {
+                    row[j] = (0..n2_rows)
+                        .map(|k| mat1[[i, k]] & mat2[[k, j]])
+                        .fold(false, |acc, v| acc ^ v)
+                }
+            });
+        Ok(result)
+    }
 }
 
 /// Gauss elimination of a matrix mat with m rows and n columns.
@@ -198,3 +224,16 @@ pub fn check_invertible_binary_matrix_inner(mat: ArrayView2<bool>) -> bool {
     let rank = compute_rank_inner(mat);
     rank == mat.nrows()
 }
+
+/// Mutate matrix ``mat`` in-place by swapping the contents of rows ``i`` and ``j``.
+pub fn swap_rows_inner(mut mat: ArrayViewMut2<bool>, i: usize, j: usize) {
+    let (mut x, mut y) = mat.multi_slice_mut((s![i, ..], s![j, ..]));
+    azip!((x in &mut x, y in &mut y) (*x, *y) = (*y, *x));
+}
+
+/// Mutate matrix ``mat`` in-place by replacing the contents of row ``i`` by ``row``.
+pub fn replace_row_inner(mut mat: ArrayViewMut2<bool>, i: usize, row: ArrayView1<bool>) {
+    let mut x = mat.slice_mut(s![i, ..]);
+    let y = row.slice(s![..]);
+    Zip::from(&mut x).and(&y).for_each(|x, &y| *x = y);
+}