How to correctly compute expectation value in QAOA?

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?

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.