quantum_phase_estimation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/quantum_algorithms_tutorials.png\" alt=\"drawing\" width=\"100\" align=\"left\"/>\n",
    "\n",
    "<h1 align=\"center\">\n",
    "\tQuantum phase estimation\n",
    "</h1>\n",
    "\n",
    "The quantum phase estimation (QPE, also known as PEA) algorithm is a quantum algorithm for estimating the phase of an eigenvector of a unitary operator [[1](https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm)]. One of the main applications of the algorithm is to estimate eigenvalues (energies) of some molecule's $H$ Hamiltonian [[2](https://www.tandfonline.com/doi/abs/10.1080/00268976.2011.552441)]. Also, phase estimation is frequently used as a subroutine in other quantum algorithms, such as Shor's algorithm and the quantum algorithm for linear systems of equations (HHL algorithm) [[1](https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm)]. \n",
    "\n",
    "The goal of the algorithm is to estimate the eigenvalue of a given Hamiltonian $\\hat{H}$. We exponentiate it and obtain $U = e^{i\\hat{H}t}$ unitary operator, where $t$ is some parameter (in our example $t = 1$). Then we apply quantum circuit that corresponds to $U = e^{i\\hat{H}t}$ on one of the eigenvectors of the $\\hat{H}$. For simplicity, we will take such $\\hat{H}$, which eigenvectors will be the computational basis states $\\left| 0 \\right\\rangle$ and $\\left| 1 \\right\\rangle$. For that purpose $\\hat{H}$ should be a diagonal matrix). In our example, we will try to find eigenvalue that corresponds to $\\left| 1 \\right\\rangle$ eigenvector. More mathematically:\n",
    "\n",
    "\\begin{align}\n",
    "\\hat{H} = \\begin{pmatrix}\n",
    "E_1 & 0\\\\\n",
    "0 & E_2\n",
    "\\end{pmatrix} \\longrightarrow\n",
    "e^{i\\hat{H}t} = \\begin{pmatrix}\n",
    "e^{i E_1 t} & 0\\\\\n",
    "0 & e^{i E_2 t} \n",
    "\\end{pmatrix} \\longrightarrow e^{i H t} \\left| 1 \\right\\rangle = e^{i E_2 t} \\left| 1 \\right\\rangle = e^{i 2 \\pi \\varphi_2} \\left| 1 \\right\\rangle,\n",
    "\\end{align}\n",
    "\n",
    "where $E_1$ and $E_2$ are the eigenvalues of the $\\hat{H}$, $\\varphi_2 = \\frac{E_2 t}{2 \\pi}$ is the phase that the algorithm is capable to estimate. Note that $U = e^{i H t}$ is a unitary operator and in the code, we will name it as ```unitary```. For simplicity, we are considering here diagonal Hermitian operator $\\hat{H}$. In general case, it is not diagonal.\n",
    "\n",
    "The code and the algorithm is similar to the our covered [IQPE](https://github.com/DavitKhach/quantum-algorithms-tutorials/blob/master/iterative_quantum_phase_estimation.ipynb) tutorial. IQPE actually is a modification of QPE algorithm. So, that's why some parts/explanations are omitted here and more details are provided about the parts of QPE that are different from IQPE.\n",
    "\n",
    "In the first code, all necessary packages are imported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from random import random\n",
    "from qiskit import *\n",
    "from qiskit.aqua.utils.controlled_circuit import get_controlled_circuit\n",
    "from qiskit.extensions.standard.barrier import Barrier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following code the quantum ```circuit``` is created with \n",
    "* $q$ quantum state register, where we encode the $\\left| 1 \\right\\rangle$ eigenvector by applying ```x``` gate.\n",
    "\n",
    "* $a$ ancillary quantum register, where, in the course of the algorithm, the bits of the phase $\\varphi_2$ will be stored.\n",
    "\n",
    "* $c$ classical register, where the measurement results will be stored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<qiskit.circuit.instructionset.InstructionSet at 0x7fb56c871780>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_bits_estimate = 10\n",
    "# For 2x2 matrix one qubit is enough\n",
    "q = QuantumRegister(1, name=\"q\")\n",
    "# In QPE we use n ancillas to estimate n bits from the phase\n",
    "a = QuantumRegister(num_bits_estimate, name=\"a\") \n",
    "# For n ancillary qubit measurment we need n cllasical bits\n",
    "c = ClassicalRegister(num_bits_estimate, name=\"c\") \n",
    "\n",
    "# Create a quantum circuit\n",
    "circuit = QuantumCircuit(q, a, c)\n",
    "\n",
    "# |1> eigenstate initialization\n",
    "circuit.x(q[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The creation of the circuit for $e^{iHt}$ unitary evolution for a random diagonal Hermitian operator $\\hat{H}$ is presented in the code below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We are going to estimate E_2 via QPE algorithm \n",
      "E_2 = 0.7294990675852743\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<qiskit.circuit.instructionset.InstructionSet at 0x7fb56c871dd8>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E_1, E_2 = (2 * np.pi * random(), 2 * np.pi * random())\n",
    "print(\"We are going to estimate E_2 via QPE algorithm \\nE_2 = {}\".format(E_2))\n",
    "\n",
    "# circuit for unitary operator exp(iHt)\n",
    "t = 1\n",
    "unitary = QuantumCircuit(q)\n",
    "\n",
    "unitary.u1(E_2 * t, q[0]) # q[0] is the only qubit in q register\n",
    "unitary.x(q[0])\n",
    "unitary.u1(E_1 * t, q[0])\n",
    "unitary.x(q[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we prepare the ancillary qubits by applying Hadamard gates to all of them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ancillary in a:\n",
    "    circuit.h(ancillary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantum state of ancillary qubits before Hadamard gates (in the text we will examine 5 ancillary qubit case, but in the code we will use 10 ancillary qubits for better results):\n",
    "\n",
    "\\begin{align*}\n",
    "\\left| \\psi_a \\right\\rangle = \\left| 00000 \\right\\rangle = \\left| 0 \\right\\rangle \\left| 0 \\right\\rangle \\left| 0 \\right\\rangle \\left| 0 \\right\\rangle \\left| 0 \\right\\rangle\n",
    "\\end{align*}\n",
    "\n",
    "The quantum state of ancillary qubits after Hadamard gates:\n",
    "\n",
    "\\begin{align*}\n",
    "\\left| \\psi_a \\right\\rangle &= \\frac{1}{2^{5/2}}\n",
    "\\left(\\left| 0 \\right\\rangle + \\left| 1 \\right\\rangle \\right)\n",
    "\\left(\\left| 0 \\right\\rangle + \\left| 1 \\right\\rangle \\right)\n",
    "\\left(\\left| 0 \\right\\rangle + \\left| 1 \\right\\rangle \\right)\n",
    "\\left(\\left| 0 \\right\\rangle + \\left| 1 \\right\\rangle \\right)\n",
    "\\left(\\left| 0 \\right\\rangle + \\left| 1 \\right\\rangle \\right) = \n",
    "\\\\\n",
    "&= \\left| + \\right\\rangle \\left| + \\right\\rangle \\left| + \\right\\rangle \\left| + \\right\\rangle \\left| + \\right\\rangle\n",
    "=\\left| +++++ \\right\\rangle ,\n",
    "\\end{align*}\n",
    "\n",
    "because $H \\left| 0 \\right\\rangle = \\frac{1}{\\sqrt{2}}\n",
    "\\left(\\left| 0 \\right\\rangle + \\left| 1 \\right\\rangle \\right) = \\left| + \\right\\rangle$. Please don't confuse Hadamard gate $H$ with Hermitian operator $\\hat{H}$, because we use the same letter for both of them.\n",
    "\n",
    "For general case we will represent $\\varphi$ in its fractional binary representation:\n",
    "\n",
    "\\begin{align*}\n",
    "\\varphi = 0.\\varphi_0 \\varphi_1 \\varphi_2 \\varphi_3... =  \\frac{\\varphi_0}{2} + \\frac{\\varphi_1}{2^2} + \\frac{\\varphi_2}{2^3} + \\frac{\\varphi_3}{2^4} + ...,\n",
    "\\end{align*}\n",
    "\n",
    "where each $\\varphi_k$ is a bit (it is equal to $0$ or $1$).\n",
    "\n",
    "Now, for each $n$th ancillary qubit, we apply the controlled version of the circuit of $U = e^{i \\hat{H} t}$ $2^{n}$ times ($U^{2^{n}}$) to $q$ quantum register controlled by $n$th ancillary qubit ($n$ starts from $0$ to $4$). Here is the circuit so far:\n",
    "\n",
    "<img src=\"images/qpe_control_5ancillary.png\" alt=\"drawing\" width=\"300\"/>\n",
    "\n",
    "This will give us a quantum state for ancillary qubits:\n",
    "\n",
    "\\begin{align*}\n",
    "\\left| \\psi_a \\right\\rangle = \\frac{1}{2^{5/2}}\n",
    "&\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi \\varphi} \\left| 1 \\right\\rangle \\right)\n",
    "\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 2\\varphi} \\left| 1 \\right\\rangle \\right)\n",
    "\\\\\n",
    "&\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 4\\varphi} \\left| 1 \\right\\rangle \\right)\n",
    "\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 8\\varphi} \\left| 1 \\right\\rangle \\right)\n",
    "\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 16\\varphi} \\left| 1 \\right\\rangle \\right).\n",
    "\\end{align*}\n",
    "\n",
    "By taking into account the fractional binary representation of $\\varphi$:\n",
    "\n",
    "\\begin{align*}\n",
    "\\left| \\psi_a \\right\\rangle = \\frac{1}{2^{5/2}}\n",
    "&\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 0.\\varphi_0 \\varphi_1 \\varphi_2 \\varphi_3 \\varphi_4} \\left| 1 \\right\\rangle \\right)\n",
    "\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 0.\\varphi_1 \\varphi_2 \\varphi_3 \\varphi_4} \\left| 1 \\right\\rangle \\right) \n",
    "\\\\\n",
    "&\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 0.\\varphi_2 \\varphi_3 \\varphi_4} \\left| 1 \\right\\rangle \\right)\n",
    "\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 0.\\varphi_3 \\varphi_4} \\left| 1 \\right\\rangle \\right)\n",
    "\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 0.\\varphi_4} \\left| 1 \\right\\rangle \\right)\n",
    "\\end{align*}.\n",
    "\n",
    "Here is the code for the controlled $U^{2^{n}}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n in range(a.size):\n",
    "    for m in range(2**n):\n",
    "        get_controlled_circuit(unitary, a[n], circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the next step, we will apply a Hadamard gate to the last $n = 5$ qubit and controlled $u1 \\left( -\\frac{2 \\pi}{2^{(n - m + 1)}} \\right)$ gates from $n$th qubit to the rest $m < n$ qubits. Here is the $u1$ gate:\n",
    "\n",
    "$$\n",
    "u1(\\lambda) = \n",
    "\\begin{pmatrix}\n",
    "1 & 0 \\\\\n",
    "0 & e^{i\\lambda}.\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "This will give us a quantum state: \n",
    "\n",
    "\\begin{align*}\n",
    "\\left| \\psi_a \\right\\rangle = \\frac{1}{2^{5/2}}\n",
    "&\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 0.\\varphi_0 \\varphi_1 \\varphi_2 \\varphi_3} \\left| 1 \\right\\rangle \\right)\n",
    "\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 0.\\varphi_1 \\varphi_2 \\varphi_3} \\left| 1 \\right\\rangle \\right) \\\\\n",
    "&\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 0.\\varphi_2 \\varphi_3} \\left| 1 \\right\\rangle \\right)\n",
    "\\left(\\left| 0 \\right\\rangle + e^{i 2 \\pi 0.\\varphi_3} \\left| 1 \\right\\rangle \\right)\n",
    "\\left| \\varphi_4 \\right\\rangle.\n",
    "\\end{align*}\n",
    "\n",
    "Here we take into account that when $\\varphi_4 = 0$ the state of the $5$th qubit is equal to $\\left| + \\right\\rangle$ state and when $\\varphi_4 = 1$ it equal to $\\left| - \\right\\rangle = \\frac{1}{\\sqrt{2}}\\left( \\left| 0 \\right\\rangle - \\left| 1 \\right\\rangle \\right) $ state. Also note that when we apply Hadamard gate on $5$th qubit, it will give us either\n",
    "\n",
    "\\begin{align*}\n",
    "H \\left| + \\right\\rangle= \\left| 0 \\right\\rangle\n",
    "\\qquad \\text{ or } \\qquad\n",
    "H \\left| - \\right\\rangle = \\left| 1 \\right\\rangle\n",
    "\\end{align*}\n",
    "\n",
    "We repeat this procedure for all qubits and at the end, the quantum state will look like: \n",
    "\n",
    "\\begin{align*}\n",
    "\\left| \\psi_a \\right\\rangle=\n",
    "\\left| \\varphi_0 \\right\\rangle\n",
    "\\left| \\varphi_1 \\right\\rangle\n",
    "\\left| \\varphi_2 \\right\\rangle\n",
    "\\left| \\varphi_3 \\right\\rangle\n",
    "\\left| \\varphi_4 \\right\\rangle.\n",
    "\\end{align*}\n",
    "\n",
    "So, by measuring each qubit we will obtain all bits of the phase. This part of the algorithm is the well-known inverse Quantum Fourier transform (inverse QFT) [[3](https://en.wikipedia.org/wiki/Quantum_Fourier_transform)] without SWAP gates. Here is the code that adds this part of the algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# inverse QFT without SWAP gates\n",
    "for n in reversed(range(a.size)):\n",
    "    circuit.h(a[n])\n",
    "    if n != 0:\n",
    "        for m in reversed(range(n)):\n",
    "            angle = -2*np.pi / (2**(n - m + 1))\n",
    "            circuit.cu1(angle, a[n], a[m])\n",
    "\n",
    "# measurements on the ancillary qubits stored in c classical register\n",
    "for n in reversed(range(a.size)):\n",
    "    circuit.measure(a[n],c[n])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantum circuit that corresponds to the whole algorithm looks like this (here we draw a circuit that has only 4 ancillary qubits):\n",
    "\n",
    "<img src=\"images/qpe_circuit_4ancillary.png\" alt=\"drawing\" width=\"800\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the final piece of code, we run the algorithm on the QC simulator and store measurement results (bits of phase) in $c$ classical register that gives us the estimated value for the phase. After it, by remembering that $E_2 = \\frac{2 \\pi \\varphi_2}{t}$, we obtain desired $E_2$ eigenvalue of the $\\hat{H}$ Hamiltonian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalue of the Hamiltonian: 0.7294990675852743\n",
      "Estimated eigenvalue of the Hamiltonian: 0.7301748550335652\n"
     ]
    }
   ],
   "source": [
    "backend = BasicAer.get_backend('qasm_simulator')\n",
    "shots = 1024  # how many time execute the algorithm\n",
    "job = execute(circuit, backend, shots=shots)\n",
    "result = job.result()\n",
    "counts = result.get_counts()\n",
    "\n",
    "phase_bits = max(counts, key=counts.get) # take the most often obtaned result\n",
    "\n",
    "phase = 0\n",
    "for index, bit in enumerate(reversed(phase_bits)):\n",
    "    phase += int(bit) / 2**(index + 1)\n",
    "    \n",
    "estimated_E_2 = 2 * np.pi * phase / t\n",
    "\n",
    "print(\"Eigenvalue of the Hamiltonian: {}\".format(E_2))\n",
    "print(\"Estimated eigenvalue of the Hamiltonian: {}\".format(estimated_E_2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Final notes**: Here we didn't implement classical textbook version of QPE algorithm. In the textbook [[4](https://www.cambridge.org/am/academic/subjects/physics/quantum-physics-quantum-information-and-quantum-computation/quantum-computation-and-quantum-information-10th-anniversary-edition?format=HB)] the algorithm is done via inverse QFT with SWAP gates. If we want to implement the algorithm with inverse QFT, the code should be changed a little bit. Firstly, the order of control ```unitary``` gates should be changed: $i$th ancillary should control $U^{2^{(n - i - 1)}}$ gate on $q$ register. Secondly, SWAP gates are added after controlled unitaries of the algorithm in order to obtain inverse QFT. Here is the textbook [[4](https://www.cambridge.org/am/academic/subjects/physics/quantum-physics-quantum-information-and-quantum-computation/quantum-computation-and-quantum-information-10th-anniversary-edition?format=HB)] implementation of the algorithm:\n",
    "\n",
    "<img src=\"images/textbook_qpe.png\" alt=\"drawing\" width=\"800\"/>\n",
    "\n",
    "where the inverse QFT is the following part of the circuit:\n",
    "\n",
    "<img src=\"images/inverse_QFT.png\" alt=\"drawing\" width=\"800\"/>\n",
    "\n",
    "But we did the algorithm's job without SWAP gates. Actually, the SWAP gates don't add anything to this algorithm and we can omit them. Not only they increase the gate number of the circuit, but also they add errors. So, that is why we are avoiding SWAP gates in our implementation.\n",
    "\n",
    "They are some details about Trotterization, the $t$ parameter in the $U = e^{iHt}$, what should we do when the eigenvalue is negative that we don't cover here. If you interested in these topics you can check our [IQPE](https://github.com/DavitKhach/quantum-algorithms-tutorials/blob/master/iterative_quantum_phase_estimation.ipynb) tutorial, that is a modification of QPE algorithm. There are two main differences between QPE and IQPE algorithms. QPE demands $n$ ancillary qubits for $n$ bit phase estimation, IQPE demands only one ancillary qubit for any $n$ bit phase estimation. The second difference is that QPE applies inverse QFT, meanwhile, IQPE applies inverse QFT iteratively ($u1$ gates are applied iteratively)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] [Quantum phase estimation algorithm (wiki)](https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm)\n",
    "\n",
    "[J.D. Whitfield,  J. Biamonte  and  A. Aspuru-Guzik,\n",
    "Molecular Physics, \"Simulation of electronic structure Hamiltonians using\n",
    "quantum computers\" (2011)](https://www.tandfonline.com/doi/abs/10.1080/00268976.2011.552441)\n",
    "\n",
    "[3] [Quantum Fourier transform (wiki)](https://en.wikipedia.org/wiki/Quantum_Fourier_transform)\n",
    "\n",
    "[4] [M.A. Nielsen, I.L. Chuang, Cambridge University Press New York, \"Quantum Computation and Quantum Information: 10th Anniversary Edition\n",
    "10th\" (2011)](https://www.cambridge.org/am/academic/subjects/physics/quantum-physics-quantum-information-and-quantum-computation/quantum-computation-and-quantum-information-10th-anniversary-edition?format=HB)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}