[BUG] Precision error in two-qubit decomposition of qml.QubitUnitary

[abhishekabhishek]

Expected behavior

The two-qubit decomposition for qml.QubitUnitary should be valid upto a high-precision.

Actual behavior

The two-qubit decomposition for qml.QubitUnitary is only valid upto around 1e-7. This leads to an error in matrix equivalence checking if the atol is not specified since default atol = 1e-8 in np.allclose.

Additional information

This seems related to the now closed issue #5308 and its corresponding fix PR #5448.

Source code

from pprint import pprint
import numpy as np
from qiskit.circuit.library import CU1Gate
from qiskit import QuantumCircuit
from pennylane_qiskit import load
from pennylane.tape import QuantumTape

op = CU1Gate(7.)
qc = QuantumCircuit(op.num_qubits)
qc.append(op, list(range(op.num_qubits)))

with QuantumTape() as exp_tape:
    load(qc)()

print(exp_tape.operations)
[QubitUnitary(array([[1.        +0.j       , 0.        +0.j       ,
        0.        +0.j       , 0.        +0.j       ],
       [0.        +0.j       , 1.        +0.j       ,
        0.        +0.j       , 0.        +0.j       ],
       [0.        +0.j       , 0.        +0.j       ,
        1.        +0.j       , 0.        +0.j       ],
       [0.        +0.j       , 0.        +0.j       ,
        0.        +0.j       , 0.75390225+0.6569866j]]), wires=[0, 1])]

U = exp_tape.operations[0].parameters[0]
print(U)
print(np.allclose(U.conj().T, np.linalg.inv(U)))
[[1.        +0.j        0.        +0.j        0.        +0.j
  0.        +0.j       ]
 [0.        +0.j        1.        +0.j        0.        +0.j
  0.        +0.j       ]
 [0.        +0.j        0.        +0.j        1.        +0.j
  0.        +0.j       ]
 [0.        +0.j        0.        +0.j        0.        +0.j
  0.75390225+0.6569866j]]
True

op_dcmp = qml.QubitUnitary(U, range(2)).decomposition()
tape = qml.tape.QuantumTape(op_dcmp, [])
pprint(tape.operations)
[[ 0.984-0.178j  0.   +0.j    -0.   +0.j     0.   -0.j   ]
 [ 0.   +0.j     0.984-0.178j -0.   -0.j    -0.   +0.j   ]
 [ 0.   +0.j    -0.   -0.j     0.984-0.178j  0.   -0.j   ]
 [-0.   -0.j     0.   +0.j     0.   -0.j     0.859+0.512j]]
True
[[ 9.83985947e-01-1.78246056e-01j  8.78875830e-25+8.65699404e-25j
  -2.07359532e-08+3.35585772e-08j  1.42132272e-17-2.38402213e-17j]
 [ 8.78875830e-25+2.85377232e-25j  9.83985947e-01-1.78246056e-01j
  -4.94732188e-18-2.73110964e-17j -1.81006909e-08+1.90109280e-08j]
 [ 3.11901081e-08+2.41523294e-08j -4.94732188e-18-2.73110964e-17j
   9.83985947e-01-1.78246056e-01j  1.12670198e-24-5.02395096e-25j]
 [-4.94732188e-18-2.73110964e-17j  1.02817984e-08+2.41523294e-08j
   7.22923602e-25-5.75538318e-25j  8.58934493e-01+5.12085477e-01j]]
[RZ(8.033185307179586, wires=[0]),
 RY(3.1415925939851483, wires=[0]),
 RZ(1.3915926535897927, wires=[0]),
 RZ(1.7499999999999996, wires=[1]),
 RY(0.0, wires=[1]),
 RZ(1.7499999999999996, wires=[1]),
 CNOT(wires=[1, 0]),
 RZ(0.35840734641020666, wires=[0]),
 RX(5.551115123125783e-17, wires=[1]),
 CNOT(wires=[1, 0]),
 RZ(9.42477796076938, wires=[0]),
 RY(3.141592611442945, wires=[0]),
 RZ(3.141592653589793, wires=[0]),
 RZ(1.5707963267948966, wires=[1]),
 RY(0.0, wires=[1]),
 RZ(1.5707963267948966, wires=[1])]

U_dcmp = qml.matrix(tape, wire_order=list(range(2)))
print(U_dcmp)
[[ 9.83985947e-01-1.78246056e-01j  8.78875830e-25+8.65699404e-25j
  -2.07359532e-08+3.35585772e-08j  1.42132272e-17-2.38402213e-17j]
 [ 8.78875830e-25+2.85377232e-25j  9.83985947e-01-1.78246056e-01j
  -4.94732188e-18-2.73110964e-17j -1.81006909e-08+1.90109280e-08j]
 [ 3.11901081e-08+2.41523294e-08j -4.94732188e-18-2.73110964e-17j
   9.83985947e-01-1.78246056e-01j  1.12670198e-24-5.02395096e-25j]
 [-4.94732188e-18-2.73110964e-17j  1.02817984e-08+2.41523294e-08j
   7.22923602e-25-5.75538318e-25j  8.58934493e-01+5.12085477e-01j]]

print(np.allclose(U_dcmp.conj().T, np.linalg.inv(U_dcmp)))
True

print(np.round(U_dcmp, decimals=3))
[[ 0.984-0.178j  0.   +0.j    -0.   +0.j     0.   -0.j   ]
 [ 0.   +0.j     0.984-0.178j -0.   -0.j    -0.   +0.j   ]
 [ 0.   +0.j    -0.   -0.j     0.984-0.178j  0.   -0.j   ]
 [-0.   -0.j     0.   +0.j     0.   -0.j     0.859+0.512j]]

product = qml.math.dot(qml.math.conj(qml.math.T(U_dcmp)), U)
product = product / product[0, 0]
print(product)
[[ 1.00000000e+00+0.00000000e+00j  8.13934100e-25-4.37463336e-25j
   2.63855706e-08-2.93250665e-08j -1.82350412e-17+2.09249910e-17j]
 [ 7.10493962e-25-1.00849220e-24j  1.00000000e+00+1.20168824e-15j
   1.21285632e-32+2.77555756e-17j  2.11994484e-08-1.54801092e-08j]
 [-2.63855706e-08-2.93250665e-08j  3.79017600e-33+2.77555756e-17j
   1.00000000e+00-1.88446565e-15j  3.26219203e-25+8.64548391e-25j]
 [ 1.82350412e-17+2.09249910e-17j -2.11994484e-08-1.54801092e-08j
   1.19820886e-24+2.93519530e-25j  1.00000000e+00-3.55044253e-16j]]

print(np.round(product, decimals=8))
[[ 1.e+00+0.e+00j  0.e+00-0.e+00j  3.e-08-3.e-08j -0.e+00+0.e+00j]
 [ 0.e+00-0.e+00j  1.e+00+0.e+00j  0.e+00+0.e+00j  2.e-08-2.e-08j]
 [-3.e-08-3.e-08j  0.e+00+0.e+00j  1.e+00-0.e+00j  0.e+00+0.e+00j]
 [ 0.e+00+0.e+00j -2.e-08-2.e-08j  0.e+00+0.e+00j  1.e+00-0.e+00j]]

qml.math.allclose(product, qml.math.eye(U.shape[0]), atol=1e-7)
True

qml.math.allclose(product, qml.math.eye(U.shape[0]), atol=1e-8)
False

Tracebacks

No response

System information

Name: PennyLane
Version: 0.36.0
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: https://github.com/PennyLaneAI/pennylane
Author: 
Author-email: 
License: Apache License 2.0
Location: /home/abhishekabhishek/anaconda3/envs/raccoon/lib/python3.11/site-packages
Requires: appdirs, autograd, autoray, cachetools, networkx, numpy, pennylane-lightning, requests, rustworkx, scipy, semantic-version, toml, typing-extensions
Required-by: PennyLane-qiskit, PennyLane_Lightning

Platform info:           Linux-6.5.0-18-generic-x86_64-with-glibc2.35
Python version:          3.11.9
Numpy version:           1.26.4
Scipy version:           1.13.0
Installed devices:
- lightning.qubit (PennyLane_Lightning-0.36.0)
- default.clifford (PennyLane-0.36.0)
- default.gaussian (PennyLane-0.36.0)
- default.mixed (PennyLane-0.36.0)
- default.qubit (PennyLane-0.36.0)
- default.qubit.autograd (PennyLane-0.36.0)
- default.qubit.jax (PennyLane-0.36.0)
- default.qubit.legacy (PennyLane-0.36.0)
- default.qubit.tf (PennyLane-0.36.0)
- default.qubit.torch (PennyLane-0.36.0)
- default.qutrit (PennyLane-0.36.0)
- default.qutrit.mixed (PennyLane-0.36.0)
- null.qubit (PennyLane-0.36.0)
- qiskit.aer (PennyLane-qiskit-0.35.1)
- qiskit.basicaer (PennyLane-qiskit-0.35.1)
- qiskit.ibmq (PennyLane-qiskit-0.35.1)
- qiskit.ibmq.circuit_runner (PennyLane-qiskit-0.35.1)
- qiskit.ibmq.sampler (PennyLane-qiskit-0.35.1)
- qiskit.remote (PennyLane-qiskit-0.35.1)

Existing GitHub issues

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

[CatalinaAlbornoz]

Hi @abhishekabhishek, thank you for opening this bug. What precision were you looking for? Is it more a matter of being aware of the precision beforehand? Or do you need a precision of 1e-9 or similar?

[abhishekabhishek]

Hi @CatalinaAlbornoz, I was expecting a precision of at least 1e-8 since that's what np.allclose uses for matrix equivalence checking by default. This is what led to some tests failing for my use case above. However, since I don't have a strict requirement on the precision, 1e-7 suffices but I think it is worth mentioning in the documentation in case other users do have a strict requirement.

[CatalinaAlbornoz]

Ah good catch @abhishekabhishek ! Yes I agree that mentioning it in the docs is important. Thanks for the feedback!

[dwierichs]

The problem seems to lie in the optimization we are using to reduce the number of CNOTs. At least the code above works out if one deactivates this optimization and always decomposes to the most general circuits containing 3 CNOTs...

[glassnotes]

I took a closer look at this (as the original author, I feel responsible for the issues here! 😅). I tried with/without the small perturbation in precision that was added in #5448, but this was not the problem, and there is precision loss in both cases.

I think what it boils down to is, the matrix used here is a special case (controlled-phase gate) whose decomposition requires only 2 CNOTs. You can check that the following circuit identity holds, up to a global phase:

Applying some simple circuit identities on the single-qubit rotations, we can see that the two_qubit_decomposition function produces a very convoluted version of this with many extra rotations, some of which are 0. The large amount of extra angles, each with some small loss of precision, contributes to the overall loss of precision.

However, I'm a bit stumped about why it works with the 3 CNOT case which has just as many angles. It could be a matter of how the rotations are derived from the matrix products computed internally by functions like _extract_su2su2_prefactors; in the 2-CNOT case for this controlled phase, many of the matrix entries are very close to 0.

[dwierichs]

Indeed, looking at all these extra rotations and setting the angles that are very close to $\pi$ or $3\pi$ to np.pi fixes the precision loss (to standard numpy precision). In particular, the middle angle of the single-qubit rotations into ZYZ seems to cause precision loss 🤔

edit: I was able to solve the issue using arccos of the diagonal entry rather than arcsin of the off-diagonal (in _zyz_get_rotation_angles). However, I'd expect that both will suffer from precision loss at extremal points, so this is not really a general fix.
Maybe we can come up with a generalized routine that smartly chooses how to extract this type of angle information from multiple (theoretically) redundant matrix entries, picking the numerically favourable one.

[glassnotes]

Perhaps a simple check like

if qml.math.isclose(np.sin(theta), 0):
    theta = np.pi

in those angle extraction functions would be useful for "resetting" things that are close-to-multiples-of-pi to pi? It could be a "good first issue". (Though TBH I could have sworn I have implemented something like that in the past, but maybe that was in somewhere in The Ionizer)

[dwierichs]

The above will no longer be differentiable at these special points (an issue I've been fighting with elsewhere (#5715) as well 😅 )

[glassnotes]

Right, of course... 🤦‍♀️

[co9olguy]

these kind of numerical-value-based checks/"fixes" can also sometimes lead to problems with jitting