Only use rust two qubit weyl coordinates functions

This commit is a follow up to Qiskit#11019 that uses the rust implementation
of computing the weyl coordinates for a unitary and the
transform_to_magic_basis functions everywhere. Previously they were only
used internally by the num_basis_gates() rust function to estimate the
number of basis gates needed for a decomposition of a given unitary.
This will likely be superseded by Qiskit#11946 when that is working, but this
is mainly an incremental step working towards finishin Qiskit#11946 that
proves the internal rust functions work in the larger decomposer code to
isolate the source of failures in Qiskit#11946.

crates/accelerate/src/two_qubit_decompose.rs

@@ -18,15 +18,17 @@
 // of real components and one of imaginary components.
 // In order to avoid copying we want to use `MatRef<c64>` or `MatMut<c64>`.
 
-use num_complex::Complex;
+use num_complex::{Complex, Complex64};
 use pyo3::prelude::*;
 use pyo3::wrap_pyfunction;
 use pyo3::Python;
 use std::f64::consts::PI;
 
 use faer::IntoFaerComplex;
+use faer::IntoNdarrayComplex;
 use faer::{mat, prelude::*, scale, Mat, MatRef};
 use faer_core::{c64, ComplexField};
+use numpy::IntoPyArray;
 use numpy::PyReadonlyArray2;
 
 use crate::utils;
@@ -44,6 +46,16 @@ const C0_5_NEG: c64 = c64 { re: -0.5, im: 0.0 };
 const C0_5_IM: c64 = c64 { re: 0.0, im: 0.5 };
 const C0_5_IM_NEG: c64 = c64 { re: 0.0, im: -0.5 };
 
+// "Magic" basis used for the Weyl decomposition. The basis and its adjoint
+// are stored individually
+// unnormalized, but such that their matrix multiplication is still the
+// identity. This is because
+// they are only used in unitary transformations (so it's safe to do so),
+// and `sqrt(0.5)` is not
+// exactly representable in floating point. Doing it this way means
+// that every element of the matrix
+// is stored exactly correctly, and the multiplication is _exactly_ the
+// identity rather than differing by 1ULP.
 const B_NON_NORMALIZED: [c64; 16] = [
     C1, C1IM, C0, C0, C0, C0, C1IM, C1, C0, C0, C1IM, C1_NEG, C1, C1_NEG_IM, C0, C0,
 ];
@@ -67,10 +79,41 @@ const B_NON_NORMALIZED_DAGGER: [c64; 16] = [
     C0,
 ];
 
-fn transform_from_magic_basis(unitary: Mat<c64>) -> Mat<c64> {
+fn transform_from_magic_basis(unitary: MatRef<c64>, reverse: bool) -> Mat<c64> {
     let _b_nonnormalized = mat::from_row_major_slice::<c64>(&B_NON_NORMALIZED, 4, 4);
     let _b_nonnormalized_dagger = mat::from_row_major_slice::<c64>(&B_NON_NORMALIZED_DAGGER, 4, 4);
-    _b_nonnormalized_dagger * unitary * _b_nonnormalized
+    if reverse {
+        _b_nonnormalized_dagger * unitary * _b_nonnormalized
+    } else {
+        _b_nonnormalized * unitary * _b_nonnormalized_dagger
+    }
+}
+/// Transform the 4-by-4 matrix ``unitary`` into the magic basis.
+///
+/// This method internally uses non-normalized versions of the basis to
+/// minimize the floating-point errors during the transformation
+///
+/// Args:
+///     unitary (np.ndarray): 4-by-4 matrix to transform.
+///     reverse (bool): Whether to do the transformation forwards
+///         (``B @ U @ B.conj().T``, the default) or backwards
+///         (``B.conj().T @ U @ B``).
+///
+/// Returns:
+///     np.ndarray: The transformed 4-by-4 matrix.
+#[pyfunction]
+#[pyo3(signature = (unitary, reverse=false))]
+pub fn transform_to_magic_basis(
+    py: Python,
+    unitary: PyReadonlyArray2<Complex64>,
+    reverse: bool,
+) -> PyObject {
+    transform_from_magic_basis(unitary.as_array().into_faer_complex(), reverse)
+        .as_ref()
+        .into_ndarray_complex()
+        .to_owned()
+        .into_pyarray(py)
+        .into()
 }
 
 // faer::c64 and num_complex::Complex<f64> are both structs
@@ -100,7 +143,7 @@ impl Arg for c64 {
 
 fn __weyl_coordinates(unitary: MatRef<c64>) -> [f64; 3] {
     let uscaled = scale(C1 / unitary.determinant().powf(0.25)) * unitary;
-    let uup = transform_from_magic_basis(uscaled);
+    let uup = transform_from_magic_basis(uscaled.as_ref(), true);
     let mut darg: Vec<_> = (uup.transpose() * &uup)
         .complex_eigenvalues()
         .into_iter()
@@ -147,6 +190,21 @@ fn __weyl_coordinates(unitary: MatRef<c64>) -> [f64; 3] {
     [cs[1], cs[0], cs[2]]
 }
 
+/// Computes the Weyl coordinates for a given two-qubit unitary matrix.
+///
+/// Args:
+///     unitary (np.ndarray): Input two-qubit unitary.
+///
+/// Returns:
+///     np.ndarray: Array of the 3 Weyl coordinates.
+#[pyfunction]
+pub fn weyl_coordinates(py: Python, unitary: PyReadonlyArray2<Complex64>) -> PyObject {
+    __weyl_coordinates(unitary.as_array().into_faer_complex())
+        .to_vec()
+        .into_pyarray(py)
+        .into()
+}
+
 #[pyfunction]
 #[pyo3(text_signature = "(basis_b, basis_fidelity, unitary, /")]
 pub fn _num_basis_gates(
@@ -195,5 +253,7 @@ fn trace_to_fid(trace: &c64) -> f64 {
 #[pymodule]
 pub fn two_qubit_decompose(_py: Python, m: &PyModule) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(_num_basis_gates))?;
+    m.add_wrapped(wrap_pyfunction!(transform_to_magic_basis))?;
+    m.add_wrapped(wrap_pyfunction!(weyl_coordinates))?;
     Ok(())
 }

qiskit/__init__.py

@@ -70,6 +70,7 @@
 sys.modules["qiskit._accelerate.convert_2q_block_matrix"] = (
     qiskit._accelerate.convert_2q_block_matrix
 )
+sys.modules["qiskit._accelerate.two_qubit_decompose"] = qiskit._accelerate.two_qubit_decompose
 
 # qiskit errors operator
 from qiskit.exceptions import QiskitError, MissingOptionalLibraryError

qiskit/synthesis/two_qubit/two_qubit_decompose.py

@@ -238,7 +238,7 @@ def __new__(cls, unitary_matrix, *, fidelity=(1.0 - 1.0e-9), _unpickling=False):
 
         # Find K1, K2 so that U = K1.A.K2, with K being product of single-qubit unitaries
         K1 = transform_to_magic_basis(Up @ P @ np.diag(np.exp(1j * d)))
-        K2 = transform_to_magic_basis(P.T)
+        K2 = transform_to_magic_basis(P.T.astype(complex, copy=False))
 
         K1l, K1r, phase_l = decompose_two_qubit_product_gate(K1)
         K2l, K2r, phase_r = decompose_two_qubit_product_gate(K2)

qiskit/synthesis/two_qubit/weyl.py

@@ -14,84 +14,4 @@
 """Routines that compute  and use the Weyl chamber coordinates.
 """
 
-from __future__ import annotations
-import numpy as np
-
-# "Magic" basis used for the Weyl decomposition. The basis and its adjoint are stored individually
-# unnormalized, but such that their matrix multiplication is still the identity.  This is because
-# they are only used in unitary transformations (so it's safe to do so), and `sqrt(0.5)` is not
-# exactly representable in floating point.  Doing it this way means that every element of the matrix
-# is stored exactly correctly, and the multiplication is _exactly_ the identity rather than
-# differing by 1ULP.
-_B_nonnormalized = np.array([[1, 1j, 0, 0], [0, 0, 1j, 1], [0, 0, 1j, -1], [1, -1j, 0, 0]])
-_B_nonnormalized_dagger = 0.5 * _B_nonnormalized.conj().T
-
-
-def transform_to_magic_basis(U: np.ndarray, reverse: bool = False) -> np.ndarray:
-    """Transform the 4-by-4 matrix ``U`` into the magic basis.
-
-    This method internally uses non-normalized versions of the basis to minimize the floating-point
-    errors that arise during the transformation.
-
-    Args:
-        U (np.ndarray): 4-by-4 matrix to transform.
-        reverse (bool): Whether to do the transformation forwards (``B @ U @ B.conj().T``, the
-        default) or backwards (``B.conj().T @ U @ B``).
-
-    Returns:
-        np.ndarray: The transformed 4-by-4 matrix.
-    """
-    if reverse:
-        return _B_nonnormalized_dagger @ U @ _B_nonnormalized
-    return _B_nonnormalized @ U @ _B_nonnormalized_dagger
-
-
-def weyl_coordinates(U: np.ndarray) -> np.ndarray:
-    """Computes the Weyl coordinates for a given two-qubit unitary matrix.
-
-    Args:
-        U (np.ndarray): Input two-qubit unitary.
-
-    Returns:
-        np.ndarray: Array of the 3 Weyl coordinates.
-    """
-    import scipy.linalg as la
-
-    pi2 = np.pi / 2
-    pi4 = np.pi / 4
-
-    U = U / la.det(U) ** (0.25)
-    Up = transform_to_magic_basis(U, reverse=True)
-    # We only need the eigenvalues of `M2 = Up.T @ Up` here, not the full diagonalization.
-    D = la.eigvals(Up.T @ Up)
-    d = -np.angle(D) / 2
-    d[3] = -d[0] - d[1] - d[2]
-    cs = np.mod((d[:3] + d[3]) / 2, 2 * np.pi)
-
-    # Reorder the eigenvalues to get in the Weyl chamber
-    cstemp = np.mod(cs, pi2)
-    np.minimum(cstemp, pi2 - cstemp, cstemp)
-    order = np.argsort(cstemp)[[1, 2, 0]]
-    cs = cs[order]
-    d[:3] = d[order]
-
-    # Flip into Weyl chamber
-    if cs[0] > pi2:
-        cs[0] -= 3 * pi2
-    if cs[1] > pi2:
-        cs[1] -= 3 * pi2
-    conjs = 0
-    if cs[0] > pi4:
-        cs[0] = pi2 - cs[0]
-        conjs += 1
-    if cs[1] > pi4:
-        cs[1] = pi2 - cs[1]
-        conjs += 1
-    if cs[2] > pi2:
-        cs[2] -= 3 * pi2
-    if conjs == 1:
-        cs[2] = pi2 - cs[2]
-    if cs[2] > pi4:
-        cs[2] -= pi2
-
-    return cs[[1, 0, 2]]
+from qiskit._accelerate.two_qubit_decompose import transform_to_magic_basis, weyl_coordinates

test/python/synthesis/test_weyl.py

@@ -32,7 +32,7 @@ class TestWeyl(QiskitTestCase):
     def test_weyl_coordinates_simple(self):
         """Check Weyl coordinates against known cases."""
         # Identity [0,0,0]
-        U = np.identity(4)
+        U = np.identity(4, dtype=complex)
         weyl = weyl_coordinates(U)
         assert_allclose(weyl, [0, 0, 0])