Failure to reproduce quantum approximate optimization algorithm results of linked paper

[or1426]

Information

• Qiskit Aer version: qiskit.__version__ =0.19.1, qiskit.providers.aer.__version__='0.10.0'
• Python version: sys.version_info(major=3, minor=10, micro=0, releaselevel='final', serial=0)
• Operating system: Linux 5.10.0-11-amd64 # 1 SMP Debian 5.10.92-1 (2022-01-18) x86_64 GNU/Linux

What is the current behavior?

Implementing the QAOA algorithm for the Max E3LIN2 as described in the paper (pages 12-14 and especially figure 4 on page 14) leads to incorrect results. Note that this is the paper which the AER extended stabilizer simulator is based on.

At gamma = pi/10 and the parameters given we should obtain QAOA expectation values around -7 (see figure 4a of the linked paper), qiskit's implementation of the extended stabilizer/metropolis simulator produces expectation values around 0.05 (values close to zero are consistent with "junk", randomly sampling bitstrings that are not related to the QAOA instance).

Steps to reproduce the problem

Code available here

https://gist.github.com/or1426/4faf96f5a71151f3aa4a3165cc8db72e

Note that the sample code takes around 10 minutes to run on my computer, to get results much faster set gamma on line 100 to (e.g.) 0.025*pi. At this gamma the correct expectation value is about -2.52, but qiskit produces values around 0.05 in my tests.

The gist above contains some code to randomly generate instances of QAOA circuits, and a fixed QAOA polynomial provided for reproducibility.

Relevant parameters (taken from linked paper):

• method='extended_stabilizer'
• shots = 40000
• extended_stabilizer_metropolis_mixing_time=10000
• extended_stabilizer_approximation_error=0.05
• extended_stabilizer_sampling_method="metropolis"

What is the expected behavior?

Producing approximately correct values as the authors claim their algorithm does in the paper.

Suggested solutions

Notes

Simulating the QAOA circuits using the extended stabilizer simulator at smaller qubit counts produces correct results. The qiskit results diverge from the true QAOA expectation values at around 31 qubits, and are completely incorrect by around 40 qubits.

Using the same code but simulating the circuits instead with the compute algorithm I am in the process of implementing (pull request #1273) produces correct results.

[or1426]

Here is a plot showing how the qiskit extended stabilizer results diverge from the correct QAOA expectation values as the number of qubits increases. For each number of qubits (x axis) I generated 10 circuits, sent them to the compute algorithm (blue dots) and the extended stabilizer subroutine in qiskit with the same parameters as I gave above (red dots) and computed the QAOA expectation value at gamma = 0.1. As you can see below about 31 qubits the extended stabilizer values are close to the true values, but then they rapidly diverge, becoming essentially random noise by 40 qubits.

[hhorii]

I have tested your scripts with small changes to run all the qubits and use all system threads, and then got following results:

python version: sys.version_info(major=3, minor=8, micro=13, releaselevel='final', serial=0)
qiskit version: 0.20.0
aer version: 0.10.4
25 -3.6083999999998557
26 -3.4318500000005057
27 -3.6478999999998107
28 -3.9636249999998316
29 -3.7865000000000215
30 -4.0264749999998735
31 -4.243624999999961
32 -4.34450000000046
33 -4.3883249999999965
34 -4.748775000000143
35 -4.711699999999726
36 -4.792400000000181
37 -4.906300000000374
38 -5.12695000000039
39 -5.191125000000185

The above looks similar to blue dots, I think. Modified QAOA() and main are as follows:

def QAOA(n, polynomial):
    time_s = 0.
    num_shots = 40000
    mixing_time = 10000
    error = 0.05
    gammas = np.array([0.1])*np.pi
    expectation_vals = np.zeros_like(gammas)
    
    for i, gamma in enumerate(gammas):
        circ = QAOA_circuit_U(n, gamma, polynomial)
        circ.measure(list(range(n)), list(range(n)))        
        extended_stabilizer_simulator = AerSimulator(method='extended_stabilizer',
                                                     extended_stabilizer_sampling_method="metropolis",
                                                     extended_stabilizer_metropolis_mixing_time=mixing_time)
        
        t_circ = transpile(circ, extended_stabilizer_simulator)        
        result = extended_stabilizer_simulator.run(t_circ, shots=num_shots,extended_stabilizer_approximation_error=error).result()
        counts = result.get_counts(0)
        for bitstring in counts:
            reversed_string = bitstring[::-1] # qiskit gives us the bitstrings "backwards" with string[0] being the outcome of the measurement on qubit[-1]
            expectation_vals[i] += C(polynomial, reversed_string) * (counts[bitstring]/float(num_shots))

        print(n, expectation_vals[i])    
        
if __name__ == "__main__":
    print("python version:", sys.version_info)
    print("qiskit version:", qiskit.__version__)
    print("aer version:", qiskit.providers.aer.__version__)
    
    for n in range(25, 40):
        D=4
        polynomial = E3LIN2(n, D)
        QAOA(n, polynomial)

[or1426]

Yes, I agree. I just tested it again and the problem is present in aer release version 10.3 but seems to be fixed in 10.4. I guess this issue should be closed then.

[jakelishman]

Sounds good, thanks for confirming!