XX and YY and ZZ Hamiltonians in vqe

There is a type cirq.PauliSum created when you add together products of Pauli operations. This type has a method expectation_from_state_vector and expectation_from_density_matrix. This is the easiest way to get the expectation values, if you're just calculating them instead of estimating them from samples taken from hardware.

Note that both methods require you to specify an index for each qubit. This is because the state vector is just a numpy array, with no information about which axes (or which bits of the index) correspond to which qubits..

import cirq

qubits = cirq.LineQubit.range(4)
pauli = cirq.Z

operator: cirq.PauliSum = sum(pauli(qubits[k - 1]) * pauli(qubits[k]) for k in range(4))
print("operator", operator)
# operator 1.000*Z(0)*Z(3)+1.000*Z(0)*Z(1)+1.000*Z(1)*Z(2)+1.000*Z(2)*Z(3)

a, b, c, d = qubits
circuit = cirq.Circuit(
    cirq.H(a),
    cirq.CNOT(a, b),
    cirq.CNOT(a, c),
    cirq.CNOT(a, d),
)
final_state = cirq.final_state_vector(circuit, qubit_order=qubits)

expectation = operator.expectation_from_state_vector(
    final_state,
    qubit_map={q: q.x for q in qubits})
print("z_expectation", expectation)
# z_expectation (3.999999761581421+0j)

If you instead want to estimate the operators based on samples, the process is more manual. For example, you could make three separate variations of the circuit. One where you measure all the qubits in the X basis, one where you measure all the qubits in the Y basis, and one where you measure all the qubits in the Z basis. You can then multiply (or xor) the individual measurement results together to get the paired measurement results and compute the average.

sampler = cirq.Simulator()  # or a hardware sampler

circuit_x = circuit + [[cirq.H(q), cirq.measure(q)] for q in qubits]
circuit_y = circuit + [[cirq.X(q)**0.5, cirq.measure(q)] for q in qubits]
circuit_z = circuit + [cirq.measure(q) for q in qubits]

import numpy as np
import pandas as pd
x_samples: pd.DataFrame = sampler.sample(circuit_x, repetitions=1000)
x_cols = [x_samples[str(q)] for q in qubits]
x_parity_bits = np.array([x_cols[k-1] ^ x_cols[k] for k in range(4)], dtype=np.int8)
x_parity_signs = 1 - 2 * x_parity_bits
x_expectation = np.mean(x_parity_signs)
print("x_expectation", x_expectation)
# x_expectation -0.011