In the last post we devised a procedure that takes a quantum algorithm and turns it into an equivalent automata based automatic one.

The trick was to replace quantum states (and operations on them) with weighted automata.

We applied the procedure to Deutsch-Jozsa algorithm, today I want to see if the same recipe works for Grover's algorithm (see also  [1] ch.6 and  [2] ch.11.4).

As usual we review the quantum algorithm first (with Python code).

Suppose that

• we have a set of \(N=2^n\) items indexed by \(x = 0 \dots N-1\)

• we want to find the items that satisfy some search criteria encoded into a black box function \(f\) that takes an index \(x\) and checks if the x-th item satisfies it

  \[f(x) = \begin{cases} 1 & \text{x-th satisfies criteria} \\ 0 & \text{otherwise} \\ \end{cases}\]

• we suppose that there is a unique \(\omega\) such that \(f(\omega) = 1\)

Classically to find the solution we would need to evaluate \(f\) on average \(N/2\) times (\(O(N)\)).

Grover's algorithm instead needs just \(O(\sqrt{N})\) evaluation thus is quadratically more efficient.

Before we delve into the algorithm itself let's introduce the oracle: instead of evaluating \(f\) directly we define the unitary operator \(U_f\)¹ (\(x_1 \dots x_n\) is the binary string corresponding to \(x\))

\[\ket{x_1 \dots x_n \; y} \xrightarrow{U_f} \ket{x_1 \dots x_n \; y \oplus f(x_1 \dots x_n)}\]

the effect of \(U_f\) is to negate \(y\) if \(x\) is a solution

\[\ket{x_1 \dots x_n \; y} \xrightarrow{U_f} \begin{cases} \ket{x_1 \dots x_n \; y} & f(x) = 0 \\ \ket{x_1 \dots x_n \; \bar{y}} & f(x) = 1 \\ \end{cases}\]

We can simplify the oracle by setting \(\ket{y} = (\ket{0} - \ket{1})/\sqrt{2}\)

\[\ket{x_1 \dots x_n} \frac{\ket{0} - \ket{1}}{\sqrt{2}} \xrightarrow{U_f} (-1)^{f(x)} \ket{x_1 \dots x_n} \frac{\ket{0} - \ket{1}}{\sqrt{2}}\]

and omitting the last qubit since it's not affected

\[\ket{x_1 \dots x_n} \xrightarrow{U_f} (-1)^{f(x)} \ket{x_1 \dots x_n}\]

also if there is a unique solution \(\omega\) we can rewrite \(U_f\) as

in Python (as usual all the code you'll need for this post is here)

import numpy as np
from wfa3 import scalar

zero, one = np.array([1.0,0]), np.array([0,1.0])

# ω = '01101'
ω = np.einsum("i,j,k,l,m->ijklm", zero, one, one, zero, one)

# Setup the oracle
U = lambda x: x - 2 * ω * scalar(ω, x) # (I - 2·|ω><ω|) |x>

Finally Grover's algorithm:

1. initialize the quantum register²

   \[\ket{\psi_0} = H^{\otimes n} \ket{0}^{\otimes n}\]

   H = (1/np.sqrt(2)) * np.array([[1, 1],
                                  [1,-1]])
   
   ψ0 = np.einsum("i,j,k,l,m->ijklm", zero, zero, zero, zero, zero)
   ψ0 = np.einsum("ij,kl,mn,op,qr,ikmoq->jlnpr", H, H, H, H, H, ψ0)
   

2. perform the Grover iteration \(r(N) = \frac{\pi}{4} \sqrt{N}\) times³

   • apply the oracle \(U_f\)
   • apply the Hadamard transform \(H^{\otimes n}\)

   • perform a conditional phase shift

     \[\ket{x_1 \dots x_n} \xrightarrow{S} -(-1)^{\delta_{x_1 0} \cdots \delta_{x_n 0}} \ket{x_1 \dots x_n}\]

     \begin{align*} S &= -\sum_{x_1 \dots x_n=0,1} (-1)^{\delta_{x,0}} \ket{x_1 \dots x_n} \bra{x_1 \dots x_n} \\ &= -\sum_{x_1 \dots x_n=0,1} (-1)^{\delta_{x,0}} \ket{x_1 \dots x_n} \bra{x_1 \dots x_n} \pm \ket{0 \dots 0}\bra{0 \dots 0} \\ &= -I + 2 \ket{0 \dots 0}\bra{0 \dots 0} \end{align*}
   • apply the Hadamard transform \(H^{\otimes n}\) again

   the full Grover iteration then amounts to

   \begin{align*} G &= H^{\otimes n} \; S \; H^{\otimes n} \; U_f \\ &= H^{\otimes n} \; (2\ket{0 \dots 0}\bra{0 \dots 0} - I) \; H^{\otimes n} \; U_f \\ &= (2\ket{\psi_0}\bra{\psi_0} - I) \; U_f \end{align*}

   at the end the registry is

   \[\ket{\psi} = G^{r(N)} \ket{\psi_0}\]

   r = int(np.round(np.pi * np.sqrt(2**5) / 4))
   
   ψ = ψ0.copy()
   for _ in range(r):
       ψ = U(ψ)                        # use the oracle
       ψ = 2 * ψ0 * scalar(ψ0, ψ) - ψ  # 2·|ψ0><ψ0| - I
   

3. perform a measurement on \(\ket{\psi}\) in the computational basis

   from itertools import product
   
   for x in product(*[[0,1]]*5):
       p = np.conj(ψ[x]) * ψ[x]
   
       if np.isclose(p, 0, atol=0.01):
           continue
   
       xs = ''.join([str(i) for i in x])
       print(f"search result: {xs} with {100*p:.1f}% probability")
   

   search result: 01101 with 99.9% probability
   

At first the whole thing seems silly: we are searching something that we already know.

The point is that the value of \(\omega\) is hidden inside the black-box and is not directly available to the search algorithm.

We are in an Alice/Bob kind of situation: Bob has some information that he makes available only through the oracle, Alice's task is to recover it in the most efficient way.

Another way to look at the oracle is that recognizing is not the same as knowing. If you are presented with an object you can easily check that it satisfies the search criteria. Actually finding it, on the other hand, is much harder (think of combinatorial search problems like the travelling salesman).

At this point we are ready to automatize Grover's algorithm: for each step we replace items on the left column with items on the right

before we start we setup the oracle

from wfa3 import norm, make_automaton_from_words, hilbert_brzozowski, normalize_wfa

n = 5

# Setup the oracle
u = np.zeros([2]*(2*n))
for x in product(*[[0,1]]*n):
    u[x + x] = -1 if ''.join([str(i) for i in x]) == '01101' else 1
u = u / norm(u)

U = hilbert_brzozowski(u)

and since we are at it we prepare the conditional shift operator

s = np.zeros([2]*(2*n))
for x in product(*[[0,1]]*n):
    s[x + x] = 1 if ''.join([str(i) for i in x]) == '00000' else -1
s = s / norm(s)

S = hilbert_brzozowski(s)

then the automatic Grover algorithm is

1. initialize the automaton

   φ0 = make_automaton_from_words('0'*n).rebase(H)
   φ0.diag('phi0')
   

   

   this automaton represents all the possible solutions to the search problem

2. perform the Grover iteration \(r(N) = \frac{\pi}{4} \sqrt{N}\) times

   • apply the oracle \(U_f\)
   • apply the Hadamard transform
   • perform a conditional phase shift
   • apply the Hadamard transform

   r = int(np.round(np.pi * np.sqrt(2**n) / 4))
   
   φ = φ0.copy()
   for i in range(r):
       φ = U.D_wfa(φ, n)        # use the oracle
   
       φ = normalize_wfa(φ, n)  # you can omit normalization but then you'll get uglier pictures
   
       if i == 0:
           φ.diag(f'phi')
   
       φ = φ.rebase(H)          # Hadamard
       φ = S.D_wfa(φ, n)        # conditional phase shift
       φ = φ.rebase(H)          # Hadamard
   
       φ = normalize_wfa(φ, n)
   

   now look at what happens when we apply the oracle \(U_f\) for the first time to \(\varphi_0\)

   

   \(U_f\) has changed the structure of the automaton in the same way as in the Deutsch-Jozsa algorithm: a new path for the solution has appeared.

   The subsequent operations do not change the structure but only update the weights: the conditional phase shift transfers weight to the word marked by the oracle

3. take a sample from \(\varphi\)

   for x in product(*[[0,1]]*n):
       xs = ''.join([str(i) for i in x])
       p = np.conj(φ(xs)) * φ(xs)
   
       if np.isclose(p, 0, atol=0.01):
           continue
   
       print(f"search result: {xs} with {100*p:.1f}% probability")
   

   search result: 01101 with 99.9% probability
   

   consistent with the quantum case.