Am I doing anything wrong when trying to calculate the expectation value in Qiskit on the real hardware?

Question

I am trying to use the method of calculating the expectation values of Pauli operators in Qiskit which I found here.

However, the results obtained via IBMQ sampling differ significantly from the exact value, even though I perform measurement error mitigation.

The state, in which I am trying to calculate the expectation values, is prepared using a simple circuit

          ┌─────────────┐                   ┌───┐                     
q_0: ─────┤ RY(-1.8018) ├──■────────■───────┤ X ├────────────────────
     ┌───┐└──────┬──────┘┌─┴─┐      │       └─┬─┘                    
q_1: ┤ X ├───────■───────┤ X ├──────┼─────────┼──────────────────────
     └───┘               └───┘┌─────┴──────┐  │                 ┌───┐
q_2: ─────────────────────────┤ RY(2.2489) ├──■─────────■───────┤ X ├
                              └────────────┘     ┌──────┴──────┐└─┬─┘
q_3: ────────────────────────────────────────────┤ RY(0.99778) ├──■──
                                                 └─────────────┘     

which after transpiling looks as follows:

For reference, for each Pauli operator I print the exact expectation value and also calculate it a few times using the QASM simulator.

Here's the code:

circuit = QuantumCircuit(4)
circuit.x(1)
circuit.cry(-1.80184863, 1, 0)
circuit.cx(0,1)
circuit.cry(2.24892942,0,2)
circuit.cx(2,0)
circuit.cry(0.9977846,2,3)
circuit.cx(3,2)
psi = CircuitStateFn( circuit )
    
paulis = [ Pauli([1,1,0,0],[1,1,0,0]),  Pauli([1,1,1,1],[1,0,0,1]) ]
shots = 8000
reps = 3   

backend_qasm = qiskit.Aer.get_backend( 'qasm_simulator' )
q_instance_qasm = QuantumInstance( backend_qasm, shots = shots )
    
load_account()
provider = get_provider( hub='ibm-q' )
backend_ibmq = least_busy( provider.backends(filters=lambda x: x.configuration().n_qubits >= 4 and not x.configuration().simulator) )
q_instance_ibmq = QuantumInstance( backend = backend_ibmq,
                                   shots = shots,
                                   measurement_error_mitigation_cls = CompleteMeasFitter,
                                   measurement_error_mitigation_shots = shots )
print(f'IBMQ backend: {backend_ibmq}.\n')     

for pauli in paulis:
    print(f'Pauli operator: {pauli}.')
    pauli = WeightedPauliOperator([[1., pauli]]).to_opflow()

    measurable_expression = StateFn( pauli, is_measurement = True ).compose( psi )
    expectation = PauliExpectation().convert( measurable_expression )

    expect_exact = psi.adjoint().compose( pauli ).compose( psi ).eval().real

    print( f'Exact expectation value: {expect_exact}.' )
    for r in range(reps):
        sampler_qasm = CircuitSampler( q_instance_qasm ).convert( expectation )
        expect_sampling_qasm = sampler_qasm.eval().real
        print( f'Exact expectation, QASM sampling: {expect_sampling_qasm}.' )
    for r in range( reps ):
        sampler_ibmq = CircuitSampler( q_instance_ibmq ).convert( expectation )
        expect_sampling_ibmq = sampler_ibmq.eval().real
        print( f'Exact expectation, IBMQ sampling: {expect_sampling_ibmq}.' )

    print()

And here's the output:

IBMQ backend: ibmq_ourense.

Pauli operator: IIYY.
WARNING - The skip Qobj validation does not work for IBMQ provider. Disable it.
Exact expectation value: -0.4201884924852.
Exact expectation, QASM sampling: -0.42275.
Exact expectation, QASM sampling: -0.4175.
Exact expectation, QASM sampling: -0.4165.
Exact expectation, IBMQ sampling: -0.19053720838213.
Exact expectation, IBMQ sampling: -0.33771371840093.
Exact expectation, IBMQ sampling: 0.14870401826006.

Pauli operator: YZZY.
Exact expectation value: 0.22895884311365.
Exact expectation, QASM sampling: 0.237.
Exact expectation, QASM sampling: 0.2385.
Exact expectation, QASM sampling: 0.2345.
Exact expectation, IBMQ sampling: 0.06862734682344.
Exact expectation, IBMQ sampling: 0.10246703115813.
Exact expectation, IBMQ sampling: 0.13078427261863.

Am I doing something conceptually wrong?

Is there an obvious way to improve the results? (except for doing more shots) Or is it what I should expect to get, given the device's gate fidelities?

Any thoughts/suggestions/corrections greatly appreciated.

Answer 1

I tried to run your code with the same backend as you, ibmq_ourense, and also got the same kind of bad results. Although, I also tried on other backends, first the ibmq_qasm_simulator and I got the exact expectation value, so I assume there is no bug on your code since it is right with the ideal machine. I also tried with ibmq_vigo, which has a better quantum volume than ibmq_ourense(16 vs 8), and I got much better results, closer to the exact expected value.

You could try the "obvious ways" to get better results as you mentioned in your question, maybe looking at the different levels of optimization in the transpile function might help, see the documentation and a tutorial from Qiskit to check how you can play with this!
Finally running your code on a device with a higher quantum volume leading to less noise thus better results might be another way around the errors!

Hope this will help you, feel free to ask any other thing :)