simulate_expectation_values() produces incorrect expectation values

[veikopalge]

Description: Sometimes simulate_expectation_values() produces incorrect expectation values. Noticed while trying to compute H2 ground state energy with the 6-31g basis.

Have confirmed this behavior for different installation instances.

Hope the script below that reproduces the behavior uses cirq and qsim routines correctly---can you check this? Thanks.

Versions:
OS: linux (Ubuntu 22.04)
Python: 3.10.4
openfermion==1.5.1
qsimcirq==0.15.0
cirq==1.1.0

How to reproduce: use the python program below. The python script is a toy example. It computes the average value of an observable in exact simulation with two different methods: statevector + expectation_from_state_vector() and simulate_expectation_values(). The script shows results for different observables (the observable is formed as a single PauliSum by adding incrementally strings from a pre-defined set HAMILTONIAN) in order to demonstrate that the behavior occurs for some but not all observables. The expectation values computed with two different methods agree most of the time, and then suddenly start to diverge. This behavior can be seen in a larger computation too. Unfortunately, for the large computation at hand the final result is 300% off when computed with simulate_expectation_values(). Since this affects both CPU and GPU, we can't use the GPU for simulations.

Why think simulate_expectation_values() gives wrong values while statevector + expectation_from_state_vector() is correct? This is based on actual values for H2 ground state energy. So far, the latter routine has produced correct values for all observables we tested.

Also, by trial and error we noticed that when identity is taken out of the observable, then the two methods always agree. To see this, comment in/out the first, identity term below in HAMILTONIAN. However, since it's just an observation based on a limited set of instances, there's no guarantee this work-around will work for an arbitrary observable.

The behavior is present in both CPU and GPU simulations and larger circuits (the one in the toy example is a trivial circuit).

Python script for reproducing the issue:

from openfermion import QubitOperator
from openfermion import qubit_operator_to_pauli_sum
import cirq
import qsimcirq
import numpy as np

# Comment in/out the identity term and observe agreement/disagreement between
# expectation values
HAMILTONIAN = {
'': -0.09706626861762856,
'X0 X1 Y2 Y3': -0.045302615508689345,
'X0 Y1 Y2 X3': 0.045302615508689345,
'Y0 X1 X2 Y3': 0.045302615508689345,
'Y0 Y1 X2 X3': -0.045302615508689345,
'Z0': 0.17141282639402444,
'Z0 Z1': 0.16868898168693283,
'Z0 Z2': 0.12062523481381827,
'Z0 Z3': 0.1659278503225076,
'Z1': 0.17141282639402455,
'Z1 Z2': 0.1659278503225076,
'Z1 Z3': 0.12062523481381827,
'Y2 Z3 Z4 Z5 Y6': 1.0,
'Z2': -0.22343153674663868,
'Z3': -0.22343153674663868}


def get_id_qubit_map(qubits: 'Qid',
                     count: int):
    """Returns the id qubit map."""

    return {qubits[index]: index for index in range(count)}


# main

# tolerance
tolerance = 1e-6

# force no of qubits
no_qubits = 8

# create a trivial circuit
qubits = cirq.LineQubit.range(no_qubits)
circuit = cirq.Circuit([cirq.Z(q) for q in qubits])

# print circuit
print("Will use this circuit with {} qubits:".format(no_qubits))
print(circuit)

# simulator
simulator = qsimcirq.QSimSimulator()

print("\nHamiltonian has {} terms.".format(len(HAMILTONIAN)))

print("\nNOTE: The two methods agree when the Hamiltonian DOES NOT contain the Id term. "
      "They may disagree when the Hamiltonian contains the Id term.")
print("\nThe following two methods should give the same average: \n")

# loop over hamiltonian, each time take one more term, run sim and compare
# the two methods of calculating avg value
qubit_map = get_id_qubit_map(qubits, no_qubits)
avg_1 = []
avg_2 = []
observable_qop = QubitOperator()
no_terms = 0
for k, v in HAMILTONIAN.items():
    observable_qop += QubitOperator(k, v)
    observable_psum = qubit_operator_to_pauli_sum(observable_qop)

    # method 1
    result = simulator.simulate(program=circuit)
    state_vec = result.final_state_vector
    avg_1.append(np.real(
                        observable_psum.expectation_from_state_vector(
                            state_vec, qubit_map)))
    # method 2
    avg_list = simulator.simulate_expectation_values(
                        program=circuit,
                        observables=observable_psum)
    avg_2.append(np.real(avg_list[0]))

    print("{}: Method (1 | 2) averages = ({} | {})".format(no_terms, avg_1[no_terms],
                                                       avg_2[no_terms]))
    if (abs(avg_1[no_terms] - avg_2[no_terms])) > tolerance:
        print("\n +++ ATTENTION +++! Methods diverge at term no {}\n".format(no_terms))
    no_terms += 1

    print("HAMILTONIAN = {}".format(observable_qop))

    print("--- --- --- ---")

[sergeisakov]

I think this issue was fixed in #588.

[veikopalge]

Yes thanks, qsim-v0.16.3 works: the above script does not show discrepancies.