diff --git a/crates/accelerate/src/commutation_analysis.rs b/crates/accelerate/src/commutation_analysis.rs index 9c23458d6e61..d497a3e26fcc 100644 --- a/crates/accelerate/src/commutation_analysis.rs +++ b/crates/accelerate/src/commutation_analysis.rs @@ -54,6 +54,7 @@ pub(crate) fn analyze_commutations_inner( py: Python, dag: &mut DAGCircuit, commutation_checker: &mut CommutationChecker, + approximation_degree: f64, ) -> PyResult<(CommutationSet, NodeIndices)> { let mut commutation_set: CommutationSet = HashMap::new(); let mut node_indices: NodeIndices = HashMap::new(); @@ -102,6 +103,7 @@ pub(crate) fn analyze_commutations_inner( qargs2, cargs2, MAX_NUM_QUBITS, + approximation_degree, )?; if !all_commute { break; @@ -132,17 +134,19 @@ pub(crate) fn analyze_commutations_inner( } #[pyfunction] -#[pyo3(signature = (dag, commutation_checker))] +#[pyo3(signature = (dag, commutation_checker, approximation_degree=1.))] pub(crate) fn analyze_commutations( py: Python, dag: &mut DAGCircuit, commutation_checker: &mut CommutationChecker, + approximation_degree: f64, ) -> PyResult> { // This returns two HashMaps: // * The commuting nodes per wire: {wire: [commuting_nodes_1, commuting_nodes_2, ...]} // * The index in which commutation set a given node is located on a wire: {(node, wire): index} // The Python dict will store both of these dictionaries in one. - let (commutation_set, node_indices) = analyze_commutations_inner(py, dag, commutation_checker)?; + let (commutation_set, node_indices) = + analyze_commutations_inner(py, dag, commutation_checker, approximation_degree)?; let out_dict = PyDict::new(py); diff --git a/crates/accelerate/src/commutation_cancellation.rs b/crates/accelerate/src/commutation_cancellation.rs index 7bd13c662a2e..45856f132632 100644 --- a/crates/accelerate/src/commutation_cancellation.rs +++ b/crates/accelerate/src/commutation_cancellation.rs @@ -57,12 +57,13 @@ struct CancellationSetKey { } #[pyfunction] -#[pyo3(signature = (dag, commutation_checker, basis_gates=None))] +#[pyo3(signature = (dag, commutation_checker, basis_gates=None, approximation_degree=1.))] pub(crate) fn cancel_commutations( py: Python, dag: &mut DAGCircuit, commutation_checker: &mut CommutationChecker, basis_gates: Option>, + approximation_degree: f64, ) -> PyResult<()> { let basis: HashSet = if let Some(basis) = basis_gates { basis @@ -97,7 +98,8 @@ pub(crate) fn cancel_commutations( sec_commutation_set_id), the value is the list gates that share the same gate type, qubits and commutation sets. */ - let (commutation_set, node_indices) = analyze_commutations_inner(py, dag, commutation_checker)?; + let (commutation_set, node_indices) = + analyze_commutations_inner(py, dag, commutation_checker, approximation_degree)?; let mut cancellation_sets: HashMap> = HashMap::new(); (0..dag.num_qubits() as u32).for_each(|qubit| { diff --git a/crates/accelerate/src/commutation_checker.rs b/crates/accelerate/src/commutation_checker.rs index d867c582c916..c99b35a47391 100644 --- a/crates/accelerate/src/commutation_checker.rs +++ b/crates/accelerate/src/commutation_checker.rs @@ -14,6 +14,7 @@ use hashbrown::{HashMap, HashSet}; use ndarray::linalg::kron; use ndarray::Array2; use num_complex::Complex64; +use num_complex::ComplexFloat; use once_cell::sync::Lazy; use smallvec::SmallVec; @@ -34,6 +35,7 @@ use qiskit_circuit::operations::{ }; use qiskit_circuit::{BitType, Clbit, Qubit}; +use crate::gate_metrics; use crate::unitary_compose; use crate::QiskitError; @@ -54,48 +56,30 @@ static SUPPORTED_OP: Lazy> = Lazy::new(|| { // and their pi-periodicity. Here we mean a gate is n-pi periodic, if for angles that are // multiples of n*pi, the gate is equal to the identity up to a global phase. // E.g. RX is generated by X and 2-pi periodic, while CRX is generated by CX and 4-pi periodic. -static SUPPORTED_ROTATIONS: Lazy)>> = Lazy::new(|| { +static SUPPORTED_ROTATIONS: Lazy>> = Lazy::new(|| { HashMap::from([ - ( - "rx", - (2, Some(OperationRef::StandardGate(StandardGate::XGate))), - ), - ( - "ry", - (2, Some(OperationRef::StandardGate(StandardGate::YGate))), - ), - ( - "rz", - (2, Some(OperationRef::StandardGate(StandardGate::ZGate))), - ), - ( - "p", - (2, Some(OperationRef::StandardGate(StandardGate::ZGate))), - ), - ( - "u1", - (2, Some(OperationRef::StandardGate(StandardGate::ZGate))), - ), - ("rxx", (2, None)), // None means the gate is in the commutation dictionary - ("ryy", (2, None)), - ("rzx", (2, None)), - ("rzz", (2, None)), + ("rx", Some(OperationRef::StandardGate(StandardGate::XGate))), + ("ry", Some(OperationRef::StandardGate(StandardGate::YGate))), + ("rz", Some(OperationRef::StandardGate(StandardGate::ZGate))), + ("p", Some(OperationRef::StandardGate(StandardGate::ZGate))), + ("u1", Some(OperationRef::StandardGate(StandardGate::ZGate))), + ("rxx", None), // None means the gate is in the commutation dictionary + ("ryy", None), + ("rzx", None), + ("rzz", None), ( "crx", - (4, Some(OperationRef::StandardGate(StandardGate::CXGate))), + Some(OperationRef::StandardGate(StandardGate::CXGate)), ), ( "cry", - (4, Some(OperationRef::StandardGate(StandardGate::CYGate))), + Some(OperationRef::StandardGate(StandardGate::CYGate)), ), ( "crz", - (4, Some(OperationRef::StandardGate(StandardGate::CZGate))), - ), - ( - "cp", - (2, Some(OperationRef::StandardGate(StandardGate::CZGate))), + Some(OperationRef::StandardGate(StandardGate::CZGate)), ), + ("cp", Some(OperationRef::StandardGate(StandardGate::CZGate))), ]) }); @@ -155,13 +139,14 @@ impl CommutationChecker { } } - #[pyo3(signature=(op1, op2, max_num_qubits=3))] + #[pyo3(signature=(op1, op2, max_num_qubits=3, approximation_degree=1.))] fn commute_nodes( &mut self, py: Python, op1: &DAGOpNode, op2: &DAGOpNode, max_num_qubits: u32, + approximation_degree: f64, ) -> PyResult { let (qargs1, qargs2) = get_bits::( py, @@ -185,10 +170,11 @@ impl CommutationChecker { &qargs2, &cargs2, max_num_qubits, + approximation_degree, ) } - #[pyo3(signature=(op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits=3))] + #[pyo3(signature=(op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits=3, approximation_degree=1.))] #[allow(clippy::too_many_arguments)] fn commute( &mut self, @@ -200,6 +186,7 @@ impl CommutationChecker { qargs2: Option<&Bound>, cargs2: Option<&Bound>, max_num_qubits: u32, + approximation_degree: f64, ) -> PyResult { let qargs1 = qargs1.map_or_else(|| Ok(PyTuple::empty(py)), PySequenceMethods::to_tuple)?; let cargs1 = cargs1.map_or_else(|| Ok(PyTuple::empty(py)), PySequenceMethods::to_tuple)?; @@ -220,6 +207,7 @@ impl CommutationChecker { &qargs2, &cargs2, max_num_qubits, + approximation_degree, ) } @@ -288,20 +276,20 @@ impl CommutationChecker { qargs2: &[Qubit], cargs2: &[Clbit], max_num_qubits: u32, + approximation_degree: f64, ) -> PyResult { - // relative and absolute tolerance used to (1) check whether rotation gates commute - // trivially (i.e. the rotation angle is so small we assume it commutes) and (2) define - // comparison for the matrix-based commutation checks - let rtol = 1e-5; - let atol = 1e-8; + // If the average gate infidelity is below this tolerance, they commute. The tolerance + // is set to max(1e-12, 1 - approximation_degree), to account for roundoffs and for + // consistency with other places in Qiskit. + let tol = 1e-12_f64.max(1. - approximation_degree); // if we have rotation gates, we attempt to map them to their generators, for example // RX -> X or CPhase -> CZ - let (op1, params1, trivial1) = map_rotation(op1, params1, rtol); + let (op1, params1, trivial1) = map_rotation(op1, params1, tol); if trivial1 { return Ok(true); } - let (op2, params2, trivial2) = map_rotation(op2, params2, rtol); + let (op2, params2, trivial2) = map_rotation(op2, params2, tol); if trivial2 { return Ok(true); } @@ -367,8 +355,7 @@ impl CommutationChecker { second_op, second_params, second_qargs, - rtol, - atol, + tol, ); } @@ -403,8 +390,7 @@ impl CommutationChecker { second_op, second_params, second_qargs, - rtol, - atol, + tol, )?; // TODO: implement a LRU cache for this @@ -439,8 +425,7 @@ impl CommutationChecker { second_op: &OperationRef, second_params: &[Param], second_qargs: &[Qubit], - rtol: f64, - atol: f64, + tol: f64, ) -> PyResult { // Compute relative positioning of qargs of the second gate to the first gate. // Since the qargs come out the same BitData, we already know there are no accidential @@ -481,81 +466,49 @@ impl CommutationChecker { None => return Ok(false), }; - if first_qarg == second_qarg { - match first_qarg.len() { - 1 => Ok(unitary_compose::commute_1q( - &first_mat.view(), - &second_mat.view(), - rtol, - atol, - )), - 2 => Ok(unitary_compose::commute_2q( - &first_mat.view(), - &second_mat.view(), - &[Qubit(0), Qubit(1)], - rtol, - atol, - )), - _ => Ok(unitary_compose::allclose( - &second_mat.dot(&first_mat).view(), - &first_mat.dot(&second_mat).view(), - rtol, - atol, - )), - } + // TODO Optimize this bit to avoid unnecessary Kronecker products: + // 1. We currently sort the operations for the cache by operation size, putting the + // *smaller* operation first: (smaller op, larger op) + // 2. This code here expands the first op to match the second -- hence we always + // match the operator sizes. + // This whole extension logic could be avoided since we know the second one is larger. + let extra_qarg2 = num_qubits - first_qarg.len() as u32; + let first_mat = if extra_qarg2 > 0 { + let id_op = Array2::::eye(usize::pow(2, extra_qarg2)); + kron(&id_op, &first_mat) } else { - // TODO Optimize this bit to avoid unnecessary Kronecker products: - // 1. We currently sort the operations for the cache by operation size, putting the - // *smaller* operation first: (smaller op, larger op) - // 2. This code here expands the first op to match the second -- hence we always - // match the operator sizes. - // This whole extension logic could be avoided since we know the second one is larger. - let extra_qarg2 = num_qubits - first_qarg.len() as u32; - let first_mat = if extra_qarg2 > 0 { - let id_op = Array2::::eye(usize::pow(2, extra_qarg2)); - kron(&id_op, &first_mat) - } else { - first_mat - }; - - // the 1 qubit case cannot happen, since that would already have been captured - // by the previous if clause; first_qarg == second_qarg (if they overlap they must - // be the same) - if num_qubits == 2 { - return Ok(unitary_compose::commute_2q( - &first_mat.view(), - &second_mat.view(), - &second_qarg, - rtol, - atol, - )); - }; + first_mat + }; - let op12 = match unitary_compose::compose( - &first_mat.view(), - &second_mat.view(), - &second_qarg, - false, - ) { - Ok(matrix) => matrix, - Err(e) => return Err(PyRuntimeError::new_err(e)), - }; - let op21 = match unitary_compose::compose( - &first_mat.view(), - &second_mat.view(), - &second_qarg, - true, - ) { - Ok(matrix) => matrix, - Err(e) => return Err(PyRuntimeError::new_err(e)), - }; - Ok(unitary_compose::allclose( - &op12.view(), - &op21.view(), - rtol, - atol, - )) - } + // the 1 qubit case cannot happen, since that would already have been captured + // by the previous if clause; first_qarg == second_qarg (if they overlap they must + // be the same) + let op12 = match unitary_compose::compose( + &first_mat.view(), + &second_mat.view(), + &second_qarg, + false, + ) { + Ok(matrix) => matrix, + Err(e) => return Err(PyRuntimeError::new_err(e)), + }; + let op21 = match unitary_compose::compose( + &first_mat.view(), + &second_mat.view(), + &second_qarg, + true, + ) { + Ok(matrix) => matrix, + Err(e) => return Err(PyRuntimeError::new_err(e)), + }; + let (fid, phase) = gate_metrics::gate_fidelity(&op12.view(), &op21.view(), None); + + // we consider the gates as commuting if the process fidelity of + // AB (BA)^\dagger is approximately the identity and there is no global phase difference + // let dim = op12.ncols() as f64; + // let matrix_tol = tol * dim.powi(2); + let matrix_tol = tol; + Ok(phase.abs() <= tol && (1.0 - fid).abs() <= matrix_tol) } fn clear_cache(&mut self) { @@ -652,13 +605,19 @@ fn map_rotation<'a>( ) -> (&'a OperationRef<'a>, &'a [Param], bool) { let name = op.name(); - if let Some((pi_multiple, generator)) = SUPPORTED_ROTATIONS.get(name) { + if let Some(generator) = SUPPORTED_ROTATIONS.get(name) { // If the rotation angle is below the tolerance, the gate is assumed to // commute with everything, and we simply return the operation with the flag that // it commutes trivially. if let Param::Float(angle) = params[0] { - let periodicity = (*pi_multiple as f64) * ::std::f64::consts::PI; - if (angle % periodicity).abs() < tol { + let gate = op + .standard_gate() + .expect("Supported gates are standard gates"); + let (tr_over_dim, dim) = gate_metrics::rotation_trace_and_dim(gate, angle) + .expect("All rotation should be covered at this point"); + let gate_fidelity = tr_over_dim.abs().powi(2); + let process_fidelity = (dim * gate_fidelity + 1.) / (dim + 1.); + if (1. - process_fidelity).abs() <= tol { return (op, params, true); }; }; diff --git a/crates/accelerate/src/gate_metrics.rs b/crates/accelerate/src/gate_metrics.rs new file mode 100644 index 000000000000..5e39d5148159 --- /dev/null +++ b/crates/accelerate/src/gate_metrics.rs @@ -0,0 +1,75 @@ +// This code is part of Qiskit. +// +// (C) Copyright IBM 2025 +// +// 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 ndarray::ArrayView2; +use num_complex::{Complex64, ComplexFloat}; +use qiskit_circuit::{operations::StandardGate, Qubit}; + +use crate::unitary_compose; + +/// For a (controlled) rotation or phase gate, return a tuple ``(Tr(gate) / dim, dim)``. +/// Returns ``None`` if the rotation gate (specified by name) is not supported. +pub fn rotation_trace_and_dim(rotation: StandardGate, angle: f64) -> Option<(Complex64, f64)> { + let dim = match rotation { + StandardGate::RXGate + | StandardGate::RYGate + | StandardGate::RZGate + | StandardGate::PhaseGate + | StandardGate::U1Gate => 2., + _ => 4., + }; + + let trace_over_dim = match rotation { + StandardGate::RXGate + | StandardGate::RYGate + | StandardGate::RZGate + | StandardGate::RXXGate + | StandardGate::RYYGate + | StandardGate::RZZGate + | StandardGate::RZXGate => Complex64::new((angle / 2.).cos(), 0.), + StandardGate::CRXGate | StandardGate::CRYGate | StandardGate::CRZGate => { + Complex64::new(0.5 + 0.5 * (angle / 2.).cos(), 0.) + } + StandardGate::PhaseGate | StandardGate::U1Gate => { + (1. + Complex64::new(0., angle).exp()) / 2. + } + StandardGate::CPhaseGate => (3. + Complex64::new(0., angle).exp()) / 4., + _ => return None, + }; + Some((trace_over_dim, dim)) +} + +pub fn gate_fidelity( + left: &ArrayView2, + right: &ArrayView2, + qargs: Option<&[Qubit]>, +) -> (f64, f64) { + let dim = left.nrows(); + + let left = left.t().mapv(|el| el.conj()); + let product = match dim { + 2 => unitary_compose::matmul_1q(&left.view(), right), + 4 => { + unitary_compose::matmul_2q(&left.view(), right, qargs.unwrap_or(&[Qubit(0), Qubit(1)])) + } + _ => left.dot(right), + }; + let trace = product.diag().sum(); + + let dim = dim as f64; + let normalized_trace = trace / dim; + let phase = normalized_trace.arg(); // compute phase difference + + let process_fidelity = normalized_trace.abs().powi(2); + let gate_fidelity = (dim * process_fidelity + 1.) / (dim + 1.); + (gate_fidelity, phase) +} diff --git a/crates/accelerate/src/lib.rs b/crates/accelerate/src/lib.rs index bcb21ee3774e..ef90ac27041e 100644 --- a/crates/accelerate/src/lib.rs +++ b/crates/accelerate/src/lib.rs @@ -31,6 +31,7 @@ pub mod error_map; pub mod euler_one_qubit_decomposer; pub mod filter_op_nodes; pub mod gate_direction; +pub mod gate_metrics; pub mod gates_in_basis; pub mod high_level_synthesis; pub mod inverse_cancellation; diff --git a/crates/accelerate/src/remove_identity_equiv.rs b/crates/accelerate/src/remove_identity_equiv.rs index 7a65a9972019..33e9d630a0b6 100644 --- a/crates/accelerate/src/remove_identity_equiv.rs +++ b/crates/accelerate/src/remove_identity_equiv.rs @@ -14,6 +14,7 @@ use num_complex::ComplexFloat; use pyo3::prelude::*; use rustworkx_core::petgraph::stable_graph::NodeIndex; +use crate::gate_metrics::rotation_trace_and_dim; use crate::nlayout::PhysicalQubit; use crate::target_transpiler::Target; use qiskit_circuit::dag_circuit::DAGCircuit; @@ -23,6 +24,8 @@ use qiskit_circuit::operations::Param; use qiskit_circuit::operations::StandardGate; use qiskit_circuit::packed_instruction::PackedInstruction; +const MINIMUM_TOL: f64 = 1e-12; + #[pyfunction] #[pyo3(signature=(dag, approx_degree=Some(1.0), target=None))] fn remove_identity_equiv( @@ -33,12 +36,14 @@ fn remove_identity_equiv( ) { let mut remove_list: Vec = Vec::new(); let mut global_phase_update: f64 = 0.; + // Minimum threshold to compare average gate fidelity to 1. This is chosen to account + // for roundoff errors and to be consistent with other places. let get_error_cutoff = |inst: &PackedInstruction| -> f64 { match approx_degree { Some(degree) => { if degree == 1.0 { - f64::EPSILON + MINIMUM_TOL } else { match target { Some(target) => { @@ -50,10 +55,10 @@ fn remove_identity_equiv( let error_rate = target.get_error(inst.op.name(), qargs.as_slice()); match error_rate { Some(err) => err * degree, - None => f64::EPSILON.max(1. - degree), + None => MINIMUM_TOL.max(1. - degree), } } - None => f64::EPSILON.max(1. - degree), + None => MINIMUM_TOL.max(1. - degree), } } } @@ -67,10 +72,10 @@ fn remove_identity_equiv( let error_rate = target.get_error(inst.op.name(), qargs.as_slice()); match error_rate { Some(err) => err, - None => f64::EPSILON, + None => MINIMUM_TOL, } } - None => f64::EPSILON, + None => MINIMUM_TOL, }, } }; @@ -83,22 +88,23 @@ fn remove_identity_equiv( let view = inst.op.view(); match view { OperationRef::StandardGate(gate) => { - let (dim, trace) = match gate { - StandardGate::RXGate | StandardGate::RYGate | StandardGate::RZGate => { - if let Param::Float(theta) = inst.params_view()[0] { - let trace = Complex64::new((theta / 2.).cos() * 2., 0.); - (2., trace) - } else { - continue; - } - } - StandardGate::RXXGate + let (tr_over_dim, dim) = match gate { + StandardGate::RXGate + | StandardGate::RYGate + | StandardGate::RZGate + | StandardGate::PhaseGate + | StandardGate::RXXGate | StandardGate::RYYGate + | StandardGate::RZXGate | StandardGate::RZZGate - | StandardGate::RZXGate => { - if let Param::Float(theta) = inst.params_view()[0] { - let trace = Complex64::new((theta / 2.).cos() * 4., 0.); - (4., trace) + | StandardGate::CRXGate + | StandardGate::CRYGate + | StandardGate::CRZGate + | StandardGate::CPhaseGate => { + if let Param::Float(angle) = inst.params_view()[0] { + let (tr_over_dim, dim) = + rotation_trace_and_dim(gate, angle).expect("Since only supported rotation gates are given, the result is not None"); + (tr_over_dim, dim) } else { continue; } @@ -106,19 +112,19 @@ fn remove_identity_equiv( _ => { if let Some(matrix) = gate.matrix(inst.params_view()) { let dim = matrix.shape()[0] as f64; - let trace = matrix.diag().iter().sum::(); - (dim, trace) + let tr_over_dim = matrix.diag().iter().sum::() / dim; + (tr_over_dim, dim) } else { continue; } } }; let error = get_error_cutoff(inst); - let f_pro = (trace / dim).abs().powi(2); + let f_pro = tr_over_dim.abs().powi(2); let gate_fidelity = (dim * f_pro + 1.) / (dim + 1.); if (1. - gate_fidelity).abs() < error { remove_list.push(op_node); - global_phase_update += (trace / dim).arg(); + global_phase_update += tr_over_dim.arg(); } } _ => { @@ -127,12 +133,12 @@ fn remove_identity_equiv( if let Some(matrix) = matrix { let error = get_error_cutoff(inst); let dim = matrix.shape()[0] as f64; - let trace: Complex64 = matrix.diag().iter().sum(); - let f_pro = (trace / dim).abs().powi(2); + let tr_over_dim = matrix.diag().iter().sum::() / dim; + let f_pro = tr_over_dim.abs().powi(2); let gate_fidelity = (dim * f_pro + 1.) / (dim + 1.); if (1. - gate_fidelity).abs() < error { remove_list.push(op_node); - global_phase_update += (trace / dim).arg(); + global_phase_update += tr_over_dim.arg(); } } } diff --git a/crates/accelerate/src/unitary_compose.rs b/crates/accelerate/src/unitary_compose.rs index c04c8103af93..4a8805a03f95 100644 --- a/crates/accelerate/src/unitary_compose.rs +++ b/crates/accelerate/src/unitary_compose.rs @@ -12,8 +12,7 @@ use ndarray::{Array, Array2, ArrayView, ArrayView2, IxDyn}; use ndarray_einsum::*; -use num_complex::{Complex, Complex64, ComplexFloat}; -use num_traits::Zero; +use num_complex::{Complex, Complex64}; use qiskit_circuit::Qubit; static LOWERCASE: [u8; 26] = [ @@ -153,58 +152,14 @@ fn _einsum_matmul_index(qubits: &[u32], num_qubits: usize) -> String { ) } -pub fn commute_1q( - left: &ArrayView2, - right: &ArrayView2, - rtol: f64, - atol: f64, -) -> bool { - // This could allow for explicit hardcoded formulas, using less FLOPS, if we only - // consider an absolute tolerance. But for backward compatibility we now implement the full - // formula including relative tolerance handling. - for i in 0..2usize { - for j in 0..2usize { - let mut ab = Complex64::zero(); - let mut ba = Complex64::zero(); - for k in 0..2usize { - ab += left[[i, k]] * right[[k, j]]; - ba += right[[i, k]] * left[[k, j]]; - } - let sum = ab - ba; - if sum.abs() > atol + ba.abs() * rtol { - return false; - } - } - } - true -} - -pub fn commute_2q( - left: &ArrayView2, - right: &ArrayView2, - qargs: &[Qubit], - rtol: f64, - atol: f64, -) -> bool { - let rev = qargs[0].0 == 1; - for i in 0..4usize { - for j in 0..4usize { - // We compute AB and BA separately, to enable checking the relative difference - // (AB - BA)_ij > atol + rtol * BA_ij. This is due to backward compatibility and could - // maybe be changed in the future to save one complex number allocation. - let mut ab = Complex64::zero(); - let mut ba = Complex64::zero(); - for k in 0..4usize { - ab += left[[_ind(i, rev), _ind(k, rev)]] * right[[k, j]]; - ba += right[[i, k]] * left[[_ind(k, rev), _ind(j, rev)]]; - } - let sum = ab - ba; - if sum.abs() > atol + ba.abs() * rtol { - return false; - } - } - } - true +#[inline] +pub fn matmul_1q(left: &ArrayView2, right: &ArrayView2) -> Array2 { + let mut out = Array2::zeros((2, 2)); + out[[0, 0]] = left[[0, 0]] * right[[0, 0]] + left[[0, 1]] * right[[1, 0]]; + out[[0, 1]] = left[[0, 0]] * right[[0, 1]] + left[[0, 1]] * right[[1, 1]]; + out[[1, 0]] = left[[1, 0]] * right[[0, 0]] + left[[1, 1]] * right[[1, 0]]; + out[[1, 1]] = left[[1, 0]] * right[[0, 1]] + left[[1, 1]] * right[[1, 1]]; + out } #[inline] @@ -217,24 +172,21 @@ fn _ind(i: usize, reversed: bool) -> usize { } } -/// For equally sized matrices, ``left`` and ``right``, check whether all entries are close -/// by the criterion -/// -/// |left_ij - right_ij| <= atol + rtol * right_ij -/// -/// This is analogous to NumPy's ``allclose`` function. -pub fn allclose( +#[inline] +pub fn matmul_2q( left: &ArrayView2, right: &ArrayView2, - rtol: f64, - atol: f64, -) -> bool { - for i in 0..left.nrows() { - for j in 0..left.ncols() { - if (left[(i, j)] - right[(i, j)]).abs() > atol + rtol * right[(i, j)].abs() { - return false; + qargs: &[Qubit], +) -> Array2 { + let mut out = Array2::zeros((4, 4)); + + let rev = qargs[0].0 == 1; + for i in 0..4usize { + for j in 0..4usize { + for k in 0..4usize { + out[[i, j]] += left[[_ind(i, rev), _ind(k, rev)]] * right[[k, j]]; } } } - true + out } diff --git a/qiskit/circuit/commutation_checker.py b/qiskit/circuit/commutation_checker.py index 2d474c7bac4e..3e4f42fe39ff 100644 --- a/qiskit/circuit/commutation_checker.py +++ b/qiskit/circuit/commutation_checker.py @@ -12,6 +12,7 @@ """Code from commutative_analysis pass that checks commutation relations between DAG nodes.""" +from __future__ import annotations from typing import List, Union, Set, Optional from qiskit.circuit.operation import Operation @@ -19,11 +20,31 @@ class CommutationChecker: - """This code is essentially copy-pasted from commutative_analysis.py. - This code cleverly hashes commutativity and non-commutativity results between DAG nodes and seems - quite efficient for large Clifford circuits. - They may be other possible efficiency improvements: using rule-based commutativity analysis, - evicting from the cache less useful entries, etc. + r"""Check commutations of two operations. + + Two unitaries :math:`A` and :math:`B` on :math:`n` qubits commute if + + .. math:: + + \frac{2^n F_{\text{process}}(AB, BA) + 1}{2^n + 1} > 1 - \varepsilon, + + where + + .. math:: + + F_{\text{process}}(U_1, U_2) = \left|\frac{\mathrm{Tr}(U_1 U_2^\dagger)}{2^n} \right|^2, + + and we set :math:`\varepsilon` to :math:`10^{-12}` to account for round-off errors on + few-qubit systems. This metric is chosen for consistency with other closeness checks in + Qiskit. + + When possible, commutation relations are queried from a lookup table. This is the case + for standard gates without parameters (such as :class:`.XGate` or :class:`.HGate`) or + gates with free parameters (such as :class:`.RXGate` with a :class:`.ParameterExpression` as + angle). Otherwise, a matrix-based check is performed, where two operations are said to + commute, if the average gate fidelity of performing the commutation is above a certain threshold + (see ``approximation_degree``). The result of this commutation is then added to the + cached lookup table. """ def __init__( @@ -40,9 +61,10 @@ def commute_nodes( op1, op2, max_num_qubits: int = 3, + approximation_degree: float = 1.0, ) -> bool: """Checks if two DAGOpNodes commute.""" - return self.cc.commute_nodes(op1, op2, max_num_qubits) + return self.cc.commute_nodes(op1, op2, max_num_qubits, approximation_degree) def commute( self, @@ -53,6 +75,7 @@ def commute( qargs2: List, cargs2: List, max_num_qubits: int = 3, + approximation_degree: float = 1.0, ) -> bool: """ Checks if two Operations commute. The return value of `True` means that the operations @@ -69,11 +92,15 @@ def commute( cargs2: second operation's clbits. max_num_qubits: the maximum number of qubits to consider, the check may be skipped if the number of qubits for either operation exceeds this amount. + approximation_degree: If the average gate fidelity in between the two operations + is above this number (up to ``1e-12``) they are assumed to commute. Returns: bool: whether two operations commute. """ - return self.cc.commute(op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits) + return self.cc.commute( + op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits, approximation_degree + ) def num_cached_entries(self): """Returns number of cached entries""" diff --git a/qiskit/transpiler/passes/optimization/commutation_analysis.py b/qiskit/transpiler/passes/optimization/commutation_analysis.py index d801e4775937..d977f3b97514 100644 --- a/qiskit/transpiler/passes/optimization/commutation_analysis.py +++ b/qiskit/transpiler/passes/optimization/commutation_analysis.py @@ -18,10 +18,10 @@ class CommutationAnalysis(AnalysisPass): - """Analysis pass to find commutation relations between DAG nodes. + r"""Analysis pass to find commutation relations between DAG nodes. - ``property_set['commutation_set']`` is a dictionary that describes - the commutation relations on a given wire, all the gates on a wire + This sets ``property_set['commutation_set']`` to a dictionary that describes + the commutation relations on a given wire: all the gates on a wire are grouped into a set of gates that commute. """ diff --git a/qiskit/transpiler/passes/optimization/remove_identity_equiv.py b/qiskit/transpiler/passes/optimization/remove_identity_equiv.py index 17445eb5ad2a..53680475e5af 100644 --- a/qiskit/transpiler/passes/optimization/remove_identity_equiv.py +++ b/qiskit/transpiler/passes/optimization/remove_identity_equiv.py @@ -32,7 +32,7 @@ class RemoveIdentityEquivalent(TransformationPass): .. math:: - \bar{F} = \frac{1 + F_{\text{process}}}{1 + d},\ + \bar{F} = \frac{1 + d F_{\text{process}}}{1 + d},\ F_{\text{process}} = \frac{|\mathrm{Tr}(G)|^2}{d^2} @@ -47,9 +47,11 @@ def __init__( Args: approximation_degree: The degree to approximate for the equivalence check. This can be a floating point value between 0 and 1, or ``None``. If the value is 1 this does not - approximate above floating point precision. For a value < 1 this is used as a scaling - factor for the cutoff fidelity. If the value is ``None`` this approximates up to the - fidelity for the gate specified in ``target``. + approximate above the floating point precision. For a value < 1 this is used as a + scaling factor for the cutoff fidelity. If the value is ``None`` this approximates up + to the fidelity for the gate specified in ``target``. In case no ``target`` is set + we approximate up to ``16 * machine_eps`` as default to account for accumulations + on few-qubit systems. target: If ``approximation_degree`` is set to ``None`` and a :class:`.Target` is provided for this field the tolerance for determining whether an operation is equivalent to diff --git a/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml b/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml new file mode 100644 index 000000000000..e36c0c84b42b --- /dev/null +++ b/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml @@ -0,0 +1,35 @@ +--- +upgrade: + - | + Increased the minimum threshold for when gates are assumed to be the identity in + :class:`.RemoveIdentityEquivalent` from machine epsilon to ``1e-12`` to + account for round-off errors in the fidelity calculation and for consistency with the + other classes, such as :class:`.CommutationAnalysis` and :class:`.TwoQubitWeylDecomposition`. + - | + Updated the metric used to check commutations in :class:`.CommutationChecker`. + Two gates are assumed to commute if the average gate fidelity of the commutation is + above ``1 - 1e-12``. This value is chosen to account for round-off errors in the + fidelity calculation and for consistency with :class:`.RemoveIdentityEquivalent` and + :class:`.TwoQubitWeylDecomposition`. See the class docstring for more information. +fixes: + - | + Fixed an inconsistency in dealing with close-to-identity gates in the transpilation process, + where gates that are close to the identity were assumed to commute with everything, but + not removed. The underlying issue were different metrics used in :class:`.RemoveIdentityEquivalent` + and :class:`.CommutationAnalysis` (and, by extension, :class:`.CommutativeInverseCancellation`). + Both now use the average gate fidelity and the same threshold to assess whether a gate + should be treated as identity (such as a rotation gate with very small angle). See either + of the aforementioned classes' docstrings for more information. + Fixed `#13547 `__. +features_circuits: + - | + Added a new argument ``approximation_degree`` to :meth:`.CommutationChecker.commute` and + :meth:`.CommutationChecker.commute_nodes`, which allows to set the approximation threshold + for when gates are said to commute. See the docstring of :class:`.CommutationChecker` for more + detail. +features_transpiler: + - | + Added a new argument ``approximation_degree`` to :class:`.CommutationAnalysis`, which allows + setting the approximation threshold for when gates are said to commute. See the class docstring + for more information. + diff --git a/test/python/circuit/test_commutation_checker.py b/test/python/circuit/test_commutation_checker.py index 222d4b59cc6b..a2eed5d75af8 100644 --- a/test/python/circuit/test_commutation_checker.py +++ b/test/python/circuit/test_commutation_checker.py @@ -54,6 +54,7 @@ RZZGate, SGate, XGate, + YGate, ZGate, HGate, UnitaryGate, @@ -112,41 +113,51 @@ def test_simple_gates(self): different orders of gates, different orders of qubits, different sets of qubits over which gates are defined, and so on.""" - # should commute + self.assertTrue(scc.commute(HGate(), [0], [], HGate(), [0], [])) + self.assertTrue(scc.commute(ZGate(), [0], [], CXGate(), [0, 1], [])) - # should not commute self.assertFalse(scc.commute(ZGate(), [1], [], CXGate(), [0, 1], [])) - # should not commute + self.assertFalse(scc.commute(XGate(), [0], [], CXGate(), [0, 1], [])) - # should commute self.assertTrue(scc.commute(XGate(), [1], [], CXGate(), [0, 1], [])) - # should not commute self.assertFalse(scc.commute(XGate(), [1], [], CXGate(), [1, 0], [])) - # should commute self.assertTrue(scc.commute(XGate(), [0], [], CXGate(), [1, 0], [])) - # should commute self.assertTrue(scc.commute(CXGate(), [1, 0], [], XGate(), [0], [])) - # should not commute self.assertFalse(scc.commute(CXGate(), [1, 0], [], XGate(), [1], [])) - # should commute + self.assertTrue(scc.commute(CXGate(), [1, 0], [], CXGate(), [1, 0], [])) - # should not commute self.assertFalse(scc.commute(CXGate(), [1, 0], [], CXGate(), [0, 1], [])) - # should commute self.assertTrue(scc.commute(CXGate(), [1, 0], [], CXGate(), [1, 2], [])) - # should not commute self.assertFalse(scc.commute(CXGate(), [1, 0], [], CXGate(), [2, 1], [])) - # should commute self.assertTrue(scc.commute(CXGate(), [1, 0], [], CXGate(), [2, 3], [])) + self.assertTrue(scc.commute(XGate(), [2], [], CCXGate(), [0, 1, 2], [])) self.assertFalse(scc.commute(CCXGate(), [0, 1, 2], [], CCXGate(), [0, 2, 1], [])) + # these would commute up to a global phase + self.assertFalse(scc.commute(HGate(), [0], [], YGate(), [0], [])) + + def test_simple_matrices(self): + """Test simple gates but matrix-based.""" + x = UnitaryGate(XGate()) + had = UnitaryGate(HGate()) + had2 = UnitaryGate(np.kron(HGate(), HGate())) + cx = UnitaryGate(CXGate()) + + self.assertTrue(scc.commute(x, [0], [], x, [0], [])) + self.assertTrue(scc.commute(had, [0], [], had, [0], [])) + + self.assertTrue(scc.commute(had2, [0, 1], [], had2, [1, 0], [])) + self.assertFalse(scc.commute(had2, [0, 1], [], cx, [1, 0], [])) + self.assertTrue(scc.commute(cx, [0, 1], [], cx, [0, 1], [])) + + self.assertFalse(scc.commute(x, [0], [], cx, [0, 1], [])) + self.assertTrue(scc.commute(x, [1], [], cx, [0, 1], [])) + def test_passing_quantum_registers(self): """Check that passing QuantumRegisters works correctly.""" qr = QuantumRegister(4) - # should commute self.assertTrue(scc.commute(CXGate(), [qr[1], qr[0]], [], CXGate(), [qr[1], qr[2]], [])) - # should not commute self.assertFalse(scc.commute(CXGate(), [qr[0], qr[1]], [], CXGate(), [qr[1], qr[2]], [])) def test_standard_gates_commutations(self): @@ -387,21 +398,25 @@ def test_cutoff_angles(self, gate_cls): generic_gate = DCXGate() # gate that does not commute with any rotation gate - cutoff_angle = 1e-5 # this is the cutoff we use in the CommutationChecker + # the cutoff angle depends on the average gate fidelity; i.e. it is the angle + # for which the average gate fidelity is smaller than 1e-12 + if gate_cls in [CPhaseGate, CRXGate, CRYGate, CRZGate]: + cutoff_angle = 3.16e-6 + else: + cutoff_angle = 2.2e-6 for i in range(1, max_power + 1): angle = 2 ** (-i) gate = gate_cls(angle) qargs = list(range(gate.num_qubits)) - if angle < cutoff_angle: self.assertTrue(scc.commute(generic_gate, [0, 1], [], gate, qargs, [])) else: self.assertFalse(scc.commute(generic_gate, [0, 1], [], gate, qargs, [])) @idata(ROTATION_GATES) - def test_controlled_rotation_mod_4pi(self, gate_cls): - """Test the rotations modulo 2pi (4pi for controlled-rx/y/z) commute with any gate.""" + def test_rotations_pi_multiples(self, gate_cls): + """Test the rotations modulo 2pi (crx/cry/crz modulo 4pi) commute with any gate.""" generic_gate = HGate() # does not commute with any rotation gate multiples = np.arange(-6, 7) @@ -413,7 +428,13 @@ def test_controlled_rotation_mod_4pi(self, gate_cls): # compute a numeric reference, that doesn't go through any special cases and # uses a matrix-based commutation check expected = scc.commute( - generic_gate, [0], [], numeric, list(range(gate.num_qubits)), [] + generic_gate, + [0], + [], + numeric, + list(range(gate.num_qubits)), + [], + approximation_degree=1 - 1e-5, ) result = scc.commute(generic_gate, [0], [], gate, list(range(gate.num_qubits)), []) @@ -463,6 +484,20 @@ def test_2q_pauli_rot_with_non_cached(self): self.assertFalse(scc.commute(something_else, [0], [], RYYGate(np.pi), [1, 0], [])) self.assertFalse(scc.commute(something_else, [1], [], RYYGate(np.pi), [1, 0], [])) + def test_approximation_degree(self): + """Test setting the approximation degree.""" + + almost_identity = RZGate(1e-5) + other = HGate() + + self.assertFalse(scc.commute(almost_identity, [0], [], other, [0], [])) + self.assertFalse( + scc.commute(almost_identity, [0], [], other, [0], [], approximation_degree=1) + ) + self.assertTrue( + scc.commute(almost_identity, [0], [], other, [0], [], approximation_degree=1 - 1e-4) + ) + if __name__ == "__main__": unittest.main() diff --git a/test/python/transpiler/test_remove_identity_equivalent.py b/test/python/transpiler/test_remove_identity_equivalent.py index 736c02b765bd..b902c7d8bff5 100644 --- a/test/python/transpiler/test_remove_identity_equivalent.py +++ b/test/python/transpiler/test_remove_identity_equivalent.py @@ -10,7 +10,7 @@ # copyright notice, and modified files need to carry a notice indicating # that they have been altered from the originals. -"""Tests for the DropNegligible transpiler pass.""" +"""Tests for the RemoveIdentityEquivalent transpiler pass.""" import ddt import numpy as np @@ -37,10 +37,10 @@ @ddt.ddt -class TestDropNegligible(QiskitTestCase): - """Test the DropNegligible pass.""" +class TestRemoveIdentityEquivalent(QiskitTestCase): + """Test the RemoveIdentityEquivalent pass.""" - def test_drops_negligible_gates(self): + def test_remove_identity_equiv_pass(self): """Test that negligible gates are dropped.""" qubits = QuantumRegister(2) circuit = QuantumCircuit(qubits) @@ -188,6 +188,9 @@ def to_matrix(self): UnitaryGate(np.array([[np.exp(1j * np.pi / 4), 0], [0, np.exp(1j * np.pi / 4)]])), GlobalPhaseGate(0), GlobalPhaseGate(np.pi / 4), + UnitaryGate(np.exp(-0.123j) * np.eye(2)), + UnitaryGate(np.exp(-0.123j) * np.eye(4)), + UnitaryGate(np.exp(-0.123j) * np.eye(8)), ) def test_remove_identity_up_to_global_phase(self, gate): """Test that gates equivalent to identity up to a global phase are removed from the circuit,