Repository: PennyLaneAI/pennylane (python)

Issue:
[BUG] Decomposition of QubitUnitary on two wires is incorrect
### Expected behavior

The `two_qubit_decomposition` is correrct.

### Actual behavior

It is, for the most part, but here's a random matrix that it fails on:

```python
import pennylane as qml
import numpy as np
from pennylane.ops.op_math.decompositions.unitary_decompositions import two_qubit_decomposition

U = np.array(
    [[-0.40836553 - 0.28850926j, 0.23864466 - 0.4393731j,
      0.02190401 + 0.06284322j, -0.23169726 - 0.66474606j],
     [-0.48865844 - 0.10589111j, 0.49770294 + 0.04787262j,
      -0.02177958 + 0.45019526j, -0.02632749 + 0.54420285j],
     [0.49276489 - 0.08475118j, 0.47943096 - 0.14193644j,
      0.51448332 - 0.17932317j, -0.42560958 + 0.14834623j],
     [0.42041002 - 0.27065738j, -0.31807806 - 0.38578018j,
      -0.03401683 + 0.70314565j, -0.00321586 + 0.06647341j]]
)

decomp = two_qubit_decomposition(U, wires=[0, 1])
matrix = qml.matrix(qml.tape.QuantumScript(decomp), wire_order=[0, 1])
```
```pycon
>>> qml.math.allclose(matrix, U, atol=1e-7)
False
```


### Additional information

This is an old bug coming back to haunt us. See https://github.com/PennyLaneAI/pennylane/issues/5308 and https://github.com/PennyLaneAI/pennylane/issues/7339

### Source code

```shell

```

### Tracebacks

```shell

```

### System information

```shell
Name: Pennylane
Version: 0.42.0.dev31
Summary: PennyLane is a cross-platform Python library for quantum computing, quantum machine learning, and quantum chemistry. Train a quantum computer the same way as a neural network.
Home-page: 
Author: 
Author-email: 
License-Expression: Apache-2.0
Location: /Users/astral.cai/Workspace/pennylane
Requires: appdirs, autograd, autoray, cachetools, diastatic-malt, networkx, numpy, packaging, pennylane-lightning, requests, rustworkx, scipy, tomlkit, typing_extensions
Required-by: PennyLane_Lightning

Platform info:           macOS-15.5-arm64-arm-64bit
Python version:          3.10.14
Numpy version:           1.26.4
Scipy version:           1.15.3
Installed devices:
- default.clifford (Pennylane-0.42.0.dev31)
- default.gaussian (Pennylane-0.42.0.dev31)
- default.mixed (Pennylane-0.42.0.dev31)
- default.qubit (Pennylane-0.42.0.dev31)
- default.qutrit (Pennylane-0.42.0.dev31)
- default.qutrit.mixed (Pennylane-0.42.0.dev31)
- default.tensor (Pennylane-0.42.0.dev31)
- null.qubit (Pennylane-0.42.0.dev31)
- reference.qubit (Pennylane-0.42.0.dev31)
- lightning.qubit (PennyLane_Lightning-0.41.1)
```

### Existing GitHub issues

- [x] I have searched existing GitHub issues to make sure the issue does not already exist.

Pull request:
New decomposition of `QubitUnitary` based on Cartan decomposition
**Context:**
See bug report. This decomposition has been a bit fragile through time.

**Description of the Change:**
New implementation for 3CNOT case.

**Benefits:**
Correct code.

**Possible Drawbacks:**
N/A

**Related GitHub Issues:**
Fixes #7467 and #7534

[sc-91723]
[sc-92318]

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 0ccbbe2e6..85a249e20 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -437,6 +437,11 @@
 
 <h3>Breaking changes 💔</h3>
 
+* A new decomposition for two-qubit unitaries was implemented in `two_qubit_decomposition`.
+  It ensures the correctness of the decomposition in some edge cases but uses 3 CNOT gates
+  even if 2 CNOTs would suffice theoretically.
+  [(#7474)](https://github.com/PennyLaneAI/pennylane/pull/7474)
+
 * The `return_type` property of `MeasurementProcess` has been removed. Please use `isinstance` for type checking instead.
   [(#7322)](https://github.com/PennyLaneAI/pennylane/pull/7322)
 
@@ -592,6 +597,12 @@ Here's a list of deprecations made this release. For a more detailed breakdown o
 
 <h3>Bug fixes 🐛</h3>
 
+* A new decomposition was implemented for two-qubit `QubitUnitary` operators in `two_qubit_decomposition`
+  based on a type-AI Cartan decomposition. It fixes previously faulty edge cases for unitaries
+  that require 2 or 3 CNOT gates. Now, 3 CNOTs are used for both cases, using one more
+  CNOT than theoretically required in the former case.
+  [(#7474)](https://github.com/PennyLaneAI/pennylane/pull/7474)
+
 * The documentation of `qml.pulse.drive` has been updated and corrected.
   [(#7459)](https://github.com/PennyLaneAI/pennylane/pull/7459)
 
diff --git a/pennylane/math/decomposition.py b/pennylane/math/decomposition.py
index 5fc165370..f9a71a88b 100644
--- a/pennylane/math/decomposition.py
+++ b/pennylane/math/decomposition.py
@@ -159,7 +159,7 @@ def zxz_rotation_angles(U, return_global_phase=False):
 
     abs_a = math.clip(math.abs(U[..., 0, 0]), 0, 1)
     abs_b = math.clip(math.abs(U[..., 0, 1]), 0, 1)
-    theta = math.where(abs_a > abs_b, 2 * math.arccos(abs_a), 2 * math.arcsin(abs_b))
+    theta = math.where(abs_a < abs_b, 2 * math.arccos(abs_a), 2 * math.arcsin(abs_b))
 
     half_phi_plus_lam = math.angle(U[..., 1, 1] + EPS)
     half_phi_minus_lam = math.angle(1j * U[..., 1, 0] + EPS)
diff --git a/pennylane/math/multi_dispatch.py b/pennylane/math/multi_dispatch.py
index 5955eb469..3b58d2a09 100644
--- a/pennylane/math/multi_dispatch.py
+++ b/pennylane/math/multi_dispatch.py
@@ -1017,7 +1017,7 @@ def detach(tensor, like=None):
 @multi_dispatch(tensor_list=[1])
 def set_index(array, idx, val, like=None):
     """Set the value at a specified index in an array.
-    Calls ``array[idx]=val`` and returns the updated array unless JAX.
+    Calls ``array[idx]=val`` and returns the updated array unless JAX or Tensorflow.
 
     Args:
         array (tensor_like): array to be modified
@@ -1028,8 +1028,6 @@ def set_index(array, idx, val, like=None):
         a new copy of the array with the specified index updated to ``val``.
 
     Whether the original array is modified is interface-dependent.
-
-    .. note:: TensorFlow EagerTensor does not support item assignment
     """
     if like == "jax":
         from jax import numpy as jnp
@@ -1038,5 +1036,10 @@ def set_index(array, idx, val, like=None):
         jax_array = jnp.array(array)
         return jax_array.at[idx].set(val)
 
+    if like == "tensorflow":
+        import tensorflow as tf
+
+        return tf.concat([array[:idx], val[None], array[idx + 1 :]], 0)
+
     array[idx] = val
     return array
diff --git a/pennylane/ops/op_math/decompositions/unitary_decompositions.py b/pennylane/ops/op_math/decompositions/unitary_decompositions.py
index e8c7662d4..5a898968f 100644
--- a/pennylane/ops/op_math/decompositions/unitary_decompositions.py
+++ b/pennylane/ops/op_math/decompositions/unitary_decompositions.py
@@ -123,14 +123,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
+    For the 2-CNOT case, the decomposition is currently not supported and will
+    instead produce a 3-CNOT circuit like above.
 
-    .. figure:: ../../_static/two_qubit_decomposition_2_cnots.svg
-        :align: center
-        :width: 50%
-        :target: javascript:void(0);
-
-    For 1 CNOT, we have a CNOT surrounded by one :math:`SU(2)` per wire on each
+    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
     :math:`SU(2)` operations.
 
@@ -217,14 +213,13 @@ 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.
             global_phase += ops.cond(
                 num_cnots == 0,
                 _decompose_0_cnots,
                 _decompose_3_cnots,
-                elifs=[
-                    (num_cnots == 1, _decompose_1_cnot),
-                    (num_cnots == 2, _decompose_2_cnots),
-                ],
+                elifs=[(num_cnots == 1, _decompose_1_cnot)],
             )(U, wires, global_phase)
 
         if _is_jax_jit(U) or not math.allclose(global_phase, 0):
@@ -378,9 +373,6 @@ def _two_qubit_resource(**_):
         ops.CNOT: 3,
         ops.RZ: 1,
         ops.RY: 2,
-        # The three-CNOT case does not involve an RX, but an RX must be accounted
-        # for in case the two-CNOT case is chosen at runtime.
-        ops.RX: 1,
         ops.GlobalPhase: 1,
     }
 
@@ -392,14 +384,13 @@ 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.
     additional_phase = ops.cond(
         num_cnots == 0,
         _decompose_0_cnots,
         _decompose_3_cnots,
-        elifs=[
-            (num_cnots == 1, _decompose_1_cnot),
-            (num_cnots == 2, _decompose_2_cnots),
-        ],
+        elifs=[(num_cnots == 1, _decompose_1_cnot)],
     )(U, wires, initial_phase)
     total_phase = initial_phase + additional_phase
     ops.cond(math.logical_not(math.allclose(total_phase, 0)), _global_phase)(total_phase)
@@ -484,9 +475,9 @@ E = np.array([[1, 1j, 0, 0], [0, 0, 1j, 1], [0, 0, 1j, -1], [1, -1j, 0, 0]]) / n
 E_dag = E.conj().T
 
 # Helpful to have static copies of these since they are needed in a few places.
-CNOT01 = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
-CNOT10 = np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]])
 SWAP = np.array([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]])
+S_0 = np.kron(np.diag([1j, 1]), np.eye(2))
+S_0_dag = S_0.conj().T
 
 
 def _compute_num_cnots(U):
@@ -526,82 +517,6 @@ def _compute_num_cnots(U):
     )()
 
 
-def _extract_su2su2_prefactors(U, V):
-    r"""This function is used for the case of 2 CNOTs and 3 CNOTs. It does something
-    similar as the 1-CNOT case, but there is no special form for one of the
-    SO(4) operations.
-
-    Suppose U, V are SU(4) matrices for which there exists A, B, C, D such that
-    (A \otimes B) V (C \otimes D) = U. The problem is to find A, B, C, D in SU(2)
-    in an analytic and fully differentiable manner.
-
-    This decomposition is possible when U and V are in the same double coset of
-    SU(4), meaning there exists G, H in SO(4) s.t. G (Edag V E) H = (Edag U
-    E). This is guaranteed here by how V was constructed in both the
-    _decomposition_2_cnots and _decomposition_3_cnots methods.
-
-    Then, we can use the fact that E SO(4) Edag gives us something in SU(2) x
-    SU(2) to give A, B, C, D.
-    """
-
-    # A lot of the work here happens in the magic basis. Essentially, we
-    # don't look explicitly at some U = G V H, but rather at
-    #     E^\dagger U E = G E^\dagger V E H
-    # so that we can recover
-    #     U = (E G E^\dagger) V (E H E^\dagger) = (A \otimes B) V (C \otimes D).
-    # There is some math in the paper explaining how when we define U in this way,
-    # we can simultaneously diagonalize functions of U and V to ensure they are
-    # in the same coset and recover the decomposition.
-    u = math.dot(math.cast_like(E_dag, V), math.dot(U, math.cast_like(E, V)))
-    v = math.dot(math.cast_like(E_dag, V), math.dot(V, math.cast_like(E, V)))
-
-    uuT = math.dot(u, math.T(u))
-    vvT = math.dot(v, math.T(v))
-
-    # Get the p and q in SO(4) that diagonalize uuT and vvT respectively (and
-    # their eigenvalues). We are looking for a simultaneous diagonalization,
-    # which we know exists because of how U and V were constructed. Furthermore,
-    # The way we will do this is by noting that, since uuT/vvT are complex and
-    # symmetric, so both their real and imaginary parts share a set of
-    # real-valued eigenvectors, which are also eigenvectors of uuT/vvT
-    # themselves. So we can use eigh, which orders the eigenvectors, and so we
-    # are guaranteed that the p and q returned will be "in the same order".
-    _, p = math.linalg.eigh(math.real(uuT) + math.imag(uuT))
-    _, q = math.linalg.eigh(math.real(vvT) + math.imag(vvT))
-
-    # If determinant of p/q is not 1, it is in O(4) but not SO(4), and has determinant
-    # We can transform it to SO(4) by simply negating one of the columns.
-    p = math.dot(p, math.diag([1, 1, 1, math.sign(math.linalg.det(p))]))
-    q = math.dot(q, math.diag([1, 1, 1, math.sign(math.linalg.det(q))]))
-
-    # Now, we should have p, q in SO(4) such that p^T u u^T p = q^T v v^T q.
-    # TODO: no you don't, the eigenvectors returned from math.linalg.eigh are not guaranteed to
-    #       be in any particular order, especially when there is degeneracy. This means that p
-    #       and q are not necessarily "in the same order" as claimed above. This may sometimes
-    #       lead to incorrect results (see https://github.com/PennyLaneAI/pennylane/issues/5308)
-    #       The current solution is to add a small perturbation to U and V to break the potential
-    #       degeneracy. We should probably find a better algorithm at some point. [sc-89460]
-
-    # Then (v^\dag q p^T u)(v^\dag q p^T u)^T = I.
-    # So we can set G = p q^T, H = v^\dag q p^T u to obtain G v H = u.
-    G = math.dot(math.cast_like(p, 1j), math.T(q))
-    H = math.dot(math.conj(math.T(v)), math.dot(math.T(G), u))
-
-    # These are still in SO(4) though - we want to convert things into SU(2) x SU(2)
-    # so use the entangler. Since u = E^\dagger U E and v = E^\dagger V E where U, V
-    # are the target matrices, we can reshuffle as in the docstring above,
-    #     U = (E G E^\dagger) V (E H E^\dagger) = (A \otimes B) V (C \otimes D)
-    # where A, B, C, D are in SU(2) x SU(2).
-    AB = math.dot(math.cast_like(E, G), math.dot(G, math.cast_like(E_dag, G)))
-    CD = math.dot(math.cast_like(E, H), math.dot(H, math.cast_like(E_dag, H)))
-
-    # Now, we just need to extract the constituent tensor products.
-    A, B = math.decomposition.su2su2_to_tensor_products(AB)
-    C, D = math.decomposition.su2su2_to_tensor_products(CD)
-
-    return A, B, C, D
-
-
 def _decompose_0_cnots(U, wires, initial_phase):
     r"""If there are no CNOTs, this is just a tensor product of two single-qubit gates.
     We can perform that decomposition directly:
@@ -696,149 +611,262 @@ def _decompose_1_cnot(U, wires, initial_phase):
     return math.cast_like(-np.pi / 4, initial_phase)
 
 
-def _decompose_2_cnots(U, wires, initial_phase):
-    r"""If 2 CNOTs are required, we can write the circuit as
-     -╭U- = -A--╭X--RZ(d)--╭X--C-
-     -╰U- = -B--╰C--RX(p)--╰C--D-
-    We need to find the angles for the Z and X rotations such that the inner
-    part has the same spectrum as U, and then we can recover A, B, C, D.
-    """
+def _multidot(*matrices):
+    mat = matrices[0]
+    for m in matrices[1:]:
+        mat = math.dot(mat, m)
+    return mat
 
-    # Compute the rotation angles
-    u = math.dot(E_dag, math.dot(U, E))
-    gammaU = math.dot(u, math.T(u))
-    evs = math.linalg.eigvals(gammaU)
 
-    # These choices are based on Proposition III.3 of
-    # https://arxiv.org/abs/quant-ph/0308045
+def _real_imag_split_eigh(A, factor):
+    _, basis = math.linalg.eigh(math.real(A) / factor + factor * math.imag(A))
+    eigvals = _multidot(math.transpose(basis), A, basis)
+    return eigvals, basis
 
-    x = math.angle(evs[0])
 
-    # If it was the conjugate, grab a different eigenvalue.
-    y = math.cond(
-        math.allclose(evs[0], math.conj(evs[1])),
-        lambda: math.angle(evs[2]),
-        lambda: math.angle(evs[1]),
+def _ai_kak(U):
+    """Compute a type-AI Cartan decomposition of a unitary ``U`` in the standard basis/representation.
+    This is done in the following steps (see App.E of https://arxiv.org/abs/2503.19014, case AI):
+    - compute a real-valued eigenbasis O_1 (orthogonal matrix) and the eigenvalues d^2 of Δ = U U^T,
+    - make O_1 *special* orthogonal by adjusting the sign of its first column,
+    - take the square root of d^2 to obtain d,
+    - compute O_2 = d* O_1^T U,
+    - adjust the sign of the first row of O_2 so that O_2 is special orthogonal, and multiply
+      the first entry of d with the same sign.
+    """
+
+    # Delta is symmetric (Delta^T=Delta) and unitary (because U is) but not real-valued.
+    Delta = math.dot(U, math.transpose(U))
+    # Denote the real and imaginary parts of Delta as R and I, respectively.
+    # Delta^T=Delta ==> R^T=R and I^T=I.  (1)
+    # Delta unitary ==> (R+iI)(R^T-iI^T)=id ==> RR^T+II^T=id and IR^T-RI^T = 0 (2)
+    # Combine (1) and second equation in (2): 0=IR-RI=[I,R] ==> I and R share an eigenbasis
+    # ==> I, R and Delta=R+iI share an eigenbasis.
+    #
+    # We need to make sure that for degenerate real or imaginary parts
+    # we actually find an eigenbasis that also is one for Delta. In a first step, we make this
+    # likely by weighting the real and imaginary part by 1/pi and pi. In a second step,
+    # we check whether this basis diagonalized Delta, and recompute with a new weighting (by 10)
+    # if it did not.
+    d_squared, o1 = _real_imag_split_eigh(Delta, np.pi)
+    d_squared, o1 = math.cond(
+        math.allclose(math.diag(math.diag(d_squared)), d_squared, atol=1e-7),
+        lambda: (d_squared, o1),
+        lambda: _real_imag_split_eigh(Delta, 10.0),
         (),
     )
 
-    delta = (x + y) / 2
-    phi = (x - y) / 2
-
-    # need to perturb x by 5 precision to avoid a discontinuity at a special case.
-    # see https://github.com/PennyLaneAI/pennylane/issues/5308
-    precision = math.finfo(delta.dtype).eps
-    RZd = ops.RZ.compute_matrix(math.cast_like(delta + 5 * precision, 1j))
-    RXp = ops.RX.compute_matrix(math.cast_like(phi + 5 * precision, 1j))
-    inner_u = math.kron(RZd, RXp)
+    # This implements o1[:, 0] *= det(o1) to ensure det(o1) = 1 afterwards
+    # No need to modify the eigenvalues or d because this change will be absorbed in o2
+    o1 = math.transpose(math.set_index(math.transpose(o1), 0, math.linalg.det(o1) * o1[:, 0]))
+
+    d = math.diag(math.sqrt(math.diag(d_squared)))
+    o2 = _multidot(math.conj(d), math.transpose(o1), U)
+
+    # Instead of guaranteeing the correct determinant while taking the square root,
+    # we correct it after the fact
+    # This implements o2[0] *= det(o2) to ensure det(o2) = 1 afterwards
+    # Here we need to also adapt d because o1 already is fixed, so it can not absorb det(o2)
+    det_o2 = math.linalg.det(o2)
+    o2 = math.set_index(o2, 0, det_o2 * o2[0])
+    d = math.set_index(d, 0, det_o2 * d[0])
+
+    return o1, d, o2
+
+
+def _extract_abde(A):
+    """Extract the parameters for the central part of a 3-CNOT circuit as well
+    as a global phase. See documentation of _decompose_3_cnots for details.
+
+    The math.cond calls are necessary in order to decide from which matrix entries to extract
+    the angles.
+    The input matrix is expected to be of the form (called C1 in _decompose_3_cnots)
+
+    exp(-i d') cos(a')           0                  0          -exp(-i d') sin(a')
+              0         exp(i e') sin(b')   exp(i e') cos(b')            0
+              0         exp(i e') cos(b')  -exp(i e') sin(b')            0
+    exp(-i d') sin(a')           0                  0           exp(-i d') cos(a')
+
+    and this function will return
+
+    a = a' + b'
+    b = a' - b'
+    d = d' + e'
+    e = (d' - e') / 2.
+
+    The performed computation steps are the following:
+    1. Compute a' from A_{00} and A_{30}
+        a. If cos(d')sin(a')≠0 or cos(d')cos(a')≠0, use
+           atan2(cos(d')sin(a'), cos(d')cos(a')) = atan2(sin(a'), cos(a')) = a'
+        b. If cos(d')sin(a')=cos(d')cos(a')=0, use
+           atan2(sin(d')sin(a'), sin(d')cos(a')) = atan2(sin(a'), cos(a')) = a'
+        Note that if cos(d')sin(a')=cos(d')cos(a')=0, we know that cos(d')=0 and sin(d')≠0
+    2. Compute b' from A_{11} and A_{21}
+        a. If cos(e')sin(b')≠0 or cos(e')cos(b')≠0, use
+           atan2(cos(e')sin(b'), cos(e')cos(b')) = atan2(sin(b'), cos(b')) = b'
+        b. If cos(e')sin(b')=cos(e')cos(b')=0, use
+           atan2(sin(e')sin(b'), sin(e')cos(b')) = atan2(sin(b'), cos(b')) = b'
+        Note that if cos(e')sin(b')=cos(e')cos(b')=0, we know that cos(e')=0 and sin(e')≠0
+    3. Compute a = a' + b' and b = a' - b'
+    4. Compute exp(-id') from A_{00} or A_{30}
+        a. If A_{00}≠0, use it and compute exp(-i d')cos(a') / cos(a')
+        b. If A_{00}=0, us A_{30} and compute exp(-i d')sin(a') / sin(a')
+    5. Compute exp(ie') from A_{11} or A_{21}
+        a. If A_{21}≠0, use it and compute exp(i e')cos(b') / cos(b')
+        b. If A_{21}=0, us A_{11} and compute exp(i e')sin(b') / sin(b')
+    6. Compute d = d' + e' as arg(exp(ie') exp(-id')*)
+    7. Compute e = (d' - e')/2 as -arg(exp(ie') exp(-id/2))
+    """
+    a_plus_b_half = math.cond(
+        math.allclose(math.real(A[::3, 0]), math.zeros_like(math.real(A[::3, 0]))),
+        lambda: math.arctan2(math.imag(A[3, 0]), math.imag(A[0, 0])),
+        lambda: math.arctan2(math.real(A[3, 0]), math.real(A[0, 0])),
+        (),
+    )
+    a_minus_b_half = math.cond(
+        math.allclose(math.real(A[1:3, 1]), math.zeros_like(math.real(A[1:3, 1]))),
+        lambda: math.arctan2(math.imag(A[1, 1]), math.imag(A[2, 1])),
+        lambda: math.arctan2(math.real(A[1, 1]), math.real(A[2, 1])),
+        (),
+    )
+    a = a_plus_b_half + a_minus_b_half
+    b = a_plus_b_half - a_minus_b_half
 
-    # We need the matrix representation of this interior part, V, in order to
-    # decompose U = (A \otimes B) V (C \otimes D)
-    V = math.dot(math.cast_like(CNOT10, U), math.dot(inner_u, math.cast_like(CNOT10, U)))
+    apb_frac = math.cond(
+        math.isclose(A[0, 0], math.zeros_like(A[0, 0])),
+        lambda: A[3, 0] / math.cast_like(math.sin(a_plus_b_half), A[3, 0]),
+        lambda: A[0, 0] / math.cast_like(math.cos(a_plus_b_half), A[0, 0]),
+        (),
+    )
+    amb_frac = math.cond(
+        math.isclose(A[2, 1], math.zeros_like(A[2, 1])),
+        lambda: A[1, 1] / math.cast_like(math.sin(a_minus_b_half), A[1, 1]),
+        lambda: A[2, 1] / math.cast_like(math.cos(a_minus_b_half), A[2, 1]),
+        (),
+    )
 
-    # Now we find the A, B, C, D in SU(2), and return the decomposition
-    A, B, C, D = _extract_su2su2_prefactors(U, V)
+    d = math.angle(amb_frac * math.conj(apb_frac))
+    e = -math.angle(amb_frac * math.exp(-1j * math.cast_like(d / 2, 1j)))
+    return a, b, d, e
 
-    ops.QubitUnitary(C, wires[0])
-    ops.QubitUnitary(D, wires[1])
-    ops.CNOT(wires=[wires[1], wires[0]])
-    ops.RZ(delta, wires=wires[0])
-    ops.RX(phi, wires=wires[1])
-    ops.CNOT(wires=[wires[1], wires[0]])
-    ops.QubitUnitary(A, wires[0])
-    ops.QubitUnitary(B, wires[1])
 
-    return math.cast_like(0, initial_phase)
+def _central_circuit(a, b, d, wires):
+    """Central part of 3-CNOT circuit."""
+    ops.CNOT([wires[1], wires[0]])
+    ops.RZ(d, wires[0])
+    ops.RY(b, wires[1])
+    ops.CNOT(wires)
+    ops.RY(a, wires[1])
+    ops.CNOT([wires[1], wires[0]])
 
 
 def _decompose_3_cnots(U, wires, initial_phase):
-    r"""The most general form of this decomposition is U = (A \otimes B) V (C \otimes D),
-    where V is as depicted in the circuit below:
-     -╭U- = -C--╭X--RZ(d)--╭C---------╭X--A-
-     -╰U- = -D--╰C--RY(b)--╰X--RY(a)--╰C--B-
+    """Decompose a unitary 4x4 matrix into a 3-CNOT circuit.
+
+    From a mathematical perspective, this decomposition mainly is from one matrix into
+    a product of three matrices, which is an instance of a Cartan, or KAK, decomposition.
+    The Cartan decomposition is of type AI, decomposing a (special) unitary matrix into
+    two special orthogonal matrices and a matrix from some representation of U(1)^r, with r=4
+    for unitary input (used here) and r=3 for special unitary input.
+    See e.g. App. E, paragraph on AI in https://arxiv.org/abs/2503.19014 for details.
+    This Cartan decomposition is implemented in _ai_kak. Here we take care
+    of translating the input and output matrices of _ai_kak into the right representations
+    of SO(4) and U(1)^4.
+    The representations that we want are given by the fixed circuit structure that we are after:
+    An arbitrary special unitary single-qubit operation on each qubit necessarily is from
+    SU(2) x SU(2), which is isomorphic to SO(4) (a so-called accidental or exceptional
+    isomorphism). The central circuit part
+
+    0: ─╭X──RZ(d)─╭●────────╭X──GlobalPhase(e)─┤
+    1: ─╰●──RY(b)─╰X──RY(a)─╰●──GlobalPhase(e)─┤
+
+    forms a representation of U(1)^4, given by the matrix (call it C1)
+
+    exp(-i d') cos(a')           0                  0          -exp(-i d') sin(a')
+              0         exp(i e') sin(b')   exp(i e') cos(b')            0
+              0         exp(i e') cos(b')  -exp(i e') sin(b')            0
+    exp(-i d') sin(a')           0                  0           exp(-i d') cos(a')
+
+    where
+    a' = a/2 + b/2
+    b' = a/2 - b/2
+    d' = d/2 + e
+    e' = d/2 - e.
+
+    Now, as is used throughout the two-qubit decompositions in this file, the transformation
+    between the canonical representation of SO(4) (real matrices with OO^T = 1 and determinant 1)
+    and the representation as single-qubit unitaries on both qubits is given by the so-called
+    magic basis E. That is, for O in the canonical representation, E O E† is of the form A⊗B with
+    A and B 2x2 special unitary matrices.
+    Simultaneously, E transforms diagonal unitary matrices into matrices of the form
+
+      exp(i t) cos(r)         0                  0         i exp(i t) sin(r)
+              0         exp(i u) cos(s)  i exp(i u) sin(s)         0
+              0       i exp(i u) sin(s)    exp(i u) cos(s)         0
+    i exp(i t) sin(r)         0                  0           exp(i t) cos(r)
+
+    generated by the "Cartan coordinate" generators X⊗X, Y⊗Y, Z⊗Z (and I⊗I). Call this matrix C2.
+
+    Finally, note that C2 can be transformed into C1 via the following static gates:
+
+    C1 = SWAP S† C2 S.
+
+    (S is just `qml.S`)
+    Now, we "just" need to combine all of these basis changes with the type-AI Cartan
+    decomposition (_ai_kak) and a function that extracts the parameters a, b, d, e from a matrix
+    of the form C1 (_extract_abde). For this, let's compute (not necessarily obvious to come up
+    with but easy to verify).
+
+    V := S SWAP U S†    (bookkeeping un-transformation)
+    W := E† V E         (magic basis un-rotation)
+    W =: K_1 A_K K_2    (computed by _ai_kak; K_i∈SO(4), A_K diagonal)
+
+    L_1 := E K_1 E†             (L_1 is of form A⊗B)
+    L_2 := E K_2 E†             (L_2 is of form A⊗B)
+    A_L := E A_K E†             (A_L is of form C2)
+
+    M_1 := SWAP S† L_1 S SWAP   (M_1 is still of the form A⊗B, with new A, B)
+    M_2 := S† L_2 S             (M_2 is still of the form A⊗B, with new A, B)
+    A_M := SWAP S† A_L S        (A_M is of the form C1)
+
+    Now we can extract a, b, d, e from A_M because it is of form C1, implemented by the
+    central circuit part. Also, we can decompose M_1 and M_2 into two single-qubit
+    unitaries via _decompose_0_cnots.
+    To verify correctness, compute
+
+    M_1 A_M M_2
+    = (SWAP S† L_1 S SWAP) (SWAP S A_L S†) (S† L_2 S)
+    = (SWAP S† L_1 A_L L_2 S)
+    = (SWAP S† (E K_1 E†) (E A_K E†) (E K_2 E†) S)
+    = (SWAP S† E (K_1 A_K K_2) E† S)
+    = (SWAP S† E W E† S)
+    = (SWAP S† E (E† V E) E† S)
+    = (SWAP S† (S SWAP U S†) S)
+    = U
+
+    So we actually implemented U!
     """
+    W = _multidot(E_dag, S_0_dag, SWAP, U, S_0, E)
+    K_1, A_K, K_2 = _ai_kak(W)
 
-    # First we add a SWAP as per v1 of arXiv:0308033, which helps with some
-    # rearranging of gates in the decomposition (it will cancel out the fact
-    # that we need to add a SWAP to fix the determinant in another part later).
-    swap_U = np.exp(1j * np.pi / 4) * math.dot(math.cast_like(SWAP, U), U)
+    L_1 = _multidot(E, K_1, E_dag)
+    A_L = _multidot(E, A_K, E_dag)
+    L_2 = _multidot(E, K_2, E_dag)
 
-    # Choose the rotation angles of RZ, RY in the two-qubit decomposition.
-    # They are chosen as per Proposition V.1 in quant-ph/0308033 and are based
-    # on the phases of the eigenvalues of :math:`E^\dagger \gamma(U) E`, where
-    #    \gamma(U) = (E^\dag U E) (E^\dag U E)^T.
-    # The rotation angles can be computed as follows (any three eigenvalues can be used)
-    u = math.dot(E_dag, math.dot(swap_U, E))
-    gammaU = math.dot(u, math.T(u))
-    # !Note: [sc-89460]
-    evs, _ = math.linalg.eig(gammaU)
-
-    x = math.angle(evs[0])
-    y = math.angle(evs[1])
-    z = math.angle(evs[2])
-
-    # Compute functions of the eigenvalues; there are different options in v1
-    # vs. v3 of the paper, I'm not entirely sure why. This is the version from v3.
-    alpha = (x + y) / 2
-    beta = (x + z) / 2
-    delta = (z + y) / 2
-
-    # We need the matrix representation of this interior part, V, in order to
-    # decompose U = (A \otimes B) V (C \otimes D)
-    #
-    # Looking at the decomposition above, V has determinant -1 (because there
-    # are 3 CNOTs, each with determinant -1). The relationship between U and V
-    # requires that both are in SU(4), so we add a SWAP after to V. We will see
-    # how this gets fixed later.
-    #
-    # -╭V- = -╭X--RZ(d)--╭C---------╭X--╭SWAP-
-    # -╰V- = -╰C--RY(b)--╰X--RY(a)--╰C--╰SWAP-
-
-    EPS = math.finfo(delta.dtype).eps
-    RZd = ops.RZ.compute_matrix(math.cast_like(delta + 5 * EPS, 1j))
-    RYb = ops.RY.compute_matrix(beta + 5 * EPS)
-    RYa = ops.RY.compute_matrix(alpha)
-
-    V_mats = [
-        CNOT10,
-        math.kron(RZd, RYb),
-        CNOT01,
-        math.kron(math.eye(2), RYa),
-        CNOT10,
-        SWAP,
-    ]
-    V = math.convert_like(math.eye(4), U)
-    for mat in V_mats:
-        V = math.dot(math.cast_like(mat, U), V)
-
-    # Now we need to find the four SU(2) operations A, B, C, D
-    A, B, C, D = _extract_su2su2_prefactors(swap_U, V)
-
-    # At this point, we have the following:
-    # -╭U-╭SWAP- = --C--╭X-RZ(d)-╭C-------╭X-╭SWAP--A
-    # -╰U-╰SWAP- = --D--╰C-RZ(b)-╰X-RY(a)-╰C-╰SWAP--B
-    #
-    # Using the relationship that SWAP(A \otimes B) SWAP = B \otimes A,
-    # -╭U-╭SWAP- = --C--╭X-RZ(d)-╭C-------╭X--B--╭SWAP-
-    # -╰U-╰SWAP- = --D--╰C-RZ(b)-╰X-RY(a)-╰C--A--╰SWAP-
-    #
-    # Now the SWAPs cancel, giving us the desired decomposition
-    # (up to a global phase).
-    # -╭U- = --C--╭X-RZ(d)-╭C-------╭X--B--
-    # -╰U- = --D--╰C-RZ(b)-╰X-RY(a)-╰C--A--
-
-    ops.QubitUnitary(C, wires[0])
-    ops.QubitUnitary(D, wires[1])
-    ops.CNOT(wires=[wires[1], wires[0]])
-    ops.RZ(delta, wires=wires[0])
-    ops.RY(beta, wires=wires[1])
-    ops.CNOT(wires=wires)
-    ops.RY(alpha, wires=wires[1])
-    ops.CNOT(wires=[wires[1], wires[0]])
-    ops.QubitUnitary(A, wires[1])
-    ops.QubitUnitary(B, wires[0])
+    M_1 = _multidot(SWAP, S_0, L_1, S_0_dag, SWAP)
+    M_2 = _multidot(S_0, L_2, S_0_dag)
+    A_M = _multidot(SWAP, S_0, A_L, S_0_dag)
 
-    return math.cast_like(-np.pi / 4, initial_phase)
+    a, b, d, e = _extract_abde(A_M)
+
+    _decompose_0_cnots(M_2, wires, 0.0)
+    _central_circuit(a, b, d, wires)
+    _decompose_0_cnots(M_1, wires, 0.0)
+    # global phases here are zero because we are guaranteed that M_1 and M_2 have unit determinant
+
+    # Return the global phase obtained from A_M. It will be combined with initial_phase
+    # in `two_qubit_decomposition`
+    return math.cast_like(-e, initial_phase)
 
 
 def _compute_udv(a, b):
@@ -911,6 +939,8 @@ def _cossin_decomposition(U, p):
 
 
 def _global_phase(phase):
+    """Calls the GlobalPhase with a negative sign. Used in a ``cond`` call which requires
+    no variable to be returned."""
     ops.GlobalPhase(-phase)
 
 
diff --git a/pennylane/ops/qubit/matrix_ops.py b/pennylane/ops/qubit/matrix_ops.py
index 9019d4c74..c6849bdf0 100644
--- a/pennylane/ops/qubit/matrix_ops.py
+++ b/pennylane/ops/qubit/matrix_ops.py
@@ -302,7 +302,8 @@ class QubitUnitary(Operation):
     # pylint: disable=arguments-renamed, invalid-overridden-method
     @property
     def has_decomposition(self) -> bool:
-        return len(self.wires) < 3 if self.has_matrix else len(self.wires) == 1
+        # We are unable to decompose sparse matrices larger than 1 qubit.
+        return self.has_matrix or len(self.wires) == 1
 
     def adjoint(self) -> "QubitUnitary":
         if self.has_matrix:


Test patch:
diff --git a/tests/ops/op_math/test_controlled_decompositions.py b/tests/ops/op_math/test_controlled_decompositions.py
index ddb51ede9..0793d0740 100644
--- a/tests/ops/op_math/test_controlled_decompositions.py
+++ b/tests/ops/op_math/test_controlled_decompositions.py
@@ -1041,21 +1041,3 @@ class TestMCXDecomposition:
             _ = _decompose_mcx_with_two_workers(
                 control_wires, target_wire, work_wires, work_wire_type="clean"
             )
-
-
-def test_ControlledQubitUnitary_has_decomposition_correct():
-    """Test that ControlledQubitUnitary reports has_decomposition=False if it is False"""
-    U = qml.Toffoli(wires=[0, 1, 2]).matrix()
-    op = qml.ControlledQubitUnitary(U, wires=[0, 1, 2, 3])
-
-    assert not op.has_decomposition
-    with pytest.raises(qml.operation.DecompositionUndefinedError):
-        op.decomposition()
-
-
-def test_ControlledQubitUnitary_has_decomposition_super_False(mocker):
-    """Test that has_decomposition returns False if super() returns False"""
-    spy = mocker.spy(qml.QueuingManager, "stop_recording")
-    op = qml.ControlledQubitUnitary(np.diag((1.0,) * 8), wires=[0, 1, 2, 3, 4])
-    assert not op.has_decomposition
-    spy.assert_not_called()
diff --git a/tests/ops/op_math/test_decompositions.py b/tests/ops/op_math/test_decompositions.py
index 2ac083a78..0dc191bd4 100644
--- a/tests/ops/op_math/test_decompositions.py
+++ b/tests/ops/op_math/test_decompositions.py
@@ -630,6 +630,32 @@ samples_3_cnots = [
             -0.4168424055846296 + 0.3576772276508062j,
         ],
     ],
+    [  # this matrix caused issues before: https://github.com/PennyLaneAI/pennylane/issues/7467
+        [
+            -0.40836553 - 0.28850926j,
+            0.23864466 - 0.4393731j,
+            0.02190401 + 0.06284322j,
+            -0.23169726 - 0.66474606j,
+        ],
+        [
+            -0.48865844 - 0.10589111j,
+            0.49770294 + 0.04787262j,
+            -0.02177958 + 0.45019526j,
+            -0.02632749 + 0.54420285j,
+        ],
+        [
+            0.49276489 - 0.08475118j,
+            0.47943096 - 0.14193644j,
+            0.51448332 - 0.17932317j,
+            -0.42560958 + 0.14834623j,
+        ],
+        [
+            0.42041002 - 0.27065738j,
+            -0.31807806 - 0.38578018j,
+            -0.03401683 + 0.70314565j,
+            -0.00321586 + 0.06647341j,
+        ],
+    ],
 ]
 
 samples_2_cnots = [
@@ -987,13 +1013,14 @@ class TestTwoQubitUnitaryDecomposition:
     @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) == 8
+        assert len(obtained_decomposition) == 11  # 8 # 8 would be the count with 2-CNOT circuit
 
         tape = qml.tape.QuantumScript(obtained_decomposition)
         obtained_matrix = qml.matrix(tape, wire_order=wires)
@@ -1545,13 +1572,7 @@ class TestQubitUnitaryDecompositionGraph:
     @pytest.mark.parametrize(
         "U, n_wires",
         [
-            pytest.param(
-                qml.matrix(qml.CRX(0.123, [0, 2]) @ qml.CRY(0.456, [1, 3])),
-                4,
-                marks=pytest.mark.xfail(
-                    strict=False, reason="https://github.com/PennyLaneAI/pennylane/issues/7467"
-                ),
-            ),
+            (qml.matrix(qml.CRX(0.123, [0, 2]) @ qml.CRY(0.456, [1, 3])), 4),
             (qml.QFT.compute_matrix(5), 5),
             (qml.GroverOperator.compute_matrix(6, []), 6),
         ],
diff --git a/tests/templates/test_subroutines/test_qmc.py b/tests/templates/test_subroutines/test_qmc.py
index 4a68de9b0..c1898280e 100644
--- a/tests/templates/test_subroutines/test_qmc.py
+++ b/tests/templates/test_subroutines/test_qmc.py
@@ -302,8 +302,8 @@ class TestQuantumMonteCarlo:
         # Do expansion in two steps to avoid also decomposing the first QubitUnitary
         queue_before_qpe = tape.operations[:2]
 
-        # 2-qubit decomposition has 11 operations, and after is a 3-qubit gate so start at 11
-        queue_after_qpe = tape.expand().operations[11:]
+        # Build a new tape from all operations following the two QubitUnitary ops and expand it
+        queue_after_qpe = qml.tape.QuantumScript(tape.operations[2:]).expand().operations
 
         A = probs_to_unitary(p)
         R = func_to_unitary(self.func, 4)
diff --git a/tests/transforms/test_cliffordt_transform.py b/tests/transforms/test_cliffordt_transform.py
index eb7019556..0fe38aa4c 100644
--- a/tests/transforms/test_cliffordt_transform.py
+++ b/tests/transforms/test_cliffordt_transform.py
@@ -410,7 +410,10 @@ class TestCliffordCompile:
     def test_raise_with_cliffordt_decomposition(self):
         """Test that exception is correctly raise when decomposing gates without any decomposition"""
 
-        tape = qml.tape.QuantumScript([qml.QubitUnitary(qml.math.eye(8), wires=[0, 1, 2])])
+        class CustomOp(qml.operation.Operation):  # pylint: disable=too-few-public-methods
+            pass
+
+        tape = qml.tape.QuantumScript([CustomOp(wires=[0, 1, 2])])
 
         with pytest.raises(ValueError, match="Cannot unroll"):
             clifford_t_decomposition(tape)
diff --git a/tests/transforms/test_unitary_to_rot.py b/tests/transforms/test_unitary_to_rot.py
index 06143bea1..2e0fbf3bd 100644
--- a/tests/transforms/test_unitary_to_rot.py
+++ b/tests/transforms/test_unitary_to_rot.py
@@ -18,6 +18,7 @@ from itertools import product
 
 import pytest
 from gate_data import H, I, S, T, X, Z
+from scipy.stats import unitary_group
 from test_optimization.utils import check_matrix_equivalence
 
 import pennylane as qml
@@ -501,15 +502,18 @@ test_two_qubit_unitaries = [
 
 
 @pytest.mark.parametrize("num_reps", [1, 2, 3, 4, 5])
-def test_unitary_to_rot_multiple_two_qubit(num_reps):
+def test_unitary_to_rot_multiple_two_qubit(num_reps, seed):
     """Test that numerous two-qubit unitaries can be decomposed sequentially."""
 
     dev = qml.device("default.qubit", wires=2)
 
-    U = np.array(test_two_qubit_unitaries[1], dtype=np.complex128)
+    Us = unitary_group.rvs(4, size=num_reps, random_state=seed)
+    if num_reps == 1:
+        # Batch dim is squeezed away by rvs() for size = 1
+        Us = [Us]
 
     def my_circuit():
-        for _ in range(num_reps):
+        for U in Us:
             qml.QubitUnitary(U, wires=[0, 1])
         return qml.expval(qml.PauliZ(0))