Oxidize two qubit local invariance functions

crates/accelerate/src/two_qubit_decompose.rs

@@ -31,8 +31,8 @@ use faer_ext::{IntoFaer, IntoFaerComplex, IntoNdarray, IntoNdarrayComplex};
 use ndarray::linalg::kron;
 use ndarray::prelude::*;
 use ndarray::Zip;
-use numpy::PyReadonlyArray2;
 use numpy::{IntoPyArray, ToPyArray};
+use numpy::{PyReadonlyArray1, PyReadonlyArray2};
 
 use pyo3::exceptions::PyValueError;
 use pyo3::prelude::*;
@@ -1883,9 +1883,110 @@ impl TwoQubitBasisDecomposer {
     }
 }
 
+static MAGIC: GateArray2Q = [
+    [
+        c64(FRAC_1_SQRT_2, 0.),
+        C_ZERO,
+        C_ZERO,
+        c64(0., FRAC_1_SQRT_2),
+    ],
+    [
+        C_ZERO,
+        c64(0., FRAC_1_SQRT_2),
+        c64(FRAC_1_SQRT_2, 0.),
+        C_ZERO,
+    ],
+    [
+        C_ZERO,
+        c64(0., FRAC_1_SQRT_2),
+        c64(-FRAC_1_SQRT_2, 0.),
+        C_ZERO,
+    ],
+    [
+        c64(FRAC_1_SQRT_2, 0.),
+        C_ZERO,
+        C_ZERO,
+        c64(0., -FRAC_1_SQRT_2),
+    ],
+];
+
+static MAGIC_DAGGER: GateArray2Q = [
+    [
+        c64(FRAC_1_SQRT_2, 0.),
+        C_ZERO,
+        C_ZERO,
+        c64(FRAC_1_SQRT_2, 0.),
+    ],
+    [
+        C_ZERO,
+        c64(0., -FRAC_1_SQRT_2),
+        c64(0., -FRAC_1_SQRT_2),
+        C_ZERO,
+    ],
+    [
+        C_ZERO,
+        c64(FRAC_1_SQRT_2, 0.),
+        c64(-FRAC_1_SQRT_2, 0.),
+        C_ZERO,
+    ],
+    [
+        c64(0., -FRAC_1_SQRT_2),
+        C_ZERO,
+        C_ZERO,
+        c64(0., FRAC_1_SQRT_2),
+    ],
+];
+
+/// Computes the local invariants for a two-qubit unitary.
+///
+/// Based on:
+///
+/// Y. Makhlin, Quant. Info. Proc. 1, 243-252 (2002).
+///
+/// Zhang et al., Phys Rev A. 67, 042313 (2003).
+#[pyfunction]
+pub fn two_qubit_local_invariants(unitary: PyReadonlyArray2<Complex64>) -> [f64; 3] {
+    let mat = unitary.as_array();
+    // Transform to bell basis
+    let bell_basis_unitary = aview2(&MAGIC_DAGGER).dot(&mat.dot(&aview2(&MAGIC)));
+    // Get determinate since +- one is allowed.
+    let det_bell_basis = bell_basis_unitary
+        .view()
+        .into_faer_complex()
+        .determinant()
+        .to_num_complex();
+    let m = bell_basis_unitary.t().dot(&bell_basis_unitary);
+    let mut m_tr2 = m.diag().sum();
+    m_tr2 *= m_tr2;
+    // Table II of Ref. 1 or Eq. 28 of Ref. 2.
+    let g1 = m_tr2 / (16. * det_bell_basis);
+    let g2 = (m_tr2 - m.dot(&m).diag().sum()) / (4. * det_bell_basis);
+
+    // Here we split the real and imag pieces of G1 into two so as
+    // to better equate to the Weyl chamber coordinates (c0,c1,c2)
+    // and explore the parameter space.
+    // Also do a FP trick -0.0 + 0.0 = 0.0
+    [g1.re + 0., g1.im + 0., g2.re + 0.]
+}
+
+#[pyfunction]
+pub fn local_equivalence(weyl: PyReadonlyArray1<f64>) -> PyResult<[f64; 3]> {
+    let weyl = weyl.as_array();
+    let weyl_2_cos_squared_product: f64 = weyl.iter().map(|x| (x * 2.).cos().powi(2)).product();
+    let weyl_2_sin_squared_product: f64 = weyl.iter().map(|x| (x * 2.).sin().powi(2)).product();
+    let g0_equiv = weyl_2_cos_squared_product - weyl_2_sin_squared_product;
+    let g1_equiv = weyl.iter().map(|x| (x * 4.).sin()).product::<f64>() / 4.;
+    let g2_equiv = 4. * weyl_2_cos_squared_product
+        - 4. * weyl_2_sin_squared_product
+        - weyl.iter().map(|x| (4. * x).cos()).product::<f64>();
+    Ok([g0_equiv + 0., g1_equiv + 0., g2_equiv + 0.])
+}
+
 #[pymodule]
 pub fn two_qubit_decompose(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(_num_basis_gates))?;
+    m.add_wrapped(wrap_pyfunction!(two_qubit_local_invariants))?;
+    m.add_wrapped(wrap_pyfunction!(local_equivalence))?;
     m.add_class::<TwoQubitGateSequence>()?;
     m.add_class::<TwoQubitWeylDecomposition>()?;
     m.add_class::<Specialization>()?;

qiskit/synthesis/two_qubit/local_invariance.py

@@ -15,17 +15,9 @@
 of two-qubit unitary operators.
 """
 from __future__ import annotations
-from math import sqrt
 import numpy as np
-
-INVARIANT_TOL = 1e-12
-
-# Bell "Magic" basis
-MAGIC = (
-    1.0
-    / sqrt(2)
-    * np.array([[1, 0, 0, 1j], [0, 1j, 1, 0], [0, 1j, -1, 0], [1, 0, 0, -1j]], dtype=complex)
-)
+from qiskit._accelerate.two_qubit_decompose import two_qubit_local_invariants as tqli_rs
+from qiskit._accelerate.two_qubit_decompose import local_equivalence as le_rs
 
 
 def two_qubit_local_invariants(U: np.ndarray) -> np.ndarray:
@@ -44,28 +36,11 @@ def two_qubit_local_invariants(U: np.ndarray) -> np.ndarray:
         Y. Makhlin, Quant. Info. Proc. 1, 243-252 (2002).
         Zhang et al., Phys Rev A. 67, 042313 (2003).
     """
-    U = np.asarray(U)
+    U = np.asarray(U, dtype=complex)
     if U.shape != (4, 4):
         raise ValueError("Unitary must correspond to a two-qubit gate.")
-
-    # Transform to bell basis
-    Um = MAGIC.conj().T.dot(U.dot(MAGIC))
-    # Get determinate since +- one is allowed.
-    det_um = np.linalg.det(Um)
-    M = Um.T.dot(Um)
-    # trace(M)**2
-    m_tr2 = M.trace()
-    m_tr2 *= m_tr2
-
-    # Table II of Ref. 1 or Eq. 28 of Ref. 2.
-    G1 = m_tr2 / (16 * det_um)
-    G2 = (m_tr2 - np.trace(M.dot(M))) / (4 * det_um)
-
-    # Here we split the real and imag pieces of G1 into two so as
-    # to better equate to the Weyl chamber coordinates (c0,c1,c2)
-    # and explore the parameter space.
-    # Also do a FP trick -0.0 + 0.0 = 0.0
-    return np.round([G1.real, G1.imag, G2.real], 12) + 0.0
+    (a, b, c) = tqli_rs(U)
+    return np.array([round(a, 12), round(b, 12), round(c, 12)])
 
 
 def local_equivalence(weyl: np.ndarray) -> np.ndarray:
@@ -83,11 +58,6 @@ def local_equivalence(weyl: np.ndarray) -> np.ndarray:
         but we multiply weyl coordinates by 2 since we are
         working in the reduced chamber.
     """
-    g0_equiv = np.prod(np.cos(2 * weyl) ** 2) - np.prod(np.sin(2 * weyl) ** 2)
-    g1_equiv = np.prod(np.sin(4 * weyl)) / 4
-    g2_equiv = (
-        4 * np.prod(np.cos(2 * weyl) ** 2)
-        - 4 * np.prod(np.sin(2 * weyl) ** 2)
-        - np.prod(np.cos(4 * weyl))
-    )
-    return np.round([g0_equiv, g1_equiv, g2_equiv], 12) + 0.0
+    mat = np.asarray(weyl, dtype=float)
+    (a, b, c) = le_rs(mat)
+    return np.array([round(a, 12), round(b, 12), round(c, 12)])