Following the first quantum algorithms we have already reviewed (Deutsch’s algorithm and the DeutschJozsa algorithm — it's better to read these articles first, otherwise all the things mentioned here will be quite vague), let’s analyze Grover’s algorithm for the unstructured quantum search. An American mathematician Lov Grover formulated the algorithm in 1996 (it was long after the model of quantum computations had become popular). The algorithm uses a feature of quantum interference in order to solve an extremely demanding task of searching the value of some parameter, at which a defined function returns certain results.
The given algorithm does not show the exponential advantage of the task in contrast to the classical computational model, but the advantage (quadratic) is quite essential for big values. However, it’s a general algorithm for solving quite a generalized problem. It has been proved that we cannot get a better result within the limits of the model of quantum computations. It is possible in more particular algorithms.
It is another post in the series of articles about the model of quantum computations. That’s why you’d better read the previous articles, as they will help you to understand all the things here. If you came across my notes for the first time, you should refer to the following articles:

A Few Words About Reversible Computing

How to Implement Deutsch’s Algorithm in Haskell

Quantum Circuit Design: Methods and Techniques

How to Implement the DeutschJozsa Algorithm in Haskell
If interested, you are most welcome under the cut.
Semiformal Description of the Algorithm
The task is formulated the following way. Let’s assume that there’s a binary function from n binary arguments, which accepts the value of 1 at one of them only (as for the other 2^{n} – 1 values, it accepts 0). It is required to find the value of input arguments, when there’s one function only, and it is provided in the form of an oracle.
In the classical variant, we should consider 2^{n}/2 variants of input values. At best, we will be able to find the “lucky number” on the first try. At the worst, we will have to go through all 2^{n}/2 variants. Grover’s algorithm allows to do it in π/4 √(2^{n}) calls to the oracle. It’s obvious that when the n number of input qubits growth, the difference in performance becomes principal. But as we have mentioned before, there’s no superpolynomial performance growth.
By definition, it is the task of the unstructured search. Suppose there is an unordered data set (a «haystack»), and it is required to find one element satisfying a specific requirement (there’s just one such element — a “needle”). Since the classical search is an alternative here, Grover’s algorithm will come in handy. For example, if we take a look at the analogy with the database of motorists containing their names that are ordered alphabetically (let there be an unambiguous correspondence between the name and the car number), the ordered search is the search of a car number by name. If there’s no index, we can do it by means of dichotomy. The unordered search is the inverse problem: the search of the name by the car number. In the given analogy, the oracle is the function (transformed accordingly) that returns the car number by name. Thus, this task comes down to Grover’s task with the help of coding names into the binary representation. The function returns the answer to the following question: “Is N motorist the owner of X car?”, where N is the input parameter of the function. X is the algorithm parameter, as though it is “sewed” inside the oracle when the latter is being built. So, this oracle returns 1 (“yes”) value only for the name that has the sought car number next to it.
Grover’s algorithm contains the following steps:

Initialize the initial state. It is necessary to prepare the equally probable superposition of states of all input qubits. We can do it by applying a corresponding Hadamard gate that is equal to the tensor product of n unary Hadamard gates multiplied by each other.

Apply Grover’s iteration. The given iteration lies in the consistent application of two gates, the oracle and the socalled Grover diffusion operator. We will review both of them below. We perform this iteration √(2^{n}) times.

Perform the measurement. After we have performed Grover’s iteration the sufficient number of times, it is necessary to measure the input register of qubits. With a really high probability, the measured value will point to the sought parameter. If it’s necessary to increase the answer authenticity, the algorithm is run several times and the composite probability of the correct answer will be calculated.
It goes without saying that the quantum circuit of an algorithm depends on the size of input data, since this depends directly on the number of applications of Grover’s iterations. We can depict the circuit, like in the following picture:
Let’s pay attention to the oracle and the diffusion gate, as they are the very heart of the algorithm. The oracle should be built in not a standard way. Remembering the process of transforming a classical binary function into an oracle, we can understand that accepting two quantum registers (x, y) as an input, the latter should return (x, y ⊕ f(x)) pair as an output. But in case with Grover’s algorithm, the oracle should invert the phase of the quantum state, that corresponds to the sought x value. (Reminding you that we are searching x by the defined f(x) value. i.e. we are solving an inverse problem.) This means that the oracle should return (1)^{f(x)}x> as an output. The reason why phase coefficient is (1)^{f(x)} is because f(x) = 1 if only if the function accepts the sought x value at the input. In this case, the phase coefficient will become equal to 1. In other words, U_{w} oracle functions the following way:

U_{w}w> = w>

U_{w}x> = x>, x ≠ w
The diffusion operator, in its turn, represents a combination of three gates, two of which are standard and they are multiqubit Hadamard gates. There’s a special gate between them that performs the flip of qubits with regard to the average value. Its representation in the form of (20^{n}><0^{n} — I_{n}) analytic formula makes the main point a bit obscure. Its matrix representation is simple: there’s value of 1 in the topleft corner of the matrix, and there are 0 in other places.
But we can achieve the same by changing 0^{n}> to +^{n}> in the mentioned above formula. In this case, we don’t have to apply Hadamard gates before and after applying this special gate. Thus, Grover diffusion operator will look like this:
(2+^{n}><+^{n} — I_{n})
To put it differently, with such a form of the diffusion operator, the quantum circuit of Grover’s algorithm will look like this:
Implementation
So that all of this would not sound unsubstantiated, we can review this algorithm in its implementation in Haskell programming language. We will use the same set of functions that has been applied more than once in this series of articles. This set will allow us to look inside the depths of the algorithm. (By the way, Quipper would not allow us to do the same, as it provides really highlevel means of generating quantum circuits).
Since we should design an oracle in the form of a gate, we should deal with the quantum circuit schematics before writing the code. For a change (and in order to solidify our knowledge of quantum system technologies), we can take a look at a function of three variables. So, the oracle will accept 3 qubits as an input (of course, there will be 3 qubits as an output). Let’s consider the following task:
f(x_{1}, x_{2}, x_{3}) = x_{1} & x_{2} & x_{3}
This function accepts value of 1 only in one of the eight variants of input values. That's exactly what we need for Grover’s algorithm. There’s no difference, what function of the eight possible ones we will take a look at, but it is simpler from the technical point of view. As usual, we can use the following table to build the oracle:
X_{1}  X_{2}  X_{3}  f(X)  (1)^{f(X)} 
0  0  0  0  1 
0  0  1  0  1 
0  1  0  0  1 
0  1  1  0  1 
1  0  0  0  1 
1  0  1  0  1 
1  1  0  0  1 
1  1  1  1  1 
The following 8×8 matrix corresponds to the table:
Since we are going to implement more serious things than those in the previous articles, we should add some service functions that are useful for building oracles and gates. Particularly, we should add a function that has been already used in the tensor product operator:
groups :: Int > [a] > [[a]]
groups i s  null s = []
 otherwise = let (h, t) = splitAt i s
in h : groups i t
Now, let’s get down to defining new gates and operators above them. According to the description of Grover’s algorithm, we are going to need the operator for subtracting matrices one from another, and also a function for creating gates for several qubits on the basis of a gate for one qubit. We will need this function at least for the tensor product of I and H gates, so that we could apply them to quantum registers containing several qubits.
To begin with, let’s define a simple operator for obtaining the difference of matrices:
(<>) :: Num a => Matrix a > Matrix a > Matrix a
m1 <> m2 = zipWith (zipWith ()) m1 m2
Let’s also define the same operator for adding matrices:
(<+>) :: Num a => Matrix a > Matrix a > Matrix a
m1 <+> m2 = zipWith (zipWith (+)) m1 m2
We should take into account the fact that the developer should control the dimensions of matrices and vectors.
It’s time to get down to the functions for generating gates. We have implemented functions for representing gates that process one or two qubits. With the help of the tensor product operator and entangle function, we can create functions for representing gates that can process an arbitrary number of qubits. For example, that’s how a generalized function that transforms the given onequbit gate processing a defined number of qubits, looks like:
gateN :: Matrix (Complex Double)
> Int > Matrix (Complex Double)
gateN g n = foldl1 (<++>) $ replicate n g
It’s really smart and simple. With the help of a standard replicate function, we can create a list from a specified number of onequbit gates. The developer should control the process, so that onequbit gates would be passed here. Nevertheless, the function will operate with multiqubit gates as well. Then, we should fold the list by means of <++> tensor product operator. That's how this function is used:
gateIn :: Int > Matrix (Complex Double)
gateIn = gateN gateI
gateHn :: Int > Matrix (Complex Double)
gateHn = gateN gateH
As you can see, the first function forms an identical transformation for the defined number of qubits, while the second one forms the Hadamard transformation for the defined number of qubits as well. The two multiqubit gates are really important in the model of quantum computations. We often use them in various quantum circuits.
Finally, we can implement Grover’s algorithm itself. Let’s start with the oracle:
oracle :: Matrix (Complex Double)
oracle = matrixToComplex [[1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 1]]
There’s nothing difficult about it. We should just code the matrix. A cautious reader has definitely extended the set of functions and defined the higher order function for generating such oracles (and many others). But we will do it in the old way.
Let’s implement a function for representing Grover diffusion operator. Considering our considerations as for the basis change (use a qubit in +> quantum state instead of 0>), its definition will be really simple:
diffusion :: Matrix (Complex Double)
diffusion = 2 <*> (qubitPlus3 >< qubitPlus3)
<> gateIn 3
where
qubitPlus3 = toVector $
foldl1 entangle $
replicate 3 qubitPlus
Take three qubits in +> (qubitPlus) state, make a list from them and fold the list with the help of entangle function. Then, transform the obtained quantum resister into the vector representation. That's how we get the local qubitPlus3 definition.
After that, multiply qubitPlus3 quantum register by itself. As a result, we will obtain +><+ matrix and multiply it by 2. Subtract an identity matrix, prepared for three qubits, from the result (i.e. it’s of 8×8 size). Using the implemented above operators, writing such functions becomes a real joy.
It’s high time to implement Grover’s algorithm itself. Here’s the definition of the function for three qubits:
grover :: Matrix (Complex Double) > IO String
grover f = initial > gateHn 3
> f > diffusion
> f > diffusion
>>> (measure . fromVector 3)
where
initial = toVector $
foldr entangle qubitZero $
replicate 2 qubitZero
In the local initial, we will prepare the initial state that is equal to 000>. After that, we will direct it to the Hadamard gate for three qubits. As a result, we will get an equally probable superposition of eight quantum states that can be on the three qubits. Run Grover’s cycle twice, as 2 ≈ √23. Finally, carry out the measurement.
But what will happen, if we add the third Grover’s cycle to this quantum circuit? It’s really simple. The results will become worse. Why? Think about it yourself.
Since Grover’s algorithm is probabilistic, it gives out a correct answer only with a really high probability. This means that when running grover function, we will sometimes get a wrong answer. Therefore, it is suggested to estimate the probability of getting the correct answer. Let’s implement such function:
main f n = do l < replicateM n $ grover f
return $
map (length &&& head) $ group $ sort l
It applies Grover’s algorithm for the defined oracle a specified number of times. Its operation result is a histogram of the algorithm results (it’s a list of pairs of (occurrence rate, result)) form). Then, we ran this function a million times and built the following chart:
The correct answer was obtained in about 94.5 % of cases. Other results were of about 0.78 % rate. It’s quite enough to run Grover’s algorithm three times and choose the result that has been repeated at least twice.
A cautious reader might have wondered, what would happen, if the oracle returned 1 phase at a several input data, not at only one. Strange as it may seem, Grover’s algorithm operates in this case as well (looking at the order of matrix multiplication, there’s nothing strange about it at all). But to find the only one correct answer out of many, we need much more iterations.
Suppose l is the number of values of input parameters, at which the function accepts the value of 1 (the oracle returns 1 phase). In that case, to find the correct answer, we will need √2n/l Grover’s iterations. We can demonstrate it with the help of the following code:
oracle' :: Matrix (Complex Double)
oracle' = matrixToComplex [[1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 1]]
If we run function main with this oracle and allow it to build a histogram at a million of runs, we will see the following chart:
What’s that? Have the correct answers got the lowest frequency? Oh, it’s all right, as grover function, called from main function, executes two Grover’s iterations, and we need just one iteration for the oracle with three correct answers. As soon as more iterations are executed, the situation turns upside down (as the flip with regard to the average). But in this very case, one iteration is not enough either, as frequency probabilities of correct and wrong answers will be quite close to each other (it’s quite reasonable, as just one iteration has been executed).
According to this example, when using Grover’s algorithm, the developer should always control the number of iterations and not allow the situation turn the wrong way.
That's when we can stop describing Grover’s algorithm. The way it operates, and the ways we can implement it in the programming language should be quite clear by now.
Summary
An interested reader can refer to the source code of the module. You are also most welcome to discuss the algorithm and the model of quantum computations in your comments below the post. I am going to dwell on more interesting algorithms in my future articles.
{# OPTIONS_HADDOCK prune, ignoreexports #}
{}
{  This module contains functions, which help to implement Grover's quantum
search algorithm
Author: Roman Dushkin
Project: Quantum Computations and Functional Programming
}
{}
module Grover
(
grover,
main
)
where
{[ IMPORT SECTION ]}
import Control.Arrow ((&&&))
import Control.Monad (replicateM)
import Data.Complex (Complex, realPart)
import Data.Function (on)
import Data.List (sort, group)
import Circuit
import Gate
import Qubit
{[ FUNCTIONS ]}
  Prepared oracle to demonstrate Grover's algorithm for
 three qubits.
oracle :: Matrix (Complex Double)
oracle = matrixToComplex [[1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 1]]
  One more oracle, which corresponds a function returning the value of 1
 for multiple (three) values.
oracle' :: Matrix (Complex Double)
oracle' = matrixToComplex [[1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 1]]
  Function which implements the diffusion gate.
diffusion :: Matrix (Complex Double)
diffusion = 2 <*> (qubitPlus3 >< qubitPlus3) <> gateIn 3
where
qubitPlus3 = toVector $ foldl1 entangle $ replicate 3 qubitPlus
  The main function of the module, which demonstrates Grover's algorithm for
 three qubits.
grover :: Matrix (Complex Double) > IO String
grover f = initial > gateHn 3
> f > diffusion
> f > diffusion
>>> (measure . fromVector 3)
where
initial = toVector $ foldr entangle qubitZero $ replicate 2 qubitZero
  The main function of the module, which builds the histogram of the results of
 measurement of the quantum register by running Grover's algorithm
 predefined number of times.
main f n = do l < replicateM n $ grover f
return $ map (length &&& head) $ group $ sort l
{[ END OF MODULE ]}