From b5a141260bf3d5b3c8b3dfa11e9becd32ca1ce99 Mon Sep 17 00:00:00 2001 From: Matthew Treinish Date: Tue, 22 Nov 2022 14:36:26 -0500 Subject: [PATCH] Oxidize parameter calculation for OneQubitEulerDecomposer This commit ports the per basis parameter calculation from python to rust. In looking at the runtime performance regression caused by #8917 the majority is just that we're doing more work synthesizing to more available basis to potentially produce better quality results. Profiling the transpiler pass shows we're spending a non-trivial amount of time in numpy/scipy (depending on whether it's before or after #9179) computing the determinant of the unitary. This is likely because those determinant functions are designed to work with an arbitrarily large square matrix while for the 1 qubit decomposer we're only ever working with a 2x2. To remove this overhead this commit writes a dedicated rust function to compute the determinant of a 2x2 complex matrix and then also adds dedicated functions to calculate the angles for given basis to rust as we can easily just return the end result from rust. Related #8774 --- qiskit/__init__.py | 3 + .../synthesis/one_qubit_decompose.py | 47 +------ src/euler_one_qubit_decomposer.rs | 118 ++++++++++++++++++ src/lib.rs | 4 + 4 files changed, 130 insertions(+), 42 deletions(-) create mode 100644 src/euler_one_qubit_decomposer.rs diff --git a/qiskit/__init__.py b/qiskit/__init__.py index 6006688fdffe..c88a9149eb7e 100644 --- a/qiskit/__init__.py +++ b/qiskit/__init__.py @@ -45,6 +45,9 @@ sys.modules["qiskit._accelerate.sampled_exp_val"] = qiskit._accelerate.sampled_exp_val sys.modules["qiskit._accelerate.vf2_layout"] = qiskit._accelerate.vf2_layout sys.modules["qiskit._accelerate.error_map"] = qiskit._accelerate.error_map +sys.modules[ + "qiskit._accelerate.euler_one_qubit_decomposer" +] = qiskit._accelerate.euler_one_qubit_decomposer # Extend namespace for backwards compat diff --git a/qiskit/quantum_info/synthesis/one_qubit_decompose.py b/qiskit/quantum_info/synthesis/one_qubit_decompose.py index 4c6fbbba0711..81ebb0d906eb 100644 --- a/qiskit/quantum_info/synthesis/one_qubit_decompose.py +++ b/qiskit/quantum_info/synthesis/one_qubit_decompose.py @@ -35,6 +35,7 @@ ) from qiskit.exceptions import QiskitError from qiskit.quantum_info.operators.predicates import is_unitary_matrix +from qiskit._accelerate import euler_one_qubit_decomposer DEFAULT_ATOL = 1e-12 @@ -220,59 +221,21 @@ def angles_and_phase(self, unitary): @staticmethod def _params_zyz(mat): """Return the Euler angles and phase for the ZYZ basis.""" - # We rescale the input matrix to be special unitary (det(U) = 1) - # This ensures that the quaternion representation is real - coeff = np.linalg.det(mat) ** (-0.5) - phase = -cmath.phase(coeff) - su_mat = coeff * mat # U in SU(2) - # OpenQASM SU(2) parameterization: - # U[0, 0] = exp(-i(phi+lambda)/2) * cos(theta/2) - # U[0, 1] = -exp(-i(phi-lambda)/2) * sin(theta/2) - # U[1, 0] = exp(i(phi-lambda)/2) * sin(theta/2) - # U[1, 1] = exp(i(phi+lambda)/2) * cos(theta/2) - theta = 2 * math.atan2(abs(su_mat[1, 0]), abs(su_mat[0, 0])) - phiplambda2 = cmath.phase(su_mat[1, 1]) - phimlambda2 = cmath.phase(su_mat[1, 0]) - phi = phiplambda2 + phimlambda2 - lam = phiplambda2 - phimlambda2 - return theta, phi, lam, phase + return euler_one_qubit_decomposer.params_zyz(mat) @staticmethod def _params_zxz(mat): """Return the Euler angles and phase for the ZXZ basis.""" - theta, phi, lam, phase = OneQubitEulerDecomposer._params_zyz(mat) - return theta, phi + np.pi / 2, lam - np.pi / 2, phase + return euler_one_qubit_decomposer.params_zxz(mat) @staticmethod def _params_xyx(mat): """Return the Euler angles and phase for the XYX basis.""" - # We use the fact that - # Rx(a).Ry(b).Rx(c) = H.Rz(a).Ry(-b).Rz(c).H - mat_zyz = 0.5 * np.array( - [ - [ - mat[0, 0] + mat[0, 1] + mat[1, 0] + mat[1, 1], - mat[0, 0] - mat[0, 1] + mat[1, 0] - mat[1, 1], - ], - [ - mat[0, 0] + mat[0, 1] - mat[1, 0] - mat[1, 1], - mat[0, 0] - mat[0, 1] - mat[1, 0] + mat[1, 1], - ], - ], - dtype=complex, - ) - theta, phi, lam, phase = OneQubitEulerDecomposer._params_zyz(mat_zyz) - newphi, newlam = _mod_2pi(phi + np.pi), _mod_2pi(lam + np.pi) - return theta, newphi, newlam, phase + (newphi + newlam - phi - lam) / 2 + return euler_one_qubit_decomposer.params_xyx(mat) @staticmethod def _params_xzx(umat): - det = np.linalg.det(umat) - phase = (-1j * np.log(det)).real / 2 - mat = umat / np.sqrt(det) - mat_zxz = _h_conjugate(mat) - theta, phi, lam, phase_zxz = OneQubitEulerDecomposer._params_zxz(mat_zxz) - return theta, phi, lam, phase + phase_zxz + return euler_one_qubit_decomposer.params_xzx(umat) @staticmethod def _params_u3(mat): diff --git a/src/euler_one_qubit_decomposer.rs b/src/euler_one_qubit_decomposer.rs new file mode 100644 index 000000000000..07a84c585959 --- /dev/null +++ b/src/euler_one_qubit_decomposer.rs @@ -0,0 +1,118 @@ +// This code is part of Qiskit. +// +// (C) Copyright IBM 2022 +// +// 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 num_complex::{Complex64, ComplexFloat}; +use pyo3::prelude::*; +use pyo3::wrap_pyfunction; +use pyo3::Python; +use std::f64::consts::PI; + +use ndarray::prelude::*; +use numpy::PyReadonlyArray2; + +#[inline] +fn det_one_qubit(mat: ArrayView2) -> Complex64 { + mat[[0, 0]] * mat[[1, 1]] - mat[[0, 1]] * mat[[1, 0]] +} + +#[inline] +fn complex_phase(x: Complex64) -> f64 { + x.im.atan2(x.re) +} + +#[inline] +fn mod_2pi(angle: f64) -> f64 { + (angle + PI) % (2. * PI) - PI +} + +fn params_zyz_inner(mat: ArrayView2) -> [f64; 4] { + let coeff: Complex64 = 1. / det_one_qubit(mat).sqrt(); + let phase = -complex_phase(coeff); + let tmp_1_0 = (coeff * mat[[1, 0]]).abs(); + let tmp_0_0 = (coeff * mat[[0, 0]]).abs(); + let theta = 2. * tmp_1_0.atan2(tmp_0_0); + let phiplambda2 = complex_phase(coeff * mat[[1, 1]]); + let phimlambda2 = complex_phase(coeff * mat[[1, 0]]); + let phi = phiplambda2 + phimlambda2; + let lam = phiplambda2 - phimlambda2; + [theta, phi, lam, phase] +} + +fn params_zxz_inner(mat: ArrayView2) -> [f64; 4] { + let [theta, phi, lam, phase] = params_zyz_inner(mat); + [theta, phi + PI / 2., lam - PI / 2., phase] +} + +#[pyfunction] +pub fn params_zyz(unitary: PyReadonlyArray2) -> [f64; 4] { + let mat = unitary.as_array(); + params_zyz_inner(mat) +} + +#[pyfunction] +pub fn params_xyx(unitary: PyReadonlyArray2) -> [f64; 4] { + let mat = unitary.as_array(); + let mat_zyz = arr2(&[ + [ + 0.5 * (mat[[0, 0]] + mat[[0, 1]] + mat[[1, 0]] + mat[[1, 1]]), + 0.5 * (mat[[0, 0]] - mat[[0, 1]] + mat[[1, 0]] - mat[[1, 1]]), + ], + [ + 0.5 * (mat[[0, 0]] + mat[[0, 1]] - mat[[1, 0]] - mat[[1, 1]]), + 0.5 * (mat[[0, 0]] - mat[[0, 1]] - mat[[1, 0]] + mat[[1, 1]]), + ], + ]); + let [theta, phi, lam, phase] = params_zyz_inner(mat_zyz.view()); + let new_phi = mod_2pi(phi + PI); + let new_lam = mod_2pi(lam + PI); + [ + theta, + new_phi, + new_lam, + phase + (new_phi + new_lam - phi - lam) / 2., + ] +} + +#[pyfunction] +pub fn params_xzx(unitary: PyReadonlyArray2) -> [f64; 4] { + let umat = unitary.as_array(); + let det = det_one_qubit(umat); + let phase = (Complex64::new(0., -1.) * det.ln()).re / 2.; + let sqrt_det = det.sqrt(); + let mat_zyz = arr2(&[ + [ + Complex64::new((umat[[0, 0]] / sqrt_det).re, (umat[[1, 0]] / sqrt_det).im), + Complex64::new((umat[[1, 0]] / sqrt_det).re, (umat[[0, 0]] / sqrt_det).im), + ], + [ + Complex64::new(-(umat[[1, 0]] / sqrt_det).re, (umat[[0, 0]] / sqrt_det).im), + Complex64::new((umat[[0, 0]] / sqrt_det).re, -(umat[[1, 0]] / sqrt_det).im), + ], + ]); + let [theta, phi, lam, phase_zxz] = params_zxz_inner(mat_zyz.view()); + [theta, phi, lam, phase + phase_zxz] +} + +#[pyfunction] +pub fn params_zxz(unitary: PyReadonlyArray2) -> [f64; 4] { + let mat = unitary.as_array(); + params_zxz_inner(mat) +} + +#[pymodule] +pub fn euler_one_qubit_decomposer(_py: Python, m: &PyModule) -> PyResult<()> { + m.add_wrapped(wrap_pyfunction!(params_zyz))?; + m.add_wrapped(wrap_pyfunction!(params_xyx))?; + m.add_wrapped(wrap_pyfunction!(params_xzx))?; + m.add_wrapped(wrap_pyfunction!(params_zxz))?; + Ok(()) +} diff --git a/src/lib.rs b/src/lib.rs index 3f5c5dbea612..b161b8833ce3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -19,6 +19,7 @@ use pyo3::Python; mod dense_layout; mod edge_collections; mod error_map; +mod euler_one_qubit_decomposer; mod nlayout; mod optimize_1q_gates; mod pauli_exp_val; @@ -55,5 +56,8 @@ fn _accelerate(_py: Python<'_>, m: &PyModule) -> PyResult<()> { m.add_wrapped(wrap_pymodule!(optimize_1q_gates::optimize_1q_gates))?; m.add_wrapped(wrap_pymodule!(sampled_exp_val::sampled_exp_val))?; m.add_wrapped(wrap_pymodule!(vf2_layout::vf2_layout))?; + m.add_wrapped(wrap_pymodule!( + euler_one_qubit_decomposer::euler_one_qubit_decomposer + ))?; Ok(()) }