[BUG] `two_qubit_decomposition` can return a non-equivalent local circuit for a slightly entangling two-qubit unitary

[YujinSong-hep]

Expected behavior

qml.ops.two_qubit_decomposition documents that it accepts a 4 x 4 unitary matrix and returns a list of operations representing the decomposition of that matrix.

For any valid two-qubit unitary input, recomposing the returned operations should recover the input unitary up to global phase.

In particular, a unitary with small but nonzero entangling content should not be silently replaced by a purely local tensor-product circuit unless the API explicitly documents approximate decomposition with an error budget.

Actual behavior

For the valid two-qubit unitary constructed below, qml.ops.two_qubit_decomposition returns only two single-qubit QubitUnitary operations plus a GlobalPhase:

returned ops: ['QubitUnitary', 'QubitUnitary', 'GlobalPhase']
num CNOTs: 0
phase-aligned error: 0.0023561653135340064

The returned circuit is therefore not equivalent to the input unitary, even after global phase alignment.

A nearby sanity case with larger entangling angles reconstructs to near machine precision under the same recomposition check, so the failure is not caused by the test oracle or wire ordering.

Additional information

The failure appears to come from the CNOT-count classifier in pennylane/ops/op_math/decompositions/unitary_decompositions.py::_compute_num_cnots. The code uses a fixed tolerance,

math.allclose(trace, 4, atol=1e-7) | math.allclose(trace, -4, atol=1e-7)

to classify the 0-CNOT/tensor-product case. In this example, that tolerance routes a weakly but genuinely entangling unitary into the 0-CNOT branch.

Source code

import numpy as np
import pennylane as qml


def random_su2(rng):
    X = rng.normal(size=(2, 2)) + 1j * rng.normal(size=(2, 2))
    Q, R = np.linalg.qr(X)
    d = np.diag(R)
    Q = Q * (d / np.abs(d))
    Q /= np.linalg.det(Q) ** 0.5
    return Q


def compose_ops(ops):
    U = np.eye(4, dtype=complex)
    for op in ops:
        U = qml.matrix(op, wire_order=(0, 1)) @ U
    return U


def phase_align_error(U, V):
    overlap = np.vdot(V.reshape(-1), U.reshape(-1))
    phase = overlap / abs(overlap)
    return np.linalg.norm(U - phase * V)


rng = np.random.default_rng(123)
A, B, C, D = [random_su2(rng) for _ in range(4)]

ops = [
    qml.QubitUnitary(A, 0),
    qml.QubitUnitary(B, 1),
    qml.CNOT((1, 0)),
    qml.RZ(8.471910266820082e-04, 0),
    qml.RX(-1.3556395803158893e-03, 1),
    qml.CNOT((1, 0)),
    qml.QubitUnitary(C, 0),
    qml.QubitUnitary(D, 1),
]

U = compose_ops(ops)
decomp = qml.ops.two_qubit_decomposition(U, wires=(0, 1))
U2 = compose_ops(decomp)

print("returned ops:", [op.name for op in decomp])
print("num CNOTs:", sum(op.name == "CNOT" for op in decomp))
print("phase-aligned error:", phase_align_error(U, U2))

Tracebacks

Output：

returned ops: ['QubitUnitary', 'QubitUnitary', 'GlobalPhase']
num CNOTs: 0
phase-aligned error: 0.0023561653135340064

System information

Version: 0.45.0
Platform info: macOS
Python version: 3.13.13

Existing GitHub issues

• I have searched existing GitHub issues to make sure the issue does not already exist.

[astralcai]

I'd say that for a 4x4 unitary, a Frobenius distance of ~0.0023 is sufficiently small that we can treat them as virtually identical.

[astralcai]

Although I agree that 1e-7 is a rather arbitrary tolerance here. We should use math.finfo(U.dtype).eps here instead.