Estimating Pi Using Qubits

Image Credit: Patrick Pallagi

Quantum Phase Estimation

The Quantum Phase Estimation (QPE) algorithm is a quantum algorithm that can be used to estimate the phase of an eigenvector of a unitary operator. It was first introduced by Kitaev in 1995.

In more detail, the QPE algorithm takes as input a unitary operator U and an eigenvector |u⟩ of U with corresponding eigenvalue e^(2πiθ), where θ is the unknown phase. The goal of the algorithm is to estimate θ.

The QPE algorithm uses a two-register quantum circuit. The first register is an n-qubit register initialized to the state |0⟩, and the second register is initialized to the state |u⟩. The algorithm works by applying a sequence of controlled-U operations to the second register, where the number of controls increases from 1 to n. After this sequence of operations, the first register contains an estimate of θ as a binary fraction.

The QPE algorithm has applications in quantum simulation, quantum chemistry, and cryptography. It is also an important subroutine in many other quantum algorithms, such as Shor’s algorithm for factoring large integers.

The Discovery of Quantum Phase Estimation

In 1995, Alexei Kitaev introduced the Quantum Phase Estimation (QPE) algorithm in his paper “Quantum measurements and the Abelian stabilizer problem”.

In his paper, Kitaev showed that the QPE algorithm can be used to solve the Abelian stabilizer problem, which is a problem in quantum error correction. The problem involves determining the stabilizer group of a quantum error-correcting code, which is a group of Pauli operators that leave the code subspace invariant.

Kitaev observed that the stabilizer group can be determined by finding the eigenvalues of certain stabilizer operators, which are unitary operators that commute with the stabilizer group. The eigenvalues of these operators are phases that can be estimated using the QPE algorithm.

Kitaev’s paper also introduced the notion of a quantum circuit with “quantum Fourier transform” (QFT) gates, which are used in the QPE algorithm to implement the quantum Fourier transform on the first register of qubits. This idea was later generalized to the more powerful notion of the quantum Fourier transform (QFT) algorithm, which can be used to efficiently implement the Fourier transform on a quantum computer.

Overall, Kitaev’s work on the QPE algorithm was an important early contribution to the field of quantum algorithms and has had a significant impact on the development of quantum computing.

Estimations Using Quantum Phase Estimation

As I mentioned before, the algorithm works by first preparing a two-register quantum state, where the first register contains n qubits initialized to the state |0⟩, and the second register is initialized to the state |u⟩. The two registers are then entangled using a controlled-U operation, where the control qubits are the qubits in the first register, and the target qubits are the qubits in the second register.

Specifically, the controlled-U operation is defined as follows:

C-U|j⟩|u⟩ = |j⟩U^j|u⟩

where |j⟩ is a basis state of the first register, and U^j is the jth power of the unitary operator U.

After applying a sequence of controlled-U operations with increasing numbers of control qubits, the state of the first register is then measured in the computational basis. The resulting measurement outcome is a binary fraction that can be used to estimate θ.

To see how this works, let’s assume that the measured outcome is x_1x_2…x_n, where x_i is either 0 or 1. Then, the estimated value of θ is given by:

θ ≈ x_1/2 + x_2/2² + … + x_n/2^n

This binary fraction approximation of θ has a precision of 1/2^n, which can be improved by increasing the number of qubits in the first register.

Overall, the QPE algorithm has a runtime that scales polynomially with the number of qubits used, making it efficient for certain types of problems. It is also a key subroutine in many other quantum algorithms, such as Shor’s algorithm for factoring large integers.

This Time We’re Estimating Pi

Estimation of Pi Using IBM Qiskit

Estimating Pi takes around 10 qubits. It is yet another proof that quantum computers can grapple with the irrational.

Thank you for reading!

Code for Estimating Pi Using Qubits

import numpy as np
import matplotlib.pyplot as plotter
from qiskit import Aer, QuantumCircuit, transpile, assemble
from qiskit.visualization import plot_histogram
from qiskit.tools.monitor import job_monitor
from qiskit_ibm_provider import IBMProvider
from numpy import pi  # Import pi from numpy


# Inverse Quantum Fourier Transform
def qft_dagger(circ, n_qubits):
    for qubit in range(int(n_qubits/2)):
        circ.swap(qubit, n_qubits-qubit-1)
    for j in range(0, n_qubits):
        for m in range(j):
            circ.cp(-np.pi/float(2**(j-m)), m, j)
        circ.h(j)

# Initial state for Quantum Phase Estimation
def qpe_pre(circ, n_qubits):
    circ.h(range(n_qubits))
    circ.x(n_qubits)

    for x in reversed(range(n_qubits)):
        for _ in range(2**(n_qubits-1-x)):
            circ.cp(1, n_qubits-1-x, n_qubits)

# Run a Qiskit job on either hardware or simulators
def run_job(circ, backend, shots=1000, optimization_level=0):
    t_circ = transpile(circ, backend, optimization_level=optimization_level)
    job = backend.run(t_circ, shots=shots)
    job_monitor(job)
    return job.result().get_counts()


# Load your IBMQ account
my_provider = IBMProvider()

simulator_cloud = my_provider.get_backend('ibmq_qasm_simulator')

# Use an available quantum device, for example, 'ibmq_belem'
device = my_provider.get_backend('ibmq_belem')

# Get Aer simulator
simulator = Aer.get_backend('aer_simulator')

# Function to estimate pi
def get_pi_estimate(n_qubits):
    circ = QuantumCircuit(n_qubits + 1, n_qubits)
    qpe_pre(circ, n_qubits)
    circ.barrier()
    qft_dagger(circ, n_qubits)
    circ.barrier()
    circ.measure(range(n_qubits), range(n_qubits))

    counts = run_job(circ, backend=simulator, shots=10000, optimization_level=3)
    max_counts_result = max(counts, key=counts.get)
    max_counts_result = int(max_counts_result, 2)
    theta = max_counts_result/2**n_qubits
    return (1./(2*theta))

# Estimate pi using different numbers of qubits
nqs = list(range(2,12+1))
pi_estimates = []
for nq in nqs:
    thisnq_pi_estimate = get_pi_estimate(nq)
    pi_estimates.append(thisnq_pi_estimate)
    print(f"{nq} qubits, pi ≈ {thisnq_pi_estimate}")

# Plot the results

plotter.plot(nqs, [pi]*len(nqs), '--r')
plotter.plot(nqs, pi_estimates, '.-', markersize=12)
plotter.xlim([1.5, 12.5])
plotter.ylim([1.5, 4.5])
plotter.legend(['$\pi$', 'estimate of $\pi$'])
plotter.xlabel('Number of qubits', fontdict={'size':20})
plotter.ylabel('$\pi$ and estimate of $\pi$', fontdict={'size':20})
plotter.tick_params(axis='x', labelsize=12)
plotter.tick_params(axis='y', labelsize=12)
plotter.show()