Energy distribution from qml.expval not matching qml.sample

[JohnKerf]

Hello,

• I have a Hermitian Hamiltonian matrix that I am using within a VQE to find the groundstate energy.
• I have been using qml.expval with a finite number of shots to calculate the expectation value
• I have noticed that at a small number of shots this does not seem to sample energies as I would expect and seems to not sample lower than the groundstate energy?

Below is a simple example using only 2 shots where we would expect a large variance in the energy distribution. We take a generic circuit (RY on all qubits) and calculate the expectation energy using a random set of parameters each time. I have then plotted these energies to see the distribution of energies returned.

import pennylane as qml
import numpy as np

num_qubits = 3
shots = 2
num_energy_samples = 2500

dev = qml.device("default.qubit", wires=num_qubits, shots=shots, seed=42)

# Hamiltonian matrix
H = np.array([[2.375     , 2.47487373, 1.76776695, 0.8660254 , 0.        ,
        0.        , 0.        , 0.        ],
       [2.47487373, 5.875     , 5.        , 3.06186218, 0.        ,
        0.        , 0.        , 0.        ],
       [1.76776695, 5.        , 9.375     , 5.51135192, 0.        ,
        0.        , 0.        , 0.        ],
       [0.8660254 , 3.06186218, 5.51135192, 5.875     , 0.        ,
        0.        , 0.        , 0.        ],
       [0.        , 0.        , 0.        , 0.        , 1.375     ,
        1.06066017, 1.76776695, 0.8660254 ],
       [0.        , 0.        , 0.        , 0.        , 1.06066017,
        4.875     , 3.        , 3.06186218],
       [0.        , 0.        , 0.        , 0.        , 1.76776695,
        3.        , 8.375     , 3.06186218],
       [0.        , 0.        , 0.        , 0.        , 0.8660254 ,
        3.06186218, 3.06186218, 4.875     ]])

H_decomp = qml.pauli_decompose(H, wire_order=range(num_qubits))
paulis = H_decomp.ops
coeffs = H_decomp.coeffs

min_eigenvalue = np.sort(np.linalg.eig(H)[0])[0]

@qml.qnode(dev)
def circuit(params):
    
    param_index=0
    for i in range(num_qubits):
        qml.RY(params[param_index], wires=i)
        param_index += 1
        
    return qml.expval(qml.Hermitian(H, wires=range(num_qubits)))

#sample energies from the circuit with a random set of parameters each time
energies = []
for i in range(num_energy_samples):

    params = np.random.rand(num_qubits)*2*np.pi
    e = circuit(params)
    energies.append(e)

As you can see it does not seem to return a normal distribution and seems to be lower bounded at the groundstate energy eigenvalue and infact the majority of samples are close to this boundary.

@qml.qnode(dev)
def circuit(params):
    param_index = 0
    for i in range(num_qubits):
        qml.RY(params[param_index], wires=i)
        param_index += 1

    return [qml.sample(op) for op in paulis] #sample from each individual pauli string


energies = []
for i in range(num_energy_samples):
 
    params = np.random.rand(num_qubits) * 2 * np.pi
    samples = circuit(params)  
    samples = np.array(samples).T 
    expectation_val = np.mean(samples @ coeffs)

    energies.append(expectation_val)

(Since I am a new user it will not let me add the histogram for this)

If I now calculate the expectation value using direct sampling from each pauli string I now see what I would expect - a normal distribution with values extending lower than the groundstate due to the variance from such a small number of shots.

I am just wondering what qml.expval is doing differently for me to see this sort of distribution in energies and why does there appear to be a bunching close to the groundstate?

Thank you in advance.

Name: PennyLane Version: 0.40.0
Platform info: Windows-11-10.0.26100-SP0
Python version: 3.13.1
Numpy version: 2.0.2
Scipy version: 1.15.1
Installed devices:

• default.clifford (PennyLane-0.40.0)
• default.gaussian (PennyLane-0.40.0)
• default.mixed (PennyLane-0.40.0)
• default.qubit (PennyLane-0.40.0)
• default.qutrit (PennyLane-0.40.0)
• default.qutrit.mixed (PennyLane-0.40.0)
• default.tensor (PennyLane-0.40.0)
• null.qubit (PennyLane-0.40.0)
• reference.qubit (PennyLane-0.40.0)
• lightning.qubit (PennyLane_Lightning-0.40.0)
• qiskit.aer (PennyLane-qiskit-0.40.0)
• qiskit.basicaer (PennyLane-qiskit-0.40.0)
• qiskit.basicsim (PennyLane-qiskit-0.40.0)
• qiskit.remote (PennyLane-qiskit-0.40.0)

[CatalinaAlbornoz]

Hi @JohnKerf ,

We’re taking a look at your question and will get back to you with an answer in the next couple of days!

[CatalinaAlbornoz]

Hi @JohnKerf ,

We’ve had a high volume of questions lately but we should get back to you very soon on this

[daniela.murcillo]

Hi!
Thanks for your question.
I think comparing these two approaches does not offer a fair comparison.
Let’s take a step back and see what we expect for the expectation value of the energy of the system. The eigenvalues of that Hamiltonian are all positive. If you measure the energy many times and collect the statistics, ideally, you should obtain a positive expectation value since none of the eigenvalues are negative to begin with. This is what PennyLane is doing: it measures the energy n times (given by the number of shots) and averages those results.
I modified your code a bit to use qml.sample() with the Hamiltonian directly, like your second approach but not with the Pauli decomposition of the Hamiltonian. This modified code works as expected.

@qml.qnode(dev)
def circuit(params):
    param_index = 0
    for i in range(num_qubits):
        qml.RY(params[param_index], wires=i)
        param_index += 1

    #return [qml.sample(op) for op in paulis] #sample from each individual pauli string
    return qml.sample(qml.Hermitian(H, wires=range(num_qubits)))

energies = []
for i in range(num_energy_samples):
 
    params = np.random.rand(num_qubits) * 2 * np.pi
    samples = circuit(params)  
    samples = np.array(samples).T 
    #expectation_val = np.mean(samples @ coeffs)
    expectation_val = np.mean(samples)

    energies.append(expectation_val)

I believe that when you collect samples for the Paulis and then reconstruct the expectation value, there will be cases such that they give a negative value, but that is because of the randomness of the sampling involved.
I want to emphasize on the fact that the measurement of the energy doesn’t include any noise, it is a sharp measurement that can only retrieve one of the eigenvalues as outcome.

Let me know if this helped.

[JohnKerf]

Hi Daniela,

Thank you for your reply.

That now makes sense why it is bounded at the groundstate energy however can this method be used with a real quantum device? I was under the impression we need to decompose the Hamiltonian into pauli strings in order to measure the energy on a real device?

I am also curious why we have such a large grouping at the groundstate if we are sampling from the eigenstates? Surely if we are sampling 2 eigenstates and averaging it should not cluster near the groundstate like we see in the histogram for 2 shots?

I ask this because I find when I run my VQE at low shots (around 2 to 1024) it can return the exact groundstate energy. However, when I run with a higher number of shots, say 10000, I get worse results and the optimiser can not achieve the groundstate energy. Even with shots=None I can not achieve the groundstate energy. This seems counter-intuitive to what I would expect. I would expect that the more shots I use, the better the results, but this is not the case. I suspect it is because of this grouping at the groundstate we observe at lower shots?

[daniela.murcillo]

Hi!

I will check more regarding the first question, but my first impression is that it can be used on a real device but internally it will do the appropriate decompositions to be implemented on hardware.

For this code in particular, we are preparing a particular state. Even if the parameters of the Y rotations are random, a particular state is prepared and that’s why it shows such a distribution for the expectation value. For example, try to change to an X or Z rotation instead and see how the distribution changes. I think the biased is being introduced by us by choosing that particular circuit that gives a state with higher amplitudes for lower energy eigenstates.

Variational algorithms are tricky and not guaranteed to converge. Have you tried any other ansantz?Or another optimizer, it might be getting stuck at a local minimum or at a plateau.

Maybe this module from the Codebook about Variational quantum algorithms might be useful.

Let me know if this helped.

[JohnKerf]

Hi Daniela,

Thank you again for your reply.

My concern is that if I use qml.expval with finite shots and the default.qubit device then I am not correctly simulating what would be done if the device was switched to a real device?

I have tried changing the ansatz to use strongly entangling layers and I encounter the same problem.

@qml.qnode(dev)
def circuit(params):

    num_layers=1
    params_shape = (num_layers, num_qubits, 3)
    params = pnp.tensor(params.reshape(params_shape), requires_grad=True)
    qml.StronglyEntanglingLayers(weights=params, wires=range(num_qubits))
    
    return qml.expval(qml.Hermitian(H, wires=range(num_qubits)))

energies = []
for i in range(num_energy_samples):

    params = np.random.rand(3*num_qubits)*2*np.pi
    e = circuit(params)
    energies.append(e)


plt.figure()
sns.histplot(energies, kde=True)
plt.axvline(x=min_eigenvalue, color='red', linestyle='--', linewidth=1.0, label='groundstate eigenvalue')
plt.xlabel("Energy")
plt.title(f"Histogram of Sampled Energies (shots={shots})")
plt.legend()
plt.show()

As you can see, at 10000 shots it does not sample near the groundstate. However, at 2 shots we see a bias towards sampling near the groundstate. (I cant add the histogram for 2 shots as I am a new user but it is a similar distribution to before)

If we are sampling from the exact eigenvalue spectrum from the Hamiltonian, why is it at 10000 shots it does not return the groundstate but at 2 shots it is more likely to return values at or close to the groundstate?

If the ansatz used can’t return the groundstate eigenvalue for 10,000 shots how is it possible that using only 2 shots returns the groundstate eigenvalue?

[daschaich]

@JohnKerf

If the ansatz used can’t return the groundstate eigenvalue for 10,000 shots how is it possible that using only 2 shots returns the groundstate eigenvalue?

If I understand Daniela’s explanation correctly, each “shot” is really an exact energy eigenvalue evaluation, with some probability of returning the ground-state energy E0, some probability of returning the first-excited-state energy E1, etc. Obtaining E0 after 2 shots then just means that E0 is returned two times in a row (which is likely) but obtaining E0 after 10k shots would require E0 being returned 10k times in a row (which is unlikely). At some point in those 10k evaluations, an excited-state energy is coming along and increasing the average away from E0.

I am new to PennyLane and more familiar with Qiskit, where I’m used to each “shot” evaluating eigenvalues for each Pauli string in the hamiltonian / cost function, which seems to correspond to qml.sample

[JohnKerf]

Thank you @daschaich that makes it more clear why we see the difference in distributions as we increase shots.

@daniela.murcillo My main question now is whether using qml.expval with a finite number of shots is suitable to use if my aim is to simulate a real quantum device with finite shots? I don’t believe sampling from the exact eigenvalue distribution is what would be done on a real device?

[daniela.murcillo]

Hi!

I think `qml.expval` is perfectly fine for VQE, as shown by this demo or the Codebook reference I previously mentioned.
However, if what you want is to simulate a real device, this is not what is happening here. In a real device, every time you measure the energy, you are still expected to get one energy eigenvalue, but since it is a real device with noise, you could end up getting any other values due to that noise. Am I understanding correctly what you mean by simulating a real quantum device?

(I also want that sampling from the Paulis, like what you were doing in the first entry of this post, ends up converging to the same distribution as qml.expval as you increase the number of samples)

VQE doesn’t require you to simulate a real device, but if that is something you want to do, you can take a look at default.mixed

We can also take a look at your optimization code, cost function, if you think that might be helpful.

[JohnKerf]

Thanks Daniela and apologies if I am not being clear.

My main concern is that the logic behind qml.expval (sampling from the exact energy eigenvaules) is not what would be done for a real qubit device. For a real device I believe the Hamiltonian would be decomposed into pauli strings and each pauli string is evaluated individually (similar to what I showed in my original post).

This means the two methods don’t really allign and therefore I shouldn’t really be using qml.expval to simulate a device with finite shots if my end goal is to run the circuit on say an IBM qubit device. I appreciate there is also additional noise from gate errors etc when using a real device, however, in terms of variance solely from finite shots, I believe the two methods dont allign.

[daniela.murcillo]

Yeah, when we simulate, we just sample from the energy spectrum, very much like what theory says. Meanwhile, in hardware you need to decompose the Hamiltonian and calculate the expectation values with respect to the Paulis.
And yes, for sure, we can already see that the two methods don’t align for a small number of samples. They are different things that eventually converge.
This is just a thought: I want to point that you would probably would want to run this, even on real hardware, with a large number or shots and that you don’t need a lot of them to see how the Pauli-measured distribution approaches the other one.

Thanks for your question.