Why is Grover's algorithm not converging monotonically? (Python implementation)

Question

I have tried to implement Grover's algorithm for three qubits in python/numpy and the first two iterations work like a charm but the third one starts to diverge. Is this expected, or is there a bug in the code? I expected the inversion around the mean to blow up the coefficient of the marked state in each iteration. The code follows Nielsen closely, with a silly Oracle that hard-codes the negation of the basis state. The Python code requires a little bit of set-up but the essential idea is:

• Define a phase shift operator of -1 on all basis states except zero.
• Define a reflection around the mean operator as Hadamard ⊗ Phase Shift ⊗ Hadamard.
• Define the full Grover step as the Oracle followed by the reflection.
• Start with an equally balanced state.
• Repeatedly apply the Grover step.

Source code:

import numpy as np

def dagger(m):
    return np.transpose(np.conjugate(m))

def proj(m):
    return m * dagger(m)

# identity matrix for 3 qubits = 8x8 matrix
id3 = np.identity(2**3)

# hadamard matrix for 1, 2 and 3 qubits
H1 = np.matrix([[1.0, 1.0], [1.0, -1.0]], dtype=np.complex256) / np.sqrt(2)
H2 = np.kron(H1, H1)
H3 = np.kron(H2, H1)

# 3 qubit zero vector |000>
zero3 = np.array([[1],[0],[0],[0],[0],[0],[0],[0]], dtype=np.complex256)

# phase shift operator  2*|0><0| - I  for 3 qubits
PS3 = 2 * proj(zero3) - id3

# reflection around the mean
R = H3 * PS3 * H3

# 3 qbit oracle, marking/negating state |101> = column vector (0 0 0 0 0 1 0 0)
O = id3
O[5,5] = -1

# grover operator
G = R * O

# start state |000>
x0 = H3 * zero3

# apply grover step three times
x1 = G * x0
print x1

x2 = G * x1
print x2

x3 = G * x2
print x3

The output of the program is shown below. The coefficient (driving the probability) for the state to search for is 0.88 after one iteration, then 0.97 but then falls back to 0.57. Am I missing any essential step in the algorithm?

[[ 0.1767767+0.0j]
 [ 0.1767767+0.0j]
 [ 0.1767767+0.0j]
 [ 0.1767767+0.0j]
 [ 0.1767767+0.0j]
 [ 0.88388348+0.0j]
 [ 0.1767767+0.0j]
 [ 0.1767767+0.0j]]
[[-0.088388348+0.0j]
 [-0.088388348+0.0j]
 [-0.088388348+0.0j]
 [-0.088388348+0.0j]
 [-0.088388348+0.0j]
 [ 0.97227182+0.0j]
 [-0.088388348+0.0j]
 [-0.088388348+0.0j]]
[[-0.30935922+0.0j]
 [-0.30935922+0.0j]
 [-0.30935922+0.0j]
 [-0.30935922+0.0j]
 [-0.30935922+0.0j]
 [ 0.57452426+0.0j]
 [-0.30935922+0.0j]
 [-0.30935922+0.0j]]

Answer 1

This seems normal by applying the definition of the inversion about average operator which transforms the amplitudes $\alpha_i $ by the formula :

$$ - \alpha_i + 2 \langle\alpha\rangle\,, $$

Apply this to your example (using the above formula for each step to verify your intermediary steps) and you should retrieve this numbers. This operator is periodic. After the maximum number of operations giving you the highest probability of the desired state to be measured, you will reset the amplitudes if you apply it another time, meaning you have to reapply again to amplify its amplitude the same way.

Answer 2

Your algorithm is working exactly as it should. Grover's algorithm does not monotonically converge. Instead, the better intuition is to think about a point moving around a circle at constant speed. If you choose the right moment, you will be close to a particular point, but leave it longer and it'll go past, although it will eventually come back again.

Indeed, the amplitude of the marked item after $n$ iterations should be $$ \sin((2n+1)\theta) $$ where $\sin\theta=\frac{1}{\sqrt{N}}$ and $N$ is the number of items being search over. In this case, $N=8$, so $\sin\theta=\frac{1}{\sqrt{8}}$ and the first few values of the marked amplitude will be $$ 0.883883, 0.972272, 0.574524, -0.110485,\ldots $$ Thus, you stop after 2 iterations and have extremely high probability of finding the marked item, but you continue to the third iteration and it gets worse again. The point is that since you know $\sin\theta$, you can predict the value of $n$ such that $\sin((2n+1)\theta)\approx 1$, and know when to stop. (Note that this is a better and better approximation as $N$ gets larger.) You measure the outcome, and get some answer $x$. You then input that to the oracle again and test if it is the correct answer. If not, you start the whole algorithm again, and repeat until success.