How to correctly compute expectation value in QAOA?

Question

In QAOA after prepating $|\psi(\gamma, \beta) \rangle$, expectation value $\langle \psi(\gamma, \beta)|H_c|\psi(\gamma, \beta) \rangle$ is computed. In tutorials, I see two approaches to this calculation. The first is to calculate the expectation of an observable $H_c$, as in this demo, which I imagine adds some CNOT gates, because measurement in $Z \otimes Z$ is accomplished by CNOT, as shown here. My code for this approach is as follows, using Pennylane (I am using the Ising model):

def construct_cost_hamiltonian(h, J):
    coeffs = []
    obs = []

    # Add terms for individual variables
    for k, v in h.items():
        coeffs.append(v)
        obs.append(qml.PauliZ(k))

    # Add terms for interactions between variables
    for k, v in J.items():
        coeffs.append(v)
        obs.append(qml.PauliZ(k[0]) @ qml.PauliZ(k[1]))
    return qml.Hamiltonian(coeffs, obs)

cost_h = construct_cost_hamiltonian(h, J)

@qml.qnode(dev)
def cost_function(params):
    if len(params.shape) < 2:
        params = np.split(params, 2)
    qaoa_circuit(params[0], params[1], h, J)
    return qml.expval(cost_h)

The second approach is to measure bitstrings, map them to $-1,1$ and compute the average of $x^T H_c x$, which makes sense since we get the final result just by measuring, and $x^T H_c x$ is our objective. This approach was used in this demo. My implementation as follows:

@qml.qnode(dev)
def samples_circuit(params):
    if len(params.shape) < 2:
        params = np.split(params, 2)
    qaoa_circuit(params[0], params[1], h, J)
    return qml.sample()

samples = samples_circuit(params)

H = np.zeros((number_of_variables, number_of_variables))
for i in range(0, number_of_variables):
    H[i, i] = h[(i,)]
    for j in range(i+1, number_of_variables):
        H[i, j] = J[i, j]

expect_val = 0
for sample in samples:
    sample_ising = -2*sample + 1
    expect_val += sample_ising.T @ H @ sample_ising
expect_val /= samples.shape[0]

And these two approaches (in my code) give drastically different results with the same parameters. So what approach is correct? Or what is my mistake?

Answer 1

My mistake is that I compute $x^T H_c x$ instead of $xh + x^T J x$, since if I put $h$ terms on the diagonal of $H_c$, all coefficients of $h_i$ would be equal to one regardless of $x_i=1$ or $x_i=-1$, because of squaring. Corrected code:

h_vector = np.zeros((number_of_variables))
J_matrix = np.zeros((number_of_variables, number_of_variables))
for i in range(0, number_of_variables):
    h_vector[i] = h[(i,)]
    for j in range(i+1, number_of_variables):
        J_matrix[i, j] = J[i, j]

expect_val = 0
for sample in samples:
    sample_ising = -2*sample + 1
    expect_val += h_vector @ sample_ising + sample_ising @ J_matrix @ sample_ising
expect_val /= samples.shape[0]

Now both approaches give similar numerical results. But I'm still not 100% sure that they are equivalent theoretically.