Grover's algorithm returns skewed probability distribution

Question

I wrote an implementation of Grover's algorithm that looks like this:

from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit, Aer, execute

# Initialize circuit
m_qubit = QuantumRegister(1)
search_register = QuantumRegister(4)
result_register = ClassicalRegister(4)
ancillaries = QuantumRegister(3)
circuit = QuantumCircuit(search_register, result_register, m_qubit, ancillaries)

# Put M qubit into 1-superposition
circuit.x(m_qubit)
circuit.h(m_qubit)

# Put search qubits into superposition
circuit.h(search_register)

for _ in range(2):

    # Encode S1 * !S2 * S3
    circuit.x( search_register[2] )
    circuit.ccx( search_register[1], search_register[2], ancillaries[0] )
    circuit.ccx( search_register[3], ancillaries[0], ancillaries[1] )
    circuit.x( search_register[2] )

    # Encode S0 * S1
    circuit.ccx( search_register[0], search_register[1], ancillaries[2] )

    # Encode oracle ((S0 * S1) + (S1 * !S2 * S3))
    circuit.x(ancillaries)
    circuit.ccx( ancillaries[1], ancillaries[2], m_qubit[0] )
    circuit.x(ancillaries)
    circuit.x(m_qubit)

    # Reset ancillaries to be used later
    circuit.reset(ancillaries)

    # Do rotation about the average
    circuit.h(search_register)
    circuit.x(search_register)
    circuit.ccx( search_register[0], search_register[1], ancillaries[0] )
    circuit.ccx( search_register[2], ancillaries[0], ancillaries[1] )
    circuit.ccx( search_register[3], ancillaries[1], m_qubit[0] )
    circuit.x(search_register)
    circuit.x(m_qubit)
    circuit.h(search_register)

    # Reset ancillaries for use later
    circuit.reset(ancillaries)

circuit.measure(search_register, result_register)

# Run the circuit with a given number of shots
backend_sim = Aer.get_backend('qasm_simulator')
job_sim = execute(circuit, backend_sim, shots = 1024)
result_sim = job_sim.result()

# get_counts returns a dictionary with the bit-strings as keys
# and the number of times the string resulted as the value
print(result_sim.get_counts(circuit))

Essentially, this circuit should search for any value that matches $X_1 X_2 X_3 + X_0 X_1$. When I do the math for this circuit I get the following:

$$\frac{1}{1024 \sqrt{2}} \left( -304 \left( \sum_{x \in \{00, 01, 11\}} \left| x \right> \left( \left|00\right> + \left|01\right> + \left|10\right> \right) + \left|1000\right> + \left|1001\right> \right) + 80 \left( \sum_{x = 00}^{11} \left|x\right> \left|11\right> + \left|1010\right> \right) \right) \left( \left|0\right> - \left|1\right> \right)$$

These amplitudes correspond to probabilities of $-\frac{{304}^{2}}{2({1024}^{2})} = -\frac{{19}^{2}}{2({64}^{2})}$ and $\frac{{80}^{2}}{2({1024}^{2})} = \frac{{5}^{2}}{2({64}^{2})}$ for the unsuccessful and successful probabilities, respectively.

What I don't understand from this is how the successful outputs would appear over the unsuccessful ones, because the magnitude of the probability for the unsuccessful results is much higher. Can someone explain how the normalization is done in this case?

Furthermore, these results do not obtain when I actually run the circuit. Instead of a step-function for the probabilities, when I run 1024 shots I tend to get between 30 and 50 for the unsuccessful results and the successful results are spread between 200 and 460. While it's obvious that the successful results have been selected, the spread is concerning because the math seems to indicate that they should be fairly close together. Moreover, the successful results always fall in the same order, which is also counter-intuitive given that results are probabilistic. This implies one of two things to me:

1. The superposition isn't actually equal (i.e. there's some bias towards one value versus the other). However, this doesn't make sense as these results were obtained with a simulator so they should be very close to optimal.
2. I have inadvertently entangled the qubits in such a way that the result is affected.

Any insight into this would be greatly appreciated.

Answer 1

It turns out the issue with this circuit was in the use of the reset operator in conjunction with entanglement and superposition.

In the oracle code I submitted, I use a Toffoli gate several times, which creates an entanglement between both input qubits and the output qubit. Since the input qubits are in superposition, that superposition state then "spreads" (for lack of a better term) to the output qubit. So, when the output qubit is reset using the reset operator, its superposition state and entanglements with the input qubits causes certain states to collapse, which changes the output of the circuit as a whole.

This had confused me earlier because using reset with either superposition or entanglement causes no issue, because both can collapse without change to the input; because, in a sense, the output and input qubits are separated. It is only, therefore, when they are used in conjunction that reset causes this behavior.

The solution then, is quite simple: reverse the circuit to return the output qubit to a $\left|0\right>$ state after it is no longer needed:

# Encode S1 * !S2 * S3
circuit.x( search_register[2] )
circuit.ccx( search_register[1], search_register[2], ancillaries[0] )
circuit.ccx( search_register[3], ancillaries[0], ancillaries[1] )

# Encode S0 * S1
circuit.ccx( search_register[0], search_register[1], ancillaries[2] )

# Encode oracle ((S0 * S1) + (S1 * !S2 * S3))
circuit.x(ancillaries)
circuit.ccx( ancillaries[1], ancillaries[2], m_qubit[0] )
circuit.x(m_qubit)

# Return ancillaries to 0s so they can be used later
circuit.x(ancillaries)
circuit.ccx( search_register[0], search_register[1], ancillaries[2] )
circuit.ccx( search_register[3], ancillaries[0], ancillaries[1] )
circuit.ccx( search_register[1], search_register[2], ancillaries[0] )
circuit.x( search_register[2] )

# Do rotation about the average
circuit.h(search_register)
circuit.x(search_register)
circuit.ccx( search_register[0], search_register[1], ancillaries[0] )
circuit.ccx( search_register[2], ancillaries[0], ancillaries[1] )
circuit.ccx( search_register[3], ancillaries[1], m_qubit[0] )
circuit.x(search_register)
circuit.x(m_qubit)

# Return ancillaries to 0s for use later
circuit.ccx( search_register[2], ancillaries[0], ancillaries[1] )
circuit.ccx( search_register[0], search_register[1], ancillaries[0] )
circuit.h(search_register)

Doing this instead produced the desired result.