Get probabilites of fock outputs from averages of samples

[snicolau]

Issue description

This issue is a follow-up from a discussion in the #strawberryfields chat in Xanadu's Slack channel with @josh146 and @mariaschuld.

Description of the issue.

• Expected behavior: My main goal would be to get the probabilites of the different outputs in the fock basis, which are the ones obtained via the qml.probs(wires), taking into account the number of shots used in the device. Then, I want to be able to use this probabilities in a cost function to update different parameters of the quantum circuit I am using via TensorFlow. The main goal of this is would be to simulate the behaviour of a real device where you have a limited amount of shots.

• Actual behavior: The parameter shots in the device does not affect the result when using the qml.probs(wires) operation, the probabilities obtained are always theoretical and not the result of averaging the samples of x shots.
  Also, I don't think there is a real way of getting samples in the fock basis from the device, as you would do in StrawberryFields with MeasureFock() | q. The most similar thing is using qml.NumberOperator() operation. When using this operation the result you get changes with the number of shots, but it looks like that it's not from averaging more samples but just having lower noise when increasing number of shots. Using this circuit:

dev = qml.device("strawberryfields.fock", wires=2, cutoff_dim=2, shots=10000)
@qml.qnode(dev, diff_method="finite-diff", interface="tf")
def circuit(weights):
    qml.FockState(1, wires=0)
    qml.Beamsplitter(weights[2], weights[3], wires=[0, 1])
    return [qml.expval(qml.NumberOperator(i)) for i in range(2)]

This is what you get for 1, 10, 1000, 10000 shots:

tf.Tensor([0.39695311 0.51878537], shape=(2,), dtype=float64)
tf.Tensor([0.26621335 0.70714082], shape=(2,), dtype=float64)
tf.Tensor([0.21179583 0.78553958], shape=(2,), dtype=float64)
tf.Tensor([0.20766148 0.79149591], shape=(2,), dtype=float64)

For 1 shot you should get either [1,0] or [0,1], but you get these noisy results that don't add to 1, and which get more accurate when increasing shot number. Moreover, for a 2 wire and 1 photon circuit, the qml.NumberOperator can be used to get the probabilities of the outputs, but for 4 wires and 2 photons, getting the average number of photons in each mode won't allow you to get the different probabilities of each output, so this operation won't be useful for my purpose.

• Reproduces how often: Always

System information

• PennyLane-SF version:
  0.15.0

• Python version:
  3.8.5

• NumPy and SciPy versions:
  NumPy 1.19.5, SciPy1.6.1

• Installation method:
  Via pip

Additional information

The circuit I am working on has 4 modes/wires and I am working with dual-rail encoded qubits. So, the outputs I'd be interested to get the probabilities of are [1,0,1,0], [1,0,0,1], [0,1,1,0] and [0,1,0,1].

[josh146]

Thanks for documenting this @snicolau! Behind the scenes, the reason for this is that the PennyLane-SF plugin is currently always computing the analytic expectation value, and then using the Central limit theorem to sample from a normal distribution with the same mean, and var = np.sqrt(mean/shots).

We should modify this logic to match other plugins. For finite shots, the plugin should instead compute samples directly, and then average these samples to determine the mean.

[josh146]

Another comment: this behaviour also leads to a discrepancy between hardware devices (which do currently average samples) and simulators. The above QNode will in fact return discrete samples if executed on the X8 chip.

[snicolau]

Hello! I just wanted to know if there are any new updates or comments regarding this issue, and if it is something which it is expected to be solved some time in the future @josh146

[josh146]

Hi @snicolau! Unfortunately I don't have any updates to report on at the moment -- in order to update the behavior here, this requires a significant refactor of both the plugin and, potentially, Strawberry Fields. We would very much like to resolve it, but I cannot provide a firm date at the moment.

Is this a feature that is blocking your work/research at the moment?