Merge branch 'master' into master

Author: jk20342

doc/releases/changelog-dev.md

@@ -30,6 +30,10 @@
 
 <h3>Improvements 🛠</h3>
 
+* Added a new decomposition, `_decompose_2_cnots`, for the two-qubit decomposition for `QubitUnitary`.
+  It supports the analytical decomposition a two-qubit unitary known to require exactly 2 CNOTs.
+  [(#8666)](https://github.com/PennyLaneAI/pennylane/issues/8666)
+  
 * Arithmetic dunder methods (`__add__`, `__mul__`, `__rmul__`) have been added to 
   :class:`~.transforms.core.TransformDispatcher`, :class:`~.transforms.core.TransformContainer`, 
   and :class:`~.transforms.core.TransformProgram` to enable intuitive composition of transform 
@@ -614,6 +618,7 @@ Mudit Pandey,
 Shuli Shu,
 Jay Soni,
 nate stemen,
+Theodoros Trochatos,
 David Wierichs,
 Hongsheng Zheng,
 Zinan Zhou

pennylane/ops/op_math/decompositions/unitary_decompositions.py

@@ -15,6 +15,7 @@
 """This module defines decomposition functions for unitary matrices."""
 
 import warnings
+from itertools import product
 
 import numpy as np
 from scipy import sparse
@@ -125,8 +126,10 @@ def two_qubit_decomposition(U, wires):
     where :math:`A, B, C, D` are :math:`SU(2)` operations, and the rotation angles are
     computed based on features of the input unitary :math:`U`.
 
-    For the 2-CNOT case, the decomposition is currently not supported and will
-    instead produce a 3-CNOT circuit like above.
+    For the 2-CNOT case, the decomposition is based on the
+    real-trace criterion of Proposition III.3 in reference (2).
+    Whenever :math:`trace(\gamma(U))` has real coefficients (equivalently :math:`trace(\gamma(U)`) ∈ R),
+    the decomposition uses exactly two CNOT gates.
 
     For a single CNOT, we have a CNOT surrounded by one :math:`SU(2)` per wire on each
     side.  The special case of no CNOTs simply returns a tensor product of two
@@ -207,13 +210,20 @@ def two_qubit_decomposition(U, wires):
             global_phase += _decompose_3_cnots(U, wires, global_phase)
         else:
             num_cnots = _compute_num_cnots(U)
-            # Use the 3-CNOT case for num_cnots=2 as well because we do not have a reliably
-            # correct implementation of the 2-CNOT case right now.
+
+            elifs = [(num_cnots == 1, _decompose_1_cnot)]
+
+            # The 2-CNOT decomposition relies on sorting eigenvalues, which is not supported
+            # with abstract tracers when capture is enabled. In that case, we fall back
+            # to the 3-CNOT decomposition.
+            if not capture.enabled():
+                elifs.append((num_cnots == 2, _decompose_2_cnots))
+
             global_phase += ops.cond(
                 num_cnots == 0,
                 _decompose_0_cnots,
                 _decompose_3_cnots,
-                elifs=[(num_cnots == 1, _decompose_1_cnot)],
+                elifs=elifs,
             )(U, wires, global_phase)
 
         if _is_jax_jit(U) or not math.allclose(global_phase, 0):
@@ -377,13 +387,20 @@ def two_qubit_decomp_rule(U, wires, **__):
 
     U, initial_phase = math.convert_to_su4(U, return_global_phase=True)
     num_cnots = _compute_num_cnots(U)
-    # Use the 3-CNOT case for num_cnots=2 as well because we do not have a reliably
-    # correct implementation of the 2-CNOT case right now.
+
+    elifs = [(num_cnots == 1, _decompose_1_cnot)]
+
+    # The 2-CNOT decomposition relies on sorting eigenvalues, which is not supported
+    # with abstract tracers when capture is enabled. In that case, we fall back
+    # to the 3-CNOT decomposition.
+    if not capture.enabled():
+        elifs.append((num_cnots == 2, _decompose_2_cnots))
+
     additional_phase = ops.cond(
         num_cnots == 0,
         _decompose_0_cnots,
         _decompose_3_cnots,
-        elifs=[(num_cnots == 1, _decompose_1_cnot)],
+        elifs=elifs,
     )(U, wires, initial_phase)
     total_phase = initial_phase + additional_phase
     ops.cond(math.logical_not(math.allclose(total_phase, 0)), ops.GlobalPhase)(-total_phase)
@@ -474,15 +491,15 @@ def multi_qubit_decomp_rule(U, wires, **__):
 
 
 def _compute_num_cnots(U):
-    r"""Compute the number of CNOTs required to implement a U in SU(4).
+    r"""
+    Compute the number of CNOTs required to implement a U in SU(4).
     This is based on the trace of
 
     .. math::
 
         \gamma(U) = (E^\dag U E) (E^\dag U E)^T,
 
     and follows the arguments of this paper: https://arxiv.org/abs/quant-ph/0308045.
-
     """
 
     U = math.dot(E_dag, math.dot(U, E))
@@ -604,6 +621,190 @@ def _decompose_1_cnot(U, wires, initial_phase):
     return math.cast_like(-np.pi / 4, initial_phase)
 
 
+def _get_basis_and_eigenvalues(M):
+    r"""
+    Helper to diagonalize M, extract diagonal eigenvalues, and sort canonically.
+    Returns eigenvalues and basis O such that :math:`D = O^T M O` is diagonal with
+    sorted eigenvalues.
+    """
+    # pylint: disable=protected-access
+    # Use split_eigh to get a basis (ignoring its mixed eigenvalues)
+    _, O = _real_imag_split_eigh(M, 1.0)
+
+    # Compute true eigenvalues: D = O.T @ M @ O
+    d_mat = math.dot(math.transpose(O), math.dot(M, O))
+    eigvals = math.diag(d_mat)
+
+    # Canonical Sort: Real part descending, then Imaginary part descending
+    r = np.round(math.real(eigvals), 6)
+    i = np.round(math.imag(eigvals), 6)
+
+    # Sort by imag then real to ensure deterministic aligning of U and V
+    sort_indices = np.lexsort((i, r))
+
+    # Reorder
+    eigvals_sorted = eigvals[sort_indices]
+    O_sorted = O[:, sort_indices]
+
+    # Enforce determinant 1 (SO(4))
+    det = math.linalg.det(O_sorted)
+    if math.real(det) < 0:
+        O_sorted = math.set_index(O_sorted, (slice(None), 3), -O_sorted[:, 3])
+
+    return eigvals_sorted, O_sorted
+
+
+def _find_so4_decomposition(U, u_mag, O_u, candidates):
+    r"""
+    Performs the exhaustive search for alpha, beta, and signs
+    to ensure Real :math:`SO(4)` correction gates.
+    Returns the best found parameters along with the basis O_v and
+    v_mag for the kernel V.
+    """
+    CNOT10_np = np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]], dtype=complex)
+    CNOT10 = math.cast_like(CNOT10_np, U)
+    best_result = None
+    min_error = np.inf
+
+    # Search Parameter Permutations
+    for alpha, beta in candidates:
+        # Construct Kernel V
+        rza = ops.RZ.compute_matrix(alpha)
+        rxb = ops.RX.compute_matrix(beta)
+        kernel_inner = math.kron(rza, rxb)
+        V = _multidot(CNOT10, kernel_inner, CNOT10)
+
+        # Get basis for V
+        v_mag = _multidot(math.cast_like(E_dag, U), V, math.cast_like(E, U))
+        gamma_v = math.dot(v_mag, math.T(v_mag))
+        _, O_v = _get_basis_and_eigenvalues(gamma_v)
+
+        # P = v^dag O_v, Q = O_u^T u
+        P = math.dot(math.conj(math.T(v_mag)), O_v)
+        Q = math.dot(math.T(O_u), u_mag)
+
+        # Sign Search to fix Gauge Freedom
+        for signs in product([1.0, -1.0], repeat=4):
+            # Enforce determinant +1 to stay in SO(4)
+            if np.prod(signs) < 0:
+                continue
+
+            S_diag = math.cast_like(signs, P)
+            # Broadcasting P * S is faster than matrix mult
+            # Metric: sum of absolute imaginary parts (should be 0 for valid decomposition)
+            R_trial = math.dot(P * S_diag, Q)
+            error = math.sum(math.abs(math.imag(R_trial)))
+
+            if error < min_error:
+                min_error = error
+                best_result = (alpha, beta, signs, O_v, v_mag)
+
+    # FALLBACK CHECK:
+    # If the best error we found is still "large" (e.g., > 1e-5),
+    # then we failed to find a valid Real-valued decomposition.
+    # We return None to signal that 2-CNOT decomposition is likely impossible/unsafe.
+    if min_error > 1e-5:
+        return None
+
+    return best_result
+
+
+def _decompose_2_cnots(U, wires, initial_phase):
+    r"""Decompose a two-qubit unitary known to require exactly 2 CNOTs.
+
+    The resulting circuit has the following canonical form:
+        0: ──A──╭X──RZ(a)──╭X──C──┤
+        1: ──B──╰●──RX(b)──╰●──D──┤
+
+    where A, B, C, D are single-qubit unitaries (:math:`SU(2)` gates) and a, b
+    are rotation angles determined by the entanglement properties of the unitary.
+
+    This implementation is based on the work by Shende, Bullock, and Markov and
+    this is done following the methods in https://arxiv.org/abs/quant-ph/0308045.
+
+    The decomposition relies on the Magic Basis, where local unitaries correspond to
+    orthogonal matrices (:math:`SO(4)`), and the entangling power of an operator is captured by
+    the spectrum (eigenvalues) of the symmetric invariant matrix: :math:`\gamma(U) = (E^\dagger U E) (E^\dagger U E)^T`.
+
+    The algorithm proceeds in four main steps:
+
+    1.  The invariant matrix :math:`\gamma(U)` is computed. Its eigenvalues
+        are extracted. These eigenvalues come in conjugate pairs,
+        and their phases directly determine the rotation parameters a, and b
+        needed for the circuit's core.
+
+    2.  A reference operator V is built using the
+        calculated parameters a, b. This operator represents the ideal
+        2-CNOT circuit core: :math:`CNOT \cdot (RZ(\alpha) \otimes RX(\beta)) \cdot CNOT`.
+        By construction, V is isospectral to U in the Magic Basis.
+
+    3.  To transform the input U into V using only
+        local gates, we align their eigenbases. This involves diagonalizing :math:`\gamma(U)` and
+        :math:`\gamma(V)` and sorting their eigenvectors canonically to match corresponding subspaces.
+
+    4.  Since eigenvectors are defined only up to a sign (parity), there are
+        :math:`2^4 = 16` possible alignments. The algorithm exhaustively searches these combinations
+        to find the specific parity that results in a valid, real-valued local transformation
+        (a matrix in :math:`SO(4)`). This transformation determines the local gates A, B, C, D.
+
+    The final circuit is then constructed by applying these local gates around the V.
+    """
+    # pylint: disable=too-many-locals
+    # 1. Compute gamma(U)
+    u_mag = _multidot(math.cast_like(E_dag, U), U, math.cast_like(E, U))
+    gamma_u = math.dot(u_mag, math.T(u_mag))
+
+    # 2. Extract interaction parameters
+    eig_u, O_u = _get_basis_and_eigenvalues(gamma_u)
+    # Extract phases and sort to group conjugate pairs
+    abs_angles = np.sort(np.abs(math.angle(eig_u)))
+    # Pick distinct representatives
+    theta1 = abs_angles[3]
+    theta2 = abs_angles[1]
+    # Map to circuit parameters
+    a_calc = (theta1 + theta2) / 2
+    b_calc = (theta1 - theta2) / 2
+    candidates = [(a_calc, b_calc), (b_calc, a_calc)]
+
+    # 3. Perform Search (Delegated to helper)
+    result = _find_so4_decomposition(U, u_mag, O_u, candidates)
+
+    # SAFETY FALLBACK to _decompose_3_cnots:
+    # If the 2-CNOT search failed (result is None), it means U is likely
+    # not a 2-CNOT gate (despite trace invariants) or numerical noise is too high.
+    # We fall back to the generic 3-CNOT decomposition to guarantee correctness.
+    if result is None:
+        return _decompose_3_cnots(U, wires, initial_phase)
+
+    alpha_f, beta_f, signs_f, O_v, v_mag = result
+
+    # 4. Compute Local Gates L (Left) and R (Right) in SO(4)
+    S_mat = math.diag(signs_f)
+    L = _multidot(O_u, S_mat, math.T(O_v))
+    R = _multidot(math.conj(math.T(v_mag)), math.T(L), u_mag)
+
+    # 5. Convert to local gates
+    AB = _multidot(math.cast_like(E, U), L, math.cast_like(E_dag, U))
+    CD = _multidot(math.cast_like(E, U), R, math.cast_like(E_dag, U))
+
+    A, B = math.decomposition.su2su2_to_tensor_products(AB)
+    C, D = math.decomposition.su2su2_to_tensor_products(CD)
+
+    # 6. Queue Circuit
+    ops.QubitUnitary(C, wires=wires[0])
+    ops.QubitUnitary(D, wires=wires[1])
+
+    ops.CNOT(wires=[wires[1], wires[0]])
+    ops.RZ(alpha_f, wires=wires[0])
+    ops.RX(beta_f, wires=wires[1])
+    ops.CNOT(wires=[wires[1], wires[0]])
+
+    ops.QubitUnitary(A, wires=wires[0])
+    ops.QubitUnitary(B, wires=wires[1])
+
+    return math.cast_like(0.0, initial_phase)
+
+
 def _multidot(*matrices):
     mat = matrices[0]
     for m in matrices[1:]:

tests/ops/op_math/test_decompositions.py

@@ -982,6 +982,28 @@ def test_rot_decomposition_jax(self, U, expected_gates, expected_params):
 class TestTwoQubitUnitaryDecomposition:
     """Test that two-qubit unitary operations are correctly decomposed."""
 
+    @pytest.mark.unit
+    @pytest.mark.parametrize("U", samples_2_cnots)
+    def test_compute_num_cnots_identifies_2_cnots(self, U):
+        """Test that the new Shende–Bullock–Markov criterion correctly
+        classifies 2-CNOT unitaries."""
+        U = qml.math.convert_to_su4(np.array(U))
+        assert _compute_num_cnots(U) == 2
+
+    @pytest.mark.unit
+    @pytest.mark.parametrize("U", samples_2_cnots)
+    def test_two_qubit_decomposition_2_cnots_gate_count(self, U):
+        """Test that the dispatcher selects the new 2-CNOT decomposition
+        and that the resulting circuit actually contains exactly 2 CNOTs."""
+        U = qml.math.convert_to_su4(np.array(U))
+
+        ops = two_qubit_decomposition(U, wires=[0, 1])
+
+        # Extract only CNOTs
+        cnot_ops = [op for op in ops if isinstance(op, qml.CNOT)]
+
+        assert len(cnot_ops) == 2, "The 2-CNOT decomposition must emit exactly 2 CNOT gates."
+
     @pytest.mark.parametrize("U_pair", samples_su2_su2)
     def test_su2su2_to_tensor_products(self, U_pair):
         """Test SU(2) x SU(2) can be correctly factored into tensor products."""
@@ -1013,14 +1035,13 @@ def test_two_qubit_decomposition_3_cnots(self, U, wires):
     @pytest.mark.parametrize("U", samples_2_cnots)
     def test_two_qubit_decomposition_2_cnots(self, U, wires):
         """Test that a two-qubit matrix using 2 CNOTs isolation is correctly decomposed."""
-        # NOTE: Currently, we defer to the 3-CNOTs function for the 2-CNOTs case.
 
         U = qml.math.convert_to_su4(np.array(U))
 
         assert _compute_num_cnots(U) == 2
 
         obtained_decomposition = two_qubit_decomposition(U, wires=wires)
-        assert len(obtained_decomposition) == 11  # 8 # 8 would be the count with 2-CNOT circuit
+        assert len(obtained_decomposition) == 8
 
         tape = qml.tape.QuantumScript(obtained_decomposition)
         obtained_matrix = qml.matrix(tape, wire_order=wires)