[BUG] Unitary Matrix Mismatch in OpenQASM cu Gate Implementation

[MattePalte]

Expected behavior

The unitary matrix generated by the cu gate in the OpenQASM code should match the expected unitary matrix.

Actual behavior

The unitary matrix generated by the cu gate in the OpenQASM code is significantly different from the expected one, leading to an incorrect implementation.

Additional information

This issue can be reproduced consistently. The discrepancy between the expected and actual unitary matrices suggests a potential bug in the conversion process.

A circuit with a cu gate, with two parameters and two qubits, is converted to OpenQASM via qml.make_qscript and qs.to_openqasm. You can run this code to verify that they implement different unitaries:

Source code

qasm_str = """
include "qelib1.inc";
qreg q[2];
cu(2.13431562633248,5.350646294405996,2.4859039172485513,4.866465765080696) q[0],q[1];
"""

import pennylane as qml
from pennylane.tape import make_qscript, QuantumScript

qc = qml.from_qasm(qasm_str)
print(qml.draw(qc)())

# 0: ─╭U(M0)─┤
# 1: ─╰U(M0)─┤
# M0 =
# [[ 1.        +0.j          0.        +0.j          0.        +0.j
#    0.        +0.j        ]
#  [ 0.        +0.j          0.07406597-0.47689822j  0.        +0.j
#   -0.42113483-0.76793736j]
#  [ 0.        +0.j          0.        +0.j          1.        +0.j
#    0.        +0.j        ]
#  [ 0.        +0.j         -0.61499667-0.62358754j  0.        +0.j
#    0.47811677+0.06574213j]]

qs = make_qscript(qc)()
qasm_str_pennylane = qs.to_openqasm(measure_all=False)
print(qasm_str_pennylane)

# OUTPUT:
# OPENQASM 2.0;
# include "qelib1.inc";
# qreg q[2];
# creg c[2];
# rz(11.05676742575686) q[0];
# ry(1.7679655525118294) q[0];
# rz(1.1945736265076667) q[0];
# rz(11.648070962090944) q[1];
# ry(1.5707963267948966) q[1];
# rz(9.42477796076938) q[1];
# cx q[1],q[0];
# rz(0.0) q[0];
# rx(1.219386317488504) q[1];
# cx q[1],q[0];
# rz(2.7952601622872653) q[0];
# ry(2.7485135499089113) q[0];
# rz(12.05550548208478) q[0];
# rz(4.440892098500626e-16) q[1];
# ry(1.5707963267948974) q[1];
# rz(3.4198552159966096) q[1];

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from qiskit.quantum_info import Operator
from qiskit.qasm2 import load, LEGACY_CUSTOM_INSTRUCTIONS

# store to files: qasm_start, qasm_pennylane
with open('qasm_start.qasm', 'w') as f:
    f.write(qasm_str)
with open('qasm_pennylane.qasm', 'w') as f:
    f.write(qasm_str_pennylane)

qc_start = load('qasm_start.qasm',
                custom_instructions=LEGACY_CUSTOM_INSTRUCTIONS)
qc_pennylane = load('qasm_pennylane.qasm',
                    custom_instructions=LEGACY_CUSTOM_INSTRUCTIONS)

op_start = Operator(qc_start)
op_pennylane = Operator(qc_pennylane)

print("Start")
print(op_start.data)

print("Pennylane Exported")
print(op_pennylane.data)

unitary_equiv = op_start.equiv(op_pennylane)
print("Are the unitaries equivalent?", unitary_equiv)

# OUTPUT:
# Start
# [[ 1.        +0.j          0.        +0.j          0.        +0.j
#    0.        +0.j        ]
#  [ 0.        +0.j          0.07406597-0.47689822j  0.        +0.j
#   -0.42113483-0.76793736j]
#  [ 0.        +0.j          0.        +0.j          1.        +0.j
#    0.        +0.j        ]
#  [ 0.        +0.j         -0.61499667-0.62358754j  0.        +0.j
#    0.47811677+0.06574213j]]
# Pennylane Exported
# [[ 9.49229584e-01+3.14584164e-01j  2.35657764e-16-6.64055839e-16j
#    1.11035105e-16-4.07918502e-16j -1.87032964e-17-3.33697653e-16j]
#  [ 1.53206331e-16-7.04242249e-16j  9.49229584e-01+3.14584164e-01j
#   -9.67733676e-17-4.90923424e-16j  1.31663013e-16-6.04771645e-17j]
#  [-5.42048940e-16-1.31959181e-16j -1.60034140e-16-2.41188979e-17j
#    2.20330240e-01-4.29385913e-01j -1.58172704e-01-8.61431204e-01j]
#  [-1.77434893e-16+3.43976794e-18j -3.84823169e-16-5.38891098e-17j
#   -3.87602263e-01-7.85395957e-01j  4.33161154e-01+2.12812338e-01j]]
# Are the unitaries equivalent? False

Tracebacks

No error tracebacks are generated, but the output matrices indicate a discrepancy.

System information

Name: PennyLane
Version: 0.38.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: .../lib/python3.10/site-packages
Requires: appdirs, autograd, autoray, cachetools, networkx, numpy, packaging, pennylane-lightning, requests, rustworkx, scipy, toml, typing-extensions
Required-by: PennyLane-qiskit, PennyLane_Lightning

Platform info:           Linux-5.15.0-130-generic-x86_64-with-glibc2.35
Python version:          3.10.14
Numpy version:           1.26.4
Scipy version:           1.14.1
Installed devices:
- qiskit.aer (PennyLane-qiskit-0.38.1)
- qiskit.basicaer (PennyLane-qiskit-0.38.1)
- qiskit.basicsim (PennyLane-qiskit-0.38.1)
- qiskit.remote (PennyLane-qiskit-0.38.1)
- default.clifford (PennyLane-0.38.0)
- default.gaussian (PennyLane-0.38.0)
- default.mixed (PennyLane-0.38.0)
- default.qubit (PennyLane-0.38.0)
- default.qubit.autograd (PennyLane-0.38.0)
- default.qubit.jax (PennyLane-0.38.0)
- default.qubit.legacy (PennyLane-0.38.0)
- default.qubit.tf (PennyLane-0.38.0)
- default.qubit.torch (PennyLane-0.38.0)
- default.qutrit (PennyLane-0.38.0)
- default.qutrit.mixed (PennyLane-0.38.0)
- default.tensor (PennyLane-0.38.0)
- null.qubit (PennyLane-0.38.0)
- lightning.qubit (PennyLane_Lightning-0.38.0)

Existing GitHub issues

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

[CatalinaAlbornoz]

Thanks for opening this bug report @MattePalte. I can replicate with PennyLane v0.40.

I'll bring this to the team and we'll come back to you with more info.

[CatalinaAlbornoz]

Hi @MattePalte,

My colleague Isaac looked into this and it looks like the problem might be caused by a wire order or lexicographical order issue.

dev = qml.device("default.qubit")

@qml.qnode(dev)
def circuit(option=1):
    qml.Hadamard(0)
    qml.Hadamard(1)

    if option == 1:
        qml.QubitUnitary(op_start.data, wires=range(2))
    elif option == 2:
        qml.QubitUnitary(op_pennylane.data, wires=range(2))

    return qml.probs(wires=range(2))

print(circuit(option=1)) # [0.25      0.4175181 0.25      0.0824819]
print(circuit(option=2)) # [0.25      0.25      0.4175181 0.0824819]

Unfortunately I don't know how long it will take the team to fix this. Is this blocking work? Or are you able to look for a workaround?