diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index d9f451a43682..8813eb486503 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/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` or `MatMut`. -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) -> Mat { +fn transform_from_magic_basis(unitary: MatRef, reverse: bool) -> Mat { let _b_nonnormalized = mat::from_row_major_slice::(&B_NON_NORMALIZED, 4, 4); let _b_nonnormalized_dagger = mat::from_row_major_slice::(&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, + 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 are both structs @@ -100,7 +143,7 @@ impl Arg for c64 { fn __weyl_coordinates(unitary: MatRef) -> [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) -> [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) -> 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(()) } diff --git a/qiskit/__init__.py b/qiskit/__init__.py index 782d8aaabfbd..4d91f6d6e86a 100644 --- a/qiskit/__init__.py +++ b/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 diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index 9985d7cf50c6..1d5ccb2cbeef 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/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) diff --git a/qiskit/synthesis/two_qubit/weyl.py b/qiskit/synthesis/two_qubit/weyl.py index b533b69a3696..a209ba16e9f2 100644 --- a/qiskit/synthesis/two_qubit/weyl.py +++ b/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 diff --git a/test/python/synthesis/test_weyl.py b/test/python/synthesis/test_weyl.py index f0773af205ef..514145f88834 100644 --- a/test/python/synthesis/test_weyl.py +++ b/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])