Quantum phase estimation is one of the most important subroutines in quantum computation. It serves as a central building block for many quantum algorithms. The objective of the algorithm is the following:

Given a unitary operator $U$, the algorithm estimates $\theta$ in $U\vert\psi \rangle =e^{\boldsymbol{2\pi i} \theta }|\psi \rangle$. Here $|\psi\rangle$ is an eigenvector and $e^{\boldsymbol{2\pi i}\theta}$ is the corresponding eigenvalue. Since $U$ is unitary, all of its eigenvalues have a norm of 1.

The general quantum circuit for phase estimation is shown below. The top register contains $t$ 'counting' qubits, and the bottom contains qubits in the state $|\psi\rangle$:

### 1.1 Intuition

The quantum phase estimation algorithm uses phase kickback to write the phase of $U$ (in the Fourier basis) to the $t$ qubits in the counting register. We then use the inverse QFT to translate this from the Fourier basis into the computational basis, which we can measure.

We remember (from the QFT chapter) that in the Fourier basis the topmost qubit completes one full rotation when counting between $0$ and $2^t$. To count to a number, $x$ between $0$ and $2^t$, we rotate this qubit by $\tfrac{x}{2^t}$ around the z-axis. For the next qubit we rotate by $\tfrac{2x}{2^t}$, then $\tfrac{4x}{2^t}$ for the third qubit.

When we use a qubit to control the $U$-gate, the qubit will turn (due to kickback) proportionally to the phase $e^{2i\pi\theta}$. We can use successive $CU$-gates to repeat this rotation an appropriate number of times until we have encoded the phase theta as a number between $0$ and $2^t$ in the Fourier basis.

Then we simply use $QFT^\dagger$ to convert this into the computational basis.

### 1.2 Mathematical Basis

As mentioned above, this circuit estimates the phase of a unitary operator $U$. It estimates $\theta$ in $U\vert\psi \rangle =e^{\boldsymbol{2\pi i} \theta }|\psi \rangle$, where $|\psi\rangle$ is an eigenvector and $e^{\boldsymbol{2\pi i}\theta}$ is the corresponding eigenvalue. The circuit operates in the following steps:

**Setup**: $\vert\psi\rangle$ is in one set of qubit registers. An additional set of $n$ qubits form the counting register on which we will store the value $2^n\theta$:

**Superposition**: Apply a $n$-bit Hadamard gate operation $H^{\otimes n}$ on the counting register:

**Controlled Unitary Operations**: We need to introduce the controlled unitary $C-U$ that applies the unitary operator $U$ on the target register only if its corresponding control bit is $|1\rangle$. Since $U$ is a unitary operatory with eigenvector $|\psi\rangle$ such that $U|\psi \rangle =e^{\boldsymbol{2\pi i} \theta }|\psi \rangle$, this means:

Applying all the $n$ controlled operations $C − U^{2^j}$ with $0\leq j\leq n-1$, and using the relation $|0\rangle \otimes |\psi \rangle +|1\rangle \otimes e^{2\pi i\theta }|\psi \rangle =\left(|0\rangle +e^{2\pi i\theta }|1\rangle \right)\otimes |\psi \rangle$:

\begin{aligned} \psi_{2} & =\frac {1}{2^{\frac {n}{2}}} \left(|0\rangle+{e^{\boldsymbol{2\pi i} \theta 2^{n-1}}}|1\rangle \right) \otimes \cdots \otimes \left(|0\rangle+{e^{\boldsymbol{2\pi i} \theta 2^{1}}}\vert1\rangle \right) \otimes \left(|0\rangle+{e^{\boldsymbol{2\pi i} \theta 2^{0}}}\vert1\rangle \right) \otimes |\psi\rangle\\\\ & = \frac{1}{2^{\frac {n}{2}}}\sum _{k=0}^{2^{n}-1}e^{\boldsymbol{2\pi i} \theta k}|k\rangle \otimes \vert\psi\rangle \end{aligned}where $k$ denotes the integer representation of n-bit binary numbers.

**Inverse Fourier Transform**: Notice that the above expression is exactly the result of applying a quantum Fourier transform as we derived in the notebook on Quantum Fourier Transform and its Qiskit Implementation. Recall that QFT maps an n-qubit input state $\vert x\rangle$ into an output as

Replacing $x$ by $2^n\theta$ in the above expression gives exactly the expression derived in step 2 above. Therefore, to recover the state $\vert2^n\theta\rangle$, apply an inverse Fourier transform on the ancilla register. Doing so, we find

$$ \vert\psi_3\rangle = \frac {1}{2^{\frac {n}{2}}}\sum _{k=0}^{2^{n}-1}e^{\boldsymbol{2\pi i} \theta k}|k\rangle \otimes | \psi \rangle \xrightarrow{\mathcal{QFT}_n^{-1}} \frac {1}{2^n}\sum _{x=0}^{2^{n}-1}\sum _{k=0}^{2^{n}-1} e^{-\frac{2\pi i k}{2^n}(x - 2^n \theta)} |x\rangle \otimes |\psi\rangle $$**Measurement**: The above expression peaks near $x = 2^n\theta$. For the case when $2^n\theta$ is an integer, measuring in the computational basis gives the phase in the ancilla register with high probability:

For the case when $2^n\theta$ is not an integer, it can be shown that the above expression still peaks near $x = 2^n\theta$ with probability better than $4/\pi^2 \approx 40\%$ [1].

Let’s take a gate we know well, the $T$-gate, and use Quantum Phase Estimation to estimate its phase. You will remember that the $T$-gate adds a phase of $e^\frac{i\pi}{4}$ to the state $|1\rangle$:

$$ T|1\rangle = \begin{bmatrix} 1 & 0\\ 0 & e^\frac{i\pi}{4}\\ \end{bmatrix} \begin{bmatrix} 0\\ 1\\ \end{bmatrix} = e^\frac{i\pi}{4}|1\rangle $$Since QPE will give us $\theta$ where:

$$ T|1\rangle = e^{2i\pi\theta}|1\rangle $$We expect to find:

$$\theta = \frac{1}{8}$$In this example we will use three qubits and obtain an *exact* result (not an estimation!)

```
#initialization
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format = 'svg' # Makes the images look nice
import numpy as np
import math
# importing Qiskit
from qiskit import IBMQ, Aer
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister, execute
# import basic plot tools
from qiskit.visualization import plot_histogram
```

Now, set up the quantum circuit. We will use four qubits -- qubits 0 to 2 as counting qubits, and qubit 3 as the eigenstate of the unitary operator ($T$).

We initialize $\vert\psi\rangle = \vert1\rangle$ by applying an $X$ gate:

```
qpe = QuantumCircuit(4, 3)
qpe.x(3)
qpe.draw(output='mpl')
```

Next, we apply Hadamard gates to the counting qubits:

```
for qubit in range(3):
qpe.h(qubit)
qpe.draw(output='mpl')
```

Next we perform the controlled unitary operations:

```
repetitions = 2**2
for counting_qubit in range(3):
for i in range(repetitions):
qpe.cu1(math.pi/4, counting_qubit, 3); # This is C-U
repetitions //= 2
qpe.draw(output='mpl')
```

We apply the inverse quantum Fourier transformation to convert the state of the counting register. Here we provide the code for $QFT^\dagger$:

```
def qft_dagger(circ, n):
"""n-qubit QFTdagger the first n qubits in circ"""
# Don't forget the Swaps!
for qubit in range(int(n/2)):
circ.swap(qubit, n-qubit-1)
for j in range(n,0,-1):
k = n - j
for m in range(k):
circ.cu1(-math.pi/float(2**(k-m)), n-m-1, n-k-1)
circ.h(n-k-1)
```

We then measure the counting register. At the moment our qubits are in reverse order (a common problem in quantum computing!) We measure to the classical bits in reverse order to fix this:

```
# Apply inverse QFT
qft_dagger(qpe, 3)
# We measure in reverse order to correct issues later
qpe.measure(0,2)
qpe.measure(1,1)
qpe.measure(2,0)
```

```
qpe.draw(output="mpl")
```

```
backend = Aer.get_backend('qasm_simulator')
shots = 2048
results = execute(qpe, backend=backend, shots=shots).result()
answer = results.get_counts()
plot_histogram(answer)
```

We see we get one result (`001`

) with certainty, which translates to the decimal: `1`

. We now need to divide our result (`1`

) by $2^n$ to get $\theta$:

This is exactly the result we expected!

```
# Create and set up circuit
qpe2 = QuantumCircuit(4, 3)
# Apply H-Gates to counting qubits:
for qubit in range(3):
qpe2.h(qubit)
# Prepare our eigenstate |psi>:
qpe2.x(3)
# Do the controlled-U operations:
angle = 2*math.pi/3
repetitions = 2**2
for counting_qubit in range(3):
for i in range(repetitions):
qpe2.cu1(angle, counting_qubit, 3);
repetitions //= 2
# Do the inverse QFT:
qft_dagger(qpe2, 3)
# Measure of course!
qpe2.measure(0,2)
qpe2.measure(1,1)
qpe2.measure(2,0)
qpe2.draw(output='mpl')
```

```
# Let's see the results!
backend = Aer.get_backend('qasm_simulator')
shots = 4096
results = execute(qpe2, backend=backend, shots=shots).result()
answer = results.get_counts()
plot_histogram(answer)
```

We are expecting the result $\theta = 0.3333\dots$, and we see our most likely results are `010(bin) = 2(dec)`

and `011(bin) = 3(dec)`

. These two results would tell us that $\theta = 0.25$ (off by 25%) and $\theta = 0.375$ (off by 13%) respectively. The true value of $\theta$ lies between the values we can get from our counting bits, and this gives us uncertainty and imprecision.

### 3.2 The Solution

To get more precision we simply add more counting qubits. We are going to add two more counting qubits:

```
# Create and set up circuit
qpe3 = QuantumCircuit(6, 5)
# Apply H-Gates to counting qubits:
for qubit in range(5):
qpe3.h(qubit)
# Prepare our eigenstate |psi>:
qpe3.x(5)
# Do the controlled-U operations:
angle = 2*math.pi/3
repetitions = 2**4
for counting_qubit in range(5):
for i in range(repetitions):
qpe3.cu1(angle, counting_qubit, 5);
repetitions //= 2
# Do the inverse QFT:
qft_dagger(qpe3, 5)
# Measure of course!
qpe3.measure(0,4)
qpe3.measure(1,3)
qpe3.measure(2,2)
qpe3.measure(3,1)
qpe3.measure(4,0)
qpe3.draw(output='mpl')
```

```
# Let's see the results!
backend = Aer.get_backend('qasm_simulator')
shots = 4096
results = execute(qpe3, backend=backend, shots=shots).result()
answer = results.get_counts()
plot_histogram(answer)
```

The two most likely measurements are now `01011`

(decimal 11) and `01010`

(decimal 10). Measuring these results would tell us $\theta$ is:

These two results differ from $\frac{1}{3}$ by 3% and 6% respectively. A much better precision!

```
qpe.draw(output='mpl')
```

```
# Load our saved IBMQ accounts and get the least busy backend device with less than or equal to n qubits
IBMQ.load_account()
from qiskit.providers.ibmq import least_busy
from qiskit.tools.monitor import job_monitor
provider = IBMQ.get_provider(hub='ibm-q')
backend = least_busy(provider.backends(filters=lambda x: x.configuration().n_qubits >= 4 and not x.configuration().simulator and x.status().operational==True))
print("least busy backend: ", backend)
# Run with 2048 shots
shots = 2048
job_exp = execute(qpe, backend=backend, shots=shots, optimization_level=3)
job_monitor(job_exp)
```

```
# get the results from the computation
results = job_exp.result()
answer = results.get_counts(qpe)
plot_histogram(answer)
```

We can hopefully see that the most likely result is `001`

which is the result we would expect from the simulator. Unlike the simulator, there is a probability of measuring something other than `001`

, this is due to noise and gate errors in the quantum computer.

## 6. Looking Forward

The quantum phase estimation algorithm may seem pointless, since we have to know $\theta$ to perform the controlled-$U$ operations on our quantum computer. We will see in later chapters that it is possible to create circuits for which we don’t know $\theta$, and for which learning theta can tell us something very useful (most famously how to factor a number!)

```
import qiskit
qiskit.__qiskit_version__
```