116. Quantum Fourier Transform

We have seen algorithms that work on shallow circuit quantum computers or noisy quantum computers. In other words, quantum computers, which don’t have a long coherence time. They have implementation difficulties. In this module, we will merely talk about algorithms that assume a perfect quantum computer and that require very high quality qubits and that allow you to run arbitrarily long algorithms. These are the kind of algorithms that can give you polynomial or exponential speed ups compared to best known classical algorithms. So these are what we call coherent protocols. Coherent not because they don’t necessarily interact with the environment. We might do some measurement as part of the protocol, but what we mean by coherent here is that we assume this level of perfection of the quantum computer that I mentioned. And what’s very important to most of these algorithms is that their inputs and outputs are all quantum states. This is slightly different than what we have seen so far because when you think about quantum annealing, you are actually preparing something classical. And then you get some classical output. If you think about QAOA, it is the same story. Now here, it’s very important that your inputs and outputs are quantum states. Because if you use state preparation, or you want to understand what the final output is and you want to do tomography, that can kill all the speed up that you get from the protocol itself. So it’s a slightly different paradigm. Before we look at one of the most important building blocks in quantum field transform, let’s take a look at what its classical variant does. So classically, we have some vector, which, for instance, describes a time series. And by applying Fourier transformation rate, we get, say, the frequency domain decomposition of the signal. So you go from the time domain to the frequency domain. And the way you do it is by applying this exponential operation for every single element on xj. And this reveals the frequency structure. As you are going to see in quantum field estimation, this is what we used to go from, say, an amplitude encoding to business encoding. So it’s very useful for a whole bunch of different things when you put it in the right context. Now, the quantum Fourier transformation achieves exactly this operation. So now you have the quantum state encoding x in amplitude encoding, so the elements of the x factor encoding in the amplitude. And you create a new vector y by applying the quantum Fourier transformation, where the elements or the amplitudes in this new vector are exactly these elements that you would do in a classical Fourier transformation. To keep things simple, I assumed that the number of elements in this vector, capital N, is some power of 2. It’s 2 to the n. So the actual transformation is surprisingly simple. So we are looking for some unitary that performs this transformation, right? And we can create this binary expansion of a number. So the binary expansion is just basically its digits expanded in the binary basis. And if you do that, then they can see the actual structure of how numbers transform, and that’s exactly what we are going to do. So this would be the Fourier transformation that we are looking for. So I just rewrote whatever was here, here. So this is what we want to get at the end of it. And we are going to rearrange things here, just some basic algebra, to see what is the actual operation we have to apply. So now we expand the k basis factor, or the label of the basis factor, in this binary expansion, right? So basically in this form. So now we go from 0 to 1 for all of these elements. So I’m expanding this sum into some n bits. And here are the n bits describing k. And I’m also expanding a little bit these parts, or k over 2n, in this binary form in this form. So I have two things going on in this step, these two. Then, when you look at just this form here, so whatever assumed exponential, that we can change it into a product. So now we can look at the product going from l equals 1 to n. And basically, I’m moving this summation out and make it into this product. And the next thing I’m going to do is I’m going to switch this product in front of the sum. So that’s what’s here. So now you have the product in front. You have the submission. kl, now here I got an index. That’s why you don’t have them explicitly written down. So kl goes from 0 to 1. And you still have this binary explanation here, with the basis vectors. And that didn’t change this, so this is the same, what you saw here. And now if we expand this very large sum, which is just 0 and 1, then we can read out what is the actual logarithm. Because what you have here is now for the 0 ket, 0 1 up for each of these kl, you have kl plus 0 here, which means that this whole phase is going to be 1. So it does not have a– well, its amplitude is 1. And then for the excited states, the 1 ket, we will have the actual corresponding value with this set to 1, which is just 2 pi i times x times 2 to the minus l, l going from 1 to m. And so this is equivalent to what we started with, and we can read the circuit from this one because this is just a couple of conditional rotations. So these will be the rotations that we want. And as you can see, the typical structure of a conditional operation, where you separate the two parts into two separate sections. So that’s it. In principle, it’s very simple, but it’s very resource intensive. You need lots of qubits, and eventually, you end up with a large number of gates, even though it is a simple structure. And therefore, you can’t really do it on this shallow circuit, shallow circuit, quantum computers that we have talked about so far. But you definitely need something that’s robust in noise and that’s large scale.

• State tomography is the procedure to acquire a complete classical description of what a quantum state is. Why is it hard to? Remember our metaphor from the very beginning of the course: a quantum state is a probability distribution with some special properties.

– Measuring the state only gives one sample, it does not reveal the entire probability distribution.

– |Since a measurement changes the state, you have to repeat the entire calculation again to obtain another sample.

• Which encoding does the quantum Fourier transformation use?

– Amplitude encoding.

• The derivation of the quantum algorithm is long, but eventually we find the U unitary that performs the desired transformation in the form $$\frac{1}{\sqrt{N}}(\mid0\rangle+e^{2\pi i0.x_{n}}\mid1\rangle)\otimes\ldots\otimes(\mid0\rangle+e^{2\pi i0.x_{1}.x_{2}\ldots x_{n}.x_{n-1}}\mid1\rangle$$. Why is this form advantageous?

– It reduces dependency between qubits: only the last qubit depends on the value of all the others, the remaining qubits depend less and less on the input qubits. #F

117. Introduction

The quantum Fourier transform and quantum phase estimation provide the foundation for many quantum algorithms, including the quantum matrix inversion, which is extensively used in quantum machine learning. It is therefore worthwhile developing a good understanding of these building blocks before moving on to more complex algorithms.

We must emphasize that starting with this notebook, the algorithms presented are coherent quantum protocols. By that, we mean that the input and output of an algorithm is a quantum state that we do not have classical information about. The protocol itself might use measurements: in this sense, they are not fully coherent, since we gain some, but incomplete classical information about the quantum system. We might also perform post-selection, which means that a gate is controlled by the classical output of a measurement. In some cases, we entirely discard a calculation based on a measurement output.

Why does it matter that we begin and end up with quantum states? Can’t we just use state preparation starting from classical data and then perform tomography on the final state? We could do that, but state preparation and tomography are resource-intensive, and they are likely to destroy any quantum advantage.

An additional problem is that the quantum Fourier transformation and other quantum algorithms similar in complexity require a very large number of gates on a large number of high-quality qubits. This is why the practical relevance of these algorithms is not immediate, but since they are the core of many quantum machine learning protocols, it is essential that we take a look at them.

import numpy as np
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit import execute
from qiskit import BasicAer as Aer
π = np.pi

118. Quantum Fourier Transform

The quantum Fourier transform is a quantum algorithm for the discrete Fourier transform over the amplitudes of a wavefunction. The exposition here follows the introduction in [1]. A similar approach can be found in the Qiskit tutorials.

The classical discrete Fourier transform acts on a vector $\vec{x}=\begin{bmatrix}x_0\\ \vdots\\ x_{N-1}\end{bmatrix}$ and maps it to the vector $\vec{y}=\begin{bmatrix}y_0\\ \vdots\\ y_{N-1}\end{bmatrix}$, where $y_k = \frac{1}{\sqrt{N}}\sum_{j=0}^{N-1}x_je^{ \boldsymbol{2\pi i} \frac{jk}{N}}$.

The quantum Fourier transform acts on an amplitude-encoded variant of this vector, the quantum state $|x\rangle=\sum_{i=0}^{N-1} x_i |i \rangle$ and maps it to the quantum state $|y\rangle=\sum_{k=0}^{N-1} y_k |k \rangle$, where $y_k = \frac{1}{\sqrt{N}}\sum_{j=0}^{N-1}x_je^{\boldsymbol{2\pi i} \frac{jk}{N}}$. Since the transformed state is also in the superposition of the computational basis, in essence, only the amplitudes are transformed.

To derive a circuit for a power of two $N=2^n$, consider the transform acts on the state $| x \rangle = | x_1...x_n \rangle$ where $x_1$ is the most significant bit, unlike the usual convention followed in the course. We will rewrite a number $y$ in the fractional binary notation as $j = 0.j_1...j_n = \sum_{k=1}^n j_k/2^k$. For example $0.8125_d = 0.1101_b = \sum_{k=1}^4 j_k/2^k = 1/2 + 1/4 + 0/8 + 1/16$, where the subscripts $d$ stands for decimal and $b$ for binary. The action of the unitary $U$ describing the transform can be expanded as \begin{aligned} U |x \rangle = U |x_1 x_2 \cdots x_n \rangle& = \frac{1}{\sqrt{N}} \sum_{k=0}^{2^n-1} e^{\boldsymbol{2\pi i} xk / 2^n} |k \rangle \ & = \frac{1}{\sqrt{N}} \sum_{k_1=0}^{1}\ldots\sum_{k_n=0}^{1} e^{\boldsymbol{2\pi i} x\left(\sum_{l=1}^n k_l2^{-l}\right) } \vert k_1 … k_n \rangle \ & = \frac{1}{\sqrt{N}} \sum_{k_1=0}^{1}\ldots\sum_{k_n=0}^{1} \bigotimes_{l=1}^n e^{\boldsymbol{2\pi i} x k_l2^{-l}} | k_1 … k_n \rangle \ & = \frac{1}{\sqrt{N}} \bigotimes_{l=1}^n \sum_{k_l=0}^{1} e^{\boldsymbol{2\pi i} x k_l2^{-l}} | k_1 … k_n \rangle \ & = \frac{1}{\sqrt{N}} \bigotimes_{l=1}^n \left(|0\rangle + e^{\boldsymbol{2\pi i} x 2^{-l} } |1\rangle \right) \ & = \frac{1}{\sqrt{N}} \left(|0\rangle + e^{\boldsymbol{2\pi i} 0.x_n} |1\rangle\right) \otimes…\otimes \left(\vert0\rangle + e^{\boldsymbol{2\pi i} 0.x_1.x_2…x_{n-1}.x_n} |1\rangle\right) \end{aligned}

This form of the QFT is useful for deriving a circuit, since only the last qubit depends on the the values of all the other input qubits. The remaining qubits depend less and less on the input qubits. The simple structure also allows to decompose the unitary as Hadamard gates and rotations. On three qubits, we can define the circuit as follows:

q = QuantumRegister(3, 'q')
c = ClassicalRegister(1, 'c')
qft = QuantumCircuit(q, c)
qft.h(q[0])
qft.cu1(π/2, q[1], q[0])
qft.h(q[1])
qft.cu1(π/4, q[2], q[0])
qft.cu1(π/2, q[2], q[1])
qft.h(q[2])

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[2], line 5
      3 qft = QuantumCircuit(q, c)
      4 qft.h(q[0])
----> 5 qft.cu1(π/2, q[1], q[0])
      6 qft.h(q[1])
      7 qft.cu1(π/4, q[2], q[0])

AttributeError: 'QuantumCircuit' object has no attribute 'cu1'

We can plot it to make the qubit dependencies more apparent:

from qiskit.tools.visualization import circuit_drawer
circuit_drawer(qft)

        ┌───┐                                     
q_0: |0>┤ H ├─■────────■──────────────────────────
        └───┘ │1.5708  │       ┌───┐              
q_1: |0>──────■────────┼───────┤ H ├─■────────────
                       │0.7854 └───┘ │1.5708 ┌───┐
q_2: |0>───────────────■─────────────■───────┤ H ├
                                             └───┘
 c_0: 0 ══════════════════════════════════════════

The conditional rotations dominate the complexity, which scales as $O(N^2)$.

119. Even more Quantum Phase Estimation

The quantum Fourier transformation is a very important algorithmic primitive in quantum computing. It is the foundation of many other quantum protocols, including quantum phase estimation, in which we want to estimate the eigenvalues of some unitary operator. So what happens is that given a unitary operator and one of its eigenvectors, then this is going to be one of its eigenvalues. But since this is a unitary operator, all of its eigenvalues are going to be on the complex unit circle. In other words, we can write it in this form. So this is e to the 2 pi i times some theta angle around the unit circle– the complex unit circle. And it is this phase– this theta phase is what we are going to estimate to some finite precision. And the way it works– it has two parts. One is an inverse quantum Fourier transformation. That’s easy. What comes before that– that’s a little bit challenging to understand. So imagine that you have two sets of registers. You have n registers which define the precision of your estimation. So this is where the estimated eigenvalue’s going to be. And then you have the actual eigenvector coming in some m registers. And what we are going to do is we are going to apply the unitary operator over and over again in these controlled operations. But in each of these applications, the unitary will be applied for a different duration. So in the first application, we apply just the unitary itself. This is the unitary to the power of 2 to the 0. And then it’s the unitary to the power of 2 the first, and so on all the way to the unitary applied to 2 to the n minus 1 times. So we create the equal superposition of the 0 cat in this ancilla register where the eigenvalue’s going to be written. And each of these controlled unitaries is going to be applied in superposition so that at the end of this procedure, this is the state that you’re going to get. This is going to be your first qubit. So here, you have the unitary applied 2 to the n minus 1 times– so it’s this one– and all the way to the very last n-th qubit in your ancilla where you are estimating your eigenvalue. And if you look at this– look at the structure of this tensor product, you can actually rewrite it in this form. It’s just a very simple sum. And what’s interesting here is that– OK, so you have your basis. You are expanding in the canonical basis. And that you pick up a global phase. So, sure, you also have your original eigenvector here, but this phrase is global to both of them. So this is called the phase kickback. And it’s a common trick to use it for something meaningful. In this case, when you look at this form, it looks like something that you would get after the Fourier transformation of some state. And now, if we apply the inverse quantum Fourier transformation, then we can write back this global phase into the register, and hence estimating the eigenvalue that we were interested in the first place.

• Why can we write the eigenvalue of a unitary mtarix as a phase $\theta$ in $$e^{2i\pi\theta}$$

– The eigenvalues of a unitary matrix always lie on the complex unit circle.

• Why do we write $\mid\hat{\theta}\rangle$ instead of $\mid\theta\rangle$?

– $\mid\hat{\theta}\rangle$ is a finite-bit representation of a continuous valued number.

• What do we achieve by the series of controlled unitary applications before the inverse quantum Fourier transformation?

– We create a superposition of U as it is applied for different durations.

120. Quantum phase estimation

The goal of a quantum phase estimation algorithm is, given a unitary operator $U$ and an eigenvector $|\psi\rangle$ of $U$, to estimate $\theta$ in $U|\psi \rangle =e^{2 i \pi \theta}|\psi \rangle$. Since $U$ is unitary, all of its eigenvalues have an absolute value of 1. By convention, $\theta$ is taken to be in $[0,1]$ and is called the phase of $U$ associated to $|\psi\rangle$.

The eigenvector $|\psi\rangle$ is encoded in one set of quantum registers. An additional set of $n$ qubits forms an ancilla register. At the end of the procedure, this ancilla register should contain an approximation of the binary fraction associated to $\theta$, with n-bits precision. A critical element is the ability to perform the controlled unitary $C-U^{2^k}$ – it is usually assumed to be provided to the phase estimation protocol.

First, the uniform superposition is prepared in the ancilla register via the application of Hadamard gates $H$. These qubits will act as controls for the unitary operators at different time steps. Our goal is to create a superposition of $U$ as the unitary is applied for different durations. Since the eigenvalues are always situated on the complex unit circle, these differently evolved components in the superposition help reveal the eigenstructure. Given that the ancilla register we have a superposition of all possible time steps between $0$ and $2^{n-1}$, we will end up with a superposition of all possible evolutions to encode binary representations of the eigenvalues. At the end of this procedure, we have the state \begin{aligned} & \frac {1}{2^{\frac {n}{2}}} (|0\rangle+{e^{2 i \pi \theta \cdot 2^{n-1}}}|1\rangle ) \otimes \cdots (|0\rangle+{e^{2 i \pi \theta \cdot 2^{1}}}|1\rangle ) \otimes (|0\rangle+{e^{2i \pi \theta \cdot 2^{0}}}|1\rangle ) = \ &\frac {1}{2^{\frac {n}{2}}}\sum _{k=0}^{2^{n}-1}e^{2 i \pi \theta k}|k\rangle \end{aligned} in the ancilla. To write the ancilla in this form, we exploit that the controlled unitary operations when applied, introduce a global phase, and it is this global phase that we see in the ancilla. This phenomenon is also known as the phase kickback.

As a final step, we apply an inverse Fourier transform on the ancilla. Measuring out in the computational basis, we get the phase in the ancilla register:

\begin{align} \frac {1}{2^{\frac {n}{2}}}\sum _{k=0}^{2^{n}-1}e^{2i \pi \theta k}|k\rangle \otimes | \psi \rangle \xrightarrow{\mathcal{QFT}_n^{-1}} | \tilde{\theta} \rangle \otimes | \psi \rangle \end{align}

where $\tilde{\theta}$ is the n-bits approximation of the binary fraction representing $\theta$.

The circuit for phase estimation is the following:

Quantum phase estimation

As a toy example, let’s take the $2\times 2$ unitary matrix $\begin{bmatrix}e^{0} & 0 \\0 & e^{i \pi}\end{bmatrix}=\begin{bmatrix}1 & 0 \\0 & -1\end{bmatrix}$, which has the eigenvectors $|0\rangle$ and $|1\rangle$, and phases $\theta_0=0$ and $\theta_1=\frac{1}{2}$. Therefore, the $C-U^{2^k}$ gate is a controlled-$Z$ gate for $k=0$ and the identity for $k\geq 1$.

Starting with $|\psi\rangle=|0\rangle$ in the main register, we prepare the superposition in the ancilla:

q = QuantumRegister(3, 'q')
c = ClassicalRegister(2, 'c')

qpe = QuantumCircuit(q, c)
qpe.h(q[0])
qpe.h(q[1])

<qiskit.extensions.standard.h.HGate at 0x7f3c8c3f1e80>

Next we perform the controlled unitary operations:

# Controlled-U0
qpe.cz(q[1], q[2])
# Controlled-U1
# nothing: identity

<qiskit.extensions.standard.cz.CzGate at 0x7f3c8c3f1860>

We apply quantum inverse Fourier transformation to write the phase to the ancilla register:

qpe.swap(q[0], q[1])
qpe.h(q[1])
qpe.cu1(-π / 2, q[0], q[1])
qpe.h(q[0])
qpe.swap(q[0], q[1])

<qiskit.extensions.standard.swap.SwapGate at 0x7f3c8c3b5550>

We will get the result from the two first registers

qpe.measure(q[0], c[0])
qpe.measure(q[1], c[1])

<qiskit.circuit.measure.Measure at 0x7f3c8c3f1be0>

We can plot the circuit:

circuit_drawer(qpe)

        ┌───┐                     ┌───┐   ┌─┐   
q_0: |0>┤ H ├────X───────■────────┤ H ├─X─┤M├───
        ├───┤    │ ┌───┐ │-1.5708 └───┘ │ └╥┘┌─┐
q_1: |0>┤ H ├─■──X─┤ H ├─■──────────────X──╫─┤M├
        └───┘ │    └───┘                   ║ └╥┘
q_2: |0>──────■────────────────────────────╫──╫─
                                           ║  ║ 
 c_0: 0 ═══════════════════════════════════╩══╬═
                                              ║ 
 c_1: 0 ══════════════════════════════════════╩═

Let’s now test our circuit:

backend = Aer.get_backend('qasm_simulator')
job = execute(qpe, backend, shots=1000)
result = job.result()
result.get_counts(qpe)

{'00': 1000}

As expected the result is $|2 \cdot \theta_0\rangle=|2\cdot 0\rangle=|00\rangle$

Let’s now run the circuit for the eigenvector $|1\rangle$:

qpe = QuantumCircuit(q, c)
qpe.h(q[0])
qpe.h(q[1])
qpe.x(q[2]) # create |1> in the main register

qpe.cz(q[1], q[2])

qpe.swap(q[0], q[1])
qpe.h(q[1])
qpe.cu1(-π / 2, q[0], q[1])
qpe.h(q[0])
qpe.swap(q[0], q[1])

qpe.measure(q[0], c[0])
qpe.measure(q[1], c[1])

<qiskit.circuit.measure.Measure at 0x7f3c8c3f1e10>

backend = Aer.get_backend('qasm_simulator')
job = execute(qpe, backend, shots=1000)
result = job.result()
result.get_counts(qpe)

{'10': 1000}

The result should be $|10\rangle$. Indeed, $10 \rightarrow 1\cdot 2^{-1} + 0 \cdot 2^0=\frac{1}{2}=\theta_1$

121. References

[1] M. Nielsen, I. Chuang. (2000). Quantum Computation and Quantum Information. Cambridge University Press.