General Grover's search showing wrong results on a particular instance

Question

I am trying to implement a general Grover's search for n qubits searching for m marked elements. I happened to observe that for the case of multiple solutions, I could stack together the oracles for each marked state and the result would be phase flip of each marked state. For example, for marked states 00 and 10,I would design an oracle for both of them and stack them together. Same goes for higher qubits and more marked elements. Then I designed a general diffusion circuit. I'm then iterating the oracle and diffusion circuits T times, where T is given by $$ T = \Bigg\lfloor{\frac{\pi\cdot k}{4\arcsin\sqrt{m/N}}}\Bigg\rfloor $$

The problem is I'm getting correct results for almost every other state except n=3 qubits and m=4 solutions. I am taking marked states to be marked = [ '111', '010', '011', '000'] and the results I am getting from running the circuit through qasm simulator is not what I expect:

This is the code I wrote:

#Define and initialize the circuit

#Number of qubits
n = 3

#Define and initialize circuit
q = QuantumRegister(n, 'q')
c = ClassicalRegister(n ,'c')
ckt = QuantumCircuit(q, c)

#Put all initial states into superposition
ckt.h(q)

#Define the marked item

marked = [ '111', '010', '011', '000']
m = len(marked)

#Implement the oracle function for given marked item
#Checked, correct

def oracle(ckt, marked):

    for j in marked:
    
    #Put a X-gate for 0 qubits
        for i in range(len(j)):
            if j[i] == '0':
                ckt.x(q[i])

        #Implement a multi-controlled Z-gate
        ckt.h(q[n-1])
        ckt.mcx([q[i] for i in range(len(j)-1)], q[n-1])
        ckt.h(q[n-1])

        #Put X gates on those states that were acted upon by X-gates before
        for i in range(len(j)):
            if j[i] == '0':
                ckt.x(q[i])


#Implement the diffuser function

def diffuser(ckt):
    #Hadamard on all qubits
    ckt.h(q)

    #X-gates on all qubits
    ckt.x(q)

    #MCZ gate
    ckt.h(q[n-1])
    ckt.mcx([q[i] for i in range(n-1)], q[n-1])
    ckt.h(q[n-1])

    #X-gate on all qubits again
    ckt.x(q)

    #H-gates on all qubits again
    ckt.h(q)

#Set number of iterations

import math
import numpy as np

N = 2**n

def T(k):
    return math.floor((np.pi*k)/(4*np.arcsin(math.sqrt(m/N))))

k = 1

while True:
    if T(k) > 0:
        break
    else:
        k = k + 2
iterations = T(k)

#Iterate the circuit for 'iterations' times

for i in range(iterations):
    oracle(ckt, marked)
    diffuser(ckt)

I would really appreciate someone telling me the reason for this malfunction.

Answer 1

Grover's Algorithm doesn't work when the number of valid solutions is exactly half the number of possibilities, similar behavior should be seen for 2 qubits and 2 solutions and 4 qubits and 8 solutions and so on. Given an oracle $O$ of this sort that is guaranteed to leave each computational basis state either unchanged or phase flipped, a quantum state $|\phi\rangle$ can be written in the form $|\phi\rangle = \cos(\theta)|s'\rangle + \sin(\theta)|\omega\rangle$ where $|s'\rangle$ and $|\omega\rangle$ are two unit vectors such that $|s'\rangle$ is unaffected by the oracle and $|\omega\rangle$ is flipped. If $B$ are the computational basis states unaffected by the oracle and $P$ the ones flipped, then $|s'\rangle = \sum_{k \in B}\frac{1}{\sqrt{2^n - m}}|k\rangle$ and $|\omega\rangle = \sum_{k \in P}\frac{1}{\sqrt{m}}|k\rangle$ allows the uniform superposition $|s\rangle = \sum_{k = 0}^{2^n - 1}\frac{1}{\sqrt{2^n}}|k\rangle$ to be written as $|s\rangle = \sqrt{\frac{2^n - m}{2^n}}|s'\rangle + \sqrt{m/2^n}|\omega\rangle$.

The relevance of this is that, for the net Grover iteration operator $G = (2|s\rangle\langle s| - I)O$, performing $q$ iterations on $|s\rangle$ gives $G^q|s\rangle = \cos((2q + 1)\theta)|s'\rangle + \sin((2q + 1)\theta)|\omega\rangle$ (the oracle flips across $|s'\rangle$ and the diffusion across $|s\rangle$, leading to a net rotation double-opposite of that of rotating from $|s\rangle$ to $|s'\rangle$ see this diagram for the base idea) where $\theta$ is such that $|s\rangle = \cos(\theta)|s'\rangle + \sin(\theta)|\omega\rangle$. The value of $\theta$ for $n$ qubits and $m$ solutions is $\arcsin(\sqrt{m/2^n})$. If $m = 2^{n -1}$, this implies $\theta = \frac{\pi}{4}$, and thus $q$ Grover iterations will lead to the coefficient of $|\omega\rangle$ being $\sin(\frac{(2q + 1)\pi}{4})$, which will always give a probability of $\frac{1}{2}$.