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)