VQE from scratch, what have I got wrong?

Question

I've been working on a code to run VQE with a grouped measurement. For some reason, my expectation values are slightly off from those computed by pennylane, the software I decided to use for this little project. I've been staring at my expectation value function exp_val, and I can't tell why it's wrong.

This is the code:

import pennylane as qml
from pennylane import qchem
import numpy as np

symbols = ["H", "H"]
coordinates = np.array([0.0, 0.0, -0.6614, 0.0, 0.0, 0.6614])
h2_ham, n_qubits = qchem.molecular_hamiltonian(symbols, coordinates)
n_shots = 1000

dev_noisy = qml.device('default.qubit', wires = n_qubits, shots = n_shots)

@qml.qnode(dev_noisy)
def VQE_circuit(params, group = None, n_qubits = None):
        qml.StronglyEntanglingLayers(params, wires = range(n_qubits))
        rotations = qml.grouping.diagonalize_qwc_pauli_words(group)[0]
        return qml.probs(wires=range(n_qubits))


import autograd.numpy as agnp
drawer = qml.draw(VQE_circuit)

def exp_val(results, coeffs, groupings):
        E = 0
        for i, result in enumerate(results):
                #Process each list of counts (probs)
                ops = groupings[i]
                coeff_list = coeffs[i]
                #print(drawer(params, group = groupings[i], n_qubits = n_qubits))
                for op_idx, op in enumerate(ops):
                        ##calculate expval for each operator in group
                        if op.name == 'Identity':
                                E += coeff_list[op_idx]
                        else:
                                exp_val = 0
                                for c_idx, count in enumerate(result):
                                        #process bitstring in result
                                        idxs = op.wires.toarray()
                                        bits = format(c_idx, "b").zfill(n_qubits)
                                        sub_bits = [bits[i] for i in idxs]
                                        par = sub_bits.count('1')%2
                                        sign = (-1)**par
                                        exp_val += sign*count
                                exp_val *= coeff_list[op_idx]
                                E += exp_val
        return E
print("\n", h2_ham, "\n")
groupings, coeffs = qml.grouping.group_observables(h2_ham.terms[1], h2_ham.terms[0], grouping_type = 'qwc', method = 'rlf')
param_shape = qml.templates.StronglyEntanglingLayers.shape(n_layers=3, n_wires=n_qubits)
params = np.random.normal(scale=0.1, size=param_shape)
results = [VQE_circuit(params, group = group, n_qubits = n_qubits) for group in groupings]
print(exp_val(results, coeffs, groupings))


##VQE execution:
def naive_cost(params):
    results = [VQE_circuit(params, group = group, n_qubits = n_qubits) for group in groupings]
    return exp_val(results, coeffs, groupings)    
sparse = qml.utils.sparse_hamiltonian(h2_ham).toarray()
gs_E = np.linalg.eigvalsh(sparse)[0]
print("GSE: ", gs_E)
import scipy.optimize as opt
param_shape = qml.templates.StronglyEntanglingLayers.shape(n_layers=5, n_wires=n_qubits)
params = np.random.normal(scale=0.1, size=param_shape)
opt = qml.AdagradOptimizer(0.05)

max_iteration = 100
conv_tol = 1e-6
energy = [naive_cost(params)]
for n in range(max_iteration):
        params, prev_E = opt.step_and_cost(naive_cost, params)
        energy.append(naive_cost(params))
        conv = np.abs(energy[-1]-prev_E)
        if n % 2 == 0:
                print(f"Step = {n},  Energy = {energy[-1]:.8f} Ha")

        if conv <= conv_tol:
                break
#print(energy)
print("FOUND GROUND STATE E: ", np.min(energy))
print("FINAL params: ", params)
print("REAL GROUND STATE E: ", gs_E)

Am I doing something obviously wrong? I've checked to ensure the endianness of the bits is correct, I made sure to calculate parities based on the active qubit indices for the desired observable, I made sure not to include a parity calculation when computing expected value of identity operators. I've triple checked that the coefficients are being applied to the correct operators, that the circuits are appended with the correct measurement circuits for the given operator in the group. Amongst these checks, I can't seem to find the flaw in my programming logic. Having an extra pair of eyes on this would be extremely helpful.

All the best,

cuhrazatee

PS:

Here's the code using Pennylane's built in exp_val.

from openfermion.ops.operators.qubit_operator import QubitOperator
import pennylane as qml
from pennylane import qchem
import numpy as np
from functools import partial
from pennylane.measure import state
from pennylane.ops.qubit import observables
from pennylane.templates import UCCSD
import matplotlib.pyplot as plt

symbols = ["H", "H"]
coordinates = np.array([0.0, 0.0, -0.6614, 0.0, 0.0, 0.6614])
h2_ham, n_qubits = qchem.molecular_hamiltonian(symbols, coordinates)
n_electrons = 2
singles, doubles = qchem.excitations(n_electrons, n_qubits)
s_wires, d_wires = qchem.excitations_to_wires(singles, doubles)
ref_state = qchem.hf_state(n_electrons, n_qubits)
ansatz = partial(UCCSD, init_state = ref_state, s_wires = s_wires, d_wires = d_wires)
groupings, coeffs = qml.grouping.group_observables(h2_ham.terms[1], h2_ham.terms[0], grouping_type = 'qwc', method = 'rlf')

n_shots = 10000
dev_noisy = qml.device('default.qubit', wires = n_qubits, shots = n_shots)

sparse = qml.utils.sparse_hamiltonian(h2_ham).toarray()
gs_E = np.linalg.eigvalsh(sparse)[0]

print("GSE: ", gs_E)
param_shape = qml.templates.StronglyEntanglingLayers.shape(n_layers=5, n_wires=n_qubits)
params = np.random.normal(scale=0.1, size=param_shape)
optimal_params = [[[0.15300575748799206, 0.0802866250748122, 0.6612327808749161],
  [-0.012197292985330403, 1.0054708209216188, -0.5031298708940922],
  [0.4871116388974964, 0.8755791036972337, 0.15300349217856668],
  [0.2125674641016197, -0.3988028820299284, -0.5903810690276766]],
 [[-0.26227065702483116, 0.7762008921102848, 0.20421958031876591],
  [-0.11210374173720475, -0.7158851608015426, -0.945090563307313],
  [-0.20041447281024863, -0.43875447105339715, -0.10552844761324888],
  [0.017868138476782234, -0.4067704016345291, 0.19911933547123295]],
 [[-0.1170885274951583, -0.40203947157121894, -0.44851762637470327],
  [-0.2715291337140317, 0.6888494094283752, 0.5389027752311034],
  [-0.5191082830999312, 0.4426962606005164, -0.25932474764548114],
  [0.8138172470220708, -0.54678942509227, 0.4102491578027457]],
 [[0.7676993827422776, -0.47721469081406376, 0.6337393057184456],
  [-0.6553402727229024, 0.8955375499127577, 0.5789282160827474],
  [0.6570560582613835, -0.8518341967262695, -0.6023881439081624],
  [0.040387149769954125, 0.3040045252649316, -0.38193967606295326]],
 [[0.42336490227621815, -0.3869902716443922, 0.27822533007353994],
  [1.2170710775127433, 0.6959673154584948, -0.11538272838636159],
  [0.33436903989516936, -0.7268804688737179, 0.49919014014531526],
  [0.818911996077618, -0.6280017753881122, 0.553169606547251]]]

cost = qml.ExpvalCost(qml.StronglyEntanglingLayers, h2_ham, dev_noisy, optimize=True)

opt = qml.AdagradOptimizer(0.06)



max_iteration = 100
conv_tol = 1e-6
energy = [cost(params)]
for n in range(max_iteration):
        params, prev_E = opt.step_and_cost(cost, params)
        energy.append(cost(params))
        conv = np.abs(energy[-1]-prev_E)
        if n % 2 == 0:
                print(f"Step = {n},  Energy = {energy[-1]:.8f} Ha")

        if conv <= conv_tol:
                break
#print(energy)
print("FOUND GROUND STATE E: ", np.min(energy))
print("REAL GROUND STATE E: ", gs_E)

```

Answer 1

Having tested your code locally, I actually don't believe there are any errors. But, there are a few reasons that you are not getting the exact expectation value.

The first is just that the device is set to use 1000 shots, so there will always be some fluctuation due to the probabilistic nature of measurement.

The second has to do with the number of optimization steps and convergence criteria. With 1000 shots, the fluctuation between the cost at each optimization step is going to be significantly larger than the 1e-6 convergence threshold; here are a couple examples from my local run:

Step = 66,  Energy = -0.82500835 Ha
Step = 68,  Energy = -0.83819441 Ha
Step = 70,  Energy = -0.84052676 Ha
Step = 72,  Energy = -0.85156650 Ha
Step = 74,  Energy = -0.89118294 Ha
Step = 76,  Energy = -0.88678496 Ha

The convergence criteria will never be met, and so instead the optimization will end when it reaches the maximum number of iterations. At 100 steps, the energy is floating around -1.0, so it actually just needs more time to approach the true value (~250 steps). A good way to sanity check is to use the analytic mode by not setting a shot number for the device; this will give a clearer picture of how the optimization will progress, and how many iterations might be needed (and, additionally, if the parameters of the Optimizer need to be tweaked).