MaximumLikelihoodAmplitudeEstimation returning wrong answer from PhaseOracle

[frankharkins]

Environment

• Qiskit Terra version: 0.18.3
• Python version: 3.9.8
• Operating system: OSX

What is happening?

For some reason, the MaximumLikelihoodAmplitudeEstimation is giving a wrong answer when the problem uses a PhaseOracle created from a DIMACS file. Seems to work correctly when using a similar, hand-madeQuantumCircuit as an oracle.

How can we reproduce the issue?

from qiskit import QuantumCircuit
from qiskit.algorithms import EstimationProblem
from qiskit.circuit.library import PhaseOracle, GroverOperator
from qiskit.algorithms import (AmplitudeEstimation,
                               MaximumLikelihoodAmplitudeEstimation)
from qiskit.providers.aer import QasmSimulator
backend = QasmSimulator()
n = 3

# create state preparation operator
state_prep = QuantumCircuit(n)
state_prep.h(range(n))

# create Grover operator from problem file
oracle = PhaseOracle.from_dimacs_file("examples/3sat.dimacs")
grover_op = GroverOperator(oracle)

problem = EstimationProblem(state_prep,
                            [*range(n)],
                            grover_operator=grover_op)

# Correct result
estimator = AmplitudeEstimation(7, quantum_instance=backend)
result = estimator.estimate(problem)
result.estimation
>>> 0.1295244

# Incorrect result!
estimator = MaximumLikelihoodAmplitudeEstimation(7, quantum_instance=backend)
result = estimator.estimate(problem)
result.estimation
>>> 0.4902973909096974

I've attached the DIMACS file (with a .txt extension so github will accept it): 3sat.dimacs.txt

Replacing oracle with:

# create Grover operator by hand
oracle = QuantumCircuit(3)
oracle.h(2)
oracle.x(2)
oracle.ccx(0,1,2)
oracle.x(2)
oracle.h(2)
grover_op = GroverOperator(oracle)

behaves correctly.

What should happen?

The final output should be close to 0.125.

Any suggestions?

No response

[Cryoris]

This is a strange behavior indeed! But it doesn't seem to be related to the PhaseOracle itself. Here is a more minimal example to recreate the same bug

from qiskit import QuantumCircuit
from qiskit.algorithms import EstimationProblem
from qiskit.circuit.library import PhaseOracle, GroverOperator
from qiskit.algorithms import MaximumLikelihoodAmplitudeEstimation, IterativeAmplitudeEstimation
from qiskit.providers.aer import StatevectorSimulator

backend = StatevectorSimulator()
n = 2

# create state preparation operator
state_prep = QuantumCircuit(n)
state_prep.h(range(n))

# create Grover operator from problem file
oracle = QuantumCircuit(n)
oracle.x(1)  # w/o the two X gates around the target qubit, MLAE gives the correct result
oracle.cz(0,1)
oracle.x(1)
grover_op = GroverOperator(oracle)

problem = EstimationProblem(state_prep,
                            [*range(n)],
                            grover_operator=grover_op)

# Correct result
estimator = IterativeAmplitudeEstimation(0.01, 0.05, quantum_instance=backend)
result1 = estimator.estimate(problem)
print(result1.estimation)

# Incorrect result!
estimator = MaximumLikelihoodAmplitudeEstimation(3, quantum_instance=backend)
result2 = estimator.estimate(problem)
print(result2.estimation)

If we take a look at the log-likelood functions the MLAE minimizes, it does find the correct minimum, so the issue doesn't seem to be in the (volatile) minimization process. Rather the bug might be before, already in the calculation of the good and bad counts?

[Cryoris]

Ok, so going through the math I'm getting the right result. Your state preparation operator is

A = H x H x H

where x is a tensorproduct, and you identify a good state as all 3 qubits being one: |111>. Hence the probability of measuring a good state, which is the output of QAE, is 1/(2^3) = 0.125 -- as you mentioned.

Now, you could just define the EstimationProblem with this information: the A operator and the index of the objective qubits:

problem = EstimationProblem(state_prep, [*range(n)])

However you also want to specify the Grover operator. This operator implements

Q = A S_0 A^\dagger S_{bad}

where S_0 is a reflection about 0 and S_{bad} a reflection about the bad state (i.e. all states except the |111> here). You can construct it, by providing it the S_{bad} operator, which would implement

S_{bad}|x> = -|111> if x = 111, else |x>

This operation is implemented e.g. by the CC-Z gate. In your code however, you're providing the S_{bad} as CC-(XZX) gate! This doesn't implement the right reflection and thus doesn't match the good state you're looking for. If we run the code with this oracle, it does give the right result:

from qiskit import QuantumCircuit
from qiskit.algorithms import EstimationProblem
from qiskit.circuit.library import PhaseOracle, GroverOperator
from qiskit.algorithms import MaximumLikelihoodAmplitudeEstimation, IterativeAmplitudeEstimation
from qiskit.providers.aer import StatevectorSimulator

backend = StatevectorSimulator()
n = 3

# create state preparation operator
state_prep = QuantumCircuit(n)
state_prep.h(range(n))

# create Grover operator from problem file
# oracle = PhaseOracle.from_dimacs_file("3sat.dimacs")

oracle = QuantumCircuit(n)
# oracle.x(2)  # no X here!
oracle.h(2)
oracle.ccx(0,1,2)
oracle.h(2)
# oracle.x(2)  # no X here!
grover_op = GroverOperator(oracle, state_preparation=state_prep)

problem = EstimationProblem(state_prep, objective_qubits=list(range(n)), grover_operator=grover_op)

# Correct result
estimator = IterativeAmplitudeEstimation(0.01, 0.05, quantum_instance=backend)
result1 = estimator.estimate(problem)
print(result1.estimation)

# Now, also correct result!
estimator = MaximumLikelihoodAmplitudeEstimation(3, quantum_instance=backend)
result2 = estimator.estimate(problem)
print(result2.estimation)

Now your definition is coming from a dimacs file, where does that definition come from?

[frankharkins]

@Cryoris thanks for looking into this! I see now I added some useless X-gates, but there still should be only one solution. If I run your second example I get:

0.12499999999999994
0.002667906599823475

which is still wrong :/

I took the dimacs file from the textbook and added a couple of lines to eliminate the other solutions.

[Cryoris]

Hi Frank -- my bad! I posted exactly the wrong oracle. The oracle should be only a CCZ, so

oracle = QuantumCircuit(n)
oracle.h(2)
oracle.ccx(0,1,2)
oracle.h(2)

I'll adjust the snippet above.

[woodsp-ibm]

I believe this issue was addressed above and as such I am closing it.