diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 568206925c2..ac2cc1d2e50 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/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) -> [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) -> 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::() / 4.; + let g2_equiv = 4. * weyl_2_cos_squared_product + - 4. * weyl_2_sin_squared_product + - weyl.iter().map(|x| (4. * x).cos()).product::(); + Ok([g0_equiv + 0., g1_equiv + 0., g2_equiv + 0.]) +} + #[pymodule] pub fn two_qubit_decompose(m: &Bound) -> 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::()?; m.add_class::()?; m.add_class::()?; diff --git a/qiskit/synthesis/two_qubit/local_invariance.py b/qiskit/synthesis/two_qubit/local_invariance.py index 346febdddf9..7b20c2e5655 100644 --- a/qiskit/synthesis/two_qubit/local_invariance.py +++ b/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)])