iterative_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",
    "Iterative quantum phase estimation\n",
    "</h1>\n",
    "\n",
    "The iterative quantum phase estimation (IQPE, also known as IPEA) algorithm is a quantum algorithm for estimating the phase (or eigenvalue) of an eigenvector of a unitary operator [[1](https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm), [2](https://www.tandfonline.com/doi/abs/10.1080/00268976.2011.552441)]. One of the main applications of the algorithm is to estimate eigenvalues (energies) of some molecule's $H$ Hamiltonian. Because $H$ is a Hermitian operator, not unitary, (the algorithm works only with unitary operators) we can't estimate directly its eigenvalues. Instead, we create some unitary operator from $H$ and estimate not the eigenvalues of $H$, but something different (the phase). From the estimated phase one can obtain the corresponding eigenvalue of $H$. So, in the end, we are not only estimating the phase but, what is more important, the desired eigenvalue. Here are the main steps of the algorithm:\n",
    "\n",
    "1) Create unitary operator $U$ from given $H$: \n",
    "$$U = e^{iHt}.$$ \n",
    "Here $t$ is a parameter.\n",
    "   \n",
    "2) Apply $U$ on the eigenstate $\\left| \\psi_k \\right\\rangle$ of $H$. For simplicity, we will assume that the eigenstate is given: \n",
    "\n",
    "$$U \\left| \\psi_k \\right\\rangle = e^{iHt} \\left| \\psi_k \\right\\rangle$$ \n",
    "\n",
    "From the Taylor series for the exponent we have:\n",
    "\n",
    "$$e^{iHt} \\left| \\psi_k \\right\\rangle = i t H \\left| \\psi_k \\right\\rangle - \\frac{t^2}{2} H^2 \\left| \\psi_k \\right\\rangle + \\frac{t^4}{4} H^4 \\left| \\psi_k \\right\\rangle - ...$$\n",
    "\n",
    "And by using the Schrödinger equation $H \\left| \\psi_k \\right\\rangle = E_{k} \\left| \\psi_k \\right\\rangle$, we will obtain: \n",
    "\n",
    "$$e^{iHt} \\left|\\psi_k \\right\\rangle = e^{iE_{k}t} \\left|\\psi_k \\right\\rangle,$$ \n",
    "\n",
    "where $E_k$ is the eigenvalue of the corresponding $\\left| \\psi_k \\right\\rangle$ eigenstate. The goal of the algorithm is to find $E_k$s.\n",
    "\n",
    "3) It easy to see that:\n",
    "\n",
    "$$e^{iHt}\\left|\\psi_k \\right\\rangle = e^{iE_{k}t} \\left|\\psi_k \\right\\rangle  = e^{i 2 \\pi \\varphi_k} \\left|\\psi_k \\right\\rangle,$$\n",
    "\n",
    "where $\\varphi_k = E_{k}t / 2 \\pi$ is the phase that the algorithm is capable to estimate. After estimating the phase it will be easy to estimate the corresponding eigenvalue:\n",
    "\n",
    "$$ E_{k} = 2 \\pi \\varphi_k / t$$\n",
    "\n",
    "Now let's see how the algorithm actually works. Firstly, we should import all the packages that we will use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we are choosing a quantum computer simulator that we will use (namely ```'qasm_simulator'```). ```a``` and ```q``` are quantum registers. ```q``` register is the main register where $\\left| \\psi_k \\right\\rangle$ eigenstate will be encoded. $a$ register contains one ancillary qubit, where, after IPEA algorithm, the phase will be stored. In the end, only the ancillary qubit is measured and the measurement result will be stored in ```c``` classical register."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "backend = BasicAer.get_backend('qasm_simulator')\n",
    "q = QuantumRegister(1)\n",
    "a = QuantumRegister(1)\n",
    "c = ClassicalRegister(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are going to choose two random eigenvalues for a diagonal Hamiltonian. We should choose them such that $E_1, E_2 < 2 \\pi$ (later you will see why). The Hamiltonian in matrix form looks like this:\n",
    "\n",
    "$$H = \\begin{pmatrix}\n",
    "E_1 & 0\\\\\n",
    "0 & E_2\n",
    "\\end{pmatrix} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We are going to estimate E_2 via IQPE algorithm \n",
      "E_2 = 4.299125626432997\n"
     ]
    }
   ],
   "source": [
    "E_1, E_2 = (2 * np.pi * random(), 2 * np.pi * random())\n",
    "print(\"We are going to estimate E_2 via IQPE algorithm \\nE_2 = {}\".format(E_2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next step, the Hamiltonian should be exponentiated and applied on $\\left| 1 \\right\\rangle$ eigenstate in order to find the $E_2$ eigenvalue (if you want you can change $\\left| 1 \\right\\rangle$ to $\\left| 0 \\right\\rangle$ and find the $E_1$ eigenvalue):\n",
    "\n",
    "$$\n",
    "e^{iHt}  \\left|1 \\right\\rangle = \\begin{pmatrix}\n",
    "e^{i E_1 t} & 0\\\\\n",
    "0 & e^{i E_2 t}\n",
    "\\end{pmatrix} \n",
    "\\begin{pmatrix}\n",
    "0\\\\\n",
    "1\n",
    "\\end{pmatrix}= e^{i E_2 t}\n",
    "\\begin{pmatrix}\n",
    "0\\\\\n",
    "1\n",
    "\\end{pmatrix}.$$\n",
    "\n",
    "Note that $e^{ix}$ is an periodic function, so we should choose $t$ parameter such that $E_1 t < 2 \\pi$ and $E_2 t < 2 \\pi$ (otherwise, we will lose some information about the eigenvalues). Because we already have chosen $E_1, E_2 < 2 \\pi$, $t = 1$ will work for us. In the code, we should construct a circuit that corresponds to $U = e^{iHt}$ ```unitary``` operator. The circuit looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<qiskit.circuit.instructionset.InstructionSet at 0x7f938f379780>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 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": [
    "Where \n",
    "\n",
    "$$\n",
    "u1(\\lambda) = \\begin{pmatrix}\n",
    "1 & 0\\\\\n",
    "0 & e^{i\\lambda}\n",
    "\\end{pmatrix},\n",
    "\\qquad\n",
    "X = \\begin{pmatrix}\n",
    "0 & 1\\\\\n",
    "1 & 0\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "are quantum gates. You can check for yourself that for an arbitrary one qubit quantum state $\\left| \\psi \\right\\rangle = c_1 \\left| 0 \\right\\rangle + c_2 \\left| 1 \\right\\rangle$ the gate sequence presented in the above code will produce the desired  $e^{iHt} \\left| \\psi \\right\\rangle$ evolution:\n",
    "\n",
    "$$\n",
    "e^{iHt} \\left| \\psi \\right\\rangle = c_1 e^{iE_1t} \\left| 0 \\right\\rangle + c_2 e^{iE_2t} \\left| 1 \\right\\rangle\n",
    "$$\n",
    "\n",
    "The next code creates a controlled version of the ```unitary``` that we will use in the algorithm. Now we will try to give you an intuition of why we use it. First of all, let's see what happens when we apply controlled ```unitary``` ($CU$) on a two-qubit arbitrary quantum state:\n",
    "\n",
    "$$\n",
    "\\left| \\psi \\right\\rangle = c_1 \\left| 00 \\right\\rangle + c_2 \\left| 01 \\right\\rangle + c_3 \\left| 10 \\right\\rangle + c_4 \\left| 11 \\right\\rangle\n",
    "$$\n",
    "\n",
    "The ```unitary``` will be applied only if the first qubit is in $\\left| 1 \\right\\rangle$, otherwise, it will not do anything. The quantum state after $CU$:\n",
    "\n",
    "$$\n",
    "\\left| \\psi \\right\\rangle = c_1 \\left| 00 \\right\\rangle + c_2 \\left| 01 \\right\\rangle + c_3 e^{iE_1t} \\left| 10 \\right\\rangle + c_4 e^{iE_2t} \\left| 11 \\right\\rangle.\n",
    "$$\n",
    "\n",
    "And if you apply $CU$ some $n$ times you will obtain:\n",
    "\n",
    "$$\n",
    "\\left| \\psi \\right\\rangle = c_1 \\left| 00 \\right\\rangle + c_2 \\left| 01 \\right\\rangle + c_3 e^{i n E_1t} \\left| 10 \\right\\rangle + c_4 e^{i n E_2t} \\left| 11 \\right\\rangle.\n",
    "$$\n",
    "\n",
    "We represent the phase in the fractional binary representation (it's up to us how to represent a number):\n",
    "\n",
    "$$\\varphi_2 = 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",
    "\n",
    "where each $\\varphi_k$ is a bit (it is equal to $0$ or $1$). The algorithm at each iteration will estimate one $\\varphi_k$ bit from the fractional binary representation of the phase.\n",
    "\n",
    "At each iteration we are going to prepare our ancillary qubit in $\\left| + \\right\\rangle = \\frac{1}{\\sqrt{2}} \\left( \\left| 0 \\right\\rangle_a + \\left| 1 \\right\\rangle_a \\right)$ state. For that we apply Hadamard $H$ gate to the ancillary qubit:\n",
    "\n",
    "$$\n",
    "H = \\frac{1}{\\sqrt{2}}\\begin{pmatrix}\n",
    "1 & 1\\\\\n",
    "1 & -1\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "Then the qubit in the $q$ register is prepared in the $\\left| 1 \\right\\rangle$ eigenstate of our diagonal Hamiltonian, which corresponding eigenvalue we want to estimate. For that, we apply $X$ gate to the qubit in the $q$ register. Here is the obtained quantum state after Hadamard and $X$ gates:\n",
    "\n",
    "$$\n",
    "\\left|+ \\right\\rangle_a \\left| 1 \\right\\rangle_q = \\frac{1}{\\sqrt{2}} \\left( \\left| 0 \\right\\rangle_a + \\left| 1 \\right\\rangle_a \\right) \\left| 1 \\right\\rangle_q = \\frac{1}{\\sqrt{2}} \\left( \\left| 01 \\right\\rangle +  \\left| 11 \\right\\rangle \\right)\n",
    "$$\n",
    "\n",
    "If we apply controlled ```unitary``` circuit it will produce the following quantum state:\n",
    "\n",
    "$$\n",
    "\\left| + \\right\\rangle_a \\left| 1 \\right\\rangle_q \\xrightarrow{\\text{CU}}\n",
    "\\frac{1}{\\sqrt{2}} \\left( \\left| 01 \\right\\rangle + e^{-i 2 \\pi E_2 t} \\left| 11 \\right\\rangle \\right) = \n",
    "\\left( \\left| 0 \\right\\rangle_a + e^{-i 2 \\pi \\varphi_2} \\left| 1 \\right\\rangle_a \\right) \\left| 1 \\right\\rangle_q,\n",
    "$$\n",
    "\n",
    "because $\\varphi_2 = E_2 t$. If we apply controlled ```unitary``` circuit $2^{k}$ times the state of the qubits will look like this:\n",
    "\n",
    "$$\n",
    "\\left| + \\right\\rangle_a \\left| 1 \\right\\rangle_q \\xrightarrow{\\text{CU}}\n",
    "\\left( \\left| 0 \\right\\rangle_a + e^{-i 2 \\pi \\left( 2^{n} \\varphi_2 \\right)} \\left| 1 \\right\\rangle_a \\right) \\left| 1 \\right\\rangle_q = \n",
    "\\left( \\left| 0 \\right\\rangle_a + e^{-i 2 \\pi 0.\\varphi_{k+1}...\\varphi_{n}} \\left| 1 \\right\\rangle_a \\right) \\left| 1 \\right\\rangle_q,\n",
    "$$\n",
    "\n",
    "where $n$ is the index of the last bit (for $5$ qubits $n=4$ ). This will let us store different portions of the phase in the ancillary qubit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "control_u = get_controlled_circuit(unitary, a[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the first iteration, we estimate the last bit $\\varphi_n$. For that, we are applying controlled ```unitary``` $2^{n}$ times. After this we will have the following quantum state:\n",
    "\n",
    "$$\n",
    "\\left( \\left| 0 \\right\\rangle_a + e^{-i 2 \\pi 0.\\varphi_{n}} \\left| 1 \\right\\rangle_a \\right) \\left| 1 \\right\\rangle_q\n",
    "$$\n",
    "\n",
    "Estimating $\\varphi_{n}$ from this state is straightforward. Depending on $\\varphi_{n}$ there are two possible quantum states for the ancillary qubit:\n",
    "\n",
    "$$\n",
    "\\left| + \\right\\rangle_a = \\frac{1}{\\sqrt{2}} \\left( \\left| 0 \\right\\rangle_a + \\left| 1 \\right\\rangle_a \\right)\n",
    "\\qquad \\text{ or } \\qquad\n",
    "\\left| - \\right\\rangle_a = \\frac{1}{\\sqrt{2}} \\left( \\left| 0 \\right\\rangle_a - \\left| 1 \\right\\rangle_a \\right)\n",
    "$$\n",
    "\n",
    "The next step is applying $H$ gate. Remember that $H \\left| + \\right\\rangle_a = \\left| 0 \\right\\rangle$ and $H \\left| - \\right\\rangle_a = \\left| 1 \\right\\rangle$ and after the measurement, we will have either $0$ or $1$ outcome. If the measurmet yields $0$, then $\\varphi_n = 0$ and if it yields $1$, then $\\varphi_n = 1$.\n",
    "\n",
    "After this step, the same procedure is done for $n-1$ bit of the $\\varphi$. This time we apply controlled ```unitary``` $2^{n-1}$ times. The quantum state:\n",
    "\n",
    "$$\n",
    "\\left( \\left| 0 \\right\\rangle_a + e^{-i 2 \\pi 0.\\varphi_{n-1}\\varphi_{n}} \\left| 1 \\right\\rangle_a \\right) \\left| 1 \\right\\rangle_q\n",
    "$$\n",
    "\n",
    "Before the Hadamard gate, we apply ```u1``` gate that subtracts from the phase already measured/known $\\varphi_{n}$. After ```u1``` gate:\n",
    "\n",
    "$$\n",
    "\\left( \\left| 0 \\right\\rangle_a + e^{-i 2 \\pi 0.\\varphi_{n-1}} \\left| 1 \\right\\rangle_a \\right) \\left| 1 \\right\\rangle_q\n",
    "$$\n",
    "\n",
    "With this wavefunction, it is easy to estimate the $\\varphi_{n-1}$. For $k$th iteration the circuit looks like this:\n",
    "\n",
    "<img src=\"images/iqpe_circuit.png\" alt=\"drawing\" width=\"400\"/>\n",
    "\n",
    "where $\\varphi^{'} = (0.0\\varphi_{k + 1}\\varphi_{k + 2}\\varphi_n)$.\n",
    "\n",
    "The following ```for``` loop iteratively creates circuits for IQPE iterations. Every iteration estimates one bit of the $\\varphi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_bits_estimate = 10\n",
    "phase = 0\n",
    "for k_precision in reversed(range(num_bits_estimate)):\n",
    "    # Create a Quantum Circuit acting on the q register\n",
    "    k_circ = QuantumCircuit(q, a, c)\n",
    "\n",
    "    # (1) |1> eigenstate initialization\n",
    "    k_circ.x(q[0])\n",
    "\n",
    "    # (2) Initial Hadamard gate applied on ancillary qubit.\n",
    "    k_circ.h(a[0])\n",
    "\n",
    "    # (3) The control Hamiltonian applied on the qubits where control qubit is the ancillary qubit.\n",
    "    for order in range(2 ** k_precision):\n",
    "        k_circ += control_u\n",
    "\n",
    "    # (4) The phase gate and final Hadamard gate on ancillary qubit.\n",
    "    phase_shift = 2 * np.pi * phase * 2 ** k_precision\n",
    "    k_circ.u1(-phase_shift, a[0])\n",
    "    k_circ.h(a[0])\n",
    "\n",
    "    # (5) Measurement of ancillary qubit (findig the bit)\n",
    "    k_circ.measure(a[0], c[0])\n",
    "\n",
    "    # (6) executing on Quantum Computer and finding a bit from the phase\n",
    "    job = execute(k_circ, backend, shots=1000)\n",
    "    result = job.result()\n",
    "    counts = result.get_counts()\n",
    "    value = int(max(counts, key=counts.get))\n",
    "\n",
    "    # (7) phase after this iteration\n",
    "    phase += value / 2 ** (k_precision + 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we are estimating the eigenvalue (energy) of the $\\left| 1 \\right\\rangle$ eigenstate and printing the result alongside the true eigenvalue $E_2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalue of the Hamiltonian: 4.299125626432997\n",
      "Estimated eigenvalue of the Hamiltonian: 4.301282129231338\n"
     ]
    }
   ],
   "source": [
    "eigenvalue = 2 * np.pi * phase / t\n",
    "print(\"Eigenvalue of the Hamiltonian: {}\".format(E_2))\n",
    "print(\"Estimated eigenvalue of the Hamiltonian: {}\".format(eigenvalue))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We omitted one step of the algorithm. That is the Trotter decomposition step. For some Hamiltonian $H = h_1 + h_2 + h_3$ we need to exponentiate it to obtain $U = e^{i H t}$ unitary operator. Actually with this algorithm (namely when we are constructing a circuit with ```evolve``` method) we obtain $\\tilde{U} = e^{i h_1 t} \\cdot e^{i h_2 t} \\cdot e^{i h_3 t}$. Actually, $e^{iHt} \\ne  e^{i h_1 t} \\cdot e^{i h_2 t} \\cdot e^{i h_3 t}$. To obtain more or less valid $U$ operator we should do Trotter decomposition:\n",
    "\n",
    "$$ U = (e^{i h_1 \\frac{t}{N}} \\cdot e^{i h_2 \\frac{t}{N}} \\cdot e^{i h_3 \\frac{t}{N}})^N + O(\\frac{t^2}{N}),$$\n",
    "\n",
    "where $N$ is the Trotter number [[2](https://www.tandfonline.com/doi/abs/10.1080/00268976.2011.552441)]. By increasing $N$ it is possible to obtain better representation for $U$ operator. Because our example of the Hamiltonian was simple (2x2 diagonal matrix) a reasonable result was obtained without doing this step.\n",
    "\n",
    "One more thing. If $E_2 < 0$ then the final step of the algorithm is changed a little bit. In that case, note that $e^{iE_2t} = e^{-i \\left|E_2 \\right| t } = e^{i \\left( 2 \\pi-\\left|E_2 \\right| t \\right) } = e^{i \\varphi_2^* }$. Because of this, the algorithm will estimate $\\varphi_2^* = 2 \\pi-\\left|E_2 \\right| t$ and in the final step instead of $E_2 = 2  \\pi  \\varphi_2^* / t$ we should write $E_2 = -2  \\pi  (1 - \\varphi_2^*) / t$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] [Quantum phase estimation algorithm](https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm)\n",
    "\n",
    "[2]  [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)"
   ]
  }
 ],
 "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
}