# Phase Estimation¶

**Introduction**

**Problem Description**

The quantum phase estimation algorithm is a quantum algorithm used to estimate the phase (or eigenvalue) of an eigenvector of a unitary operator. More precisely, given a unitary matrix \(U\) and a quantum state \(\vert\psi\rangle\) such that \(U\vert\psi\rangle=e^{2\pi i\theta}\vert\psi\rangle\) (that is, \(\vert\psi\rangle\) is an eigenstate of \(U\)) and \(\theta\in[0,1)\), the algorithm estimates the value of \(\theta\) with high probability within additive error \(\varepsilon\), using \(O(1/\varepsilon)\) controlled-\(U\) operations.

*Note:* We can write the eigenvalue of \(U\) in the form
\(e^{2\pi i\theta}\) because \(U\) is a unitary operator over a
complex vector space, so its eigenvalues must be complex numbers with
absolute value \(1\).

**Blackbox (aka. Oracles)**

To perform the estimation we assume that we have available blackbox capable of preparing the state \(\vert\psi\rangle\) and performing the controlled-\(U^{2^j}\) operation for suitable \(j\in\mathbb{Z}^+\).

**Input**

The input consists of two quantum registers:

First, the upper \(n\) qubits comprise the first register. This register is used to store the value of \(\left\lfloor 2^n\theta \right\rfloor\). If we express \(\theta\) in the binary form: \(\theta=0. \theta_0\theta_1\cdots\), then \(\left\lfloor 2^n\theta \right\rfloor = \theta_0\theta_1\cdots\theta_{n-1}\). The \(n\) qubits stores the \(n\) binaries \(\theta_i\). The choice of \(n\) depends on two things: the number of digits of accuracy we wish to have for the estimation of \(\theta\), and with what probability we wish the estimation algorithm to be successful. There is a tradeoff between accuracy and probability.

Second, the lower \(m\) qubits are the second register, which stores the input state \(\vert\psi\rangle\). It keeps unchanged throughout the algorithm.

**The Circuit**

The quantum circuit for phase estimation is shown as follows (excerpted from Wikipedia)

We analyze the circuit step by step. We use \(\vert\Psi_t\rangle\) to represent the state of the system after \(t\) steps. The initial system state is given by

*Step 1. Create a superposition.*

As most quantum algorithms do, we put the initial state of the first
register into a superposition. After *Step 1*, the system state becomes

Here we assume \(t\) is expressed in its binary form.

*Step 2. Apply controlled unitary operations.*

It is easy to verify that

Assume the binary expansion of \(t = t_0\cdots t_{n-1}\), each \(t_i\in\{0,1\}\), we have

After *Step 2*, the system state becomes

As we can see, by applying controlled unitaries, we map \(\theta\) to the amplitude of the state.

*Step 3. Apply inverse Quantum Fourier transform.*

Further analyzing \(\vert\Psi_2\rangle\), we find that it can be obtained from state \(\vert\theta\rangle\vert\psi\rangle\) by performing a quantum Fourier transform on the first quantum register (with some error). Now we perform a inverse quantum Fourier transform on \(\vert\Psi_2\rangle\) (particularly, on state \(\vert t\rangle\)), we obtain

Without loss of generality, we can assume \(2^n\theta = a + 2^n\delta\), where \(a\) is the nearest integer to \(2^n\theta\). The difference \(2^n\delta\) satisfies \(0 \leq 2^n\delta \leq 1/2\) (guaranteed by the fact that \(a\) is the nearest integer to \(2^n\theta\)). With this expansion, \(\vert\Psi_3\rangle\) can be rewritten as

*Step 4. Measurement.*

Now we perform a measurement in the standard basis on the first register. We obtain outcome \(\vert a\rangle\) with probability

**Case 1. :math:`delta = 0`.**In this case, the approximation is precise, we always measure the accurate value of the phase: \(\theta = a/2^n\). The state of the system after the measurement is \(\vert a \rangle\otimes\vert \psi \rangle\).**Case 2. :math:`delta neq 0`.**In this case, we have to lower bound the probability \(\operatorname{Pr}\left(\vert a \rangle\right)\). Since \(2^n\delta\leq 1/2\), the algorithm yields the correct result with probability\[\operatorname{Pr}\left(\vert a \rangle\right) \geq \frac{4}{\pi^2} \approx 0.405\]This result shows that we can measure the best \(n\)-bit estimate of \(\theta\) with high probability. By increasing the amount of qubits \(n\) by \(O(\log(1/\epsilon))\) and ignoring those last qubits we can increase the probability to \(1-\epsilon\).

**HiQ Implementation**

Based on the **Introduction**, we need to be able to implement the
following operations:

First, controlled unitary operations. This can be done on the HiQ using
the `Control`

meta.

Second, unitary operation raised to some power. This can be done on the
HiQ using the `Loop`

meta.

Third, inversed QFT. This can be done on the HiQ by calling
`get_inverse(QFT)`

, as `QFT`

is already a built-in operation.

We first write a function to calculate the number of bits required to estimate \(\theta\) accurate to n bits with probability of success at least \(1 - \epsilon\). The analytic expression can be found in Nielsen Book, Eq. (5.35):

1 2 3 4 5 | ```
def _bits_required_to_achieve_accuracy(n, epsilon):
tmp = math.log2(2 + 1/(2*epsilon))
tmp = math.ceil(tmp)
n_updated = int(n + tmp)
return n_updated
``` |

Now we write the main routine used to implement *Steps 1-3*. We remark
that *Step 4* is not necessary if we use phase estimation as a
subroutine for other algorithms, and this is why we do not count it in
the main routine.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 | ```
def run_phase_estimation(eng, U, state, m, n):
"""
Phase estimation algorithm on unitary U.
Args:
eng (MainEngine): Main compiler engine to run the phase_estimation algorithm
U (np.matrix): the unitary that is to be estimated
state (qureg): the input state, which has the form |psi> \otimes |0>.
Similar to the Introduction, state is composed of two parts:
1. The first m qubits are the input state |psi>
2. The last n qubits are initialized to |0>, which is used to
store the binary digits of the estimated phase
m (int): number of qubits for input state psi
n (int): number of qubits used to store the phase to given accuracy
"""
# The beginning index for the phase qubits will have an offset
OFFSET = m
# Step 1 and 2. Create a superposition and perform the controlled U^{2^j} operations
for k in range(n):
H | state[OFFSET+k]
with Control(eng, state[OFFSET+k]):
# number of loops required to perform the controlled U^{2^j} operation
num = int(math.pow(2, k))
with Loop(eng, num):
U | state[:OFFSET]
# Step 3. Perform the inverse QFT operation
# Step 4. Swap the qubits
for k in range(math.floor(n/2)):
Swap | (state[OFFSET+k], state[OFFSET+n-k-1])
# Step 5. Perform the original inverse QFT
get_inverse(QFT) | state[OFFSET:]
``` |

*Remark:* The default HiQ implementation of `QFT`

ignores the `SWAP`

operations required (see the caption of Figure 5.1 in Nielsen Book).
Here we must explicitly call the `SWAP`

operations to swap the qubits
and then perform the inversed QFT.

After the `run_phase_estimation`

function, the estimated phase is
stored in the last n qubits of the `state`

variable. We can either
measure it (This is *Step 4* in introduction), or we can use it as the
input to other routines.

**References:**