Hi,

Python functions that define circuits which return sampling are not differentiable. This is quite understandable. But doesn’t `qml.expVal()`

use sampling to estimate the energy? Why circuits that return `qml.expVal()`

can be differentiated?

# Differentiable circuit with sampling

Hey @Einar_Gabbassov!

But doesn’t

`qml.expVal()`

use sampling to estimate the energy?

Expectation values can be calculated either using sampling or analytically with access to the underlying state (only compatible with simulators). Analytic mode is specified by setting `shots=None`

when loading your PennyLane device and is the default.

Why circuits that return

`qml.expVal()`

can be differentiated?

While it doesn’t make sense to differentiate a given sample from a random process (e.g. a quantum circuit), it is ok to differentiate a statistic defined on the random process (such as the expectation value).

The expectation value is a function \mu = f(\boldsymbol{\theta}) that maps from the circuit parameters \boldsymbol{\theta} to the corresponding value. With an exact simulator, we know \mu exactly while on an approximate simulator or quantum hardware we find an estimate \hat{\mu}.

The gradient of the expectation value g(\boldsymbol{\theta}) := \nabla_{\boldsymbol{\theta}} f(\boldsymbol{\theta}) is also well defined. Again, we can find g(\boldsymbol{\theta}) exactly with an exact simulator or estimate it when we have a limited number of shots. Check out Sec. 4 of this paper to get a deeper understanding.

In other words, having a finite number of shots just means that the gradient of the expectation value is approximate.

hi @Tom_Bromley

Thank you for the detailed reply!

This brings me to the question which I wanted to ask but decided to differ until I get the answer for the first question.

In practice, one might want to use bias estimates. As far as I know `expVal()`

gives unbiased estimate. This is not flexible.

Ideally, a circuit like this would be very useful

```
ansatz(params)
samples = qml.sample(Operator)
return compute_biased_average(samples)
```

The above could be differentiated. However, as far as I know we can’t do this for two reasons:

- Pennylane gives warning that circuits which use sampling can’t be differentiated.
- According to the Docs, circuit functions must not contain arbitrary code, except multiplication operations.

Unless, I miss some technicality (which is very likely), the above puzzles me a lot. We can use samples to compute expVal and do differentiation, but we can’t use samples to compute other estimates such that the circuit is still differentiable.

I hope I’m wrong, otherwise it makes my day so much harder

Hi @Einar_Gabbassov!

I believe what you want to do is possible, using PennyLane’s concept of ‘shot-batching’.

For example, consider the following QNode:

```
import pennylane as qml
from pennylane import numpy as np
dev = qml.device("default.qubit", wires=2, shots=[(1, 100)])
@qml.qnode(dev)
def circuit(params):
for i in range(2):
qml.RX(params[i, 0], wires=i)
qml.RY(params[i, 1], wires=i)
qml.CNOT(wires=[0, 1])
for i in range(2):
qml.RX(params[2 + i, 0], wires=i)
qml.RY(params[2 + i, 1], wires=i)
qml.CNOT(wires=[0, 1])
return qml.expval(qml.PauliZ(0))
```

We have specified that the device has shots `[(1, 100)]`

, which corresponds to `[(size_of_batch, number of batches)]`

. This means that the QNode will be executed for 100 shots on the remote hardware/simulator, however, *these 100 shots will be post-processed in batches of 1*. In essence, we are computing single-shot expectation values:

```
>>> params = np.ones([4, 2], dtype=np.float64)
>>> circuit(params)
[ 1. 1. 1. -1. 1. -1. -1. -1. -1. 1. 1. -1. -1. -1. -1. 1. 1. -1.
1. 1. -1. 1. 1. 1. -1. -1. 1. -1. -1. 1. -1. -1. -1. -1. -1. -1.
-1. -1. -1. -1. -1. 1. -1. 1. 1. 1. -1. -1. -1. -1. 1. -1. 1. -1.
-1. -1. -1. -1. -1. -1. -1. 1. -1. 1. -1. 1. -1. -1. 1. -1. -1. -1.
1. 1. -1. 1. -1. 1. -1. 1. -1. 1. 1. -1. -1. 1. -1. -1. -1. -1.
1. 1. -1. 1. -1. 1. -1. 1. 1. 1.]
```

Note that we are not using `sample`

, we are using `expval`

, so this is differentiable!

Let’s create a cost function that processes this to compute the population variance:

```
def cost(params):
samples = circuit(params)
return np.var(samples, ddof=1)
```

This can be differentiated:

```
>>> cost(params)
0.9418181818181818
>>> qml.grad(cost)(params)
[[-0.03636364 0.14343434]
[ 0.13131313 0.04848485]
[-0.02626263 -0.05252525]
[ 0.07474747 -0.03030303]]
```

More details on shot-batching is available here: https://pennylane.readthedocs.io/en/stable/introduction/circuits.html#shot-batches

Hi, thanks everyone for the nice discussion.

I am trying to do the same with the strawberry.fock backend but I get a single output, even if I specify shots = [(1,100)].

More specifically, I am using

`dev = qml.device("strawberryfields.fock", wires=1, cutoff_dim=10, shots=[(1,100)])`

and the same photonic quantum neural network as in the tutorial (with `expval(qml.X(0))`

Does anyone have an idea why i do not get 100 output like in the answer from Josh? Thanks a lot

Hi @Oriel, unfortunately I believe that this is because the `strawberryfields.fock`

device has not yet been updated to support shot batching.

Thanks for catching this! I will create a GitHub issue to alert the team.

Here is the issue:

Hi,

I’ve been trying the proposed approach but it messes up outputs from the `qml.ExpvalCost()`

.

For example, I have a device with shots = 100.

I call `qml.ExpvalCost()`

, this returns a single value, ok!

Next, I modify my device so that it does batch shots, `dev.shots=[(1,100)]`

and I run my custom cost function which is now differentiable because I use batches with `qml.expval()`

. Then, I try to set the device back to the default behaviour `dev.shots=100`

.

But this breaks `qml.ExpvalCost()`

as it does not output a single value but an average energy for each Hamiltonian term.

I would post an error, but it is pretty hard to pin point a meaningful error trace.

How I could temporarily set batch shots and then revert back to default single batch behaviour without `qml.ExpvalCost()`

doing funny stuff?

ah, I just realized that I have a single objective function which consists of two functions. One function tries to make use of the default device and the other function is trying to use the device with batch shots.

Trying manipulate the device from inside the functions doesn’t work.

Hi,

On the same note, i would like to efficiently sample from circuits with shifted parameters, where the shifts correspond to the shift in the shift rule. I mean, if I have a circuit with a single qubit rotation with angle theta, I would like to sample from the circuit with angle theta + pi/2 and theta - pi/2.

I tried to do it myself but I am lacking efficiency since I basically do a for loop over all parameters. Does anyone know if this is already implemented or if there is a more efficient way to do so?

Thanks a lot

Hey @Oriel, there is a manual implementation of the shift rule available in the backpropagation tutorial:

```
def parameter_shift_term(qnode, params, i):
shifted = params.copy()
shifted[i] += np.pi/2
forward = qnode(shifted) # forward evaluation
shifted[i] -= np.pi
backward = qnode(shifted) # backward evaluation
return 0.5 * (forward - backward)
def parameter_shift(qnode, params):
gradients = np.zeros([len(params)])
for i in range(len(params)):
gradients[i] = parameter_shift_term(qnode, params, i)
return gradients
```

This assumes that the QNode only takes a single parameter, and this single parameter is a 1-dimensional NumPy array

hi @josh

Sorry for taking so long, been busy.

I attached the snippet of the code at the end of this post. Basically, the code has a Hamiltonian with 4 terms and the device is configured to use batch shots i.e, `shots=[1, num_batches]`

. Strangely, the output from the `expValCost`

function are 4 values. This doesn’t make sense, because based on the default mechanics, `expValCost`

should output expectations of energy of the **entire** Hamiltonian and not expectations of each Hamiltonian term.

I was able to rectify this by writing my own `expValCost`

which works with batch shots. However, the fundamental problem arises when trying to do gradient descent optimization.

From what I understand, suppose we use `MyExpValCost()`

that produces Hamiltonian energy samples (each energy sample is evaluated from a single shot but we have many batches just as you recommended).

This would yield an array of Hamiltonian energies:

```
energy_samples = MyExpValCost(params)
print("energy samples:", energy_samples)
```

The output for 10 batches is 10 expected energies of the Hamiltonian H. All 10 expectations are estimated from 1 shot.

That is, we have,

\{<H>_1, <H>_2, ..., <H>_{10}\}.

Next, I create a cost function like so

cost = \sum_{i=1}^{10} <H>_i.

As you claimed this is differentiable.

The fundamental problem with this is that it takes orders of magnitude more time to evaluate. Let’s say I do 100000 batches (each batch is one shot). Then the autograd should differentiate a sum which contains 100000 differentiable terms.

Just to conclude, there are two problems:

- ExpValCost has strange output when batch shots are used.
- Even if 1. is not the case, the differentiation with batch shots is extremely slow compared to just differentiation of a cost function which uses
`expValCost`

with NO batch shots.

I still need a differentiable circuit which can do sampling then process sampling for producing some biased estimates.

```
ansatz(params)
samples = qml.sample(Operator)
return compute_biased_average(samples)
```

Here is the code which shows that `expValCost`

produces 4 terms instead of Hamiltonian energies:

```
objective_hamiltonian = qml.Hamiltonian([1,1,1,1], [qml.PauliZ(0), qml.PauliZ(1), qml.PauliZ(2), qml.PauliZ(3)])
def ansatz(params, wires, depth, h_objective_coeffs, **kwargs):
params = params.reshape((depth, 2))
# Initial state
qml.BasisState(np.array([0]*num_qubits, requires_grad=False), wires=wires)
# Add Hadamard gates
for i in wires:
qml.Hadamard(wires=i)
num_qubits = 4
num_shots = 10000
# Use batch shots for differentiability
dev = qml.device('default.qubit', wires=num_qubits, shots=[1,num_shots])
depth = 1
penalty = 2
ansatz_settings = {
"h_objective_coeffs": objective_hamiltonian.coeffs,
"wires": dev.wires,
"depth": depth
}
ansatz = partial(ansatz, **ansatz_settings)
objective_cost_fn = qml.ExpvalCost(ansatz, objective_hamiltonian, dev, optimize=True)
params = np.zeros((1,2))
objective_cost_fn(params)
```