From 653b5b66f5d856290acf26d40a9dac29ccfa21d5 Mon Sep 17 00:00:00 2001 From: Shelly Garion <46566946+ShellyGarion@users.noreply.github.com> Date: Fri, 13 Sep 2024 17:35:03 +0300 Subject: [PATCH] Clean two_qubit_decompose code (#13093) * add pyfunction for decompose_two_qubit_product_gate * replace python code by rust for decompose_two_qubit_product_gate * remove unused code from two_qubit_decompose.py * updates following review comments * port trace_to_fid to rust * port Ud to rust * fix lint errors * fix lint errors * format * remove PyResult following review --- crates/accelerate/src/two_qubit_decompose.rs | 73 ++++++++++++++++++- .../two_qubit/two_qubit_decompose.py | 67 +---------------- test/python/synthesis/test_synthesis.py | 2 +- test/randomized/test_synthesis.py | 2 +- 4 files changed, 78 insertions(+), 66 deletions(-) diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs index 92ad4724682f..76b92acd3faf 100644 --- a/crates/accelerate/src/two_qubit_decompose.rs +++ b/crates/accelerate/src/two_qubit_decompose.rs @@ -32,7 +32,7 @@ use ndarray::linalg::kron; use ndarray::prelude::*; use ndarray::Zip; use numpy::{IntoPyArray, ToPyArray}; -use numpy::{PyReadonlyArray1, PyReadonlyArray2}; +use numpy::{PyArray2, PyArrayLike2, PyReadonlyArray1, PyReadonlyArray2}; use pyo3::exceptions::PyValueError; use pyo3::prelude::*; @@ -154,6 +154,15 @@ impl TraceToFidelity for c64 { } } +#[pyfunction] +#[pyo3(name = "trace_to_fid")] +/// Average gate fidelity is :math:`Fbar = (d + |Tr (Utarget \\cdot U^dag)|^2) / d(d+1)` +/// M. Horodecki, P. Horodecki and R. Horodecki, PRA 60, 1888 (1999) +fn py_trace_to_fid(trace: Complex64) -> PyResult { + let fid = trace.trace_to_fid(); + Ok(fid) +} + fn decompose_two_qubit_product_gate( special_unitary: ArrayView2, ) -> PyResult<(Array2, Array2, f64)> { @@ -182,9 +191,31 @@ fn decompose_two_qubit_product_gate( } l.mapv_inplace(|x| x / det_l.sqrt()); let phase = det_l.arg() / 2.; + Ok((l, r, phase)) } +#[pyfunction] +#[pyo3(name = "decompose_two_qubit_product_gate")] +/// Decompose :math:`U = U_l \otimes U_r` where :math:`U \in SU(4)`, +/// and :math:`U_l,~U_r \in SU(2)`. +/// Args: +/// special_unitary_matrix: special unitary matrix to decompose +/// Raises: +/// QiskitError: if decomposition isn't possible. +fn py_decompose_two_qubit_product_gate( + py: Python, + special_unitary: PyArrayLike2, +) -> PyResult<(PyObject, PyObject, f64)> { + let view = special_unitary.as_array(); + let (l, r, phase) = decompose_two_qubit_product_gate(view)?; + Ok(( + l.into_pyarray_bound(py).unbind().into(), + r.into_pyarray_bound(py).unbind().into(), + phase, + )) +} + fn __weyl_coordinates(unitary: MatRef) -> [f64; 3] { let uscaled = scale(C1 / unitary.determinant().powf(0.25)) * unitary; let uup = transform_from_magic_basis(uscaled); @@ -301,6 +332,43 @@ fn rz_matrix(theta: f64) -> Array2 { array![[(-ilam2).exp(), C_ZERO], [C_ZERO, ilam2.exp()]] } +/// Generates the array :math:`e^{(i a XX + i b YY + i c ZZ)}` +fn ud(a: f64, b: f64, c: f64) -> Array2 { + array![ + [ + (IM * c).exp() * (a - b).cos(), + C_ZERO, + C_ZERO, + IM * (IM * c).exp() * (a - b).sin() + ], + [ + C_ZERO, + (M_IM * c).exp() * (a + b).cos(), + IM * (M_IM * c).exp() * (a + b).sin(), + C_ZERO + ], + [ + C_ZERO, + IM * (M_IM * c).exp() * (a + b).sin(), + (M_IM * c).exp() * (a + b).cos(), + C_ZERO + ], + [ + IM * (IM * c).exp() * (a - b).sin(), + C_ZERO, + C_ZERO, + (IM * c).exp() * (a - b).cos() + ] + ] +} + +#[pyfunction] +#[pyo3(name = "Ud")] +fn py_ud(py: Python, a: f64, b: f64, c: f64) -> Py> { + let ud_mat = ud(a, b, c); + ud_mat.into_pyarray_bound(py).unbind() +} + fn compute_unitary(sequence: &TwoQubitSequenceVec, global_phase: f64) -> Array2 { let identity = aview2(&ONE_QUBIT_IDENTITY); let phase = c64(0., global_phase).exp(); @@ -2278,9 +2346,12 @@ pub fn local_equivalence(weyl: PyReadonlyArray1) -> PyResult<[f64; 3]> { pub fn two_qubit_decompose(m: &Bound) -> PyResult<()> { m.add_wrapped(wrap_pyfunction!(_num_basis_gates))?; + m.add_wrapped(wrap_pyfunction!(py_decompose_two_qubit_product_gate))?; m.add_wrapped(wrap_pyfunction!(two_qubit_decompose_up_to_diagonal))?; m.add_wrapped(wrap_pyfunction!(two_qubit_local_invariants))?; m.add_wrapped(wrap_pyfunction!(local_equivalence))?; + m.add_wrapped(wrap_pyfunction!(py_trace_to_fid))?; + m.add_wrapped(wrap_pyfunction!(py_ud))?; m.add_class::()?; m.add_class::()?; m.add_class::()?; diff --git a/qiskit/synthesis/two_qubit/two_qubit_decompose.py b/qiskit/synthesis/two_qubit/two_qubit_decompose.py index 0e3427773387..d4c7702da35a 100644 --- a/qiskit/synthesis/two_qubit/two_qubit_decompose.py +++ b/qiskit/synthesis/two_qubit/two_qubit_decompose.py @@ -24,8 +24,6 @@ arXiv:1811.12926 [quant-ph] (2018). """ from __future__ import annotations -import cmath -import math import io import base64 import warnings @@ -91,41 +89,18 @@ def decompose_two_qubit_product_gate(special_unitary_matrix: np.ndarray): QiskitError: if decomposition isn't possible. """ special_unitary_matrix = np.asarray(special_unitary_matrix, dtype=complex) - # extract the right component - R = special_unitary_matrix[:2, :2].copy() - detR = R[0, 0] * R[1, 1] - R[0, 1] * R[1, 0] - if abs(detR) < 0.1: - R = special_unitary_matrix[2:, :2].copy() - detR = R[0, 0] * R[1, 1] - R[0, 1] * R[1, 0] - if abs(detR) < 0.1: - raise QiskitError("decompose_two_qubit_product_gate: unable to decompose: detR < 0.1") - R /= np.sqrt(detR) - - # extract the left component - temp = np.kron(np.eye(2), R.T.conj()) - temp = special_unitary_matrix.dot(temp) - L = temp[::2, ::2] - detL = L[0, 0] * L[1, 1] - L[0, 1] * L[1, 0] - if abs(detL) < 0.9: - raise QiskitError("decompose_two_qubit_product_gate: unable to decompose: detL < 0.9") - L /= np.sqrt(detL) - phase = cmath.phase(detL) / 2 + (L, R, phase) = two_qubit_decompose.decompose_two_qubit_product_gate(special_unitary_matrix) temp = np.kron(L, R) deviation = abs(abs(temp.conj().T.dot(special_unitary_matrix).trace()) - 4) + if deviation > 1.0e-13: raise QiskitError( "decompose_two_qubit_product_gate: decomposition failed: " f"deviation too large: {deviation}" ) - return L, R, phase - - -_ipx = np.array([[0, 1j], [1j, 0]], dtype=complex) -_ipy = np.array([[0, 1], [-1, 0]], dtype=complex) -_ipz = np.array([[1j, 0], [0, -1j]], dtype=complex) -_id = np.array([[1, 0], [0, 1]], dtype=complex) + return (L, R, phase) class TwoQubitWeylDecomposition: @@ -239,7 +214,7 @@ def actual_fidelity(self, **kwargs) -> float: """Calculates the actual fidelity of the decomposed circuit to the input unitary.""" circ = self.circuit(**kwargs) trace = np.trace(Operator(circ).data.T.conj() @ self.unitary_matrix) - return trace_to_fid(trace) + return two_qubit_decompose.trace_to_fid(trace) def __repr__(self): """Represent with enough precision to allow copy-paste debugging of all corner cases""" @@ -460,40 +435,6 @@ def _weyl_gate(self, circ: QuantumCircuit, atol=1.0e-13): return circ -def Ud(a, b, c): - r"""Generates the array :math:`e^{(i a XX + i b YY + i c ZZ)}`""" - return np.array( - [ - [cmath.exp(1j * c) * math.cos(a - b), 0, 0, 1j * cmath.exp(1j * c) * math.sin(a - b)], - [0, cmath.exp(-1j * c) * math.cos(a + b), 1j * cmath.exp(-1j * c) * math.sin(a + b), 0], - [0, 1j * cmath.exp(-1j * c) * math.sin(a + b), cmath.exp(-1j * c) * math.cos(a + b), 0], - [1j * cmath.exp(1j * c) * math.sin(a - b), 0, 0, cmath.exp(1j * c) * math.cos(a - b)], - ], - dtype=complex, - ) - - -def trace_to_fid(trace): - r"""Average gate fidelity is - - .. math:: - - \bar{F} = \frac{d + |\mathrm{Tr} (U_\text{target} \cdot U^{\dag})|^2}{d(d+1)} - - M. Horodecki, P. Horodecki and R. Horodecki, PRA 60, 1888 (1999)""" - return (4 + abs(trace) ** 2) / 20 - - -def rz_array(theta): - """Return numpy array for Rz(theta). - - Rz(theta) = diag(exp(-i*theta/2),exp(i*theta/2)) - """ - return np.array( - [[cmath.exp(-1j * theta / 2.0), 0], [0, cmath.exp(1j * theta / 2.0)]], dtype=complex - ) - - class TwoQubitBasisDecomposer: """A class for decomposing 2-qubit unitaries into minimal number of uses of a 2-qubit basis gate. diff --git a/test/python/synthesis/test_synthesis.py b/test/python/synthesis/test_synthesis.py index b4fac3f427f0..b08240197211 100644 --- a/test/python/synthesis/test_synthesis.py +++ b/test/python/synthesis/test_synthesis.py @@ -59,11 +59,11 @@ two_qubit_cnot_decompose, TwoQubitBasisDecomposer, TwoQubitControlledUDecomposer, - Ud, decompose_two_qubit_product_gate, ) from qiskit._accelerate.two_qubit_decompose import two_qubit_decompose_up_to_diagonal from qiskit._accelerate.two_qubit_decompose import Specialization +from qiskit._accelerate.two_qubit_decompose import Ud from qiskit.synthesis.unitary import qsd from test import combine # pylint: disable=wrong-import-order from test import QiskitTestCase # pylint: disable=wrong-import-order diff --git a/test/randomized/test_synthesis.py b/test/randomized/test_synthesis.py index 9f0619a0c296..5e500b78dafc 100644 --- a/test/randomized/test_synthesis.py +++ b/test/randomized/test_synthesis.py @@ -23,8 +23,8 @@ from qiskit.synthesis.two_qubit.two_qubit_decompose import ( two_qubit_cnot_decompose, TwoQubitBasisDecomposer, - Ud, ) +from qiskit._accelerate.two_qubit_decompose import Ud class TestSynthesis(CheckDecompositions):