$\newcommand{\I}{\mathrm{i}} \newcommand{\E}{\mathrm{e}} \newcommand{\D}{\mathop{}\!\mathrm{d}} \DeclareMathOperator{\Tr}{Tr} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\ket}[1]{|{#1}\rangle} \newcommand{\braket}[1]{\langle{#1}\rangle} \newcommand{\bbraket}[1]{\langle\!\langle{#1}\rangle\!\rangle} \newcommand{\bm}{\boldsymbol}$

Lectures on advanced quantum mechanics

Quantum Algorithms

Since physical laws are inherently quantum, Feynman proposed, in a famous talk given in 1981, to use “quantum computers” to efficiently solve physical problems¹. A specific model of quantum computation was devised by Deutsch, who, in 1989, found a set of universal unitary gates, extending to the quantum domain the universality of the classical Turing machine². In fact, the exponential growth of the Hilbert space, prevents a quantum system to be simulated by a classic computer, with polynomial classical ressources. In Feynman’s view, a quantum computer should solve a quantum system with resources proportional to the size (number of degrees of freedom) of the physical system:

…the physical world is quantum mechanical, and therefore the proper problem is the simulation of quantum physics¹

The discovery by Shor (1994)³ of a quantum algorithm that can factor large integers in a time which is polynomial in the number of digits, triggered interest in the whole new domain of quantum computing. In this section we discuss quantum algorithms, using the circuit model of computation, which show some of the basic ingredients allowing quantum computers to solve hard classical problems.

Grover search

Searching for an item in a unstructured data base is a $O(N)$ task, if $N$ is the number of items. However, Grover (1996) found that using a quantum computer, finding an item among $N$ can be done in $O(\sqrt{N})$ steps. The square speedup results from quantum interference: the algorithm uses amplitude amplification to reinforce the probability of the searched item.

Let start with the classical search algorithm. Instead of searching for the data base item, we label each item with an integer $x=1, \ldots, N$. To certify that $x$ correspond to the label that we are searching for, $x_0$ say, we define an oracle, a boolean function which is $f(x) = 1$ if $x = x_0$ and zero otherwise. Therefore, the search can be implemented as in the (psudo)code:

input: N and f(x)
for x in range(N):
    if f(x) == 1:
        return x

it uses a loop of length $N$, meaning that we need about $N$ oracle evaluations to find $x_0$.

To search using quantum states we need two qubit registers, one addressing the data, the labels are encoded in a quantum state $\ket{x}$, and a single qubit register, the ancilla to apply the oracle. The quantum algorithm replaces the classical oracle by a unitary operator, to test if the data state corresponds to the searched item,

$$ O\ket{x} \otimes \ket{a} = \ket{x} \otimes \ket{a \oplus f(x)}, \quad \ket{x} \in \mathcal{H}, \;\ket{a} \in \mathcal{H}_A $$

($\mathcal{H}$ is the data space, and $\mathcal{H}_A$ is the ancilla space) which modifies the ancilla qubit according to the state of the data register (control) EX:

$$ O \ket{x} \ket{-} = (-1)^{f(x)} \ket{x}\ket{-}, $$

it marks the $x_0$ state with a phase $\pi$. It is just this amplitude that, in the superposition state, is amplified. (See an example of oracle below.)

The quantum amplification is effectively obtained by applying, after the oracle, the reflection operator

$$ R = (2\ket{s} \bra{s} - 1_N) \otimes 1_2 $$

where $\ket{s} = H^{\otimes N} \ket{0}$ is the uniform amplitude global superposition state:

$$ \ket{s} = \frac{1}{2^{N/2}} \sum_{n=0}^{2^N - 1} \ket{n}. $$

The Grover operator is $G = RO$ is the composition of the two reflections acting on the data register (the ancilla state do not change, it always $\ket{-}$) $O$, implements a reflection about an axis orthogonal to $\ket{x_0}$, and $R$, another reflection about the axis $\ket{s}$ of the uniform state EX. Two reflections are equivalent to a rotation: at each iteration of $G$ rotates the state to “approach” the state labeled by $x_0$; the rotation angle depends on the value of $\braket{x_0|s}$, the angle between the searched item and the uniform state. The optimal number of iterations is EX:

$$ T = \frac{\pi}{2} \sqrt{N}; $$

at this “time” the probability to get $x = x_0$ approaches 1.

Quantum algorithm:

input: $N$, $f(x)$
intialize the data and ancilla registers in the state $\ket{0}^{\otimes N} \otimes \ket{-}$
apply $H^{\otimes N}$ to the data register, the state is $\ket{\psi} = \ket{s}\ket{-}$
iterate $G$:

psi
for t in range(T): # T approx pi/2 sqrt(N)
    psi = G * psi # apply G to the total state

measure the state register to get the answer

Note that the oracle is applied $T \sim O(\sqrt{N})$, leading to a quadratic speedup with respect to the classical algorithm.⁴

Example of a circuit implementing the oracle operator to identify the item $x_0 = 6 = 110$: the output state is marked by a $-1 = \E^{\I \pi}$ phase.

Fourier transform and phase estimation

The classical fourier transform has an inverse: applying the transform to a function and then applying its inverse, lead to the identity; the fact that the inverse is related to the complex conjugate of the direct transformation, naturally leads to the observation that the fourier transformation can be seen as a unitary operator. The quantum algorithm is simply an implementation of this unitary operator acting on a space of dimension $2^N$, as a circuit of one and two qubits gates.

Let us start with the fourier transform of the basis states. We denote a basis state of the $N$ qubits hilbert space,

$$ \ket{s} \in \mathcal{H}^{\otimes N} $$

by the binary string

$$ s = s_1 \ldots s_{N} = \sum_{n=1}^{N} s_n 2^{N-n},\; s_n = \{0,1\}. $$

Therefore,

$$\ket{s} = \frac{1}{\sqrt{2^N}} \sum_{k=0}^{2^N-1} \E^{2\pi \I s k/2^N} \ket{k}$$

is the fourier transform of $\ket{s}$. Knowing the transformation of the basis vector we can compute the transformation of any superposition. The basic idea is to transform the unitary operator, here defined as a sum over the basis vectors in fourier space $\ket{k}$, into a product. This can be achieved using the decomposition of $s$ and $k$ into its binary digits:

$$ \begin{align*} \ket{s} &= \frac{1}{\sqrt{2^N}} \sum_{k=0}^{2^N-1} \exp\left( 2\pi \I \frac{s k}{2^N} \right) \ket{k} \\ & = \frac{1}{\sqrt{2^N}} \sum_{k_1=0,1} \ldots \sum_{k_{N}=0,1} \exp\left( 2\pi \I s \sum_{n=1}^{N} \frac{k_n}{2^n} \right) \ket{k_1 \ldots k_{N}} \\ & = \frac{1}{\sqrt{2^N}} \sum_{k_1=0,1} \ldots \sum_{k_{N}=0,1} \bigotimes_{n=1}^{N} \exp\left( 2\pi \I \frac{k_n}{2^n} \right) \ket{k_n} \\ & = \frac{1}{\sqrt{2^N}} \bigotimes_{n=1}^{N} \sum_{k_n=0,1} \exp\left( 2\pi \I \frac{k_n}{2^n} \right) \ket{k_n} \\ & = \frac{1}{\sqrt{2^N}} \bigotimes_{n=1}^{N} \left[ \ket{0} + \exp\left( 2\pi \I \frac{s}{2^n} \right) \ket{1} \right]. \end{align*} $$

The exponent selects the fractional part of $s/2^n$ for $n=1,\ldots,N$:

$$ \{s/2^n\} = 0.s_{N-n+1}\ldots s_N\,, $$

which allows the previous product to be written in the simple form,

$$ \begin{equation} \label{e:prod} \ket{s_1 \ldots s_N} = \frac{1}{2^{N/2}} \left(\ket{0} + \E^{2 \pi \I 0.s_N} \ket{1} \right) \otimes \ldots \otimes \left(\ket{0} + \E^{2 \pi \I 0.s_1 \ldots s_N} \ket{1} \right)\,. \end{equation} $$

Consider the last exponential, each term $s_n/2^n$ of the phase $2\pi 0.s_1\ldots s_N$,

$$ 0.s_1 \ldots s_N = \frac{s_1}{2} + \cdots + \frac{s_N}{2^N} $$

can be obtained by a rotation of an angle $2\pi/2^n$ if $s_n = 1$, otherwise it remains $0$; this operation can be implemented by the rotation matrix

$$ R_n = \begin{pmatrix} 1 & 0 \\ 0 & \E^{2 \pi \I/ 2^n} \end{pmatrix}\,, $$

controlled by $s_n$ and applied to the first qubit state (initially $\ket{s_1}$). We observe that this term depends on the values of all qubits, leading to a decomposition in at least $N$ operators (one qubit or with one control qubit). The first decimal can be simply obtained from the initial $\ket{s_1}$ by the application of an hadamard gate

$$ H \otimes^{N-1} 1_2 \ket{s_1 \ldots s_N} = \left(\ket{0} + \E^{2 \pi \I 0.s_1} \ket{1} \right) \ket{s_2 \ldots s_N}, $$

note $\E^{2 \pi \I s_1/2} = -1$ only if $s_1=1$. Now, applying $R_2$ controlled by the second qubit, we obtain

$$ \left(\ket{0} + \E^{2 \pi \I 0.s_1s_2} \ket{1} \right) \ket{s_2 \ldots s_N}. $$

Iterating this last step $N-1$ times, up to applying $R_N$, we get the last term in the product:

$$ \left(\ket{0} + \E^{2 \pi \I 0.s_1s_2 \ldots s_N} \ket{1} \right) \ket{s_2 \ldots s_N}. $$

For the other qubits we use the same recipe, $N-1$ times for $0.s_2 \ldots s_N$, $N-3$ on $0.s_3 \ldots s_N$, and so on: the total number of steps is then $N(N-1)/2$. We represent the whole algorithm in the “QFT” circuit of the figure.

Note that we need to reverse the order of the qubits to obtain the result; this can be achieved by applying $N/2$ swap gates,

$$ \mathsf{SW} = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}, $$

at the end of the QFT circuit.

The classical algorithm, the fast fourier transform, needs instead $2^N N$ steps, which is exponential in $N$, compared to the $N^2$ quantum algorithm. This speed up in the computation of the fourier transform is at the origin of the Shor algorithm performance, which allows the factorization of large integers in polynomial time.

Phase estimation

The quantum fourier transform may be used to compute the phase $\theta \in (0,1)$, related to the eigenvalue of a unitary $U$ by,

$$U \ket{\theta} = \E^{2 \pi \I \theta} \ket{\theta}$$

where $\ket{\theta}$ is a known eigenvector (of dimension $D$). The circuit (phase estimation algorithm)

implements the first part of the algorithm (Cleve et alia, 1998)⁵. One register, whose input is $\ket{\theta}$, is used to compute the iterates $U^{2^n}$ ($n=0,\ldots,N-1$) which do not modify the eigenstate, and an “index” register of $N$ ancilla qubits, used as controls of the $U$ iterates; the value of $N$ is chosen in order to get the eigenvalue $\bar{\theta} \approx \theta$ with a $N$ bits precision. The state obtained after the circuit application corresponds to the fourier transform of the state EX:

$$ \ket{2^N \theta} \ket{\theta} = \frac{1}{2^{N/2}} \sum_{k=0}^{2^N - 1} \E^{2\pi \I \theta k} \ket{k} \ket{\theta} . $$

Applying the inverse fourier transform

$$ \ket{k} = \mathcal{F^\dagger} \ket{s} = \frac{1}{2^{N/2}} \sum_{s=0}^{2^N - 1} \E^{-2\pi \I s k/2^N} \ket{s}. $$

to this state, one obtains,

$$ \frac{1}{2^N} \sum_{s=0}^{2^N - 1} \left[ \sum_{k=0}^{2^N - 1} \E^{-2\pi \I k (s - 2^N \theta)/2^N} \right] \ket{s}\ket{\theta} \approx \ket{2^N \bar{\theta}}\ket{\theta} . $$

The expression in brackets is picked at $\theta \approx s/2^N = \bar{\theta}$ (where $\bar{\theta}$ is a binary string), therefore, after this inverse transformation, measurement of the ancilla qubits gives us $s$, a $N$ bit approximation of $\theta$.

Eigenvectors and eigenvalues finding

An influential paper of Abrams and Lloyd⁶ adapted the phase estimation algorithm to calculate the eigenvector $\ket{\theta_k}$ of $U$ and its eigenvalue $\E^{2\pi\I \theta_k}$. They assumed an initial state $\ket{\bar{\theta}}$ possessing a substantial overlap with the true eigenvector $\ket{\theta}$; this means that the expansion

$$ \ket{\bar{\theta}} = \sum_k a_k \ket{k}, \quad U \ket{k} = \E^{2\pi\I \theta_k} \ket{k}, $$

contains essentially order one amplitudes

$$ a_k = \braket{\theta_k|\bar{\theta}}. $$

The idea is to define a circuit using a set of $N$ ancilla qubits (the index register) such that their measure leads with probability

$$ p_k = |\braket{\theta_k | \bar{\theta}}|^2 = \big| a_k \big|^2 $$

to the correct eigenvector $\ket{\theta_k}$ and its corresponding eigenvalue $\theta_k$. Consider the circuit,

the initial state is

$$ \ket{\psi} = \ket{0}^{\otimes N} \ket{\theta}, $$

the hadamard gates creates the superposition state,

$$ \ket{\psi} = \frac{1}{2^{N/2}} \sum_{s=0}^{2^N-1} \ket{s} \ket{\theta} = \frac{1}{2^{N/2}} \sum_{s=0}^{2^N-1} \ket{s} \sum_k a_k \ket{\theta_k} , $$

where, in the last step, we expanded $\theta$ in the $U$ basis $\ket{\theta_k}$. We apply next the control-$U$ gates: the $s$ iteration gives a factor $\exp(2\pi\I s \theta_k$ for each eigenvector,

$$ \ket{\psi} = \frac{1}{2^{N/2}} \sum_k \sum_{s=0}^{2^N-1} a_k \E^{2\pi\I s \theta_k} \ket{s} \ket{\theta_k} . $$

We observe that the sum over $s$ is the index register fourier transform of the set of eigenvalues $\theta_k$; we then apply the corresponding inverse fourier transform:

$$ \ket{\psi} = \sum_k a_k \sum_{t} \Big(\frac{1}{2^{N}} \sum_{s=0}^{2^N-1} \E^{2\pi\I s \theta_k} \E^{-2\pi\I st/2^N} \Big) \ket{t} \ket{\theta_k} = \sum_k\sum_t a_k \delta(2^N \theta_k, t)\ket{t} \ket{\theta_k}, $$

which leads, using the Konecker delta, to the state

$$ \ket{\psi} = \sum_k a_k \ket{2^N\bar{\theta}_k} \ket{\theta_k} $$

with, similarly to the phase estimation case, $\bar{\theta}_k$ the $N$-bits approximation of the actual eigenvalue $\theta_k$. Therefore, a measure of the index register gives the eigenvalue

$$ \bar{\theta_\bar{k}} \approx \theta_\bar{k} $$

for some $\bar{k}$, with probability $|a_\bar{k}|^2$ (supposed to be of order one for the initial guess $\ket{\theta}$), and the after-measure state is projected to the corresponding eigenvector $\ket{2^N \theta_{\bar{k}}}$. We may iterate the algorithm to obtain other eigenvalues and eigenvectors.

Exercises

Compute the final state after the application of the first part of the phase estimation algorithm (see the circuit figure). Investigate the performance of the algorithm.
Analyze the circuit,

where $H$ is the hadamard gate,
$$ Z_\varphi = \mathrm{diag} \left(1, \E^{\I \varphi} \right) $$
is a phase gate, and $L$ an integer, and show that it transforms the initial state $\ket{0} \otimes \ket{u_\theta}$, $\ket{u_\theta}$ is an eigenvector of $U$, into $\ket{L, \varphi, \theta} \otimes \ket{u_\theta}$,
$$ \ket{L,\varphi,\theta} = \frac{1}{2} \left( 1 + \E^{2\pi \I L \theta + \I \varphi} \right) \ket{0} + \frac{1}{2} \left( 1 - \E^{2\pi \I L \theta + \I \varphi} \right) \ket{1} $$
Compute the probabilities of measuring the ancilla qubit to be in state $0$ or $1$. (Hint $P_{L,\varphi}(0, \theta) = |\braket{0|L,\varphi, \theta}|^2$)

Show that with a choice of $\varphi$ one can accurately estimate
$$\cos(2\pi L \theta) \; \text{and} \; \sin(2\pi L \theta)$$
from which one finally gets an approximate value of the phase $\theta$.

Practicum

Quantum Fourier Transform

We follow the QFT circuit to write the corresponding algorithm. To compute the quantum Fourier transform of a bit string $x=0, \ldots, 2^{N}-1$ needs $N$ qubits, and can be written in the form of two loops because of the recursive nature of the product formula $\eqref{e:prod}$

Loop $n$ over the $N$ qubits
Apply Hadamard and loop $k \ge 2$ over the remaining $N-n+1$ qubits to apply the control phase gates of decreasing angles $2\pi/2^k$
$$ \mathsf{CPhase}(k) = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & \E^{2\I \pi/2^k} \end{pmatrix}. $$
End both loops.
Apply the swap gates to reverse the qubit order.
Choose a state product of up and down qubits and apply to it steps 1-3. Verify the result; you may compute, for instance, $\braket{y|x}$ if $x$ is the input state and $y$ the output one.

Show that the result of applying the above circuit to an arbitrary state $\ket{x}$ is equivalent to apply the unitary operator

$$ \braket{x|\mathrm{QFT}|y} = \exp\big(2\I\pi xy/2^N\big). $$

Modify the code to compute also the inverse Fourier transform IQFT; verify that applying QFT, followed by IQFT, you recover in initial state.

Phase estimation for the cluster model

We consider a three spin system, forming a triangle, whose Hamiltonian is

$$ H_C(J) = -J\big( X \otimes Z \otimes Z + Z \otimes X \otimes Z + Z \otimes Z \otimes X \big), $$

called a “cluster” Hamiltonian, used in some models of quantum computing (see the section about the cluster model).

One eigenvector $\ket{C}$ of $H_C$ can be constructed by applying the control-phase gate to each edge $e \in E$ of the triangle graph $(V,E)$:

$$ \ket{C} = \prod_{e \in E} \mathsf{CZ}(e), \quad e = \{ (0,1), (0,2), (1,2)\}, $$

(note that the terms in the product commute). Note that the parity operator $P=X \otimes X \otimes X$ commutes with the Hamiltonian. Our goal is to find the eigenvalue $\E^{2\I \theta}$ ($theta$ is proportional to $J$), of the cluster operator

$$ U(J) = \E^{-\I H_C(J)}, \quad P U(J) P = U(J), $$

using the phase estimation algorithm.

Write a python code to compute theta to four significant binary (fraction) digits; use your code for (i) the cluster state $\ket{C}$, and (ii) its parity related state $\ket{\bar{C}}$ ($X$-cluster state) obtained by application of $P$. The “QPE” algorithm is, schematically,

Create a control register with $NC$ qubits in the up state, and apply Hadamard gates to each one.
Create a target register in the initial eigenvector state (cluster, $X$-cluster, for instance).
Define a controlled $U$ and apply it $2^n$ times to the target register, for $n=0,ldots,NC$, using the corresponding qubits of the control register.
Apply the inverse Fourier transform to the control qubits, using the previous code but now applied only to the control register (the target register contains the initial eigenvector).
Measure the control register qubits to obtain the digits of $2^N\theta$ in binary form (note that you must reverse the binary sequence). Verity the results using the exact eigenvalues.

Remark that, contrary to other circuits, the target register is not explicitly given by a set of qubits, but its dimension is the dimension of the eigenvector of $U$; this implies that you must change the definition of some gates, to take into account the tensor structure of the whole circuit state.

Notes

The qiskit.org website contains complete documentation on the quantum fourier transform and its application to the phase estimation.

The diagrams of quantum circuits were drawn using the latex package quantikz written by Alastair Kay arXiv:1809.03842 v.4 (2020)

Feynman, Richard P. Simulating Physics with Computers. Inter. J. Theor. Phys. 21 467 (1982) DOI ↩↩
Deutsch, D. Quantum Theory, the Church-Turing Principle and the Universal Quantum Computer. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 400, 97-117 (1985) DOI ↩
Shor, P. W. Algorithms for Quantum Computation: Discrete Logarithms and Factoring. In Foundations of Computer Science, 1994 Proceedings., 35th Annual Symposium On, 124-34 (1994) DOI ↩
An implementation of the Grover algorithm to a real quantum computer can be found in the IBM qiskit site by John Watrous. ↩
Cleve et al. Quantum algorithms revisited, Proc. R. Soc. Lond. A 454, 339 (1998) .pdf ↩
Abrams, D., and Lloyd, S., Quantum algorithm providing exponential speed increase for finding eigenvalues and eigenvectors, Phys. Rev. Lett. 83, 5162 (1999) DOI ↩