# How linear combination of unitaries gradient work (Qiskit, PennyLane)?

I'm trying to implement linear combination of unitaries(LCU) gradient from Qiskit Gradient Framework but on PennyLane.

First, i looked through the source code in Qiskit. In Qiskit LCU gradient if we have 2 parameterized gate and 2 hamiltonian terms it will look like this

for each parameterized gate they use control unitary gate that unitary gate is same as parameterized gate (if parameterized gate is RZ it will be CZ) and calculate expectation value. So i tried to implement this on PennyLane like this

in this case i have 3 parameterized gates for each qubit. Expectation value of each circuit is gradient of each parameterized gate.

The code looks correct but the result when i use this gradient to do optimization is really bad. So i want to ask if there anything wrong in my idea of how LCU gradient work?

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.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):
# Prepare auxiliary qubit:
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

# 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...