Example 4: Quantum Phase Estimation (QPE)

Purpose

Quantum Phase Estimation (QPE) is an important part of many quantum algorithms, such as Shor's factorization algorithm. This example explains how QPE works and how it can be implemented uisng TinyQsim.

This notebook builds on the previous ones on Phase Kickback and the Quantum Fourier Transform (QFT), so it is recommended to study those first.

Imports and Definitions

We will start by importing the libraries that we need and the inverse QFT from the earlier example notebook.

from tinyqsim.qcircuit import QCircuit
from tinyqsim import plotting, quantum
from tinyqsim.utils import int_to_bits
from examples_lib import iqft
import numpy as np
from math import pi
from cmath import exp
from functools import partial

U_PI = '\u03C0'  # Unicode pi
U_PHI = '\u03d5'  # Unicode phi
U_THETA = '\u03b8'  # Unicode theta

Introduction

The aim of QPE is to find the phase associated with an eigenvector of a unitary \(U\).

Phase Kickback

We will quickly recap the main points from the Phase Kickback notebook.

Any unitary operator can be expressed in terms of the phases \(\phi_j\) of its eigenvectors:

\(\qquad U=\sum_j{e^{i\phi_j}\,\ket{u_j}\bra{u_j}}\)

where \(\lambda_j = e^{i\phi_j}\) are the corresponding eigenvalues.

The eigenvalue always has a magnitude of one because the operator is unitary.

If we apply the unitary to one of its eigenvectors \(\ket{\psi_k}\), the result will be the corresponding eigenvalue \(e^{i\phi_k}\) multiplied by the eigenstate:

\[U \ket{u_k} = e^{i\phi_k}\ket{u_k} \]

If we apply a controlled unitary to one of its eigenvectors, with the control in an equal superposition state, the phase is "kicked back" to the control qubit and the target qubit is left unchanged.

Fractional Phase

The principal value of the phase \(\phi\) can be expressed as a fraction \(\theta\) of \(2\pi\).

\[\theta=\frac{\phi}{2\pi},\quad\text{where }\theta\in[0,1) \]

This can be approximated as an \(m\)-bit binary fraction:

\[\quad\theta\approx\sum_{k=1}^m{\theta_k} 2^{-k},\quad\text{where }\theta_k\in\{0,1\} \]

which can be written using '.' for the binary point, analogous to a decimal point:

\[\quad\theta\approx 0.\theta_1\theta_2\dots\theta_m \]

Encoding the phase this way allows us to represent it using a qubit register, with the qubits corresponding to the bits of the binary fraction. The aim is to evolve the phase so that, when measured, it will give the bits of the binary fraction. For example, with a 3-qubit register, a measurement outcome of \(\ket{011}\) would correspond to the phase \(\theta=.011\) in binary, if we place a binary point before it.

Example Unitary Operator

To demonstate QPE, we will define a 2\(\times\)2 matrix which has eigenvalue phases of 0, p1, p2, p3 for eigenstates \(\ket{00}, \ket{01}, \ket{10}, \ket{11}\), respectively. The phases are expressed as multiples of \(2\pi\).

\[\begin{bmatrix} 1&0&0&0\\0&e^{i \phi_1}&0&0\\0&0&e^{i \phi_2}&0\\0&0&0&e^{i \phi_3} \end{bmatrix} \]

where \(\phi_k= 2\pi p_k,\quad p_k\in [0,1)\)

Note that the suffix \(k\) specifies which eigenvalue we are considering. This should not be confused with the suffixes on \(\theta\) which are bit indices.

The following function returns \(U^k\) where the exponent \(k\) is given as a parameter.

def create_unitary(p1, p2, p3, k):
    # Create an 2-qubit unitary matrix with eigenvalue phases: 
    # 0,p1,p2,p2 x 2PI for eigenstates |00>, |01>, |10>, |11>
    x = 2 * pi * k
    a = exp(1j * p1 * x)
    b = exp(1j * p2 * x)
    c = exp(1j * p3 * x)
    return np.array([[1,0,0,0],
                     [0,a,0,0],
                     [0,0,b,0],
                     [0,0,0,c],
                    ])

The QPE Circuit

The following function will be used to construct the QPE circuit. The parameter 'eigenvec' specifies the eigenvector of interest, in decimal. For example, eigenvec=2 sets the eigenvector to \(\ket{2}=\ket{10}\).

def create_qpe_circuit(unitary, nqc, nqt, eigenvec=1, pre_iqft=None):
    """Create Quantum Phase Estimation (QPE) circuit, excluding measurement.
    :param unitary: function returning the unitary matrix
    :param nqc: Number of control qubits
    :param nqt: Number of target qubits
    :param eigenvec: e.g. 3 for |0011> when nqt=4 (default=1)
    :param pre_iqft: Function called just before IQFT
    :return: QPE circuit with 'nqc' control qubits and 'nqt' target qubits."""

    # Registers
    regC = range(nqc)  # Control register
    regT = range(nqc, nqc + nqt)  # Target register

    # Qubit labels
    labels = [f'C{i}' for i in range(nqc)] + [f'T{i}' for i in range(nqt)]
    label_dict = dict(zip(range(nqc + nqt), labels))

    # Create the circuit
    qc = QCircuit(nqc + nqt)
    qc.qubit_labels(label_dict, numbers=False)

    # Initialialize registers
    x = int_to_bits(eigenvec,nqt)
    qc.x([regT[i] for i,s in enumerate(x) if s==1])
    qc.h(regC)
    qc.barrier('1')
 
    # Apply unitary operators
    for i, j in enumerate(reversed(regC)):
        u = unitary(2**i)
        qc.cu(u, f'$U^{{{2**i}}}$', j, *regT)
    qc.barrier('2')

    # Call function to examine state just before the IQFT
    if pre_iqft is not None:
        pre_iqft(qc, nqc)
    
    # Inverse QFT
    iqft(qc, nqc)
    return qc

The following diagram shows an example of the QPE circuit, for the case of 3 control qubits and 2 target qubits.

nqc = 3  # Control qubits
nqt = 2  # Target qubits
n = 2**nqc
u1 = partial(create_unitary, 1/n, 3/n, 5/n)
qc = create_qpe_circuit(u1, nqc, nqt, eigenvec=2)
qc.barrier('3')
m = qc.measure(*range(nqc))
print(f'\nMeasurement result = {m}\n')
qc.draw()

Measurement result = [0 1 1]

There are two quantum registers: a Control register C with \(m\) qubits and a Target register T with \(t\) qubits.
In this example \(m=3, t=2\).

The circuit has four stages, which are separated by barrier symbols in the diagram:

• Initialisation
• Phase kickback
• Inverse QFT
• Measurement

We will consider each of these stages in turn:

Initialization

The T register is initialized to the eigenvector of interest \(\ket{u}\), which has eigenvalue \(\lambda=e^{2\pi i\theta}\), where \(\theta\in [0,1)\) is the phase we want to find.

The eigenvector \(\ket{u}\) and eigenvalue \(\lambda\) are related by the eigenvalue equation:

\[\qquad\begin{align*} U \ket{u} = \lambda\ket{u}=e^{2\pi i\theta}\ket{u} \end{align*} \]

A Hadamard transform is applied to the control register, which places it in a superposition state:

\[\begin{align*} \ket{\psi_1}=(H^{\otimes m}\ket{0}^{\otimes m})\otimes\ket{u}\\ ={\small\frac{1}{\sqrt{M}}}\sum_{k=0}^{M-1}\ket{k} \otimes\ket{u} \end{align*} \]

where \(M=2^m\)

Phase Kickback

You will probably recognize that the Hadamard gates and the controlled unitaries form three instances of the Phase Kickback circuit discussed in a previous example. The final Hadamard gates of the Phase Kickback circuits have been replaced by a single inverse Quantum Fourier Transform (QFT\(^\dagger\)), but we will come back to that shortly.

The action of applying a unitary operator \(U\) to one of its eigenvectors \(\ket{u}\) is to multiply the target eigenvector by its eigenvalue. If the unitary is a controlled-unitary, the phase is kicked-back from the target to the control and the target is left unchanged:

\[\begin{align*} \ket{\psi_2}&={\small\frac{1}{\sqrt{M}}}\sum_{k=0}^{M-1}\lambda\ket{k} \otimes\ket{u}\\ &={\small\frac{1}{\sqrt{M}}}\sum_{k=0}^{M-1}e^{2\pi i k\ \theta}\ket{k} \otimes\ket{u} \end{align*} \]

Inverse QFT

The Quantum Fourier Transform (QFT) is defined by:

\[QFT: \ket{x}\mapsto{\small\frac{1}{\sqrt{M}}}\sum_{k=0}^{M-1}e^{2\pi i k\frac{x}{M}}\ket{k} \]

This is identical to the control register part of \(\ket{\psi_2}\) if we equate \(\theta\) with \(\frac{x}{M}\) :

\[\theta = \frac{x}{M} \]

So, taking the inverse QFT of the control register should result in a single basis state \(\ket{x}\), provided that \(\theta\) can be expressed as fraction with an integer numerator and denominator \(M\). Otherwise the probability distribution will peak close to this value and \(\frac{x}{M}\) is the nearest approximation to \(\theta\) expressed as a binary fraction with \(m\) bits. In this case, \(\frac{x}{M}\) is an estimate \(\hat{\theta}\):

\[\hat{\theta}=\frac{x}{M} \]

However, the phase we want to find might be a value that is not expressible as a fraction with an integer numerator and a denominator which is a power of 2, for example \(\frac{5}{26}\). In the next notebook, we will look at how Shor's factorization algorithm solves this problem using partial fraction convergents.

Measurement Result

The example circuit above was configured with the target register initialized to the eigenvector \(\ket{10}\), which corresponds to the eigenvalue with phase \(3/2^3=3/8\). The measurement result was \(011_2\) in binary. If we precede this by a binary point, it is the phase as the binary fraction \(\theta=.011_2=\frac{3}{8}\), which is correct. These phases are expressed as fractions of \(2\pi\).

Since we are using a state simulator, instead of a real quantum computer, we will omit the measurement stage from the subsequent examples. Simulation allows us to look at the phases and probabilities, which will give us a much better insight into how the algorithm works.

Plotting the Phase

It is not meaningful to plot the phase of just the control register, unless its state is separable from that of the target register. However, with phase kickback, the states of the control and target registers are separable, so we can plot just the control register.

The earlier Phase Kickback notebook showed that the result of phase kickback with one control qubit could be factored into a tensor product, showing that the control and target qubits are not entangled. In our QPE example, we simply have three instances of a similar circuit.

def plot_phase(qc, nq):
    """Plot fractional phase of first 'nq' qubits assuming 
       they are not entangled with the rest of the state."""
    state = qc.state_vector.reshape(2**nq,-1).sum(1)
    nq = quantum.n_qubits(state)
    phase = np.atan2(state.imag, state.real) / np.pi
    phase = np.where(phase<0, phase+2, phase) / 2
    plotting.plot_bars(quantum.basis_names(nq), [phase], ylabels=[f'Phase'],ylims=[(0,1)])

Examples

The next example considers a QPE circuit with 5 control qubits and 2 target qubits. The unitary matrix will be created with phases corresponding to \(\frac{1}{32},\frac{4}{32}\text{ and }\frac{5.5}{32}\), for eigenvectors \(\ket{1},\ket{2}\text{ and }\ket{3}\) respectively.

nqc = 5
nqt = 2
n = 2**nqc
u = partial(create_unitary, 1/32, 5/32, 11/64)

Example 1

For the first example, we will set the eigenvector to \(\ket{01}\), which corresponds to the eigenvalue with fractional phase \(\theta=\frac{1}{32}\).

qc = create_qpe_circuit(u, nqc, nqt, 1, plot_phase)
qc.plot_probabilities(*range(nqc), height=0.7)

The above plots show the phase of the Control register just before the inverse QFT (IQFT) and the measurement probabilities after the IQFT. The plots were generated by the state simulator. This should look familar if you have studied the preceeding QFT example notebook.

Looking at this in reverse, an output which correponds to a single basis state \(\ket{k}\), would result in a linear phase ramp with \(k\) cycles of \(2\pi\). In this case \(k=1\).

The probability distribution shows that measuring the control register would give the result 00001. If we precede this by a binary point, it is the phase as the binary fraction \(\theta=.00001_2 = \frac{1}{32}\), which is correct.

Example 2

For this example, we will set the eigenvector to \(\ket{10}\), which corresponds to the eigenvalue with fractional phase \(\theta=\frac{5}{32}\).

nqc = 5
nqt = 2
qc = create_qpe_circuit(u, nqc, nqt, 2, plot_phase)
qc.plot_probabilities(*range(nqc), height=0.7)

The result of measurement would be \(00101\). If we precede this with a binary point we get \(\theta=.00101_2 = \frac{5}{32}\), as expected.

The phase has 5 cycles in a state vector of length 32 and the result is \(\frac{5}{32}\).

Example 3

The previous examples have deliberately used fractional phases that are integer multiples of \(\large\frac{1}{2^5}\) as the control register has 5 qubits. The following example sets the eigenvector to \(\ket{11}\) which has an eigenvalue phase of \(\theta=\frac{11}{64} = \frac{5.5}{32}\).

qc = create_qpe_circuit(u, nqc, nqt, 3, plot_phase)
qc.plot_probabilities(*range(nqc), height=0.7)

There are now 5.5 phase cycles, which doesn't correspond to a single basis state. The probabilities are divided between results 5 and 6, but there is some leakage into adjacent values.

Note: If you are familiar with Fourier Transforms in signal processing, you may recognize this as spectral leakage.

In this case, we could just repeat the QPE with 6 control qbits and we would get the result \(\theta=001011_2=\frac{11}{64}\).

However, the phase we want to find might not be expressible as a fraction with an integer numerator and a denominator that is a power of 2, such as \(\frac{5}{26}\). In the next notebook we look at how Shor's factorization algorithm solves this problem using partial fraction convergents.