Oxidize parameter calculation for OneQubitEulerDecomposer (#9185)

* 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

* Eliminate python function for staticmethod definition

This commit removes one layer of function calls for the
OneQubitEulerDecomposer's staticmethods for calculating parameters.
Previously, there was a python function which called an inner rust
function, this eliminates one layer and attaches the rust function
directly as a static method to the class definition.

qiskit/__init__.py

@@ -45,6 +45,9 @@ sys.modules["qiskit._accelerate.optimize_1q_gates"] = qiskit._accelerate.optimiz
 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.circuit.library.standard_gates import (
 )
 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
 
@@ -217,62 +216,10 @@ class OneQubitEulerDecomposer:
         """
         return self._params(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
-
-    @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
-
-    @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
-
-    @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
+    _params_zyz = staticmethod(euler_one_qubit_decomposer.params_zyz)
+    _params_zxz = staticmethod(euler_one_qubit_decomposer.params_zxz)
+    _params_xyx = staticmethod(euler_one_qubit_decomposer.params_xyx)
+    _params_xzx = staticmethod(euler_one_qubit_decomposer.params_xzx)
 
     @staticmethod
     def _params_u3(mat):
@@ -593,14 +540,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(())
 }