{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"tags": [
"remove_cell"
]
},
"source": [
"# Estimating pi ($\\pi$) using Quantum Phase Estimation Algorithm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Quick overview of the [Quantum Phase Estimation Algorithm](https://qiskit.org/textbook/ch-algorithms/quantum-phase-estimation.html)\n",
"\n",
"Quantum Phase Estimation (QPE) is a quantum algorithm that forms the building block of many more complex quantum algorithms. At its core, QPE solves a fairly straightforward problem: given an operator $U$ and a quantum state $\\vert\\psi\\rangle$ that is an eigenvalue of $U$ with $U\\vert\\psi\\rangle = \\exp\\left(2 \\pi i \\theta\\right)\\vert\\psi\\rangle$, can we obtain an estimate of $\\theta$?\n",
"\n",
"The answer is yes. The QPE algorithm gives us $2^n\\theta$, where $n$ is the number of qubits we use to estimate the phase $\\theta$.\n",
"\n",
"## 2. Estimating $\\pi$\n",
"\n",
"In this demo, we choose\n",
"$$U = p(\\theta), \\vert\\psi\\rangle = \\vert1\\rangle$$\n",
"where \n",
"$$\n",
"p(\\theta) = \\begin{bmatrix}\n",
"1 & 0\\\\ 0 & \\exp(i\\theta)\n",
"\\end{bmatrix}\n",
"$$\n",
"is one of the quantum gates available in Qiskit, and\n",
"$$p(\\theta)\\vert1\\rangle = \\exp(i\\theta)\\vert1\\rangle.$$ \n",
"\n",
" By choosing the phase for our gate to be $\\theta = 1$, we can solve for $\\pi$ using the following two relations:\n",
"\n",
"1. From the output of the QPE algorithm, we measure an estimate for $2^n\\theta$. Then, $\\theta = \\text{measured} / 2^n$ \n",
"2. From the definition of the $p(\\theta)$ gate above, we know that $2\\pi\\theta = 1 \\Rightarrow \\pi = 1 / 2\\theta$\n",
"\n",
"Combining these two relations, $\\pi = 1 / \\left(2 \\times (\\text{(measured)}/2^n)\\right)$. \n",
"\n",
"For detailed understanding of the QPE algorithm, please refer to the chapter dedicated to it in the Qiskit Textbook located at [qiskit.org/textbook](https://qiskit.org/textbook/ch-algorithms/quantum-phase-estimation.html)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Time to write code"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We begin by importing the necessary libraries."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"## import the necessary tools for our work\n",
"from IPython.display import clear_output\n",
"from qiskit import *\n",
"from qiskit.visualization import plot_histogram\n",
"import numpy as np\n",
"import matplotlib.pyplot as plotter\n",
"from qiskit.tools.monitor import job_monitor\n",
"# Visualisation settings\n",
"import seaborn as sns, operator\n",
"sns.set_style(\"dark\")\n",
"\n",
"pi = np.pi"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The function `qft_dagger` computes the inverse Quantum Fourier Transform. For a detailed understanding of this algorithm, see the dedicated chapter for it in the [Qiskit Textbook](https://qiskit.org/textbook/ch-algorithms/quantum-fourier-transform.html)."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"## Code for inverse Quantum Fourier Transform\n",
"## adapted from Qiskit Textbook at\n",
"## qiskit.org/textbook\n",
"\n",
"def qft_dagger(circ_, n_qubits):\n",
" \"\"\"n-qubit QFTdagger the first n qubits in circ\"\"\"\n",
" for qubit in range(int(n_qubits/2)):\n",
" circ_.swap(qubit, n_qubits-qubit-1)\n",
" for j in range(0,n_qubits):\n",
" for m in range(j):\n",
" circ_.cp(-np.pi/float(2**(j-m)), m, j)\n",
" circ_.h(j)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The next function, `qpe_pre`, prepares the initial state for the estimation. Note that the starting state is created by applying a Hadamard gate on the all but the last qubit, and setting the last qubit to $\\vert1\\rangle$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"## Code for initial state of Quantum Phase Estimation\n",
"## adapted from Qiskit Textbook at qiskit.org/textbook\n",
"## Note that the starting state is created by applying \n",
"## H on the first n_qubits, and setting the last qubit to |psi> = |1>\n",
"\n",
"def qpe_pre(circ_, n_qubits):\n",
" circ_.h(range(n_qubits))\n",
" circ_.x(n_qubits)\n",
"\n",
" for x in reversed(range(n_qubits)):\n",
" for _ in range(2**(n_qubits-1-x)):\n",
" circ_.cp(1, n_qubits-1-x, n_qubits)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we write a quick function, `run_job`, to run a quantum circuit and return the results."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"## Run a Qiskit job on either hardware or simulators\n",
"\n",
"def run_job(circ, backend, shots=1000, optimization_level=0):\n",
" t_circ = transpile(circ, backend, optimization_level=optimization_level)\n",
" qobj = assemble(t_circ, shots=shots)\n",
" job = backend.run(qobj)\n",
" job_monitor(job)\n",
" return job.result().get_counts()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then, load your account to use the cloud simulator or real devices."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"tags": [
"uses-hardware"
]
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"ibmqfactory.load_account:WARNING:2021-02-18 11:32:59,136: Credentials are already in use. The existing account in the session will be replaced.\n"
]
}
],
"source": [
"## Load your IBMQ account if \n",
"## you'd like to use the cloud simulator or real quantum devices\n",
"my_provider = IBMQ.load_account()\n",
"simulator_cloud = my_provider.get_backend('ibmq_qasm_simulator')\n",
"device = my_provider.get_backend('ibmq_16_melbourne')"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"simulator = Aer.get_backend('qasm_simulator')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we bring everything together in a function called `get_pi_estimate` that uses `n_qubits` to get an estimate for $\\pi$."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"## Function to estimate pi\n",
"## Summary: using the notation in the Qiskit textbook (qiskit.org/textbook),\n",
"## do quantum phase estimation with the 'phase' operator U = p(theta) and |psi> = |1>\n",
"## such that p(theta)|1> = exp(2 x pi x i x theta)|1>\n",
"## By setting theta = 1 radian, we can solve for pi\n",
"## using 2^n x 1 radian = most frequently measured count = 2 x pi\n",
"\n",
"def get_pi_estimate(n_qubits):\n",
"\n",
" # create the circuit\n",
" circ = QuantumCircuit(n_qubits + 1, n_qubits)\n",
" # create the input state\n",
" qpe_pre(circ, n_qubits)\n",
" # apply a barrier\n",
" circ.barrier()\n",
" # apply the inverse fourier transform\n",
" qft_dagger(circ, n_qubits)\n",
" # apply a barrier\n",
" circ.barrier()\n",
" # measure all but the last qubits\n",
" circ.measure(range(n_qubits), range(n_qubits))\n",
"\n",
" # run the job and get the results\n",
" counts = run_job(circ, backend=simulator, shots=10000, optimization_level=0)\n",
" # print(counts) \n",
"\n",
" # get the count that occurred most frequently\n",
" max_counts_result = max(counts, key=counts.get)\n",
" max_counts_result = int(max_counts_result, 2)\n",
" \n",
" # solve for pi from the measured counts\n",
" theta = max_counts_result/2**n_qubits\n",
" return (1./(2*theta))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, run the `get_pi_estimate` function with different numbers of qubits and print the estimates."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Job Status: job has successfully run\n",
"2 qubits, pi ≈ 2.0\n",
"Job Status: job has successfully run\n",
"3 qubits, pi ≈ 4.0\n",
"Job Status: job has successfully run\n",
"4 qubits, pi ≈ 2.6666666666666665\n",
"Job Status: job has successfully run\n",
"5 qubits, pi ≈ 3.2\n",
"Job Status: job has successfully run\n",
"6 qubits, pi ≈ 3.2\n",
"Job Status: job has successfully run\n",
"7 qubits, pi ≈ 3.2\n",
"Job Status: job has successfully run\n",
"8 qubits, pi ≈ 3.1219512195121952\n",
"Job Status: job has successfully run\n",
"9 qubits, pi ≈ 3.1604938271604937\n",
"Job Status: job has successfully run\n",
"10 qubits, pi ≈ 3.1411042944785277\n",
"Job Status: job has successfully run\n",
"11 qubits, pi ≈ 3.1411042944785277\n",
"Job Status: job has successfully run\n",
"12 qubits, pi ≈ 3.1411042944785277\n"
]
}
],
"source": [
"# estimate pi using different numbers of qubits\n",
"nqs = list(range(2,12+1))\n",
"pi_estimates = []\n",
"for nq in nqs:\n",
" thisnq_pi_estimate = get_pi_estimate(nq)\n",
" pi_estimates.append(thisnq_pi_estimate)\n",
" print(f\"{nq} qubits, pi ≈ {thisnq_pi_estimate}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And plot all the results."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
"