From fea4e0dffe9afff4117d536daca9a1bfa6ffb23a Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Thu, 13 Jun 2024 04:39:47 -0400 Subject: [PATCH 1/5] Oxidize TwoQubitDecomposeUpToDiagonal MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This commit ports the TwoQubitDecomposeUpToDiagonal class from Python to rust. This internal private class is used internally by the quantum shannon decomposition code, and while not performance critical was simple to port. One difference is while the original Python implementation was a class, it acted more like a function in practice. So the new rust version is exposed as a function. Co-authored-by: Luciano Bello Co-authored-by: Elena Peña Tapia <57907331+ElePT@users.noreply.github.com> Co-authored-by: Sebastian Brandhofer <148463728+sbrandhsn@users.noreply.github.com> Co-authored-by: Jake Lishman Co-authored-by: John Lapeyre Co-authored-by: Julien Gacon Co-authored-by: Eli Arbel <46826214+eliarbel@users.noreply.github.com> Co-authored-by: Raynel Sanchez <87539502+raynelfss@users.noreply.github.com> Co-authored-by: Henry Zou <87874865+henryzou50@users.noreply.github.com> Co-authored-by: Shelly Garion <46566946+ShellyGarion@users.noreply.github.com> Co-authored-by: Alexander Ivrii --- crates/accelerate/src/two_qubit_decompose.rs | 391 +++++++++++++------ qiskit/synthesis/unitary/qsd.py | 8 +- 2 files changed, 285 insertions(+), 114 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 5e833bd86fda..b3b0b76ce351 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -550,75 +550,10 @@ impl TwoQubitWeylDecomposition { } } } -} - -static IPZ: [[Complex64; 2]; 2] = [ - [C1_IM, Complex64::new(0., 0.)], - [Complex64::new(0., 0.), Complex64::new(0., -1.)], -]; -static IPY: [[Complex64; 2]; 2] = [ - [Complex64::new(0., 0.), Complex64::new(1., 0.)], - [Complex64::new(-1., 0.), Complex64::new(0., 0.)], -]; -static IPX: [[Complex64; 2]; 2] = [ - [Complex64::new(0., 0.), C1_IM], - [C1_IM, Complex64::new(0., 0.)], -]; - -#[pymethods] -impl TwoQubitWeylDecomposition { - #[staticmethod] - fn _from_state( - angles: [f64; 4], - matrices: [PyReadonlyArray2; 5], - specialization: Specialization, - default_euler_basis: EulerBasis, - calculated_fidelity: f64, - requested_fidelity: Option, - ) -> Self { - let [a, b, c, global_phase] = angles; - Self { - a, - b, - c, - global_phase, - K1l: matrices[0].as_array().to_owned(), - K1r: matrices[1].as_array().to_owned(), - K2l: matrices[2].as_array().to_owned(), - K2r: matrices[3].as_array().to_owned(), - specialization, - default_euler_basis, - calculated_fidelity, - requested_fidelity, - unitary_matrix: matrices[4].as_array().to_owned(), - } - } - - fn __reduce__(&self, py: Python) -> PyResult> { - Ok(( - py.get_type_bound::().getattr("_from_state")?, - ( - [self.a, self.b, self.c, self.global_phase], - [ - self.K1l.to_pyarray_bound(py), - self.K1r.to_pyarray_bound(py), - self.K2l.to_pyarray_bound(py), - self.K2r.to_pyarray_bound(py), - self.unitary_matrix.to_pyarray_bound(py), - ], - self.specialization, - self.default_euler_basis, - self.calculated_fidelity, - self.requested_fidelity, - ), - ) - .into_py(py)) - } - - #[new] - #[pyo3(signature=(unitary_matrix, fidelity=DEFAULT_FIDELITY, _specialization=None))] - fn new( - unitary_matrix: PyReadonlyArray2, + /// Instantiate a new TwoQubitWeylDecomposition with rust native + /// data structures + fn new_inner( + unitary_matrix: ArrayView2, fidelity: Option, _specialization: Option, ) -> PyResult { @@ -626,8 +561,8 @@ impl TwoQubitWeylDecomposition { let ipy: ArrayView2 = aview2(&IPY); let ipx: ArrayView2 = aview2(&IPX); - let mut u = unitary_matrix.as_array().to_owned(); - let unitary_matrix = unitary_matrix.as_array().to_owned(); + let mut u = unitary_matrix.to_owned(); + let unitary_matrix = unitary_matrix.to_owned(); let det_u = u.view().into_faer_complex().determinant().to_num_complex(); let det_pow = det_u.powf(-0.25); u.mapv_inplace(|x| x * det_pow); @@ -1093,6 +1028,80 @@ impl TwoQubitWeylDecomposition { specialized.global_phase += tr.arg(); Ok(specialized) } +} + +static IPZ: [[Complex64; 2]; 2] = [ + [C1_IM, Complex64::new(0., 0.)], + [Complex64::new(0., 0.), Complex64::new(0., -1.)], +]; +static IPY: [[Complex64; 2]; 2] = [ + [Complex64::new(0., 0.), Complex64::new(1., 0.)], + [Complex64::new(-1., 0.), Complex64::new(0., 0.)], +]; +static IPX: [[Complex64; 2]; 2] = [ + [Complex64::new(0., 0.), C1_IM], + [C1_IM, Complex64::new(0., 0.)], +]; + +#[pymethods] +impl TwoQubitWeylDecomposition { + #[staticmethod] + fn _from_state( + angles: [f64; 4], + matrices: [PyReadonlyArray2; 5], + specialization: Specialization, + default_euler_basis: EulerBasis, + calculated_fidelity: f64, + requested_fidelity: Option, + ) -> Self { + let [a, b, c, global_phase] = angles; + Self { + a, + b, + c, + global_phase, + K1l: matrices[0].as_array().to_owned(), + K1r: matrices[1].as_array().to_owned(), + K2l: matrices[2].as_array().to_owned(), + K2r: matrices[3].as_array().to_owned(), + specialization, + default_euler_basis, + calculated_fidelity, + requested_fidelity, + unitary_matrix: matrices[4].as_array().to_owned(), + } + } + + fn __reduce__(&self, py: Python) -> PyResult> { + Ok(( + py.get_type_bound::().getattr("_from_state")?, + ( + [self.a, self.b, self.c, self.global_phase], + [ + self.K1l.to_pyarray_bound(py), + self.K1r.to_pyarray_bound(py), + self.K2l.to_pyarray_bound(py), + self.K2r.to_pyarray_bound(py), + self.unitary_matrix.to_pyarray_bound(py), + ], + self.specialization, + self.default_euler_basis, + self.calculated_fidelity, + self.requested_fidelity, + ), + ) + .into_py(py)) + } + + #[new] + #[pyo3(signature=(unitary_matrix, fidelity=DEFAULT_FIDELITY, _specialization=None))] + fn new( + unitary_matrix: PyReadonlyArray2, + fidelity: Option, + _specialization: Option, + ) -> PyResult { + TwoQubitWeylDecomposition::new_inner(unitary_matrix.as_array(), fidelity, _specialization) + } #[allow(non_snake_case)] #[getter] @@ -1647,55 +1656,17 @@ impl TwoQubitBasisDecomposer { } Ok(res) } -} - -static K12R_ARR: [[Complex64; 2]; 2] = [ - [ - Complex64::new(0., FRAC_1_SQRT_2), - Complex64::new(FRAC_1_SQRT_2, 0.), - ], - [ - Complex64::new(-FRAC_1_SQRT_2, 0.), - Complex64::new(0., -FRAC_1_SQRT_2), - ], -]; - -static K12L_ARR: [[Complex64; 2]; 2] = [ - [Complex64::new(0.5, 0.5), Complex64::new(0.5, 0.5)], - [Complex64::new(-0.5, 0.5), Complex64::new(0.5, -0.5)], -]; -fn decomp0_inner(target: &TwoQubitWeylDecomposition) -> SmallVec<[Array2; 8]> { - smallvec![target.K1r.dot(&target.K2r), target.K1l.dot(&target.K2l),] -} - -#[pymethods] -impl TwoQubitBasisDecomposer { - fn __getnewargs__(&self, py: Python) -> (String, PyObject, f64, &str, Option) { - ( - self.gate.clone(), - self.basis_decomposer - .unitary_matrix - .to_pyarray_bound(py) - .into(), - self.basis_fidelity, - self.euler_basis.as_str(), - self.pulse_optimize, - ) - } - - #[new] - #[pyo3(signature=(gate, gate_matrix, basis_fidelity=1.0, euler_basis="U", pulse_optimize=None))] - fn new( + fn new_inner( gate: String, - gate_matrix: PyReadonlyArray2, + gate_matrix: ArrayView2, basis_fidelity: f64, euler_basis: &str, pulse_optimize: Option, ) -> PyResult { let ipz: ArrayView2 = aview2(&IPZ); let basis_decomposer = - TwoQubitWeylDecomposition::new(gate_matrix, Some(DEFAULT_FIDELITY), None)?; + TwoQubitWeylDecomposition::new_inner(gate_matrix, Some(DEFAULT_FIDELITY), None)?; let super_controlled = relative_eq!(basis_decomposer.a, PI4, max_relative = 1e-09) && relative_eq!(basis_decomposer.c, 0.0, max_relative = 1e-09); @@ -1843,6 +1814,158 @@ impl TwoQubitBasisDecomposer { }) } + fn call_inner( + &self, + unitary: ArrayView2, + basis_fidelity: Option, + approximate: bool, + _num_basis_uses: Option, + ) -> PyResult { + let basis_fidelity = if !approximate { + 1.0 + } else { + basis_fidelity.unwrap_or(self.basis_fidelity) + }; + let target_decomposed = + TwoQubitWeylDecomposition::new_inner(unitary, Some(DEFAULT_FIDELITY), None)?; + let traces = self.traces(&target_decomposed); + let best_nbasis = traces + .into_iter() + .enumerate() + .map(|(idx, trace)| (idx, trace.trace_to_fid() * basis_fidelity.powi(idx as i32))) + .min_by(|(_idx1, fid1), (_idx2, fid2)| fid2.partial_cmp(fid1).unwrap()) + .unwrap() + .0; + let best_nbasis = _num_basis_uses.unwrap_or(best_nbasis as u8); + let decomposition = match best_nbasis { + 0 => decomp0_inner(&target_decomposed), + 1 => self.decomp1_inner(&target_decomposed), + 2 => self.decomp2_supercontrolled_inner(&target_decomposed), + 3 => self.decomp3_supercontrolled_inner(&target_decomposed), + _ => unreachable!("Invalid basis to use"), + }; + let pulse_optimize = self.pulse_optimize.unwrap_or(true); + let sequence = if pulse_optimize { + self.pulse_optimal_chooser(best_nbasis, &decomposition, &target_decomposed)? + } else { + None + }; + if let Some(seq) = sequence { + return Ok(seq); + } + let target_1q_basis_list = vec![self.euler_basis]; + let euler_decompositions: SmallVec<[Option; 8]> = decomposition + .iter() + .map(|decomp| { + unitary_to_gate_sequence_inner( + decomp.view(), + &target_1q_basis_list, + 0, + None, + true, + None, + ) + }) + .collect(); + // Worst case length is 5x 1q gates for each 1q decomposition + 1x 2q gate + // We might overallocate a bit if the euler basis is different but + // the worst case is just 16 extra elements with just a String and 2 smallvecs + // each. This is only transient though as the circuit sequences aren't long lived + // and are just used to create a QuantumCircuit or DAGCircuit when we return to + // Python space. + let mut gates = Vec::with_capacity(21); + let mut global_phase = target_decomposed.global_phase; + global_phase -= best_nbasis as f64 * self.basis_decomposer.global_phase; + if best_nbasis == 2 { + global_phase += PI; + } + for i in 0..best_nbasis as usize { + if let Some(euler_decomp) = &euler_decompositions[2 * i] { + for gate in &euler_decomp.gates { + gates.push((gate.0.clone(), gate.1.clone(), smallvec![0])); + } + global_phase += euler_decomp.global_phase + } + if let Some(euler_decomp) = &euler_decompositions[2 * i + 1] { + for gate in &euler_decomp.gates { + gates.push((gate.0.clone(), gate.1.clone(), smallvec![1])); + } + global_phase += euler_decomp.global_phase + } + gates.push((self.gate.clone(), smallvec![], smallvec![0, 1])); + } + if let Some(euler_decomp) = &euler_decompositions[2 * best_nbasis as usize] { + for gate in &euler_decomp.gates { + gates.push((gate.0.clone(), gate.1.clone(), smallvec![0])); + } + global_phase += euler_decomp.global_phase + } + if let Some(euler_decomp) = &euler_decompositions[2 * best_nbasis as usize + 1] { + for gate in &euler_decomp.gates { + gates.push((gate.0.clone(), gate.1.clone(), smallvec![1])); + } + global_phase += euler_decomp.global_phase + } + Ok(TwoQubitGateSequence { + gates, + global_phase, + }) + } +} + +static K12R_ARR: [[Complex64; 2]; 2] = [ + [ + Complex64::new(0., FRAC_1_SQRT_2), + Complex64::new(FRAC_1_SQRT_2, 0.), + ], + [ + Complex64::new(-FRAC_1_SQRT_2, 0.), + Complex64::new(0., -FRAC_1_SQRT_2), + ], +]; + +static K12L_ARR: [[Complex64; 2]; 2] = [ + [Complex64::new(0.5, 0.5), Complex64::new(0.5, 0.5)], + [Complex64::new(-0.5, 0.5), Complex64::new(0.5, -0.5)], +]; + +fn decomp0_inner(target: &TwoQubitWeylDecomposition) -> SmallVec<[Array2; 8]> { + smallvec![target.K1r.dot(&target.K2r), target.K1l.dot(&target.K2l),] +} + +#[pymethods] +impl TwoQubitBasisDecomposer { + fn __getnewargs__(&self, py: Python) -> (String, PyObject, f64, &str, Option) { + ( + self.gate.clone(), + self.basis_decomposer + .unitary_matrix + .to_pyarray_bound(py) + .into(), + self.basis_fidelity, + self.euler_basis.as_str(), + self.pulse_optimize, + ) + } + + #[new] + #[pyo3(signature=(gate, gate_matrix, basis_fidelity=1.0, euler_basis="U", pulse_optimize=None))] + fn new( + gate: String, + gate_matrix: PyReadonlyArray2, + basis_fidelity: f64, + euler_basis: &str, + pulse_optimize: Option, + ) -> PyResult { + TwoQubitBasisDecomposer::new_inner( + gate, + gate_matrix.as_array(), + basis_fidelity, + euler_basis, + pulse_optimize, + ) + } + fn traces(&self, target: &TwoQubitWeylDecomposition) -> [Complex64; 4] { [ 4. * Complex64::new( @@ -2042,9 +2165,53 @@ impl TwoQubitBasisDecomposer { } } +fn u4_to_su4(u4: ArrayView2) -> (Array2, f64) { + let det_u = u4.into_faer_complex().determinant().to_num_complex(); + let phase_factor = det_u.powf(-0.25).conj(); + let su4 = u4.mapv(|x| x / phase_factor); + (su4, phase_factor.arg()) +} + +fn real_trace_transform(mat: ArrayView2) -> Array2 { + let a1 = -mat[[1, 3]] * mat[[2, 0]] + mat[[1, 2]] * mat[[2, 1]] + mat[[1, 1]] * mat[[2, 2]] + - mat[[1, 0]] * mat[[2, 3]]; + let a2 = mat[[0, 3]] * mat[[3, 0]] - mat[[0, 2]] * mat[[3, 1]] - mat[[0, 1]] * mat[[3, 2]] + + mat[[0, 0]] * mat[[3, 3]]; + let theta = 0.; // Arbitrary! + let phi = 0.; // This is extra arbitrary! + let psi = f64::atan2(a1.im + a2.im, a1.re - a2.re) - phi; + let im = Complex64::new(0., -1.); + let temp = [ + (theta * im).exp(), + (phi * im).exp(), + (psi * im).exp(), + (-(theta + phi + psi) * im).exp(), + ]; + Array2::from_diag(&arr1(&temp)) +} + +#[pyfunction] +fn two_qubit_decompose_up_to_diagonal( + py: Python, + mat: PyReadonlyArray2, +) -> PyResult<(PyObject, TwoQubitGateSequence)> { + let mat_arr: ArrayView2 = mat.as_array(); + let (su4, phase) = u4_to_su4(mat_arr); + let mut real_map = real_trace_transform(su4.view()); + let mapped_su4 = real_map.dot(&su4.view()); + let decomp = + TwoQubitBasisDecomposer::new_inner("cx".to_string(), aview2(&CXGATE), 1.0, "U", None)?; + + let mut circ = decomp.call_inner(mapped_su4.view(), None, true, None)?; + circ.global_phase += phase; + real_map.mapv_inplace(|x| x.conj()); + Ok((real_map.into_pyarray_bound(py).into(), circ)) +} + #[pymodule] pub fn two_qubit_decompose(m: &Bound) -> PyResult<()> { m.add_wrapped(wrap_pyfunction!(_num_basis_gates))?; + m.add_wrapped(wrap_pyfunction!(two_qubit_decompose_up_to_diagonal))?; m.add_class::()?; m.add_class::()?; m.add_class::()?; diff --git a/qiskit/synthesis/unitary/qsd.py b/qiskit/synthesis/unitary/qsd.py index 80a8afc1311b..d3cde37de653 100644 --- a/qiskit/synthesis/unitary/qsd.py +++ b/qiskit/synthesis/unitary/qsd.py @@ -30,6 +30,7 @@ from qiskit.circuit.library.generalized_gates.uc_pauli_rot import UCPauliRotGate, _EPS from qiskit.circuit.library.generalized_gates.ucry import UCRYGate from qiskit.circuit.library.generalized_gates.ucrz import UCRZGate +from qiskit._accelerate.two_qubit_decompose import two_qubit_decompose_up_to_diagonal def qs_decomposition( @@ -253,7 +254,7 @@ def _apply_a2(circ): from qiskit.quantum_info import Operator from qiskit.circuit.library.generalized_gates.unitary import UnitaryGate - decomposer = two_qubit_decompose.TwoQubitDecomposeUpToDiagonal() + decomposer = two_qubit_decompose_up_to_diagonal ccirc = transpile(circ, basis_gates=["u", "cx", "qsd2q"], optimization_level=0) ind2q = [] # collect 2q instrs @@ -275,7 +276,10 @@ def _apply_a2(circ): instr2 = ccirc.data[ind2] mat2 = Operator(instr2.operation).data # rollover - dmat, qc2cx = decomposer(mat1) + dmat, qc2cx_seq = decomposer(mat1) + qc2cx = QuantumCircuit(2, global_phase=qc2cx_seq.global_phase) + for gate, param, qargs in qc2cx_seq: + getattr(qc2cx, gate)(*param, *qargs) ccirc.data[ind1] = instr1.replace(operation=qc2cx.to_gate()) mat2 = mat2 @ dmat ccirc.data[ind2] = instr2.replace(UnitaryGate(mat2)) From d8d1d899b7a6fed18a096b087a8d0998f02f0345 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Fri, 14 Jun 2024 01:26:59 -0400 Subject: [PATCH 2/5] Build circuit from rust Since #12459 recently merged we now have a mechanism to build a circuit from rust. This commit updates the synthesis function to build the circuit directly in rust instead of returning a circuit sequence and building the circuit from Python. This should speed up the construction substantially. --- crates/accelerate/src/two_qubit_decompose.rs | 34 ++++++++++++++++---- qiskit/synthesis/unitary/qsd.py | 6 ++-- 2 files changed, 30 insertions(+), 10 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index b2502b217c60..9364a1ad8ae1 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -51,8 +51,11 @@ use rand::prelude::*; use rand_distr::StandardNormal; use rand_pcg::Pcg64Mcg; +use qiskit_circuit::circuit_data::CircuitData; use qiskit_circuit::gate_matrix::{CX_GATE, H_GATE, ONE_QUBIT_IDENTITY, SX_GATE, X_GATE}; -use qiskit_circuit::SliceOrInt; +use qiskit_circuit::operations::Param; +use qiskit_circuit::operations::StandardGate; +use qiskit_circuit::{Qubit, SliceOrInt}; const PI2: f64 = PI / 2.0; const PI4: f64 = PI / 4.0; @@ -2142,16 +2145,35 @@ fn real_trace_transform(mat: ArrayView2) -> Array2 { fn two_qubit_decompose_up_to_diagonal( py: Python, mat: PyReadonlyArray2, -) -> PyResult<(PyObject, TwoQubitGateSequence)> { +) -> PyResult<(PyObject, CircuitData)> { let mat_arr: ArrayView2 = mat.as_array(); let (su4, phase) = u4_to_su4(mat_arr); let mut real_map = real_trace_transform(su4.view()); let mapped_su4 = real_map.dot(&su4.view()); let decomp = - TwoQubitBasisDecomposer::new_inner("cx".to_string(), aview2(&CXGATE), 1.0, "U", None)?; - - let mut circ = decomp.call_inner(mapped_su4.view(), None, true, None)?; - circ.global_phase += phase; + TwoQubitBasisDecomposer::new_inner("cx".to_string(), aview2(&CX_GATE), 1.0, "U", None)?; + + let circ_seq = decomp.call_inner(mapped_su4.view(), None, true, None)?; + let circ = CircuitData::from_standard_gates( + py, + 2, + circ_seq + .gates + .into_iter() + .map(|(gate, param_floats, qubit_index)| { + let params: SmallVec<[Param; 3]> = + param_floats.into_iter().map(Param::Float).collect(); + let qubits: SmallVec<[Qubit; 2]> = + qubit_index.into_iter().map(|x| Qubit(x as u32)).collect(); + let gate = if gate == "cx" { + StandardGate::CXGate + } else { + StandardGate::UGate + }; + (gate, params, qubits) + }), + Param::Float(circ_seq.global_phase + phase), + )?; real_map.mapv_inplace(|x| x.conj()); Ok((real_map.into_pyarray_bound(py).into(), circ)) } diff --git a/qiskit/synthesis/unitary/qsd.py b/qiskit/synthesis/unitary/qsd.py index d3cde37de653..5682591a2358 100644 --- a/qiskit/synthesis/unitary/qsd.py +++ b/qiskit/synthesis/unitary/qsd.py @@ -276,10 +276,8 @@ def _apply_a2(circ): instr2 = ccirc.data[ind2] mat2 = Operator(instr2.operation).data # rollover - dmat, qc2cx_seq = decomposer(mat1) - qc2cx = QuantumCircuit(2, global_phase=qc2cx_seq.global_phase) - for gate, param, qargs in qc2cx_seq: - getattr(qc2cx, gate)(*param, *qargs) + dmat, qc2cx_data = decomposer(mat1) + qc2cx = QuantumCircuit._from_circuit_data(qc2cx_data) ccirc.data[ind1] = instr1.replace(operation=qc2cx.to_gate()) mat2 = mat2 @ dmat ccirc.data[ind2] = instr2.replace(UnitaryGate(mat2)) From c56f8af531f24392859931e8cf0fc138c1058cbe Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Fri, 14 Jun 2024 01:37:48 -0400 Subject: [PATCH 3/5] Remove unused private Python class This commit removes the Python implementation of the function. This is now unused in Qiskit and was never a public class so nothing external should be depending on it. Since it's not used we should just remove it. --- .../two_qubit/two_qubit_decompose.py | 89 ----------- test/python/synthesis/test_synthesis.py | 140 +----------------- 2 files changed, 7 insertions(+), 222 deletions(-) diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index 41ba75c6b237..a189ef5cf2c2 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -691,95 +691,6 @@ def traces(self, target): return self._inner_decomposer.traces(target._inner_decomposition) -class TwoQubitDecomposeUpToDiagonal: - """ - Class to decompose two qubit unitaries into the product of a diagonal gate - and another unitary gate which can be represented by two CX gates instead of the - usual three. This can be used when neighboring gates commute with the diagonal to - potentially reduce overall CX count. - """ - - def __init__(self): - sy = np.array([[0, -1j], [1j, 0]]) - self.sysy = np.kron(sy, sy) - - def _u4_to_su4(self, u4): - phase_factor = np.conj(np.linalg.det(u4) ** (-1 / u4.shape[0])) - su4 = u4 / phase_factor - return su4, cmath.phase(phase_factor) - - def _gamma(self, mat): - """ - proposition II.1: this invariant characterizes when two operators in U(4), - say u, v, are equivalent up to single qubit gates: - - u ≡ v -> Det(γ(u)) = Det(±(γ(v))) - """ - sumat, _ = self._u4_to_su4(mat) - sysy = self.sysy - return sumat @ sysy @ sumat.T @ sysy - - def _cx0_test(self, mat): - # proposition III.1: zero cx sufficient - gamma = self._gamma(mat) - evals = np.linalg.eigvals(gamma) - return np.all(np.isclose(evals, np.ones(4))) - - def _cx1_test(self, mat): - # proposition III.2: one cx sufficient - gamma = self._gamma(mat) - evals = np.linalg.eigvals(gamma) - uvals, ucnts = np.unique(np.round(evals, 10), return_counts=True) - return ( - len(uvals) == 2 - and all(ucnts == 2) - and all((np.isclose(x, 1j)) or np.isclose(x, -1j) for x in uvals) - ) - - def _cx2_test(self, mat): - # proposition III.3: two cx sufficient - gamma = self._gamma(mat) - return np.isclose(np.trace(gamma).imag, 0) - - def _real_trace_transform(self, mat): - """ - Determine diagonal gate such that - - U3 = D U2 - - Where U3 is a general two-qubit gate which takes 3 cnots, D is a - diagonal gate, and U2 is a gate which takes 2 cnots. - """ - a1 = ( - -mat[1, 3] * mat[2, 0] - + mat[1, 2] * mat[2, 1] - + mat[1, 1] * mat[2, 2] - - mat[1, 0] * mat[2, 3] - ) - a2 = ( - mat[0, 3] * mat[3, 0] - - mat[0, 2] * mat[3, 1] - - mat[0, 1] * mat[3, 2] - + mat[0, 0] * mat[3, 3] - ) - theta = 0 # arbitrary - phi = 0 # arbitrary - psi = np.arctan2(a1.imag + a2.imag, a1.real - a2.real) - phi - diag = np.diag(np.exp(-1j * np.array([theta, phi, psi, -(theta + phi + psi)]))) - return diag - - def __call__(self, mat): - """do the decomposition""" - su4, phase = self._u4_to_su4(mat) - real_map = self._real_trace_transform(su4) - mapped_su4 = real_map @ su4 - if not self._cx2_test(mapped_su4): - warnings.warn("Unitary decomposition up to diagonal may use an additionl CX gate.") - circ = two_qubit_cnot_decompose(mapped_su4) - circ.global_phase += phase - return real_map.conj(), circ - - # This weird duplicated lazy structure is for backwards compatibility; Qiskit has historically # always made ``two_qubit_cnot_decompose`` available publicly immediately on import, but it's quite # expensive to construct, and we want to defer the obejct's creation until it's actually used. We diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index 8d953ff1b83c..6f001a6fbfda 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -62,8 +62,8 @@ TwoQubitControlledUDecomposer, Ud, decompose_two_qubit_product_gate, - TwoQubitDecomposeUpToDiagonal, ) +from qiskit._accelerate.two_qubit_decompose import two_qubit_decompose_up_to_diagonal from qiskit._accelerate.two_qubit_decompose import Specialization from qiskit.synthesis.unitary import qsd from test import combine # pylint: disable=wrong-import-order @@ -1673,149 +1673,23 @@ def test_a2_opt_single_2q(self): class TestTwoQubitDecomposeUpToDiagonal(QiskitTestCase): """test TwoQubitDecomposeUpToDiagonal class""" - def test_prop31(self): - """test proposition III.1: no CNOTs needed""" - dec = TwoQubitDecomposeUpToDiagonal() - # test identity - mat = np.identity(4) - self.assertTrue(dec._cx0_test(mat)) - - sz = np.array([[1, 0], [0, -1]]) - zz = np.kron(sz, sz) - self.assertTrue(dec._cx0_test(zz)) - - had = np.matrix([[1, 1], [1, -1]]) / np.sqrt(2) - hh = np.kron(had, had) - self.assertTrue(dec._cx0_test(hh)) - - sy = np.array([[0, -1j], [1j, 0]]) - hy = np.kron(had, sy) - self.assertTrue(dec._cx0_test(hy)) - - qc = QuantumCircuit(2) - qc.cx(0, 1) - self.assertFalse(dec._cx0_test(Operator(qc).data)) - - def test_prop32_true(self): - """test proposition III.2: 1 CNOT sufficient""" - dec = TwoQubitDecomposeUpToDiagonal() - qc = QuantumCircuit(2) - qc.ry(np.pi / 4, 0) - qc.ry(np.pi / 3, 1) - qc.cx(0, 1) - qc.ry(np.pi / 4, 0) - qc.y(1) - mat = Operator(qc).data - self.assertTrue(dec._cx1_test(mat)) - - qc = QuantumCircuit(2) - qc.ry(np.pi / 5, 0) - qc.ry(np.pi / 3, 1) - qc.cx(1, 0) - qc.ry(np.pi / 2, 0) - qc.y(1) - mat = Operator(qc).data - self.assertTrue(dec._cx1_test(mat)) - - # this SU4 is non-local - mat = scipy.stats.unitary_group.rvs(4, random_state=84) - self.assertFalse(dec._cx1_test(mat)) - - def test_prop32_false(self): - """test proposition III.2: 1 CNOT not sufficient""" - dec = TwoQubitDecomposeUpToDiagonal() - qc = QuantumCircuit(2) - qc.ry(np.pi / 4, 0) - qc.ry(np.pi / 3, 1) - qc.cx(0, 1) - qc.ry(np.pi / 4, 0) - qc.y(1) - qc.cx(0, 1) - qc.ry(np.pi / 3, 0) - qc.rx(np.pi / 2, 1) - mat = Operator(qc).data - self.assertFalse(dec._cx1_test(mat)) - - def test_prop33_true(self): - """test proposition III.3: 2 CNOT sufficient""" - dec = TwoQubitDecomposeUpToDiagonal() - qc = QuantumCircuit(2) - qc.rx(np.pi / 4, 0) - qc.ry(np.pi / 2, 1) - qc.cx(0, 1) - qc.rx(np.pi / 4, 0) - qc.ry(np.pi / 2, 1) - qc.cx(0, 1) - qc.rx(np.pi / 4, 0) - qc.y(1) - mat = Operator(qc).data - self.assertTrue(dec._cx2_test(mat)) - - def test_prop33_false(self): - """test whether circuit which requires 3 cx fails 2 cx test""" - dec = TwoQubitDecomposeUpToDiagonal() - qc = QuantumCircuit(2) - qc.u(0.1, 0.2, 0.3, 0) - qc.u(0.4, 0.5, 0.6, 1) - qc.cx(0, 1) - qc.u(0.1, 0.2, 0.3, 0) - qc.u(0.4, 0.5, 0.6, 1) - qc.cx(0, 1) - qc.u(0.5, 0.2, 0.3, 0) - qc.u(0.2, 0.4, 0.1, 1) - qc.cx(1, 0) - qc.u(0.1, 0.2, 0.3, 0) - qc.u(0.4, 0.5, 0.6, 1) - mat = Operator(qc).data - self.assertFalse(dec._cx2_test(mat)) - - def test_ortho_local_map(self): - """test map of SO4 to SU2⊗SU2""" - dec = TwoQubitDecomposeUpToDiagonal() - emap = np.array([[1, 1j, 0, 0], [0, 0, 1j, 1], [0, 0, 1j, -1], [1, -1j, 0, 0]]) / math.sqrt( - 2 - ) - so4 = scipy.stats.ortho_group.rvs(4, random_state=284) - sy = np.array([[0, -1j], [1j, 0]]) - self.assertTrue(np.allclose(-np.kron(sy, sy), emap @ emap.T)) - self.assertFalse(dec._cx0_test(so4)) - self.assertTrue(dec._cx0_test(emap @ so4 @ emap.T.conj())) - - def test_ortho_local_map2(self): - """test map of SO4 to SU2⊗SU2""" - dec = TwoQubitDecomposeUpToDiagonal() - emap = np.array([[1, 0, 0, 1j], [0, 1j, 1, 0], [0, 1j, -1, 0], [1, 0, 0, -1j]]) / math.sqrt( - 2 - ) - so4 = scipy.stats.ortho_group.rvs(4, random_state=284) - sy = np.array([[0, -1j], [1j, 0]]) - self.assertTrue(np.allclose(-np.kron(sy, sy), emap @ emap.T)) - self.assertFalse(dec._cx0_test(so4)) - self.assertTrue(dec._cx0_test(emap @ so4 @ emap.T.conj())) - - def test_real_trace_transform(self): - """test finding diagonal factor of unitary""" - dec = TwoQubitDecomposeUpToDiagonal() - u4 = scipy.stats.unitary_group.rvs(4, random_state=83) - su4, _ = dec._u4_to_su4(u4) - real_map = dec._real_trace_transform(su4) - self.assertTrue(dec._cx2_test(real_map @ su4)) - def test_call_decompose(self): """ test __call__ method to decompose """ - dec = TwoQubitDecomposeUpToDiagonal() + dec = two_qubit_decompose_up_to_diagonal u4 = scipy.stats.unitary_group.rvs(4, random_state=47) - dmat, circ2cx = dec(u4) + dmat, circ2cx_data = dec(u4) + circ2cx = QuantumCircuit._from_circuit_data(circ2cx_data) dec_diag = dmat @ Operator(circ2cx).data self.assertTrue(Operator(u4) == Operator(dec_diag)) def test_circuit_decompose(self): """test applying decomposed gates as circuit elements""" - dec = TwoQubitDecomposeUpToDiagonal() + dec = two_qubit_decompose_up_to_diagonal u4 = scipy.stats.unitary_group.rvs(4, random_state=47) - dmat, circ2cx = dec(u4) + dmat, circ2cx_data = dec(u4) + circ2cx = QuantumCircuit._from_circuit_data(circ2cx_data) qc1 = QuantumCircuit(2) qc1.append(UnitaryGate(u4), range(2)) From cf422e02d096fe397e0f808bdc6c5c6552ee2116 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Wed, 26 Jun 2024 10:17:13 -0400 Subject: [PATCH 4/5] Remove unused import --- test/python/synthesis/test_synthesis.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index 4fd8145c35da..05739db249b7 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -16,7 +16,6 @@ import unittest import contextlib import logging -import math import numpy as np import scipy import scipy.stats From 655ce1d677d0f66265b8dd9a1e0631c31f2c26f9 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Wed, 24 Jul 2024 12:38:45 -0400 Subject: [PATCH 5/5] Calculate best_nbasis in unwrap_or_else() --- crates/accelerate/src/two_qubit_decompose.rs | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 781545e0351b..dfee4c5b0bff 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -1733,14 +1733,15 @@ impl TwoQubitBasisDecomposer { let target_decomposed = TwoQubitWeylDecomposition::new_inner(unitary, Some(DEFAULT_FIDELITY), None)?; let traces = self.traces(&target_decomposed); - let best_nbasis = traces - .into_iter() - .enumerate() - .map(|(idx, trace)| (idx, trace.trace_to_fid() * basis_fidelity.powi(idx as i32))) - .min_by(|(_idx1, fid1), (_idx2, fid2)| fid2.partial_cmp(fid1).unwrap()) - .unwrap() - .0; - let best_nbasis = _num_basis_uses.unwrap_or(best_nbasis as u8); + let best_nbasis = _num_basis_uses.unwrap_or_else(|| { + traces + .into_iter() + .enumerate() + .map(|(idx, trace)| (idx, trace.trace_to_fid() * basis_fidelity.powi(idx as i32))) + .min_by(|(_idx1, fid1), (_idx2, fid2)| fid2.partial_cmp(fid1).unwrap()) + .unwrap() + .0 as u8 + }); let decomposition = match best_nbasis { 0 => decomp0_inner(&target_decomposed), 1 => self.decomp1_inner(&target_decomposed),