Skip to main content
IBM Quantum Platform

Phase estimation procedure

Next, we'll discuss the phase-estimation procedure, which is a quantum algorithm for solving the phase estimation problem.

We'll begin with a low-precision warm-up, which explains some of the basic intuition behind the method. We'll then talk about the quantum Fourier transform, which is an important quantum operation used in the phase-estimation procedure, as well as its quantum circuit implementation. Once we have the quantum Fourier transform in hand, we'll describe the phase-estimation procedure in full generality and analyze its performance.


Warm-up: approximating phases with low precision

We'll begin with a couple of simple versions of the phase-estimation procedure that provide low-precision solutions to the phase-estimation problem. This is helpful for explaining the intuition behind the general procedure that we'll see a bit later in the lesson.

Using the phase kickback

A simple approach to the phase-estimation problem, which allows us to learn something about the value θ\theta we seek, is based on the phase kick-back phenomenon. As we'll see, this is essentially a single-qubit version of the general phase-estimation procedure to be discussed later in the lesson.

As part of the input to the phase estimation problem, we have a unitary quantum circuit for the operation U.U. We can use the description of this circuit to create a circuit for a controlled-UU operation, which can be depicted as this figure suggests (with the operation U,U, viewed as a quantum gate, on the left and a controlled-UU operation on the right).

Uncontrolled and controlled versions of a unitary operation

We can create a quantum circuit for a controlled-UU operation by first adding a control qubit to the circuit for U,U, and then replacing every gate in the circuit for UU with a controlled version of that gate — so our one new control qubit effectively controls every single gate in the circuit for U.U. This requires that we have a controlled version of every gate in our circuit, but we can always build circuits for these controlled operations in case they're not included in our gate set.

Now consider the following circuit, where the input state ψ\vert\psi\rangle of all of the qubits except the top one is the quantum state eigenvector of U.U.

A single-qubit circuit for phase-estimation

The measurement outcome probabilities for this circuit depend on the eigenvalue of UU corresponding to the eigenvector ψ.\vert\psi\rangle. Let's analyze the circuit in detail to determine exactly how.

States of a single-qubit circuit for phase-estimation

The initial state of the circuit is

π0=ψ0\vert\pi_0\rangle = \vert\psi\rangle \vert 0\rangle

and the first Hadamard gate transforms this state to

π1=ψ+=12ψ0+12ψ1.\vert\pi_1\rangle = \vert\psi\rangle \vert +\rangle = \frac{1}{\sqrt{2}} \vert\psi\rangle \vert 0\rangle + \frac{1}{\sqrt{2}} \vert\psi\rangle \vert 1\rangle.

Next, the controlled-UU operation is performed, which results in the state

π2=12ψ0+12(Uψ)1.\vert\pi_2\rangle = \frac{1}{\sqrt{2}} \vert\psi\rangle \vert 0\rangle + \frac{1}{\sqrt{2}} \bigl(U \vert\psi\rangle\bigr) \vert 1\rangle.

Using the assumption that ψ\vert\psi\rangle is an eigenvector of UU having eigenvalue λ=e2πiθ,\lambda = e^{2\pi i\theta}, we can alternatively express this state as follows.

π2=12ψ0+e2πiθ2ψ1=ψ(120+e2πiθ21)\vert\pi_2\rangle = \frac{1}{\sqrt{2}} \vert\psi\rangle \vert 0\rangle + \frac{e^{2\pi i \theta}}{\sqrt{2}} \vert\psi\rangle \vert 1\rangle = \vert\psi\rangle \otimes \left( \frac{1}{\sqrt{2}} \vert 0\rangle + \frac{e^{2\pi i \theta}}{\sqrt{2}} \vert 1\rangle\right)

Here we observe the phase kickback phenomenon. It is slightly different this time than it was for Deutsch's algorithm and the Deutsch-Jozsa algorithm because we're not working with a query gate — but the idea is similar.

Finally, the second Hadamard gate is performed. After just a bit of simplification, we obtain this expression for this state.

π3=ψ(1+e2πiθ20+1e2πiθ21)\vert\pi_3\rangle = \vert\psi\rangle \otimes \left( \frac{1+ e^{2\pi i \theta}}{2} \vert 0\rangle + \frac{1 - e^{2\pi i \theta}}{2} \vert 1\rangle\right)

The measurement therefore yields the outcomes 00 and 11 with these probabilities:

p0=1+e2πiθ22=cos2(πθ)p1=1e2πiθ22=sin2(πθ).\begin{aligned} p_0 &= \left\vert \frac{1+ e^{2\pi i \theta}}{2} \right\vert^2 = \cos^2(\pi\theta)\\[1mm] p_1 &= \left\vert \frac{1- e^{2\pi i \theta}}{2} \right\vert^2 = \sin^2(\pi\theta). \end{aligned}

Here's a plot of the probabilities for the two possible outcomes, 00 and 1,1, as functions of θ.\theta.

Outcome probabilities from phase kickback

Naturally, the two probabilities always sum to 1.1. Notice that when θ=0,\theta = 0, the measurement outcome is always 0,0, and when θ=1/2,\theta = 1/2, the measurement outcome is always 1.1. So, although the measurement result doesn't reveal exactly what θ\theta is, it does provide us with some information about it — and if we were promised that either θ=0\theta = 0 or θ=1/2,\theta = 1/2, we could learn from the circuit which one is correct without error.

Intuitively speaking, we can think of the circuit's measurement outcome as being a guess for θ\theta to "one bit of accuracy." In other words, if we were to write θ\theta in binary notation and round it off to one bit, we'd have a number like this:

0.a={0a=012a=1.0.a = \begin{cases} 0 & a = 0\\ \frac{1}{2} & a = 1. \end{cases}

The measurement outcome can be viewed as a guess for the bit a.a. When θ\theta is neither 00 nor 1/2,1/2, there's a nonzero probability that the guess will be wrong — but the probability of making an error becomes smaller and smaller as we get closer to 00 or 1/2.1/2.

It's natural to ask what role the two Hadamard gates play in this procedure:

  • The first Hadamard gate sets the control qubit to a uniform superposition of 0\vert 0\rangle and 1,\vert 1\rangle, so that when the phase kickback occurs, it happens for the 1\vert 1\rangle state and not the 0\vert 0\rangle state, creating a relative phase difference that affects the measurement outcomes. If we didn't do this and the phase kickback produced a global phase, it would have no effect on the probabilities of obtaining different measurement outcomes.

  • The second Hadamard gate allows us to learn something about the number θ\theta through the phenomenon of interference. Prior to the second Hadamard gate, the state of the top qubit is

    120+e2πiθ21,\frac{1}{\sqrt{2}} \vert 0\rangle + \frac{e^{2\pi i \theta}}{\sqrt{2}} \vert 1\rangle,

    and if we were to measure this state, we would obtain 00 and 11 each with probability 1/2,1/2, telling us nothing about θ.\theta. By performing the second Hadamard gate, however, we cause the number θ\theta to affect the output probabilities.

Doubling the phase

The circuit above uses the phase kickback phenomenon to approximate θ\theta to a single bit of accuracy. One bit of accuracy may be all we need in some situations — but for factoring we're going to need a lot more accuracy than that. A natural question is, how can we learn more about θ?\theta?

One very simple thing we can do is to replace the controlled-UU operation in our circuit with two copies of this operation, like in this circuit:

Single-bit phase estimation doubled

Two copies of a controlled-UU operation is equivalent to a controlled-U2U^2 operation. If ψ\vert\psi\rangle is an eigenvector of UU having eigenvalue λ=e2πiθ,\lambda = e^{2\pi i \theta}, then this state is also an eigenvector of U2,U^2, this time having eigenvalue λ2=e2πi(2θ).\lambda^2 = e^{2\pi i (2\theta)}.

So, if we run this version of the circuit, we're effectively performing the same computation as before, except that the number θ\theta is replaced by 2θ.2\theta. Here's a plot illustrating the output probabilities as θ\theta ranges from 00 to 1.1.

Outcome probabilities from double-phase kickback

Doing this can indeed provide us with some additional information about θ.\theta. If the binary representation of θ\theta is

θ=0.a1a2a3\theta = 0.a_1 a_2 a_3\cdots

then doubling θ\theta effectively shifts the binary point one position to the right:

2θ=a1.a2a32\theta = a_1. a_2 a_3\cdots

And because we're equating θ=1\theta = 1 with θ=0\theta = 0 as we move around the unit circle, we see that the bit a1a_1 has no influence on our probabilities, and we're effectively obtaining a guess for the second bit after the binary point if we round θ\theta to two bits. For instance, if we knew in advance that θ\theta was either 00 or 1/4,1/4, then we could fully trust the measurement outcome to tell us which.

It's not immediately clear, though, how this estimation should be reconciled with what we learned from the original (non-doubled) phase kickback circuit to give us the most accurate information possible about θ.\theta. So let's take a step back and consider how to proceed.

Two-qubit phase estimation

Rather than considering the two options described above separately, let's combine them into a single circuit like so.

The initial set-up for phase estimation with two qubits

The Hadamard gates after the controlled operations have been removed and there are no measurements here yet. We'll add more to the circuit as we consider our options for learning as much as we can about θ.\theta.

If we run this circuit when ψ\vert\psi\rangle is an eigenvector of U,U, the state of the bottom qubits will remain ψ\vert\psi\rangle throughout the entire circuit, and phases will be "kicked" into the state of the top two qubits. Let's analyze the circuit carefully, by means of the following figure.

States for phase estimation with two qubits

We can write the state π1\vert\pi_1\rangle like this:

π1=ψ12a0=01a1=01a1a0.\vert\pi_1\rangle = \vert \psi\rangle \otimes \frac{1}{2} \sum_{a_0 = 0}^1 \sum_{a_1 = 0}^1 \vert a_1 a_0 \rangle.

When the first controlled-UU operation is performed, the eigenvalue λ=e2πiθ\lambda = e^{2\pi i\theta} gets kicked into the phase when a0a_0 (the top qubit) is equal to 1,1, but not when it's 0.0. So, we can express the resulting state like this:

π2=ψ12a0=01a1=01e2πia0θa1a0.\vert\pi_2\rangle = \vert\psi\rangle \otimes \frac{1}{2} \sum_{a_0=0}^1 \sum_{a_1=0}^1 e^{2 \pi i a_0 \theta} \vert a_1 a_0 \rangle.

The second and third controlled-UU gates do something similar, except for a1a_1 rather than a0,a_0, and with θ\theta replaced by 2θ.2\theta. We can express the resulting state like this:

π3=ψ12a0=01a1=01e2πi(2a1+a0)θa1a0.\vert\pi_3\rangle = \vert\psi\rangle\otimes\frac{1}{2}\sum_{a_0 = 0}^1 \sum_{a_1 = 0}^1 e^{2\pi i (2 a_1 + a_0)\theta} \vert a_1 a_0 \rangle.

If we think about the binary string a1a0a_1 a_0 as representing an integer x{0,1,2,3}x \in \{0,1,2,3\} in binary notation, which is x=2a1+a0,x = 2 a_1 + a_0, we can alternatively express this state as follows.

π3=ψ12x=03e2πixθx\vert\pi_3\rangle = \vert \psi\rangle \otimes \frac{1}{2} \sum_{x = 0}^3 e^{2\pi i x \theta} \vert x \rangle

Our goal is to extract as much information about θ\theta as we can from this state.

At this point we'll consider a special case, where we're promised that θ=y4\theta = \frac{y}{4} for some integer y{0,1,2,3}.y\in\{0,1,2,3\}. In other words, we have θ{0,1/4,1/2,3/4},\theta\in \{0, 1/4, 1/2, 3/4\}, so we can express this number exactly using binary notation with two bits, as .00,00, .01,01, .10,10, or .11.11. In general, θ\theta might not be one of these four values, but thinking about this special case will help us to figure out how to most effectively extract information about θ\theta in general.

First we'll define a two-qubit state vector for each possible value y{0,1,2,3}.y \in \{0, 1, 2, 3\}.

ϕy=12x=03e2πix(y4)x=12x=03e2πixy4x\vert \phi_y\rangle = \frac{1}{2} \sum_{x = 0}^3 e^{2\pi i x (\frac{y}{4})} \vert x \rangle = \frac{1}{2} \sum_{x = 0}^3 e^{2\pi i \frac{x y}{4}} \vert x \rangle

After simplifying the exponentials, we can write these vectors as follows.

ϕ0=120+121+122+123ϕ1=120+i21122i23ϕ2=120121+122123ϕ3=120i21122+i23\begin{aligned} \vert\phi_0\rangle & = \frac{1}{2} \vert 0 \rangle + \frac{1}{2} \vert 1 \rangle + \frac{1}{2} \vert 2 \rangle + \frac{1}{2} \vert 3 \rangle \\[3mm] \vert\phi_1\rangle & = \frac{1}{2} \vert 0 \rangle + \frac{i}{2} \vert 1 \rangle - \frac{1}{2} \vert 2 \rangle - \frac{i}{2} \vert 3 \rangle \\[3mm] \vert\phi_2\rangle & = \frac{1}{2} \vert 0 \rangle - \frac{1}{2} \vert 1 \rangle + \frac{1}{2} \vert 2 \rangle - \frac{1}{2} \vert 3 \rangle \\[3mm] \vert\phi_3\rangle & = \frac{1}{2} \vert 0 \rangle - \frac{i}{2} \vert 1 \rangle - \frac{1}{2} \vert 2 \rangle + \frac{i}{2} \vert 3 \rangle \end{aligned}

These vectors are orthogonal: if we choose any pair of them and compute their inner product, we get 0.0. Each one is also a unit vector, so {ϕ0,ϕ1,ϕ2,ϕ3}\{\vert\phi_0\rangle, \vert\phi_1\rangle, \vert\phi_2\rangle, \vert\phi_3\rangle\} is an orthonormal basis. We therefore know right away that there is a measurement that can discriminate them perfectly — meaning that, if we're given one of them but we don't know which, then we can figure out which one it is without error.

To perform such a discrimination with a quantum circuit, we can first define a unitary operation VV that transforms standard basis states into the four states listed above.

V00=ϕ0V01=ϕ1V10=ϕ2V11=ϕ3\begin{aligned} V \vert 00 \rangle & = \vert\phi_0\rangle \\ V \vert 01 \rangle & = \vert\phi_1\rangle \\ V \vert 10 \rangle & = \vert\phi_2\rangle \\ V \vert 11 \rangle & = \vert\phi_3\rangle \end{aligned}

To write down VV as a 4×44\times 4 matrix, it's just a matter of taking the columns of VV to be the states ϕ0,,ϕ3.\vert\phi_0\rangle,\ldots,\vert\phi_3\rangle.

V=12(11111i1i11111i1i)V = \frac{1}{2} \begin{pmatrix} 1 & 1 & 1 & 1\\[1mm] 1 & i & -1 & -i\\[1mm] 1 & -1 & 1 & -1\\[1mm] 1 & -i & -1 & i \end{pmatrix}

This is a special matrix, and it's likely that some readers will have encountered it before: it's the matrix associated with the 44-dimensional discrete Fourier transform. In light of this fact, let us call it by the name QFT4\mathrm{QFT}_4 rather than V.V. The name QFT\mathrm{QFT} is short for quantum Fourier transform — which is essentially just the discrete Fourier transform, viewed as a unitary operation. We'll discuss the quantum Fourier transform in greater detail and generality shortly.

QFT4=12(11111i1i11111i1i)\mathrm{QFT}_4 = \frac{1}{2} \begin{pmatrix} 1 & 1 & 1 & 1\\[1mm] 1 & i & -1 & -i\\[1mm] 1 & -1 & 1 & -1\\[1mm] 1 & -i & -1 & i \end{pmatrix}

We can perform the inverse of this operation to go the other way, to transform the states ϕ0,,ϕ3\vert\phi_0\rangle,\ldots,\vert\phi_3\rangle into the standard basis states 0,,3.\vert 0\rangle,\ldots,\vert 3\rangle. If we do this, then we can measure to learn which value y{0,1,2,3}y\in\{0,1,2,3\} describes θ\theta as θ=y/4.\theta = y/4. Here's a diagram of a quantum circuit that does this.

Phase estimation with two qubits

To summarize, if we run this circuit when θ=y/4\theta = y/4 for y{0,1,2,3},y\in\{0,1,2,3\}, the state immediately before the measurements take place will be ψy\vert \psi\rangle \vert y\rangle (for yy encoded as a two-bit binary string), so the measurements will reveal the value yy without error.

This circuit is motivated by the special case that θ{0,1/4,1/2,3/4}\theta \in \{0,1/4,1/2,3/4\} — but we can run it for any choice of UU and ψ,\vert \psi\rangle, and hence any value of θ,\theta, that we wish. Here's a plot of the output probabilities the circuit produces for arbitrary choices of θ.\theta.

Outcome probabilities from two-qubit phase estimation

This is a clear improvement over the single-qubit variant described earlier in the lesson. It's not perfect — it can give us the wrong answer — but the answer is heavily skewed toward values of yy for which y/4y/4 is close to θ.\theta. In particular, the most likely outcome always corresponds to the closest value of y/4y/4 to θ\theta (equating θ=0\theta = 0 and θ=1\theta = 1 as before), and from the plot it looks like this closest value for xx always appears with probability just above 40%.40\%. When θ\theta is exactly halfway between two such values, like θ=0.375\theta = 0.375 for instance, the two equally close values of yy are equally likely.

Preparing to generalize to many qubits

Given the improvement we've just obtained by using two control qubits rather than one, in conjunction with the inverse of the 44-dimensional quantum Fourier transform, it's natural to consider generalizing it further — by adding more control qubits. When we do this, we obtain the general phase estimation procedure. We'll see how this works shortly, but in order to describe it precisely we're going to need to discuss the quantum Fourier transform in greater generality, to see how it's defined for other dimensions and to see how we can implement it (or its inverse) with a quantum circuit.


Quantum Fourier transform

The quantum Fourier transform is a unitary operation that can be defined for any positive integer dimension N.N. In this section, we'll see how this operation is defined and how it can be implemented with a quantum circuit on mm qubits with cost O(m2)O(m^2) when N=2m.N = 2^m.

The matrices that describe the quantum Fourier transform are derived from an analogous operation on NN-dimensional vectors known as the discrete Fourier transform. This operation can be thought about in different ways. For instance, we can think about the discrete Fourier transform in purely abstract, mathematical terms as a linear mapping. Or we can think about it in computational terms, where we're given an NN-dimensional vector of complex numbers (using binary notation to encode the real and imaginary parts of the entries, let us suppose) and the goal is to calculate the NN-dimensional vector obtained by applying the discrete Fourier transform. Our focus will be on third way, which is viewing this transformation as a unitary operation that can be performed on a quantum system.

There's an efficient algorithm for computing the discrete Fourier transform on a given input vector known as the fast Fourier transform. It has applications in signal processing and many other areas, and is considered by many to be one of the most important algorithms ever discovered. As it turns out, the implementation of the quantum Fourier transform when NN is a power of 2 that we'll study is based on precisely the same underlying structure that make the fast Fourier transform possible.

Definition of the quantum Fourier transform

To define the quantum Fourier transform, we'll first define a complex number ωN,\omega_N, for each positive integer N,N, like this:

ωN=e2πiN=cos(2πN)+isin(2πN).\omega_N = e^{\frac{2\pi i}{N}} = \cos\left(\frac{2\pi}{N}\right) + i \sin\left(\frac{2\pi}{N}\right).

This is the number on the complex unit circle we obtain if we start at 11 and move counter-clockwise by an angle of 2π/N2\pi/N radians, or a fraction of 1/N1/N of the circumference of the circle. Here are a few examples:

ω1=1ω2=1ω3=12+32iω4=iω8=1+i2ω16=2+22+222iω1000.998+0.063i\begin{gathered} \omega_1 = 1\\[1mm] \omega_2 = -1\\[1mm] \omega_3 = -\frac{1}{2} + \frac{\sqrt{3}}{2} i\\[2mm] \omega_4 = i\\[1mm] \omega_8 = \frac{1+i}{\sqrt{2}}\\[3mm] \omega_{16} = \frac{\sqrt{2 + \sqrt{2}}}{2} + \frac{\sqrt{2 - \sqrt{2}}}{2} i\\[2mm] \omega_{100} \approx 0.998 + 0.063 i \end{gathered}

Now we can define the NN-dimensional quantum Fourier transform, which is described by an N×NN\times N matrix whose rows and columns are associated with the standard basis states 0,,N1.\vert 0\rangle,\ldots,\vert N-1\rangle. We're only going to need this operation for when N=2mN = 2^m is a power of 22 for phase estimation, but the operation can be defined for any positive integer N.N.

QFTN=1Nx=0N1y=0N1ωNxyxy\mathrm{QFT}_N = \frac{1}{\sqrt{N}} \sum_{x = 0}^{N-1} \sum_{y = 0}^{N-1} \omega_N^{xy} \vert x \rangle\langle y\vert

As was already stated, this is the matrix associated with the NN-dimensional discrete Fourier transform. Often the leading factor of 1/N1/\sqrt{N} is not included in the definition of this matrix, but we need to include it to obtain a unitary matrix.

Here's the quantum Fourier transform, written as a matrix, for some small values of N.N.

QFT1=(1)\mathrm{QFT}_1 = \begin{pmatrix} 1 \end{pmatrix} QFT2=12(1111)\mathrm{QFT}_2 = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1\\[1mm] 1 & -1 \end{pmatrix} QFT3=13(11111+i321i3211i321+i32)\mathrm{QFT}_3 = \frac{1}{\sqrt{3}} \begin{pmatrix} 1 & 1 & 1\\[2mm] 1 & \frac{-1 + i\sqrt{3}}{2} & \frac{-1 - i\sqrt{3}}{2}\\[2mm] 1 & \frac{-1 - i\sqrt{3}}{2} & \frac{-1 + i\sqrt{3}}{2} \end{pmatrix} QFT4=12(11111i1i11111i1i)\mathrm{QFT}_4 = \frac{1}{2} \begin{pmatrix} 1 & 1 & 1 & 1\\[1mm] 1 & i & -1 & -i\\[1mm] 1 & -1 & 1 & -1\\[1mm] 1 & -i & -1 & i \end{pmatrix} QFT8=122(1111111111+i2i1+i211i2i1i21i1i1i1i11+i2i1+i211i2i1i21111111111i2i1i211+i2i1+i21i1i1i1i11i2i1i211+i2i1+i2)\mathrm{QFT}_8 = \frac{1}{2\sqrt{2}} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\[2mm] 1 & \frac{1+i}{\sqrt{2}} & i & \frac{-1+i}{\sqrt{2}} & -1 & \frac{-1-i}{\sqrt{2}} & -i & \frac{1-i}{\sqrt{2}}\\[2mm] 1 & i & -1 & -i & 1 & i & -1 & -i\\[2mm] 1 & \frac{-1+i}{\sqrt{2}} & -i & \frac{1+i}{\sqrt{2}} & -1 & \frac{1-i}{\sqrt{2}} & i & \frac{-1-i}{\sqrt{2}}\\[2mm] 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1\\[2mm] 1 & \frac{-1-i}{\sqrt{2}} & i & \frac{1-i}{\sqrt{2}} & -1 & \frac{1+i}{\sqrt{2}} & -i & \frac{-1+i}{\sqrt{2}}\\[2mm] 1 & -i & -1 & i & 1 & -i & -1 & i\\[2mm] 1 & \frac{1-i}{\sqrt{2}} & -i & \frac{-1-i}{\sqrt{2}} & -1 & \frac{-1+i}{\sqrt{2}} & i & \frac{1+i}{\sqrt{2}}\\[2mm] \end{pmatrix}

Notice, in particular, that QFT2\mathrm{QFT}_2 is another name for a Hadamard operation.

Unitarity

Let's check that QFTN\mathrm{QFT}_N is unitary, for any selection of N.N. One way to do this is to show that its columns form an orthonormal basis. We can define a vector corresponding to column number y,y, starting from y=0y = 0 and going up to y=N1,y = N-1, like this:

ϕy=1Nx=0N1ωNxyx.\vert\phi_y\rangle = \frac{1}{\sqrt{N}} \sum_{x = 0}^{N-1} \omega_N^{xy} \vert x \rangle.

Taking the inner product between any two of these vectors gives us this expression:

ϕzϕy=1Nx=0N1ωNx(yz)\langle \phi_z \vert \phi_y \rangle = \frac{1}{N} \sum_{x = 0}^{N-1} \omega_N^{x (y - z)}

We can evaluate sums like this using the following formula for the sum of the first NN terms of a geometric series.

1+α+α2++αN1={αN1α1if α1Nif α=11 + \alpha + \alpha^2 + \cdots + \alpha^{N-1} = \begin{cases} \frac{\alpha^N - 1}{\alpha - 1} & \text{if } \alpha\neq 1\\[2mm] N & \text{if } \alpha=1 \end{cases}

Specifically, we can use this formula when α=ωNyz.\alpha = \omega_N^{y-z}. When y=z,y = z, we have α=1,\alpha = 1, so using the formula and dividing by NN gives

ϕyϕy=1.\langle \phi_y \vert \phi_y \rangle = 1.

When yz,y\neq z, we have α1,\alpha \neq 1, so the formula reveals this:

ϕzϕy=1NωNN(yz)1ωNyz1=1N11ωNyz1=0.\langle \phi_z \vert \phi_y \rangle = \frac{1}{N} \frac{\omega_N^{N(y-z)} - 1}{\omega_N^{y-z} - 1} = \frac{1}{N} \frac{1 - 1}{\omega_N^{y-z} - 1} = 0.

This happens because ωNN=e2πi=1,\omega_N^N = e^{2\pi i} = 1, so ωNN(yz)=1yz=1,\omega_N^{N(y-z)} = 1^{y-z} = 1, making numerator zero, while the denominator is nonzero because ωNyz1.\omega_N^{y-z} \neq 1. Intuitively speaking, what we're doing is summing a bunch of points that are distributed around the unit circle, and they cancel out and leave 00 when summed.

We have therefore established that {ϕ0,,ϕN1}\{\vert\phi_0\rangle,\ldots,\vert\phi_{N-1}\rangle\} is an orthonormal set,

ϕzϕy={1y=z0yz,\langle \phi_z \vert \phi_y \rangle = \begin{cases} 1 & y=z\\ 0 & y\neq z, \end{cases}

which reveals that QFTN\mathrm{QFT}_N is unitary.

Controlled-phase gates

To implement the quantum Fourier transform with a quantum circuit, we'll need to make use of controlled-phase gates. Recall that a phase operation is a single-qubit unitary operation of the form

Pα=(100eiα)P_{\alpha} = \begin{pmatrix} 1 & 0\\[1mm] 0 & e^{i\alpha} \end{pmatrix}

for any real number α.\alpha. A controlled version of this gate has the following matrix:

CPα=(100001000010000eiα)CP_{\alpha} = \begin{pmatrix} 1 & 0 & 0 & 0\\[1mm] 0 & 1 & 0 & 0\\[1mm] 0 & 0 & 1 & 0\\[1mm] 0 & 0 & 0 & e^{i\alpha} \end{pmatrix}

For this controlled gate, it doesn't actually matter which qubit is the control and which is the target because the two possibilities are equivalent. We can use any of the following symbols to represent this gate in quantum circuit diagrams.

Quantum circuit diagram representation for controlled-phase gates

For the third form, the label α\alpha is also sometimes placed on the side of the control line or under the lower control when that's convenient.

To perform the quantum Fourier transform when N=2mN = 2^m and m2,m\geq 2, we're going to need to perform an operation on mm qubits whose action on standard basis states can be described as

yaω2mayya,\vert y \rangle \vert a \rangle \mapsto \omega_{2^m}^{ay} \vert y \rangle \vert a \rangle,

where aa is a bit and y{0,,2m11}y \in \{0,\ldots,2^{m-1} - 1\} is a number encoded in binary notation as a string of m1m-1 bits. This can be done using controlled-phase gates by generalizing the following example, for which m=5.m=5.

Quantum circuit diagram for phase injection

In general, for an arbitrary choice of m2,m\geq 2, the top qubit corresponding to the bit aa can be viewed as the control, with the phase gates PαP_{\alpha} ranging from α=π/2m1\alpha = \pi/2^{m-1} on the qubit corresponding to the least significant bit of yy to α=π2\alpha = \frac{\pi}{2} on the qubit corresponding to the most significant bit of y.y. These controlled-phase gates all commute with one another and could be performed in any order.

Circuit implementation of the QFT

Now we'll see how we can implement the quantum Fourier transform with a circuit when the dimension N=2mN = 2^m is a power of 2.2. There are, in fact, multiple ways to implement the quantum Fourier transform, but this is arguably the simplest method known. Once we know how to implement the quantum Fourier transform with a quantum circuit, it's straightforward to implement its inverse: we can replace each gate with its inverse (or, equivalently, conjugate transpose) and apply the gates in the reverse order. Every quantum circuit composed of unitary gates alone can be inverted in this way.

The implementation is recursive in nature, so that's how it's most naturally described. The base case is m=1,m=1, in which case the quantum Fourier transform is a Hadamard operation.

To perform the quantum Fourier transform on mm qubits when m2,m \geq 2, we can perform the following steps, whose actions we'll describe for standard basis states of the form xa,\vert x \rangle \vert a\rangle, where x{0,,2m11}x\in\{0,\ldots,2^{m-1} - 1\} is an integer encoded as m1m-1 bits using binary notation and aa is a single bit.

  1. First apply the 2m12^{m-1}-dimensional quantum Fourier transform to the bottom/leftmost m1m-1 qubits to obtain this state:

    (QFT2m1x)a=12m1y=02m11ω2m1xyya.\Bigl(\mathrm{QFT}_{2^{m-1}} \vert x \rangle\Bigr) \vert a\rangle = \frac{1}{\sqrt{2^{m-1}}} \sum_{y = 0}^{2^{m-1} - 1} \omega_{2^{m-1}}^{xy} \vert y \rangle \vert a \rangle.

    This is done by recursively applying the method being described for one fewer qubit, using the Hadamard operation on a single qubit as the base case.

  2. Use the top/rightmost qubit as a control to inject the phase ω2my\omega_{2^m}^y for each standard basis state y\vert y\rangle of the remaining m1m-1 qubits (as is described above) to obtain this state:

    12m1y=02m11ω2m1xyω2mayya.\frac{1}{\sqrt{2^{m-1}}} \sum_{y = 0}^{2^{m-1} - 1} \omega_{2^{m-1}}^{xy} \omega_{2^m}^{ay} \vert y \rangle \vert a \rangle.
  3. Perform a Hadamard gate on the top/rightmost qubit to obtain this state:

    12my=02m11b=01(1)abω2m1xyω2mayyb.\frac{1}{\sqrt{2^{m}}} \sum_{y = 0}^{2^{m-1} - 1} \sum_{b=0}^1 (-1)^{ab} \omega_{2^{m-1}}^{xy} \omega_{2^m}^{ay} \vert y \rangle \vert b \rangle.
  4. Permute the order of the qubits so that the least significant bit becomes the most significant bit, with all others shifted up/right:

    12my=02m11b=01(1)abω2m1xyω2mayby.\frac{1}{\sqrt{2^{m}}} \sum_{y = 0}^{2^{m-1} - 1} \sum_{b=0}^1 (-1)^{ab} \omega_{2^{m-1}}^{xy} \omega_{2^m}^{ay} \vert b \rangle \vert y \rangle.

For example, here's the circuit we obtain for N=32=25.N = 32 = 2^5. In this diagram, the qubits are given names that correspond to the standard basis vectors xa\vert x\rangle \vert a\rangle (for the input) and by\vert b\rangle \vert y\rangle (for the output) for clarity.

Quantum circuit diagram for the 32-dimensional quantum Fourier transform

Analysis

The key formula we need to verify that the circuit just described implements the 2m2^m-dimensional quantum Fourier transform is this one:

(1)abω2m1xyω2may=ω2m(2x+a)(2m1b+y).(-1)^{ab} \omega_{2^{m-1}}^{xy} \omega_{2^m}^{ay} = \omega_{2^m}^{(2x+ a)(2^{m-1}b + y)}.

This formula works for any choice of integers a,a, b,b, x,x, and y,y, but we'll only need it for a,b{0,1}a,b\in\{0,1\} and x,y{0,,2m11}.x,y\in\{0,\ldots,2^{m-1}-1\}. It can be checked by expanding the product in the exponent on the right-hand side,

ω2m(2x+a)(2m1b+y)=ω2m2mxbω2m2xyω2m2m1abω2may=(1)abω2m1xyω2may, \omega_{2^m}^{(2x+ a)(2^{m-1}b + y)} = \omega_{2^m}^{2^m xb} \omega_{2^m}^{2xy} \omega_{2^m}^{2^{m-1}ab} \omega_{2^m}^{ay} = (-1)^{ab} \omega_{2^{m-1}}^{xy} \omega_{2^m}^{ay},

where the second equality makes use of the observation that

ω2m2mxb=(ω2m2m)xb=1xb=1.\omega_{2^m}^{2^m xb} = \bigl(\omega_{2^m}^{2^m}\bigr)^{xb} = 1^{xb} = 1.

The 2m2^m-dimensional quantum Fourier transform is defined as follows for every u{0,,2m1}.u\in\{0,\ldots,2^m - 1\}.

QFT2mu=12mv=02m1ω2muvv\mathrm{QFT}_{2^m} \vert u\rangle = \frac{1}{\sqrt{2^m}} \sum_{v = 0}^{2^m - 1} \omega_{2^m}^{uv} \vert v\rangle

If we write uu and vv as

u=2x+av=2m1b+y\begin{aligned} u & = 2x + a\\ v & = 2^{m-1}b + y \end{aligned}

for a,b{0,1}a,b\in\{0,1\} and x,y{0,,2m11},x,y\in\{0,\ldots,2^{m-1} - 1\}, we obtain

QFT2m2x+a=12my=02m11b=01ω2m(2x+a)(2m1b+y)b2m1+y=12my=02m11b=01(1)abω2m1xyω2mayb2m1+y.\begin{aligned} \mathrm{QFT}_{2^m} \vert 2x + a\rangle & = \frac{1}{\sqrt{2^m}} \sum_{y = 0}^{2^{m-1} - 1} \sum_{b=0}^1 \omega_{2^m}^{(2x+ a)(2^{m-1}b + y)} \vert b 2^{m-1} + y\rangle\\[2mm] & = \frac{1}{\sqrt{2^m}} \sum_{y = 0}^{2^{m-1} - 1} \sum_{b=0}^1 (-1)^{ab} \omega_{2^{m-1}}^{xy} \omega_{2^m}^{ay} \vert b 2^{m-1} + y\rangle. \end{aligned}

Finally, by thinking about the standard basis states xa\vert x \rangle \vert a\rangle and by\vert b \rangle \vert y \rangle as binary encodings of integers in the range {0,,2m1},\{0,\ldots,2^m-1\},

xa=2x+aby=2m1b+y,\begin{aligned} \vert x \rangle \vert a\rangle & = \vert 2x + a \rangle\\ \vert b \rangle \vert y \rangle & = \vert 2^{m-1}b + y\rangle, \end{aligned}

we see that the circuit above implements the required operation. If this method for performing the quantum Fourier transform seems remarkable, it's because it is: it's essentially the fast Fourier transform in the form of a quantum circuit.

Finally, let's count how many gates are used in the circuit just described. The controlled-phase gates aren't in the standard gate set that we discussed in the previous lesson, but to begin we'll ignore this and count each of them as a single gate.

Let's let sms_m denote the number of gates we need for each possible choice of m.m. If m=1,m=1, the quantum Fourier transform is just a Hadamard operation, so

s1=1.s_1 = 1.

If m2,m\geq 2, then in the circuit above we need sm1s_{m-1} gates for the quantum Fourier transform on m1m-1 qubits, plus m1m-1 controlled-phase gates, plus a Hadamard gate, plus m1m-1 swap gates, so

sm=sm1+(2m1).s_m = s_{m-1} + (2m - 1).

We can obtain a closed-form expression by summing:

sm=k=1m(2k1)=m2.s_m = \sum_{k = 1}^m (2k - 1) = m^2.

We don't actually need as many swap gates as the method describes. If we rearrange the gates just a bit, we can push all of the swap gates out to the right and reduce the number of swap gates required to m/2.\lfloor m/2\rfloor. Asymptotically speaking this isn't a major improvement: we still obtain circuits with size O(m2)O(m^2) for performing QFT2m.\mathrm{QFT}_{2^m}.

If we wish to implement the quantum Fourier transform using only gates from our standard gate set, we need to either build or approximate each of the controlled-phase gates with gates from our set. The number required depends on how much accuracy we require, but as a function of mm the total cost remains quadratic.

It is, in fact, possible to approximate the quantum Fourier transform quite closely with a sub-quadratic number of gates by using the fact that PαP_{\alpha} is very close to the identity operation when α\alpha is very small — which means that we can simply leave out most of the controlled-phase gates without suffering too much of a loss in terms of accuracy.


General procedure and analysis

Now we'll examine the phase-estimation procedure in general. The idea is to extend the two-qubit version of phase estimation that we considered above in the natural way suggested by the following diagram.

Phase estimation procedure

Notice that, for each new control qubit added on the top, we double the number of times the unitary operation UU is performed. This is indicated in the diagram by the powers on UU for each of the controlled-unitary operations.

The most straightforward way to implement a controlled-UkU^k operation for some choice of kk is simply to repeat a controlled-UU operation kk times. If this is indeed the methodology that is used, it must be recognized that the addition of control qubits contributes significantly to the size of the circuit: if we have mm control qubits, like the diagram depicts, a total of 2m12^m - 1 copies of the controlled-UU operation are required. This means that a significant computational cost is incurred as mm is increased — but as we will see, it also leads to a significantly more accurate approximation of θ.\theta.

It is important to note, however, that for some choices of UU it may be possible to create a circuit that implements the operation UkU^k for large values of kk in a more efficient way than simply repeating kk times the circuit for U.U. We'll see a specific example of this in the context of integer factorization later in the lesson, where the efficient algorithm for modular exponentiation discussed in the previous lesson comes to the rescue.

Now let us analyze the circuit just described. The state immediately prior to the inverse quantum Fourier transform looks like this:

12mx=02m1(Uxψ)x=ψ12mx=02m1e2πixθx.\frac{1}{\sqrt{2^m}} \sum_{x = 0}^{2^m - 1} \bigl( U^x \vert\psi\rangle \bigr) \vert x\rangle = \vert\psi\rangle \otimes \frac{1}{\sqrt{2^m}} \sum_{x = 0}^{2^m - 1} e^{2\pi i x\theta} \vert x\rangle.

A special case

Along similar lines to what we did in the m=2m=2 case, we'll first consider the special case that θ=y/2m\theta = y/2^m for y{0,,2m1}.y\in\{0,\ldots,2^m-1\}. In this case the state prior to the inverse quantum Fourier transform can alternatively be written like this:

ψ12mx=02m1e2πixy2mx=ψ12mx=02m1ω2mxyx=ψQFT2my.\vert\psi\rangle \otimes \frac{1}{\sqrt{2^m}} \sum_{x = 0}^{2^m - 1} e^{2\pi i \frac{xy}{2^m}} \vert x\rangle = \vert\psi\rangle \otimes \frac{1}{\sqrt{2^m}} \sum_{x = 0}^{2^m - 1} \omega_{2^m}^{xy} \vert x\rangle = \vert\psi\rangle \otimes \mathrm{QFT}_{2^m} \vert y\rangle.

So, when the inverse quantum Fourier transform is applied, the state becomes

ψy\vert\psi\rangle \vert y\rangle

and the measurements reveal yy (encoded in binary).

Bounding the probabilities

For other values of θ,\theta, meaning ones that don't take the form y/2my/2^m for an integer y,y, the measurement outcomes won't be certain, but we can prove bounds on the probabilities for different outcomes. Going forward, let's consider an arbitrary choice of θ\theta satisfying 0θ<1.0\leq \theta < 1.

After the inverse quantum Fourier transform is performed, the state of the circuit is this:

ψ12my=02m1x=02m1e2πix(θy/2m)y.\vert \psi \rangle \otimes \frac{1}{2^m} \sum_{y=0}^{2^m - 1} \sum_{x=0}^{2^m-1} e^{2\pi i x (\theta - y/2^m)} \vert y\rangle.

So, when the measurements on the top mm qubits are performed, we see each outcome yy with probability

py=12mx=02m1e2πix(θy/2m)2.p_y = \left\vert \frac{1}{2^m} \sum_{x=0}^{2^m - 1} e^{2\pi i x (\theta - y/2^m)} \right\vert^2.

To get a better handle on these probabilities, we'll make use of the same formula that we saw before, for the sum of the initial portion of a geometric series.

1+α+α2++αN1={αN1α1if α1Nif α=11 + \alpha + \alpha^2 + \cdots + \alpha^{N-1} = \begin{cases} \frac{\alpha^N - 1}{\alpha - 1} & \text{if } \alpha\neq 1\\[2mm] N & \text{if } \alpha=1 \end{cases}

We can simplify the sum appearing in the formula for pyp_y by taking α=e2πi(θy/2m).\alpha = e^{2\pi i (\theta - y/2^m)}. Here's what we obtain.

x=02m1e2πix(θy/2m)={2mθ=y/2me2π(2mθy)1e2π(θy/2m)1θy/2m\sum_{x=0}^{2^m - 1} e^{2\pi i x (\theta - y/2^m)} = \begin{cases} 2^m & \theta = y/2^m\\[2mm] \frac{e^{2\pi (2^m \theta - y)} - 1}{e^{2\pi (\theta - y/2^m)} - 1} & \theta\neq y/2^m \end{cases}

So, in the case that θ=y/2m,\theta = y/2^m, we find that py=1p_y = 1 (as we already knew from considering this special case), and in the case that θy/2m,\theta \neq y/2^m, we find that

py=122me2πi(2mθy)1e2πi(θy/2m)12.p_y = \frac{1}{2^{2m}} \left\vert \frac{e^{2\pi i (2^m \theta - y)} - 1}{e^{2\pi i (\theta - y/2^m)} - 1}\right\vert^2.

We can learn more about these probabilities by thinking about how arc lengths and chord lengths on the unit circle are related. Here's a figure that illustrates the relationships we need for any real number δ[12,12].\delta\in \bigl[ -\frac{1}{2},\frac{1}{2}\bigr].

Illustration of the relationship between arc and chord lengths

First, the chord length (drawn in blue) can't possibly be larger than the arc length (drawn in purple):

e2πiδ12πδ.\bigl\vert e^{2\pi i \delta} - 1\bigr\vert \leq 2\pi\vert\delta\vert.

Relating these lengths in the other direction, we see that the ratio of the arc length to the chord length is greatest when δ=±1/2,\delta = \pm 1/2, and in this case the ratio is half the circumference of the circle divided by the diameter, which is π/2.\pi/2. Thus, we have

2πδe2πiδ1π2,\frac{2\pi\vert\delta\vert}{\bigl\vert e^{2\pi i \delta} - 1\bigr\vert} \leq \frac{\pi}{2},

and so

e2πiδ14δ.\bigl\vert e^{2\pi i \delta} - 1\bigr\vert \geq 4\vert\delta\vert.

An analysis based on these relations reveals the following two facts.

  1. Suppose that θ\theta is a real number and y{0,,2m1}y\in \{0,\ldots,2^m-1\} satisfies

    θy2m2(m+1).\Bigl\vert \theta - \frac{y}{2^m}\Bigr\vert \leq 2^{-(m+1)}.

    This means that y/2my/2^m is either the best mm-bit approximation to θ,\theta, or it's exactly halfway between y/2my/2^m and either (y1)/2m(y-1)/2^m or (y+1)/2m,(y+1)/2^m, so it's one of the two best approximations to θ.\theta.

    We'll prove that pyp_y has to be pretty large in this case. By the assumption we're considering, it follows that 2mθy1/2,\vert 2^m \theta - y \vert \leq 1/2, so we can use the second observation above relating arc and chord lengths to conclude that

    e2πi(2mθy)142mθy=42mθy2m.\left\vert e^{2\pi i (2^m \theta - y)} - 1\right\vert \geq 4 \vert 2^m \theta - y \vert = 4 \cdot 2^m \cdot \Bigl\vert \theta - \frac{y}{2^m}\Bigr\vert.

    We can also use the first observation about arc and chord lengths to conclude that

    e2πi(θy/2m)12πθy2m.\left\vert e^{2\pi i (\theta - y/2^m)} - 1\right\vert \leq 2\pi \Bigl\vert \theta - \frac{y}{2^m}\Bigr\vert.

    Putting these two inequalities to use on pyp_y reveals

    py122m1622m4π2=4π20.405.p_y \geq \frac{1}{2^{2m}} \frac{16 \cdot 2^{2m}}{4 \pi^2} = \frac{4}{\pi^2} \approx 0.405.

    This explains our observation that the best outcome occurs with probability greater than 40%40\% in the m=2m=2 version of phase estimation discussed earlier. It's not really 40%, it's 4/π2,4/\pi^2, and in fact this bound holds for every choice of m.m.

  2. Now suppose that y{0,,2m1}y\in \{0,\ldots,2^m-1\} satisfies

    2mθy2m12.2^{-m} \leq \Bigl\vert \theta - \frac{y}{2^m}\Bigr\vert \leq \frac{1}{2}.

    This means that there's a better approximation z/2mz/2^m to θ\theta in between θ\theta and y/2m.y/2^m.

    This time we'll prove that pyp_y can't be too big. We can start with the simple observation that

    e2πi(2mθy)12,\left\vert e^{2\pi i (2^m \theta - y)} - 1\right\vert \leq 2,

    which follows from the fact that any two points on the unit circle can differ in absolute value by at most 2.2.

    We can also use the second observation about arc and chord lengths from above, this time working with the denominator of pyp_y rather than the numerator, to conclude

    e2πi(θy/2m)14θy2m42m.\left\vert e^{2\pi i (\theta - y/2^m)} - 1\right\vert \geq 4\Bigl\vert \theta - \frac{y}{2^m}\Bigr\vert \geq 4 \cdot 2^{-m}.

    Putting the two inequalities together reveals

    py122m41622m=14.p_y \leq \frac{1}{2^{2m}} \frac{4}{16 \cdot 2^{-2m}} = \frac{1}{4}.

    Note that, while this bound is good enough for our purposes, it is fairly crude — the probability is usually much lower than 1/4.1/4.

The important take-away from this analysis is that very close approximations to θ\theta are likely to occur — we'll get a best mm-bit approximation with probability greater than 40%40\% — whereas approximations off by more than 2m2^{-m} are less likely to occur, with probability upper bounded by 25%.25\%.

Given these guarantees, it is possible to boost our confidence by repeating the phase estimation procedure several times, to gather statistical evidence about θ.\theta. It is important to note that the state ψ\vert\psi\rangle of the bottom collection of qubits is unchanged by the phase estimation procedure, so it can be used to run the procedure as many times as we like. In particular, each time we run the circuit, we get a best mm-bit approximation to θ\theta with probability greater than 40%,40\%, while the probability of being off by more than 2m2^{-m} is bounded by 25%.25\%. If we run the circuit several times and take the most commonly appearing outcome of the runs, it's therefore exceedingly likely that the outcome that appears most commonly will not be one that occurs at most 25%25\% of the time. As a result, we'll be very likely to obtain an approximation y/2my/2^m that's within 1/2m1/2^m of the value θ.\theta. Indeed, the unlikely chance that we're off by more than 1/2m1/2^m decreases exponentially in the number of times the procedure is run.

Here are two plots showing the probabilities for three consecutive values for yy when m=3m = 3 and m=4m=4 as functions of θ.\theta. (Only three outcomes are shown for clarity. Probabilities for other outcomes are obtained by cyclically shifting the same underlying function.)

Plot showing outcome probabilities for three-qubit phase estimation Plot showing outcome probabilities for four-qubit phase estimation
Was this page helpful?
Report a bug or request content on GitHub.