Last time we devised a recipe to turn quantum algorithms into automata based ones and used it on the Deutsch-Jozsa algorithm. But why stop there? Let’s automatize Grover’s algorithm too.

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)1 ((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

\begin{align*} U_f &= \sum_{x_1 \dots x_n=0,1} (-1)^{{f(x)} \ket{x_1 \dots x_n} \bra{x_1 \dots x_n} \ &= \sum_{x_1 \dots x_n=0,1} (-1)}{f(x)} \ket{x_1 \dots x_n} \bra{x_1 \dots x_n} \pm \ket{\omega}\bra{\omega} \ &= I - 2 \ket{\omega}\bra{\omega} \end{align*}

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])

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

U = lambda x: x - 2 * ω * scalar(ω, x)

Finally Grover’s algorithm:

initialize the quantum register2

[\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)

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

• 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(ψ)                            ψ = 2 * ψ0 * scalar(ψ0, ψ) - ψ

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

quantum Grover	automatic Grover
quantum states	automata
Hadamard transform	alphabet remapping
(U_f) operator	string replacement
quantum measurement	sampling

before we start we setup the oracle

from wfa3 import norm, make_automaton_from_words, hilbert_brzozowski, normalize_wfa

n = 5

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

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 1.

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)
φ = normalize_wfa(φ, n)
if i == 0:
φ.diag(f'phi')

φ = φ.rebase(H)              φ = S.D_wfa(φ, n)            φ = φ.rebase(H)
φ = 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 1.

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.

---

[1]

[2]

---

Footnotes:

1

it is the same as in Deutsch-Jozsa algorithm

2

it’s the superposition of all possible solutions

3

for the complete treatment see [1] 6.1.3