I reconstructed a circuit like yours:
dev = qml.device("lightning.qubit", wires=5)
target_wires = dev.wires[:4]
H = qml.Hamiltonian([0.5]*4, [qml.PauliY(w) for w in target_wires])
@qml.qnode(dev)
def circuit(x, y, z):
[qml.Hadamard(w) for w in target_wires]
[qml.RX(_x, w) for _x, w in zip(x, target_wires)]
[qml.RY(_y, w) for _y, w in zip(y, target_wires)]
[qml.RZ(_z, w) for _z, w in zip(z, target_wires)]
return qml.expval(H)
x, y, z = pnp.random.random((3, 4), requires_grad=True)
Here I just assumed whatever terms in the Hamiltonian, these details should not matter. I believe that the bit that is missing from your circuits is the changed measurement: In addition to the auxiliary qubit being entangled with the target qubit(s), one needs to actually execute the "linear combination" in the LCU, by measuring PauliZ
on the auxiliary qubit in addition to measuring the original Hamiltonian. That means we want to measure their tensor product.
I manually coded up a function for this specific circuit that spits out the derivative circuit for a given parameter to be differentiated:
H_LCU = qml.PauliZ(4) @ H
@qml.qnode(dev)
def LCU_circuit(x, y, z, diff_idx):
[qml.Hadamard(w) for w in target_wires]
# Prepare auxiliary qubit:
qml.Hadamard(4)
qml.adjoint(qml.S)(4)
if diff_idx in range(4):
qml.CNOT([4, diff_idx])
[qml.RX(_x, w) for _x, w in zip(x, target_wires)]
if diff_idx in range(4, 8):
qml.CY([4, diff_idx%4])
[qml.RY(_y, w) for _y, w in zip(y, target_wires)]
if diff_idx in range(8, 12):
qml.CZ([4, diff_idx%4])
[qml.RZ(_z, w) for _z, w in zip(z, target_wires)]
#postprocess auxiliary qubit
qml.Hadamard(4)
# Measure altered observable
return qml.expval(H_LCU)
Here, H_LCU
is the mentioned tensor product of PauliZ
on the auxiliary qubit with the original Hamiltonian. I counted the parameters, as indexed by diff_idx
in the order "first all x, then all y, then all z parameters", which of course is a somewhat arbitrary choice again.
However, this works :)
>>> print(qml.jacobian(circuit)(x, y, z)) # Autodiff derivative
(array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.11022302e-16]),
array([-0.00346884, -0.02131515, -0.11414253, -0.06634378]),
array([0.44512374, 0.38051277, 0.38258548, 0.38570437]))
>>> print(tuple(np.array([LCU_circuit(x, y, z, i) for i in range(j*4, (j+1)*4)]) for j in range(3))) # LCU circuit derivative for all parameters
(array([0., 0., 0., 0.]),
array([-0.00346884, -0.02131515, -0.11414253, -0.06634378]),
array([0.44512374, 0.38051277, 0.38258548, 0.38570437]))
The created circuit, say for diff_idx=6
, looks like yours, up to the changed measurement:
0: ──H──RX(1.00)──RY(0.73)──RZ(0.59)───────────┤ ╭<𝓗>
1: ──H──RX(0.21)──RY(0.57)──RZ(0.24)───────────┤ ├<𝓗>
2: ──H──RX(0.59)─╭Y─────────RY(0.26)──RZ(0.63)─┤ ├<𝓗>
3: ──H──RX(0.37)─│──────────RY(0.68)──RZ(0.34)─┤ ├<𝓗>
4: ──H──S†───────╰●─────────H──────────────────┤ ╰<𝓗>
Now, can we take this further and automate the generation of these circuits?
Definitely, this technique basically is used also when you compute the metric_tensor
in PennyLane via the "Hadamard test" method (another name for the LCU method, because this trick using an auxiliary qubit is a generalized Hadamard test). Unfortunately, currently the tape transform needed for this is not implemented in PennyLane, maybe it is something you want to consider contributing...
N.B. The derivative with respect to all x
parameters is zero, because they parametrize RX
rotations acting on a |+>
state, which does not do anything...