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 Qiskit#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 Qiskit#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 Qiskit#8774

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

qiskit/quantum_info/synthesis/one_qubit_decompose.py

@@ -14,8 +14,6 @@
 Decompose a single-qubit unitary via Euler angles.
 """
 
-import math
-import cmath
 import numpy as np
 
 from qiskit.circuit.quantumcircuit import QuantumCircuit
@@ -35,6 +33,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 +219,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):
@@ -593,14 +554,3 @@ def _mod_2pi(angle: float, atol: float = 0):
     if abs(wrapped - np.pi) < atol:
         wrapped = -np.pi
     return wrapped
-
-
-def _h_conjugate(su2):
-    """Return su2 conjugated by Hadamard gate. No warning if input matrix is not in su2."""
-    return np.array(
-        [
-            [su2[0, 0].real + 1j * su2[1, 0].imag, 1j * su2[0, 0].imag + su2[1, 0].real],
-            [1j * su2[0, 0].imag - su2[1, 0].real, su2[0, 0].real - 1j * su2[1, 0].imag],
-        ],
-        dtype=complex,
-    )

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>) -> 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<Complex64>) -> [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<Complex64>) -> [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<Complex64>) -> [f64; 4] {
+    let mat = unitary.as_array();
+    params_zyz_inner(mat)
+}
+
+#[pyfunction]
+pub fn params_xyx(unitary: PyReadonlyArray2<Complex64>) -> [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<Complex64>) -> [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<Complex64>) -> [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(())
+}

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(())
 }