`ApproxTimeEvo` decomposition bug

Dear PennyLane team,
Hope all is well.
I’m having problems with using ApproxTimeEvolution.
I would like to simulate a simple Hamiltonian given by HH molecule using ApproxTimeEvolution operator. However, the resulting circuit has a number of issues:

  1. It seems the circuit contains both the operator ApproxTimeEvolution and its decomposition. Having just a container ApproxTimeEvolution doesn’t make any sense so I attempt to decompose it. However, this produces a weird result (see the figure).
  2. When I export the circuit to QASM, the PauliRotation gates are decomposed into CNOTS, and single-qubit rotation. This is undesired behavior as I did not specify that the circuit should be decomposed into single- and two-qubit gates.

Here is the code to reproduce the issue:

symbols = ["H", "H"]
coordinates = np.array([0, 0, 0, 0, 0, 0.74137727])
hamiltonian, num_qubits = qchem.molecular_hamiltonian(symbols, coordinates, method="pyscf")
print(f"Number of qubits: {num_qubits}")
print("Qubit Hamiltonian")
print(hamiltonian)

wires = range(num_qubits)
num_wires = num_qubits

dev = qml.device('default.qubit', wires=wires)

evolution_time=0.5
trotter_steps=1
@qml.qnode(dev)
def circuit(evolution_time):
    ApproxTimeEvolution(hamiltonian, evolution_time, trotter_steps).decompose()
    return [qml.expval(qml.PauliZ(wires=i)) for i in wires]

specs_func = qml.specs(circuit, expansion_strategy='device')
print(specs_func(evolution_time))
fig, ax = qml.draw_mpl(circuit)(evolution_time)

Output of spec_func:

{'gate_sizes': defaultdict(<class 'int'>, {1: 8, 2: 12, 4: 8}), 'gate_types': defaultdict(<class 'int'>, {'PauliRot': 28}), 'num_operations': 28, 'num_observables': 4, 'num_diagonalizing_gates': 0, 'num_used_wires': 4, 'depth': 19, 'num_trainable_params': 28, 'num_device_wires': 4, 'device_name': 'default.qubit.autograd', 'expansion_strategy': 'gradient', 'gradient_options': {}, 'interface': 'autograd', 'diff_method': 'best', 'gradient_fn': 'backprop'} 

Note that there is a wrong number of Pauli rotation gates (twice of what is expected). There should be 14 Pauli rotation gates because there are 14 non-tirival Hamiltonian terms.

The drawing of the circuit shows that the circuit contains both the decomposition and the original ApproxTimeEvo operator which is also incorrect.

Name: PennyLane
Version: 0.28.0
Summary: PennyLane is a Python quantum machine learning library by Xanadu Inc.
Home-page: https://github.com/XanaduAI/pennylane
Author: 
Author-email: 
License: Apache License 2.0
Requires: appdirs, autograd, autoray, cachetools, networkx, numpy, pennylane-lightning, requests, retworkx, scipy, semantic-version, toml
Required-by: PennyLane-Lightning

Platform info:           Linux-5.19.0-46-generic-x86_64-with-glibc2.35
Python version:          3.10.6
Numpy version:           1.23.5
Scipy version:           1.10.0

Thank you!

Hey @Einar_Gabbassov! Welcome back :smiley:

You’re a few versions behind. We’re currently on v0.31 :slight_smile:. That said, I don’t think your issue has anything to do with the v0.28 of pennylane, it all has to do with how operators / gates are “placed” in a circuit under the hood.

The long and short of it is that anytime an operator/gate is instantiated inside of a QNode, it gets placed in the circuit. So, when you do ApproxTimeEvolution(hamiltonian, evolution_time, trotter_steps).decompose(), two things are happening:

  1. An instance of ApproxTimeEvolution is instantiated
  2. All of the operators in its decomposition are instantiated as well

That’s why both ApproxTimeEvolution and the elements of its decomposition show up in your circuit :slight_smile:.

What you want to do is the following:

coeffs_and_time = [*hamiltonian.coeffs, evolution_time]

@qml.qnode(dev)
def circuit(evolution_time):
    qml.ApproxTimeEvolution.compute_decomposition(*coeffs_and_time, wires=wires, n=trotter_steps, hamiltonian=hamiltonian)
    return [qml.expval(qml.PauliZ(wires=i)) for i in wires]

specs_func = qml.specs(circuit, expansion_strategy='device')
print(specs_func(evolution_time))
fig, ax = qml.draw_mpl(circuit)(evolution_time)

Screenshot 2023-07-12 at 10.55.22 AM

This way, only its decomposition elements are instantiated, since compute_decomposition is a static method.

1 Like

Hi @isaacdevlugt,
Thanks for the reply! It’s been a while since I used PennyLane, but it’s always good to use domestic products :slight_smile:. I might update to a newer version at some point.

I have a minor comment. Why does compute_decomposition() takes both coeffs_and_time and hamiltonian? Seems a bit redundant as both variables contain coefficients.

Aside that, all seems to work now. However, my second issue is still pending. Getting the right circuit is only half of the job. Now, I need to export the PennyLane circuit as is into QASM.
Right now, the circuit is broken down into 1- and 2-qubit gates.

qasm_string = circuit.qtape.to_openqasm()
print(qasm_string)

The output is:

OPENQASM 2.0;
include "qelib1.inc";
qreg q[4];
creg c[4];
rz(0.23702219030799287) q[0];
rz(0.23702219030799287) q[1];
rz(-0.4605771375540184) q[2];
rz(-0.4605771375540184) q[3];
cx q[1],q[0];
rz(0.18453162325559797) q[0];
cx q[1],q[0];
rx(1.5707963267948966) q[0];
h q[1];
h q[2];
...

The command qasm_string = circuit.qtape.to_openqasm() was kindly recommended to me by someone from PennyLane team. Otherwise, I would never figure out how to export to QASM :smiley:

We haven’t really updated the openqasm method in quite some time, and its still openqasm2, not 3. So it has a fairly limited gate set. It’s on the list of things we want to improve on :grin:

FYI: I made a PR to make the relevant docstrings more clear