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)
```
exp_val
is the most important in the code above. The way the execution works is for each set of compatible observables there is a probability distribution inresults
corresponding to the statistics of the prepared wavefunction. Then I do parity averaging on the operators in that set to compute their individual expectation values. This is done for every set of compatible observables, and is added together to produce the expected value of the entire hamiltonian. $\endgroup$