diff --git a/crates/accelerate/src/synthesis/clifford/greedy_synthesis.rs b/crates/accelerate/src/synthesis/clifford/greedy_synthesis.rs new file mode 100644 index 000000000000..e53e25282005 --- /dev/null +++ b/crates/accelerate/src/synthesis/clifford/greedy_synthesis.rs @@ -0,0 +1,441 @@ +// 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 indexmap::IndexSet; +use ndarray::{s, ArrayView2}; +use smallvec::smallvec; + +use crate::synthesis::clifford::utils::CliffordGatesVec; +use crate::synthesis::clifford::utils::{adjust_final_pauli_gates, SymplecticMatrix}; +use qiskit_circuit::operations::StandardGate; +use qiskit_circuit::Qubit; + +/// Converts a pair of Paulis pauli_x and pauli_z acting on a specific qubit +/// to the corresponding index in [PauliPairsClass] or [SingleQubitGate] classes. +/// The input is given as a 4-tuple: (pauli_x stabilizer, pauli_x destabilizer, +/// pauli_z stabilizer, pauli_z destabilizer), and the output is an unsigned +/// integer from 0 to 15. +fn pauli_pair_to_index(xs: bool, xd: bool, zs: bool, zd: bool) -> usize { + ((xs as usize) << 3) | ((xd as usize) << 2) | ((zs as usize) << 1) | (zd as usize) +} + +/// The five classes of Pauli 2-qubit operators as described in the paper. +#[derive(Clone, Copy)] +enum PauliPairsClass { + ClassA, + ClassB, + ClassC, + ClassD, + ClassE, +} + +/// The 16 Pauli 2-qubit operators are divided into 5 equivalence classes +/// under the action of single-qubit Cliffords. +static PAULI_INDEX_TO_CLASS: [PauliPairsClass; 16] = [ + PauliPairsClass::ClassE, // 'II' + PauliPairsClass::ClassD, // 'IX' + PauliPairsClass::ClassD, // 'IZ' + PauliPairsClass::ClassD, // 'IY' + PauliPairsClass::ClassC, // 'XI' + PauliPairsClass::ClassB, // 'XX' + PauliPairsClass::ClassA, // 'XZ' + PauliPairsClass::ClassA, // 'XY' + PauliPairsClass::ClassC, // 'ZI' + PauliPairsClass::ClassA, // 'ZX' + PauliPairsClass::ClassB, // 'ZZ' + PauliPairsClass::ClassA, // 'ZY' + PauliPairsClass::ClassC, // 'YI' + PauliPairsClass::ClassA, // 'YX' + PauliPairsClass::ClassA, // 'YZ' + PauliPairsClass::ClassB, // 'YY' +]; + +/// Single-qubit Clifford gates modulo Paulis. +#[derive(Clone, Copy)] +enum SingleQubitGate { + GateI, + GateS, + GateH, + GateSH, + GateHS, + GateSHS, +} + +/// Maps pair of pauli operators to the single-qubit gate required +/// for the decoupling step. +static PAULI_INDEX_TO_1Q_GATE: [SingleQubitGate; 16] = [ + SingleQubitGate::GateI, // 'II' + SingleQubitGate::GateH, // 'IX' + SingleQubitGate::GateI, // 'IZ' + SingleQubitGate::GateSH, // 'IY' + SingleQubitGate::GateI, // 'XI' + SingleQubitGate::GateI, // 'XX' + SingleQubitGate::GateI, // 'XZ' + SingleQubitGate::GateSHS, // 'XY' + SingleQubitGate::GateH, // 'ZI' + SingleQubitGate::GateH, // 'ZX' + SingleQubitGate::GateH, // 'ZZ' + SingleQubitGate::GateSH, // 'ZY' + SingleQubitGate::GateS, // 'YI' + SingleQubitGate::GateHS, // 'YX' + SingleQubitGate::GateS, // 'YZ' + SingleQubitGate::GateS, // 'YY' +]; + +pub struct GreedyCliffordSynthesis<'a> { + /// The Clifford tableau to be synthesized. + tableau: ArrayView2<'a, bool>, + + /// The total number of qubits. + num_qubits: usize, + + /// Symplectic matrix being reduced. + symplectic_matrix: SymplecticMatrix, + + /// Unprocessed qubits. + unprocessed_qubits: IndexSet, +} + +impl GreedyCliffordSynthesis<'_> { + pub(crate) fn new(tableau: ArrayView2) -> Result, String> { + let tableau_shape = tableau.shape(); + if (tableau_shape[0] % 2 == 1) || (tableau_shape[1] != tableau_shape[0] + 1) { + return Err("The shape of the Clifford tableau is invalid".to_string()); + } + + let num_qubits = tableau_shape[0] / 2; + + // We are going to modify symplectic_matrix in-place until it + // becomes the identity. + let symplectic_matrix = SymplecticMatrix { + num_qubits, + smat: tableau.slice(s![.., 0..2 * num_qubits]).to_owned(), + }; + + let unprocessed_qubits: IndexSet = (0..num_qubits).collect(); + + Ok(GreedyCliffordSynthesis { + tableau, + num_qubits, + symplectic_matrix, + unprocessed_qubits, + }) + } + + /// Computes the CX cost of decoupling the symplectic matrix on the + /// given qubit. + fn compute_cost(&self, qubit: usize) -> Result { + let mut a_num = 0; + let mut b_num = 0; + let mut c_num = 0; + let mut d_num = 0; + + let mut qubit_is_in_a = false; + + for q in &self.unprocessed_qubits { + let pauli_pair_index = pauli_pair_to_index( + self.symplectic_matrix.smat[[*q, qubit + self.num_qubits]], + self.symplectic_matrix.smat[[*q + self.num_qubits, qubit + self.num_qubits]], + self.symplectic_matrix.smat[[*q, qubit]], + self.symplectic_matrix.smat[[*q + self.num_qubits, qubit]], + ); + let pauli_class = PAULI_INDEX_TO_CLASS[pauli_pair_index]; + + match pauli_class { + PauliPairsClass::ClassA => { + a_num += 1; + if *q == qubit { + qubit_is_in_a = true; + } + } + PauliPairsClass::ClassB => { + b_num += 1; + } + PauliPairsClass::ClassC => { + c_num += 1; + } + PauliPairsClass::ClassD => { + d_num += 1; + } + PauliPairsClass::ClassE => {} + } + } + + if a_num % 2 == 0 { + return Err("Symplectic Gaussian elimination failed.".to_string()); + } + + let mut cnot_cost: usize = + 3 * (a_num - 1) / 2 + (b_num + 1) * ((b_num > 0) as usize) + c_num + d_num; + + if !qubit_is_in_a { + cnot_cost += 3; + } + + Ok(cnot_cost) + } + + /// Calculate a decoupling operator D: + /// D^{-1} * Ox * D = x1 + /// D^{-1} * Oz * D = z1 + /// and reduces the clifford such that it will act trivially on min_qubit. + fn decouple_qubit( + &mut self, + gate_seq: &mut CliffordGatesVec, + min_qubit: usize, + ) -> Result<(), String> { + let mut a_qubits = IndexSet::new(); + let mut b_qubits = IndexSet::new(); + let mut c_qubits = IndexSet::new(); + let mut d_qubits = IndexSet::new(); + + for qubit in &self.unprocessed_qubits { + let pauli_pair_index = pauli_pair_to_index( + self.symplectic_matrix.smat[[*qubit, min_qubit + self.num_qubits]], + self.symplectic_matrix.smat + [[*qubit + self.num_qubits, min_qubit + self.num_qubits]], + self.symplectic_matrix.smat[[*qubit, min_qubit]], + self.symplectic_matrix.smat[[*qubit + self.num_qubits, min_qubit]], + ); + + let single_qubit_gate = PAULI_INDEX_TO_1Q_GATE[pauli_pair_index]; + match single_qubit_gate { + SingleQubitGate::GateS => { + gate_seq.push(( + StandardGate::SGate, + smallvec![], + smallvec![Qubit(*qubit as u32)], + )); + self.symplectic_matrix.prepend_s(*qubit); + } + SingleQubitGate::GateH => { + gate_seq.push(( + StandardGate::HGate, + smallvec![], + smallvec![Qubit(*qubit as u32)], + )); + self.symplectic_matrix.prepend_h(*qubit); + } + SingleQubitGate::GateSH => { + gate_seq.push(( + StandardGate::SGate, + smallvec![], + smallvec![Qubit(*qubit as u32)], + )); + gate_seq.push(( + StandardGate::HGate, + smallvec![], + smallvec![Qubit(*qubit as u32)], + )); + self.symplectic_matrix.prepend_s(*qubit); + self.symplectic_matrix.prepend_h(*qubit); + } + SingleQubitGate::GateHS => { + gate_seq.push(( + StandardGate::HGate, + smallvec![], + smallvec![Qubit(*qubit as u32)], + )); + gate_seq.push(( + StandardGate::SGate, + smallvec![], + smallvec![Qubit(*qubit as u32)], + )); + self.symplectic_matrix.prepend_h(*qubit); + self.symplectic_matrix.prepend_s(*qubit); + } + SingleQubitGate::GateSHS => { + gate_seq.push(( + StandardGate::SGate, + smallvec![], + smallvec![Qubit(*qubit as u32)], + )); + gate_seq.push(( + StandardGate::HGate, + smallvec![], + smallvec![Qubit(*qubit as u32)], + )); + gate_seq.push(( + StandardGate::SGate, + smallvec![], + smallvec![Qubit(*qubit as u32)], + )); + self.symplectic_matrix.prepend_s(*qubit); + self.symplectic_matrix.prepend_h(*qubit); + self.symplectic_matrix.prepend_s(*qubit); + } + SingleQubitGate::GateI => {} + } + + let pauli_class = PAULI_INDEX_TO_CLASS[pauli_pair_index]; + match pauli_class { + PauliPairsClass::ClassA => { + a_qubits.insert(*qubit); + } + PauliPairsClass::ClassB => { + b_qubits.insert(*qubit); + } + PauliPairsClass::ClassC => { + c_qubits.insert(*qubit); + } + PauliPairsClass::ClassD => { + d_qubits.insert(*qubit); + } + PauliPairsClass::ClassE => {} + } + } + + if a_qubits.len() % 2 != 1 { + return Err("Symplectic Gaussian elimination failed.".to_string()); + } + + if !a_qubits.contains(&min_qubit) { + let qubit_a = a_qubits[0]; + gate_seq.push(( + StandardGate::SwapGate, + smallvec![], + smallvec![Qubit(min_qubit as u32), Qubit(qubit_a as u32)], + )); + self.symplectic_matrix.prepend_swap(min_qubit, qubit_a); + + if b_qubits.contains(&min_qubit) { + b_qubits.swap_remove(&min_qubit); + b_qubits.insert(qubit_a); + } else if c_qubits.contains(&min_qubit) { + c_qubits.swap_remove(&min_qubit); + c_qubits.insert(qubit_a); + } else if d_qubits.contains(&min_qubit) { + d_qubits.swap_remove(&min_qubit); + d_qubits.insert(qubit_a); + } + + a_qubits.swap_remove(&qubit_a); + a_qubits.insert(min_qubit); + } + + for qubit in c_qubits { + gate_seq.push(( + StandardGate::CXGate, + smallvec![], + smallvec![Qubit(min_qubit as u32), Qubit(qubit as u32)], + )); + self.symplectic_matrix.prepend_cx(min_qubit, qubit); + } + + for qubit in d_qubits { + gate_seq.push(( + StandardGate::CXGate, + smallvec![], + smallvec![Qubit(qubit as u32), Qubit(min_qubit as u32)], + )); + self.symplectic_matrix.prepend_cx(qubit, min_qubit); + } + + if b_qubits.len() > 1 { + let qubit_b = b_qubits[0]; + for qubit in &b_qubits[1..] { + gate_seq.push(( + StandardGate::CXGate, + smallvec![], + smallvec![Qubit(qubit_b as u32), Qubit(*qubit as u32)], + )); + self.symplectic_matrix.prepend_cx(qubit_b, *qubit); + } + } + + if !b_qubits.is_empty() { + let qubit_b = b_qubits[0]; + gate_seq.push(( + StandardGate::CXGate, + smallvec![], + smallvec![Qubit(min_qubit as u32), Qubit(qubit_b as u32)], + )); + self.symplectic_matrix.prepend_cx(min_qubit, qubit_b); + + gate_seq.push(( + StandardGate::HGate, + smallvec![], + smallvec![Qubit(qubit_b as u32)], + )); + self.symplectic_matrix.prepend_h(qubit_b); + + gate_seq.push(( + StandardGate::CXGate, + smallvec![], + smallvec![Qubit(qubit_b as u32), Qubit(min_qubit as u32)], + )); + self.symplectic_matrix.prepend_cx(qubit_b, min_qubit); + } + + let a_len: usize = (a_qubits.len() - 1) / 2; + if a_len > 0 { + a_qubits.swap_remove(&min_qubit); + } + + for qubit in 0..a_len { + gate_seq.push(( + StandardGate::CXGate, + smallvec![], + smallvec![ + Qubit(a_qubits[2 * qubit + 1] as u32), + Qubit(a_qubits[2 * qubit] as u32) + ], + )); + self.symplectic_matrix + .prepend_cx(a_qubits[2 * qubit + 1], a_qubits[2 * qubit]); + + gate_seq.push(( + StandardGate::CXGate, + smallvec![], + smallvec![Qubit(a_qubits[2 * qubit] as u32), Qubit(min_qubit as u32)], + )); + self.symplectic_matrix + .prepend_cx(a_qubits[2 * qubit], min_qubit); + + gate_seq.push(( + StandardGate::CXGate, + smallvec![], + smallvec![ + Qubit(min_qubit as u32), + Qubit(a_qubits[2 * qubit + 1] as u32) + ], + )); + self.symplectic_matrix + .prepend_cx(min_qubit, a_qubits[2 * qubit + 1]); + } + + Ok(()) + } + + /// The main synthesis function. + pub(crate) fn run(&mut self) -> Result<(usize, CliffordGatesVec), String> { + let mut clifford_gates = CliffordGatesVec::new(); + + while !self.unprocessed_qubits.is_empty() { + let costs: Vec<(usize, usize)> = self + .unprocessed_qubits + .iter() + .map(|q| self.compute_cost(*q).map(|cost| (cost, *q))) + .collect::, _>>()?; + + let min_cost_qubit = costs.iter().min_by_key(|(cost, _)| cost).unwrap().1; + + self.decouple_qubit(&mut clifford_gates, min_cost_qubit)?; + + self.unprocessed_qubits.swap_remove(&min_cost_qubit); + } + + adjust_final_pauli_gates(&mut clifford_gates, self.tableau, self.num_qubits)?; + + Ok((self.num_qubits, clifford_gates)) + } +} diff --git a/crates/accelerate/src/synthesis/clifford/mod.rs b/crates/accelerate/src/synthesis/clifford/mod.rs new file mode 100644 index 000000000000..6772228acf88 --- /dev/null +++ b/crates/accelerate/src/synthesis/clifford/mod.rs @@ -0,0 +1,48 @@ +// 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. + +mod greedy_synthesis; +mod utils; + +use crate::synthesis::clifford::greedy_synthesis::GreedyCliffordSynthesis; +use crate::QiskitError; +use numpy::PyReadonlyArray2; +use pyo3::prelude::*; +use qiskit_circuit::circuit_data::CircuitData; +use qiskit_circuit::operations::Param; + +/// Create a circuit that synthesizes a given Clifford operator represented as a tableau. +/// +/// This is an implementation of the "greedy Clifford compiler" presented in +/// Appendix A of the paper "Clifford Circuit Optimization with Templates and Symbolic +/// Pauli Gates" by Bravyi, Shaydulin, Hu, and Maslov (2021), ``__. +/// +/// This method typically yields better CX cost compared to the Aaronson-Gottesman method. +/// +/// Note that this function only implements the greedy Clifford compiler and not the +/// templates and symbolic Pauli gates optimizations that are also described in the paper. +#[pyfunction] +#[pyo3(signature = (clifford))] +fn synth_clifford_greedy(py: Python, clifford: PyReadonlyArray2) -> PyResult { + let tableau = clifford.as_array(); + let mut greedy_synthesis = + GreedyCliffordSynthesis::new(tableau.view()).map_err(QiskitError::new_err)?; + let (num_qubits, clifford_gates) = greedy_synthesis.run().map_err(QiskitError::new_err)?; + + CircuitData::from_standard_gates(py, num_qubits as u32, clifford_gates, Param::Float(0.0)) +} + +#[pymodule] +pub fn clifford(m: &Bound) -> PyResult<()> { + m.add_function(wrap_pyfunction!(synth_clifford_greedy, m)?)?; + Ok(()) +} diff --git a/crates/accelerate/src/synthesis/clifford/utils.rs b/crates/accelerate/src/synthesis/clifford/utils.rs new file mode 100644 index 000000000000..766d84ed179d --- /dev/null +++ b/crates/accelerate/src/synthesis/clifford/utils.rs @@ -0,0 +1,289 @@ +// 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::calc_inverse_matrix_inner; +use ndarray::{azip, s, Array1, Array2, ArrayView2}; +use qiskit_circuit::operations::{Param, StandardGate}; +use qiskit_circuit::Qubit; +use smallvec::{smallvec, SmallVec}; + +/// Symplectic matrix. +/// Currently this class is internal to the synthesis library. +pub struct SymplecticMatrix { + /// Number of qubits. + pub num_qubits: usize, + /// Matrix with dimensions (2 * num_qubits) x (2 * num_qubits). + pub smat: Array2, +} + +/// Clifford. +/// Currently this class is internal to the synthesis library and +/// has a very different functionality from Qiskit's python-based +/// Clifford class. +pub struct Clifford { + /// Number of qubits. + pub num_qubits: usize, + /// Matrix with dimensions (2 * num_qubits) x (2 * num_qubits + 1). + pub tableau: Array2, +} + +impl SymplecticMatrix { + /// Modifies the matrix in-place by appending S-gate + #[allow(dead_code)] + pub fn append_s(&mut self, qubit: usize) { + let (x, mut z) = self + .smat + .multi_slice_mut((s![.., qubit], s![.., self.num_qubits + qubit])); + azip!((z in &mut z, &x in &x) *z ^= x); + } + + /// Modifies the matrix in-place by prepending S-gate + pub fn prepend_s(&mut self, qubit: usize) { + let (x, mut z) = self + .smat + .multi_slice_mut((s![self.num_qubits + qubit, ..], s![qubit, ..])); + azip!((z in &mut z, &x in &x) *z ^= x); + } + + /// Modifies the matrix in-place by appending H-gate + #[allow(dead_code)] + pub fn append_h(&mut self, qubit: usize) { + let (mut x, mut z) = self + .smat + .multi_slice_mut((s![.., qubit], s![.., self.num_qubits + qubit])); + azip!((x in &mut x, z in &mut z) (*x, *z) = (*z, *x)); + } + + /// Modifies the matrix in-place by prepending H-gate + pub fn prepend_h(&mut self, qubit: usize) { + let (mut x, mut z) = self + .smat + .multi_slice_mut((s![qubit, ..], s![self.num_qubits + qubit, ..])); + azip!((x in &mut x, z in &mut z) (*x, *z) = (*z, *x)); + } + + /// Modifies the matrix in-place by appending SWAP-gate + #[allow(dead_code)] + pub fn append_swap(&mut self, qubit0: usize, qubit1: usize) { + let (mut x0, mut z0, mut x1, mut z1) = self.smat.multi_slice_mut(( + s![.., qubit0], + s![.., self.num_qubits + qubit0], + s![.., qubit1], + s![.., self.num_qubits + qubit1], + )); + azip!((x0 in &mut x0, x1 in &mut x1) (*x0, *x1) = (*x1, *x0)); + azip!((z0 in &mut z0, z1 in &mut z1) (*z0, *z1) = (*z1, *z0)); + } + + /// Modifies the matrix in-place by prepending SWAP-gate + pub fn prepend_swap(&mut self, qubit0: usize, qubit1: usize) { + let (mut x0, mut z0, mut x1, mut z1) = self.smat.multi_slice_mut(( + s![qubit0, ..], + s![self.num_qubits + qubit0, ..], + s![qubit1, ..], + s![self.num_qubits + qubit1, ..], + )); + azip!((x0 in &mut x0, x1 in &mut x1) (*x0, *x1) = (*x1, *x0)); + azip!((z0 in &mut z0, z1 in &mut z1) (*z0, *z1) = (*z1, *z0)); + } + + /// Modifies the matrix in-place by appending CX-gate + #[allow(dead_code)] + pub fn append_cx(&mut self, qubit0: usize, qubit1: usize) { + let (x0, mut z0, mut x1, z1) = self.smat.multi_slice_mut(( + s![.., qubit0], + s![.., self.num_qubits + qubit0], + s![.., qubit1], + s![.., self.num_qubits + qubit1], + )); + azip!((x1 in &mut x1, &x0 in &x0) *x1 ^= x0); + azip!((z0 in &mut z0, &z1 in &z1) *z0 ^= z1); + } + + /// Modifies the matrix in-place by prepending CX-gate + pub fn prepend_cx(&mut self, qubit0: usize, qubit1: usize) { + let (x0, mut z0, mut x1, z1) = self.smat.multi_slice_mut(( + s![qubit1, ..], + s![self.num_qubits + qubit1, ..], + s![qubit0, ..], + s![self.num_qubits + qubit0, ..], + )); + azip!((x1 in &mut x1, &x0 in &x0) *x1 ^= x0); + azip!((z0 in &mut z0, &z1 in &z1) *z0 ^= z1); + } +} + +impl Clifford { + /// Modifies the tableau in-place by appending S-gate + pub fn append_s(&mut self, qubit: usize) { + let (x, mut z, mut p) = self.tableau.multi_slice_mut(( + s![.., qubit], + s![.., self.num_qubits + qubit], + s![.., 2 * self.num_qubits], + )); + + azip!((p in &mut p, &x in &x, &z in &z) *p ^= x & z); + azip!((z in &mut z, &x in &x) *z ^= x); + } + + /// Modifies the tableau in-place by appending Sdg-gate + #[allow(dead_code)] + pub fn append_sdg(&mut self, qubit: usize) { + let (x, mut z, mut p) = self.tableau.multi_slice_mut(( + s![.., qubit], + s![.., self.num_qubits + qubit], + s![.., 2 * self.num_qubits], + )); + + azip!((p in &mut p, &x in &x, &z in &z) *p ^= x & !z); + azip!((z in &mut z, &x in &x) *z ^= x); + } + + /// Modifies the tableau in-place by appending H-gate + pub fn append_h(&mut self, qubit: usize) { + let (mut x, mut z, mut p) = self.tableau.multi_slice_mut(( + s![.., qubit], + s![.., self.num_qubits + qubit], + s![.., 2 * self.num_qubits], + )); + + azip!((p in &mut p, &x in &x, &z in &z) *p ^= x & z); + azip!((x in &mut x, z in &mut z) (*x, *z) = (*z, *x)); + } + + /// Modifies the tableau in-place by appending SWAP-gate + pub fn append_swap(&mut self, qubit0: usize, qubit1: usize) { + let (mut x0, mut z0, mut x1, mut z1) = self.tableau.multi_slice_mut(( + s![.., qubit0], + s![.., self.num_qubits + qubit0], + s![.., qubit1], + s![.., self.num_qubits + qubit1], + )); + azip!((x0 in &mut x0, x1 in &mut x1) (*x0, *x1) = (*x1, *x0)); + azip!((z0 in &mut z0, z1 in &mut z1) (*z0, *z1) = (*z1, *z0)); + } + + /// Modifies the tableau in-place by appending CX-gate + pub fn append_cx(&mut self, qubit0: usize, qubit1: usize) { + let (x0, mut z0, mut x1, z1, mut p) = self.tableau.multi_slice_mut(( + s![.., qubit0], + s![.., self.num_qubits + qubit0], + s![.., qubit1], + s![.., self.num_qubits + qubit1], + s![.., 2 * self.num_qubits], + )); + azip!((p in &mut p, &x0 in &x0, &z0 in &z0, &x1 in &x1, &z1 in &z1) *p ^= (x1 ^ z0 ^ true) & z1 & x0); + azip!((x1 in &mut x1, &x0 in &x0) *x1 ^= x0); + azip!((z0 in &mut z0, &z1 in &z1) *z0 ^= z1); + } + + /// Creates a Clifford from a given sequence of Clifford gates. + /// In essence, starts from the identity tableau and modifies it + /// based on the gates in the sequence. + pub fn from_gate_sequence( + gate_seq: &CliffordGatesVec, + num_qubits: usize, + ) -> Result { + // create the identity + let mut clifford = Clifford { + num_qubits, + tableau: Array2::from_shape_fn((2 * num_qubits, 2 * num_qubits + 1), |(i, j)| i == j), + }; + + gate_seq + .iter() + .try_for_each(|(gate, _params, qubits)| match *gate { + StandardGate::SGate => { + clifford.append_s(qubits[0].0 as usize); + Ok(()) + } + StandardGate::HGate => { + clifford.append_h(qubits[0].0 as usize); + Ok(()) + } + StandardGate::CXGate => { + clifford.append_cx(qubits[0].0 as usize, qubits[1].0 as usize); + Ok(()) + } + StandardGate::SwapGate => { + clifford.append_swap(qubits[0].0 as usize, qubits[1].0 as usize); + Ok(()) + } + _ => Err(format!("Unsupported gate {:?}", gate)), + })?; + Ok(clifford) + } +} + +/// A sequence of Clifford gates. +/// Represents the return type of Clifford synthesis algorithms. +pub type CliffordGatesVec = Vec<(StandardGate, SmallVec<[Param; 3]>, SmallVec<[Qubit; 2]>)>; + +/// Given a sequence of Clifford gates that correctly implements the symplectic matrix +/// of the target clifford tableau, adds the Pauli gates to also match the phase of +/// the tableau. +pub fn adjust_final_pauli_gates( + gate_seq: &mut CliffordGatesVec, + target_tableau: ArrayView2, + num_qubits: usize, +) -> Result<(), String> { + // simulate the clifford circuit that we have constructed + let simulated_clifford = Clifford::from_gate_sequence(gate_seq, num_qubits)?; + + // compute the phase difference + let target_phase = target_tableau.column(2 * num_qubits); + let sim_phase = simulated_clifford.tableau.column(2 * num_qubits); + + let delta_phase: Vec = target_phase + .iter() + .zip(sim_phase.iter()) + .map(|(&a, &b)| a ^ b) + .collect(); + + // compute inverse of the symplectic matrix + let smat = target_tableau.slice(s![.., ..2 * num_qubits]); + let smat_inv = calc_inverse_matrix_inner(smat, false)?; + + // compute smat_inv * delta_phase + let arr1 = smat_inv.map(|v| *v as usize); + let vec2: Vec = delta_phase.into_iter().map(|v| v as usize).collect(); + let arr2 = Array1::from(vec2); + let delta_phase_pre = arr1.dot(&arr2).map(|v| v % 2 == 1); + + // add pauli gates + for qubit in 0..num_qubits { + if delta_phase_pre[qubit] && delta_phase_pre[qubit + num_qubits] { + // println!("=> Adding Y-gate on {}", qubit); + gate_seq.push(( + StandardGate::YGate, + smallvec![], + smallvec![Qubit(qubit as u32)], + )); + } else if delta_phase_pre[qubit] { + // println!("=> Adding Z-gate on {}", qubit); + gate_seq.push(( + StandardGate::ZGate, + smallvec![], + smallvec![Qubit(qubit as u32)], + )); + } else if delta_phase_pre[qubit + num_qubits] { + // println!("=> Adding X-gate on {}", qubit); + gate_seq.push(( + StandardGate::XGate, + smallvec![], + smallvec![Qubit(qubit as u32)], + )); + } + } + + Ok(()) +} diff --git a/crates/accelerate/src/synthesis/linear/mod.rs b/crates/accelerate/src/synthesis/linear/mod.rs index 2fa158ea761f..b184a170fa5f 100644 --- a/crates/accelerate/src/synthesis/linear/mod.rs +++ b/crates/accelerate/src/synthesis/linear/mod.rs @@ -14,7 +14,7 @@ use crate::QiskitError; use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2, PyReadwriteArray2}; use pyo3::prelude::*; -mod utils; +pub mod utils; #[pyfunction] #[pyo3(signature = (mat, ncols=None, full_elim=false))] diff --git a/crates/accelerate/src/synthesis/mod.rs b/crates/accelerate/src/synthesis/mod.rs index db28751437f6..1b9908ef80cf 100644 --- a/crates/accelerate/src/synthesis/mod.rs +++ b/crates/accelerate/src/synthesis/mod.rs @@ -10,7 +10,8 @@ // copyright notice, and modified files need to carry a notice indicating // that they have been altered from the originals. -mod linear; +mod clifford; +pub mod linear; mod permutation; use pyo3::prelude::*; @@ -18,7 +19,8 @@ use pyo3::wrap_pymodule; #[pymodule] pub fn synthesis(m: &Bound) -> PyResult<()> { - m.add_wrapped(wrap_pymodule!(permutation::permutation))?; m.add_wrapped(wrap_pymodule!(linear::linear))?; + m.add_wrapped(wrap_pymodule!(permutation::permutation))?; + m.add_wrapped(wrap_pymodule!(clifford::clifford))?; Ok(()) } diff --git a/qiskit/__init__.py b/qiskit/__init__.py index aca555da8cb8..9aaa7a76a68e 100644 --- a/qiskit/__init__.py +++ b/qiskit/__init__.py @@ -82,6 +82,7 @@ sys.modules["qiskit._accelerate.vf2_layout"] = _accelerate.vf2_layout sys.modules["qiskit._accelerate.synthesis.permutation"] = _accelerate.synthesis.permutation sys.modules["qiskit._accelerate.synthesis.linear"] = _accelerate.synthesis.linear +sys.modules["qiskit._accelerate.synthesis.clifford"] = _accelerate.synthesis.clifford from qiskit.exceptions import QiskitError, MissingOptionalLibraryError diff --git a/qiskit/synthesis/clifford/clifford_decompose_greedy.py b/qiskit/synthesis/clifford/clifford_decompose_greedy.py index 784e6706d62e..0a679a8a7a6f 100644 --- a/qiskit/synthesis/clifford/clifford_decompose_greedy.py +++ b/qiskit/synthesis/clifford/clifford_decompose_greedy.py @@ -12,22 +12,16 @@ """ Circuit synthesis for the Clifford class. """ -# pylint: disable=invalid-name # --------------------------------------------------------------------- # Synthesis based on Bravyi et. al. greedy clifford compiler # --------------------------------------------------------------------- - -import numpy as np from qiskit.circuit import QuantumCircuit -from qiskit.exceptions import QiskitError -from qiskit.quantum_info import Clifford, Pauli -from qiskit.quantum_info.operators.symplectic.clifford_circuits import ( - _append_cx, - _append_h, - _append_s, - _append_swap, +from qiskit.quantum_info import Clifford + +from qiskit._accelerate.synthesis.clifford import ( + synth_clifford_greedy as synth_clifford_greedy_inner, ) @@ -56,296 +50,8 @@ def synth_clifford_greedy(clifford: Clifford) -> QuantumCircuit: *Clifford Circuit Optimization with Templates and Symbolic Pauli Gates*, `arXiv:2105.02291 [quant-ph] `_ """ - - num_qubits = clifford.num_qubits - circ = QuantumCircuit(num_qubits, name=str(clifford)) - qubit_list = list(range(num_qubits)) - clifford_cpy = clifford.copy() - - # Reducing the original Clifford to identity - # via symplectic Gaussian elimination - while len(qubit_list) > 0: - # Calculate the adjoint of clifford_cpy without the phase - clifford_adj = clifford_cpy.copy() - tmp = clifford_adj.destab_x.copy() - clifford_adj.destab_x = clifford_adj.stab_z.T - clifford_adj.destab_z = clifford_adj.destab_z.T - clifford_adj.stab_x = clifford_adj.stab_x.T - clifford_adj.stab_z = tmp.T - - list_greedy_cost = [] - for qubit in qubit_list: - pauli_x = Pauli("I" * (num_qubits - qubit - 1) + "X" + "I" * qubit) - pauli_x = pauli_x.evolve(clifford_adj, frame="s") - - pauli_z = Pauli("I" * (num_qubits - qubit - 1) + "Z" + "I" * qubit) - pauli_z = pauli_z.evolve(clifford_adj, frame="s") - list_pairs = [] - pauli_count = 0 - - # Compute the CNOT cost in order to find the qubit with the minimal cost - for i in qubit_list: - typeq = _from_pair_paulis_to_type(pauli_x, pauli_z, i) - list_pairs.append(typeq) - pauli_count += 1 - cost = _compute_greedy_cost(list_pairs) - list_greedy_cost.append([cost, qubit]) - - _, min_qubit = (sorted(list_greedy_cost))[0] - - # Gaussian elimination step for the qubit with minimal CNOT cost - pauli_x = Pauli("I" * (num_qubits - min_qubit - 1) + "X" + "I" * min_qubit) - pauli_x = pauli_x.evolve(clifford_adj, frame="s") - - pauli_z = Pauli("I" * (num_qubits - min_qubit - 1) + "Z" + "I" * min_qubit) - pauli_z = pauli_z.evolve(clifford_adj, frame="s") - - # Compute the decoupling operator of cliff_ox and cliff_oz - decouple_circ, decouple_cliff = _calc_decoupling( - pauli_x, pauli_z, qubit_list, min_qubit, num_qubits, clifford_cpy - ) - circ = circ.compose(decouple_circ) - - # Now the clifford acts trivially on min_qubit - clifford_cpy = decouple_cliff.adjoint().compose(clifford_cpy) - qubit_list.remove(min_qubit) - - # Add the phases (Pauli gates) to the Clifford circuit - for qubit in range(num_qubits): - stab = clifford_cpy.stab_phase[qubit] - destab = clifford_cpy.destab_phase[qubit] - if destab and stab: - circ.y(qubit) - elif not destab and stab: - circ.x(qubit) - elif destab and not stab: - circ.z(qubit) - - return circ - - -# --------------------------------------------------------------------- -# Helper functions for Bravyi et. al. greedy clifford compiler -# --------------------------------------------------------------------- - -# Global arrays of the 16 pairs of Pauli operators -# divided into 5 equivalence classes under the action of single-qubit Cliffords - -# Class A - canonical representative is 'XZ' -A_class = [ - [[False, True], [True, True]], # 'XY' - [[False, True], [True, False]], # 'XZ' - [[True, True], [False, True]], # 'YX' - [[True, True], [True, False]], # 'YZ' - [[True, False], [False, True]], # 'ZX' - [[True, False], [True, True]], -] # 'ZY' - -# Class B - canonical representative is 'XX' -B_class = [ - [[True, False], [True, False]], # 'ZZ' - [[False, True], [False, True]], # 'XX' - [[True, True], [True, True]], -] # 'YY' - -# Class C - canonical representative is 'XI' -C_class = [ - [[True, False], [False, False]], # 'ZI' - [[False, True], [False, False]], # 'XI' - [[True, True], [False, False]], -] # 'YI' - -# Class D - canonical representative is 'IZ' -D_class = [ - [[False, False], [False, True]], # 'IX' - [[False, False], [True, False]], # 'IZ' - [[False, False], [True, True]], -] # 'IY' - -# Class E - only 'II' -E_class = [[[False, False], [False, False]]] # 'II' - - -def _from_pair_paulis_to_type(pauli_x, pauli_z, qubit): - """Converts a pair of Paulis pauli_x and pauli_z into a type""" - - type_x = [pauli_x.z[qubit], pauli_x.x[qubit]] - type_z = [pauli_z.z[qubit], pauli_z.x[qubit]] - return [type_x, type_z] - - -def _compute_greedy_cost(list_pairs): - """Compute the CNOT cost of one step of the algorithm""" - - A_num = 0 - B_num = 0 - C_num = 0 - D_num = 0 - - for pair in list_pairs: - if pair in A_class: - A_num += 1 - elif pair in B_class: - B_num += 1 - elif pair in C_class: - C_num += 1 - elif pair in D_class: - D_num += 1 - - if (A_num % 2) == 0: - raise QiskitError("Symplectic Gaussian elimination fails.") - - # Calculate the CNOT cost - cost = 3 * (A_num - 1) / 2 + (B_num + 1) * (B_num > 0) + C_num + D_num - if list_pairs[0] not in A_class: # additional SWAP - cost += 3 - - return cost - - -def _calc_decoupling(pauli_x, pauli_z, qubit_list, min_qubit, num_qubits, cliff): - """Calculate a decoupling operator D: - D^{-1} * Ox * D = x1 - D^{-1} * Oz * D = z1 - and reduce the clifford such that it will act trivially on min_qubit - """ - - circ = QuantumCircuit(num_qubits) - - # decouple_cliff is initialized to an identity clifford - decouple_cliff = cliff.copy() - num_qubits = decouple_cliff.num_qubits - decouple_cliff.phase = np.zeros(2 * num_qubits) - decouple_cliff.symplectic_matrix = np.eye(2 * num_qubits) - - qubit0 = min_qubit # The qubit for the symplectic Gaussian elimination - - # Reduce the pair of Paulis to a representative in the equivalence class - # ['XZ', 'XX', 'XI', 'IZ', 'II'] by adding single-qubit gates - for qubit in qubit_list: - - typeq = _from_pair_paulis_to_type(pauli_x, pauli_z, qubit) - - if typeq in [ - [[True, True], [False, False]], # 'YI' - [[True, True], [True, True]], # 'YY' - [[True, True], [True, False]], - ]: # 'YZ': - circ.s(qubit) - _append_s(decouple_cliff, qubit) - - elif typeq in [ - [[True, False], [False, False]], # 'ZI' - [[True, False], [True, False]], # 'ZZ' - [[True, False], [False, True]], # 'ZX' - [[False, False], [False, True]], - ]: # 'IX' - circ.h(qubit) - _append_h(decouple_cliff, qubit) - - elif typeq in [ - [[False, False], [True, True]], # 'IY' - [[True, False], [True, True]], - ]: # 'ZY' - circ.s(qubit) - circ.h(qubit) - _append_s(decouple_cliff, qubit) - _append_h(decouple_cliff, qubit) - - elif typeq == [[True, True], [False, True]]: # 'YX' - circ.h(qubit) - circ.s(qubit) - _append_h(decouple_cliff, qubit) - _append_s(decouple_cliff, qubit) - - elif typeq == [[False, True], [True, True]]: # 'XY' - circ.s(qubit) - circ.h(qubit) - circ.s(qubit) - _append_s(decouple_cliff, qubit) - _append_h(decouple_cliff, qubit) - _append_s(decouple_cliff, qubit) - - # Reducing each pair of Paulis (except of qubit0) to 'II' - # by adding two-qubit gates and single-qubit gates - A_qubits = [] - B_qubits = [] - C_qubits = [] - D_qubits = [] - - for qubit in qubit_list: - typeq = _from_pair_paulis_to_type(pauli_x, pauli_z, qubit) - if typeq in A_class: - A_qubits.append(qubit) - elif typeq in B_class: - B_qubits.append(qubit) - elif typeq in C_class: - C_qubits.append(qubit) - elif typeq in D_class: - D_qubits.append(qubit) - - if len(A_qubits) % 2 != 1: - raise QiskitError("Symplectic Gaussian elimination fails.") - - if qubit0 not in A_qubits: # SWAP qubit0 and qubitA - qubitA = A_qubits[0] - circ.swap(qubit0, qubitA) - _append_swap(decouple_cliff, qubit0, qubitA) - if qubit0 in B_qubits: - B_qubits.remove(qubit0) - B_qubits.append(qubitA) - A_qubits.remove(qubitA) - A_qubits.append(qubit0) - elif qubit0 in C_qubits: - C_qubits.remove(qubit0) - C_qubits.append(qubitA) - A_qubits.remove(qubitA) - A_qubits.append(qubit0) - elif qubit0 in D_qubits: - D_qubits.remove(qubit0) - D_qubits.append(qubitA) - A_qubits.remove(qubitA) - A_qubits.append(qubit0) - else: - A_qubits.remove(qubitA) - A_qubits.append(qubit0) - - # Reduce pairs in Class C to 'II' - for qubit in C_qubits: - circ.cx(qubit0, qubit) - _append_cx(decouple_cliff, qubit0, qubit) - - # Reduce pairs in Class D to 'II' - for qubit in D_qubits: - circ.cx(qubit, qubit0) - _append_cx(decouple_cliff, qubit, qubit0) - - # Reduce pairs in Class B to 'II' - if len(B_qubits) > 1: - for qubit in B_qubits[1:]: - qubitB = B_qubits[0] - circ.cx(qubitB, qubit) - _append_cx(decouple_cliff, qubitB, qubit) - - if len(B_qubits) > 0: - qubitB = B_qubits[0] - circ.cx(qubit0, qubitB) - circ.h(qubitB) - circ.cx(qubitB, qubit0) - _append_cx(decouple_cliff, qubit0, qubitB) - _append_h(decouple_cliff, qubitB) - _append_cx(decouple_cliff, qubitB, qubit0) - - # Reduce pairs in Class A (except of qubit0) to 'II' - Alen = int((len(A_qubits) - 1) / 2) - if Alen > 0: - A_qubits.remove(qubit0) - for qubit in range(Alen): - circ.cx(A_qubits[2 * qubit + 1], A_qubits[2 * qubit]) - circ.cx(A_qubits[2 * qubit], qubit0) - circ.cx(qubit0, A_qubits[2 * qubit + 1]) - _append_cx(decouple_cliff, A_qubits[2 * qubit + 1], A_qubits[2 * qubit]) - _append_cx(decouple_cliff, A_qubits[2 * qubit], qubit0) - _append_cx(decouple_cliff, qubit0, A_qubits[2 * qubit + 1]) - - return circ, decouple_cliff + circuit = QuantumCircuit._from_circuit_data( + synth_clifford_greedy_inner(clifford.tableau.astype(bool)) + ) + circuit.name = str(clifford) + return circuit diff --git a/releasenotes/notes/oxidize-synth-clifford-greedy-0739e9688bc4eedd.yaml b/releasenotes/notes/oxidize-synth-clifford-greedy-0739e9688bc4eedd.yaml new file mode 100644 index 000000000000..044fb7d44315 --- /dev/null +++ b/releasenotes/notes/oxidize-synth-clifford-greedy-0739e9688bc4eedd.yaml @@ -0,0 +1,6 @@ +--- +upgrade_synthesis: + - | + The function :func:`.synth_clifford_greedy` that synthesizes :class:`.Clifford` operators + was ported to Rust, leading to a significant increase in performance for all numbers of + qubits. For Cliffords over 50 qubits, the speedup is on the order of 1000 times. diff --git a/test/python/synthesis/test_clifford_sythesis.py b/test/python/synthesis/test_clifford_sythesis.py index 887f1af5ad99..8ca11f1ef251 100644 --- a/test/python/synthesis/test_clifford_sythesis.py +++ b/test/python/synthesis/test_clifford_sythesis.py @@ -16,6 +16,7 @@ import numpy as np from ddt import ddt from qiskit.circuit.random import random_clifford_circuit +from qiskit.quantum_info import random_clifford from qiskit.quantum_info.operators import Clifford from qiskit.synthesis.clifford import ( synth_clifford_full, @@ -99,8 +100,7 @@ def test_synth_greedy(self, num_qubits): rng = np.random.default_rng(1234) samples = 50 for _ in range(samples): - circ = random_clifford_circuit(num_qubits, 5 * num_qubits, seed=rng) - target = Clifford(circ) + target = random_clifford(num_qubits, rng) synth_circ = synth_clifford_greedy(target) value = Clifford(synth_circ) self.assertEqual(value, target)