In quantum computing, the quantum phase estimation algorithm is a quantum algorithm to estimate the phase corresponding to an eigenvalue of a given unitary operator. Because the eigenvalues of a unitary operator always have unit modulus, they are characterized by their phase, and therefore the algorithm can be equivalently described as retrieving either the phase or the eigenvalue itself. The algorithm was initially introduced by Alexei Kitaev in 1995.^{[1][2]: 246}

Phase estimation is frequently used as a subroutine in other quantum algorithms, such as Shor's algorithm,^[2]: 131 the quantum algorithm for linear systems of equations, and the quantum counting algorithm.

Formal description of the problem

Let ${\displaystyle U}$ be a unitary operator acting on an ${\displaystyle m}$-qubit register. Unitarity implies that all the eigenvalues of ${\displaystyle U}$ have unit modulus, and can therefore be characterized by their phase. Thus if ${\displaystyle |\psi \rangle }$ is an eigenvector of ${\displaystyle U}$, then ${\displaystyle U|\psi \rangle =e^{2\pi i\theta }\left|\psi \right\rangle }$ for some ${\displaystyle \theta \in \mathbb {R} }$. Due to the periodicity of the complex exponential, we can always assume ${\displaystyle 0\leq \theta <1}$.

Our goal is to find a good approximation to ${\displaystyle \theta }$ with a small number of gates and with high probability. The quantum phase estimation algorithm achieves this under the assumptions of having oracular access to ${\displaystyle U}$, and having ${\displaystyle |\psi \rangle }$ available as a quantum state.

More precisely, the algorithm returns an approximation for ${\displaystyle \theta }$, with high probability within additive error ${\displaystyle \varepsilon }$, using ${\displaystyle O(\log(1/\varepsilon ))}$ qubits (without counting the ones used to encode the eigenvector state) and ${\displaystyle O(1/\varepsilon )}$ controlled-U operations. Furthermore, we can improve the success probability to ${\displaystyle 1-\Delta }$ for any ${\displaystyle \Delta >0}$ by using a total of ${\displaystyle O(\log(1/\Delta )/\varepsilon )}$ uses of controlled-U, and this is optimal.^[3]

The algorithm

Setup

The input consists of two registers (namely, two parts): the upper ${\displaystyle n}$ qubits comprise the first register, and the lower ${\displaystyle m}$ qubits are the second register.

The initial state of the system is:

${\displaystyle |0\rangle ^{\otimes n}|\psi \rangle .}$

After applying n-bit Hadamard gate operation ${\displaystyle H^{\otimes n))$ on the first register, the state becomes:

${\displaystyle {\frac {1}{2^{\frac {n}{2))))(|0\rangle +|1\rangle )^{\otimes n}|\psi \rangle ={\frac {1}{2^{n/2))}\sum _{j=0}^{2^{n}-1}|j\rangle |\psi \rangle }$.

Let ${\displaystyle U}$ be a unitary operator with eigenvector ${\displaystyle |\psi \rangle }$ such that ${\displaystyle U|\psi \rangle =e^{2\pi i\theta }|\psi \rangle }$. Thus,

${\displaystyle U^{2^{j))|\psi \rangle =e^{2\pi i2^{j}\theta }|\psi \rangle }$.

Overall, the transformation implemented on the two registers by the controlled gates applying ${\displaystyle U,U^{2},U^{2^{2)),\ldots ,U^{2^{n-1))}$ is

$${\displaystyle |k\rangle |\psi \rangle \mapsto |k\rangle U^{k}|\psi \rangle }$$

${\displaystyle k}$

${\displaystyle k_{n-1}k_{n-2}\ldots k_{1}k_{0))$

${\displaystyle 2^{n-1}k_{n-1}+2^{n-2}k_{n-2}+\ldots +2k_{1}+k_{0))$

${\displaystyle k_{n-1},\ldots ,k_{1},k_{0}\in \{0,1\))$

${\displaystyle U^{k))$

$${\displaystyle U^{2^{n-1}k_{n-1}+\ldots +2^{2}k_{2}+2k_{1}+k_{0))=U^{2^{n-1}k_{n-1))\ldots U^{2^{2}k_{2))U^{2k_{1))U^{k_{0))}$$

${\displaystyle U^{2^{j}k_{j))}$

${\displaystyle k_{j))$

${\displaystyle 1}$

${\displaystyle |k\rangle U^{k}|\psi \rangle }$

${\displaystyle U^{2^{j))}$

${\displaystyle j}$

Therefore, the state will be transformed by the controlled ${\displaystyle U^{2^{j))}$ gates like so:

$${\displaystyle {\frac {1}{\sqrt {2^{n))))\sum _{j=0}^{2^{n}-1}|j\rangle |\psi \rangle \mapsto {\frac {1}{\sqrt {2^{n))))\sum _{j=0}^{2^{n}-1}e^{2\pi ij\theta }|j\rangle |\psi \rangle }$$

${\displaystyle |\psi \rangle }$

${\displaystyle {\frac {1}{\sqrt {2^{n))))\sum _{j=0}^{2^{n}-1}e^{2\pi ij\theta }|j\rangle }$

Apply inverse quantum Fourier transform

Applying the inverse quantum Fourier transform on

${\displaystyle {\frac {1}{2^{\frac {n}{2))))\sum _{k=0}^{2^{n}-1}e^{2\pi i\theta k}|k\rangle }$

yields

${\displaystyle {\frac {1}{2^{\frac {n}{2))))\sum _{k=0}^{2^{n}-1}e^{2\pi i\theta k}\left({\frac {1}{2^{\frac {n}{2))))\sum _{x=0}^{2^{n}-1}e^{\frac {-2\pi ikx}{2^{n))}|x\rangle \right)={\frac {1}{2^{n))}\sum _{x=0}^{2^{n}-1}\sum _{k=0}^{2^{n}-1}e^{-{\frac {2\pi ik}{2^{n))}\left(x-2^{n}\theta \right)}|x\rangle .}$

We can approximate the value of ${\displaystyle \theta \in [0,1]}$ by rounding ${\displaystyle 2^{n}\theta }$ to the nearest integer. This means that ${\displaystyle 2^{n}\theta =a+2^{n}\delta ,}$ where ${\displaystyle a}$ is the nearest integer to ${\displaystyle 2^{n}\theta ,}$ and the difference ${\displaystyle 2^{n}\delta }$ satisfies ${\displaystyle 0\leqslant |2^{n}\delta |\leqslant {\tfrac {1}{2))}$.

Using this decomposition we can rewrite the state as ${\textstyle \sum _{x=0}^{2^{n}-1}c_{x}|x\rangle ,}$ where

${\displaystyle c_{x}\equiv {\frac {1}{2^{n))}\sum _{k=0}^{2^{n}-1}e^{-{\frac {2\pi ik}{2^{n))}(x-2^{n}\theta )}={\frac {1}{2^{n))}\sum _{k=0}^{2^{n}-1}e^{-{\frac {2\pi ik}{2^{n))}\left(x-a\right)}e^{2\pi i\delta k}.}$

Measurement

Performing a measurement in the computational basis on the first register yields the outcome ${\displaystyle |y\rangle }$ with probability

$${\displaystyle \Pr(y)=|c_{y}|^{2}=\left|{\frac {1}{2^{n))}\sum _{k=0}^{2^{n}-1}e^((\frac {-2\pi ik}{2^{n))}(y-a)}e^{2\pi i\delta k}\right|^{2}.}$$

${\displaystyle \operatorname {Pr} (a)=1}$

${\displaystyle \delta =0}$

${\displaystyle \theta }$

${\displaystyle \theta =a/2^{n))$

${\displaystyle y=a}$

${\displaystyle \delta \neq 0}$

$${\displaystyle \operatorname {Pr} (a)={\frac {1}{2^{2n))}\left|\sum _{k=0}^{2^{n}-1}e^{2\pi i\delta k}\right|^{2}={\frac {1}{2^{2n))}\left|{\frac {1-{e^{2\pi i2^{n}\delta ))}{1-{e^{2\pi i\delta ))))\right|^{2}.}$$

${\displaystyle \Pr(a)\geqslant {\frac {4}{\pi ^{2))}\approx 0.405}$

${\displaystyle \delta \neq 0}$

${\displaystyle \delta }$

${\displaystyle |\delta |\leqslant {\tfrac {1}{2^{n+1))))$

$${\displaystyle {\begin{aligned}\Pr(a)&={\frac {1}{2^{2n))}\left|{\frac {1-{e^{2\pi i2^{n}\delta ))}{1-{e^{2\pi i\delta ))))\right|^{2}&&{\text{for ))\delta \neq 0\\[6pt]&={\frac {1}{2^{2n))}\left|{\frac {2\sin \left(\pi 2^{n}\delta \right)}{2\sin(\pi \delta )))\right|^{2}&&\left|1-e^{2ix}\right|^{2}=4\left|\sin(x)\right|^{2}\\[6pt]&={\frac {1}{2^{2n))}{\frac {\left|\sin \left(\pi 2^{n}\delta \right)\right|^{2)){|\sin(\pi \delta )|^{2))}\\[6pt]&\geqslant {\frac {1}{2^{2n))}{\frac {\left|\sin \left(\pi 2^{n}\delta \right)\right|^{2)){|\pi \delta |^{2))}&&|\sin(\pi \delta )|\leqslant |\pi \delta |\\[6pt]&\geqslant {\frac {1}{2^{2n))}{\frac {|2\cdot 2^{n}\delta |^{2)){|\pi \delta |^{2))}&&|2\cdot 2^{n}\delta |\leqslant |\sin(\pi 2^{n}\delta )|{\text{ for ))|\delta |\leqslant {\frac {1}{2^{n+1))}\\[6pt]&\geqslant {\frac {4}{\pi ^{2))}.\end{aligned))}$$

We conclude that the algorithm provides the best ${\displaystyle n}$-bit estimate (i.e., one that is within ${\displaystyle 1/2^{n))$ of the correct answer) of ${\displaystyle \theta }$ with probability at least ${\displaystyle 4/\pi ^{2))$. By adding a number of extra qubits on the order of ${\displaystyle O(\log(1/\epsilon ))}$ and truncating the extra qubits the probability can increase to ${\displaystyle 1-\epsilon }$.^[5]

Toy examples

Consider the simplest possible instance of the algorithm, where only ${\displaystyle n=1}$ qubit, on top of the qubits required to encode ${\displaystyle |\psi \rangle }$, is involved. Suppose the eigenvalue of ${\displaystyle |\psi \rangle }$ reads ${\displaystyle \lambda =e^{2\pi i\theta ))$. The first part of the algorithm generates the one-qubit state ${\textstyle |\phi \rangle \equiv {\frac {1}{\sqrt {2))}(|0\rangle +\lambda |1\rangle )}$. Applying the inverse QFT amounts in this case to applying a Hadamard gate. The final outcome probabilities are thus ${\displaystyle p_{\pm }=|\langle \pm |\phi \rangle |^{2))$ where ${\textstyle |\pm \rangle \equiv {\frac {1}{\sqrt {2))}(|0\rangle \pm |1\rangle )}$, or more explicitly,

$${\displaystyle p_{\pm }={\frac {|1\pm \lambda |^{2)){4))={\frac {1\pm \cos(2\pi \theta )}{2)).}$$

${\displaystyle \lambda =1}$

${\displaystyle |\phi \rangle =|+\rangle }$

${\displaystyle p_{+}=1}$

${\displaystyle p_{-}=0}$

${\displaystyle \lambda }$

${\displaystyle \lambda =-1}$

If on the other hand ${\displaystyle \lambda =e^{2\pi i/3))$, then ${\displaystyle p_{\pm }=[1\pm \cos(2\pi /3)]/2}$, that is, ${\displaystyle p_{+}=1/4}$ and ${\displaystyle p_{-}=3/4}$. In this case the result is not deterministic, but we still find the outcome ${\displaystyle |-\rangle }$ as more likely, compatibly with the fact that ${\displaystyle 2/3}$ is close to 1 than to 0.

More generally, if ${\displaystyle \lambda =e^{2\pi i\theta ))$, then ${\displaystyle p_{+}\geq 1/2}$ if and only if ${\displaystyle |\theta |\leq 1/4}$. This is consistent with the results above because in the cases ${\displaystyle \lambda =\pm 1}$, corresponding to ${\displaystyle \theta =0,1/2}$, the phase is retrieved deterministically, and the other phases are retrieved with higher accuracy the closer they are to these two.