Oxidize two qubit local invariance functions (#12739)

* Oxidize two qubit local invariance functions

This commit migrates the two functions in the private module
`qiskit.synthesis.two_qubit.local_invariance` to primarily be
implemented in rust. Since the two_qubit_local_invariants() function is
being used in #12727 premptively porting these functions to rust will
both potentially speed up the new transpiler pass in that PR and also
facilitate a future porting of that pass to rust. The two python space
functions continue to exist as a small wrapper to do input type
checking/conversion and rounding of the result (since the python API for
rounding is simpler). There is no release note included for these
functions as they are internal utilities in Qiskit and not exposed as a
public interface.

* Add docstring to rust function with citation

* Store adjoint magic array statically

* Use arraview1 instead of slice for local_equivalence()

* Fix rustfmt

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)])