Measuring an operator using construct_evaluation_circuit yields incorrect results

[kuehnste]

Informations

• Qiskit version: qiskit.version yields '0.12.0'
• Operating system: Windows, Linux
• Python version: 3.7.6 on Windows, 3.7.4 on Linux

What is the current behavior?

Using the output of the function construct_evaluation_circuit of the class WeightedPauliOperator to measure an operator yields an incorrect result in a certain case. The value obtained is twice the correct one.

Steps to reproduce the problem

# Imports
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit import execute
from qiskit import IBMQ
from qiskit.quantum_info.operators.pauli import Pauli
rom qiskit.aqua.operators import *
from qiskit.aqua.operators.op_converter import to_matrix_operator
from qiskit import BasicAer
from scipy.linalg import eig

# Define a basic circuit
q = QuantumRegister(2, name='q')
qc = QuantumCircuit(q)
qc.rx(10.9891251356965,0)
qc.rx(6.286692023269373,1)
qc.rz(7.848801398269382, 0)
qc.rz(9.42477796076938,1)
qc.cx(0,1)
qc.draw()

# Define the Hamiltonian
pauli_string = []
pauli_string.append([1.0, Pauli.from_label('XX')])
pauli_string.append([-1.0, Pauli.from_label('YY')])
pauli_string.append([-1.0, Pauli.from_label('ZZ')])
H = WeightedPauliOperator(pauli_string)

# Show the Hamiltonian and compute its spectrum with exact diagonalization
Hmatrix = to_matrix_operator(H)
eigvals = eig(Hmatrix.dense_matrix)[0]
print("Hamiltonian: " + H.print_details())
print("Eigenvalues: " + str(eigvals))

# Output obtained:
# Hamiltonian: XX	(1+0j)
# YY	(-1+0j)
# ZZ	(-1+0j)

# Eigenvalues: [ 1.+0.j -3.+0.j  1.+0.j  1.+0.j]

# Now we measure the expectation value using built-in Qiskit functionality
evaluation_circuits = H.construct_evaluation_circuit(qc, False)
backend = BasicAer.get_backend('qasm_simulator')
job = execute(evaluation_circuits, backend, shots=1024)
# Evaluate the results
expectation_value, standard_deviation = H.evaluate_with_result(job.result(), False)
print("Expectation value: " + str(expectation_value) + " +- " + str(standard_deviation))

# Output obtained:
# Expectation value: (-3+0j) +- 0j

# Define the Hamiltonian but this time we split each of the terms in a sum of twice the identical part
pauli_string = []
pauli_string.append([0.5, Pauli.from_label('XX')])
pauli_string.append([-0.5, Pauli.from_label('YY')])
pauli_string.append([-0.5, Pauli.from_label('ZZ')])
pauli_string.append([0.5, Pauli.from_label('XX')])
pauli_string.append([-0.5, Pauli.from_label('YY')])
pauli_string.append([-0.5, Pauli.from_label('ZZ')])
H2 = WeightedPauliOperator(pauli_string)

# Show the Hamiltonian and compute its spectrum with exact diagonalization
Hmatrix2 = to_matrix_operator(H2)
eigvals2 = eig(Hmatrix2.dense_matrix)[0]
print("Hamiltonian: " + H2.print_details())
print("Eigenvalues: " + str(eigvals2))

# Output obtained:
# Hamiltonian: XX	(1+0j)
# YY	(-1+0j)
# ZZ	(-1+0j)

# Eigenvalues: [ 1.+0.j -3.+0.j  1.+0.j  1.+0.j]

# Now we measure the expectation value using built-in Qiskit functionality
evaluation_circuits2 = H2.construct_evaluation_circuit(qc, False)
backend2 = BasicAer.get_backend('qasm_simulator')
job2 = execute(evaluation_circuits2, backend2, shots=1024)
# Evaluate the results
expectation_value2, standard_deviation2 = H2.evaluate_with_result(job2.result(), False)
print("Expectation value: " + str(expectation_value2) + " +- " + str(standard_deviation2))

# Output obtained:
# Expectation value: (-6+0j) +- 0j <-- twice the correct result

What is the expected behavior?

Since it is the same Hamiltonian and also the function print_details() shows that it is the same WeightedPauliOperator, the measurement should yield a value of -3 in both cases.

Suggested solutions

Measurement should yield consistent results.

[woodsp-ibm]

The measurement is an evaluation of the operator against the state as prepared by the circuit. To have the eigenvalue as you expect that state supplied by the circuit needs to be the corresponding eigenstate.

Try running VQE on H2 which uses the variational principle to find the minimum eigenvalue. The circuit as described by the variational form (ansatz) at that point will describe eigenstate. For instance the following added to the end of your file correctly evaluates to -3.0 . If you get the parameters for RY at that point you can get the circuit corresponding to that and use it in your code to see that with the right state you will get -3.0

# Get min eigval using VQE
from qiskit.aqua.algorithms import VQE
from qiskit.aqua.components.variational_forms import RY
from qiskit.aqua.components.optimizers import SLSQP

result = VQE(H2, RY(H2.num_qubits), SLSQP()).run(quantum_instance=BasicAer.get_backend('statevector_simulator'))
print('VQE result: {}'.format(result['energy']))

[kuehnste]

The circuit given in my example prepares the ground state of the Hamiltonian, which is also unique as the exact diagonalization results show. My point is not that I want to get the ground state, but rather that H and H2 are the same operator, the only difference is that H2 is constructed as H2 = 0.5 * XX - 0.5 * YY - 0.5 * ZZ + 0.5 * XX - 0.5 * YY -0.5 * ZZ, whereas H is constructed as H = XX - YY - ZZ. Hence, up to statistical fluctuations due to a finite number of shots, the qasm simulator should yield identical results in both cases if I evaluate it with the same circuit. In fact <psi|H|psi> / <psi|psi> cannot be lower than the ground state energy of the system (up to accumulating errors from measuring the Hamiltonian terms individually with a finite number of shots). So how can measuring H2 yield -6.0?

[woodsp-ibm]

Thanks it turned out to be a bug in the processing/simplification of the pauli list that is passed.