variational_quantum_eigensolver.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[<img src=\"images/quantum_algorithms_tutorials.png\" alt=\"drawing\" width=\"100\" align=\"left\"/>][7]\n",
    "\n",
    "<h1 align=\"center\">\n",
    "\tVariational quantum eigensolver\n",
    "</h1>\n",
    "\n",
    "**[[Homepage][7]]**\n",
    "**[[Open with the nbviewer][8]]**\n",
    "\n",
    "Variational quantum eigensolver (VQE) is a hybrid quantum-classical algorithm [[1]] that finds the smallest eigenvalue (and corresponding eigenvector) of a given Hamiltonian. One of the main applications of the algorithm is finding ground state energy of molecules. It has a big advantage over [IQPE][5] (iterative quantum phase estimation) and [QPE][6] (quantum phase estimation) algorithms, that also can be used for finding the ground state energy of a molecule. The main advantage is that VQE uses  much smaller circuit depths (or gates) then IQPE and QPE, what is very important for NISQ (Noisy Intermediate-Scale Quantum) era quantum computation. In the NISQ era (now!) we are working with qubits that are very noisy because they are not isolated from the environment well enough. Thus, there is small and finite time to work with qubits until they will be \"spoiled\", because of the environment, imperfect gates and etc. This restriction gives a big advantage to those algorithms (like VQE) that are using small depth circuits. \n",
    "\n",
    "The idea of the VQE algorithm is as follows. We have a Hamiltonian that can be expressed by the sum of tensor products of Pauli operators (Pauli terms):\n",
    "\n",
    "$$H = 0.4 \\cdot IX + 0.6 \\cdot IZ + 0.8 \\cdot XY.$$\n",
    "\n",
    "For a given $\\left| \\psi \\right\\rangle$ we want to measure the expectation value of the Hamiltonian:\n",
    "\n",
    "$$\\left\\langle H \\right\\rangle = \\left\\langle \\psi \\right| H \\left| \\psi \\right\\rangle = 0.4 \\cdot \\left\\langle \\psi \\right| IX \\left| \\psi \\right\\rangle + 0.6 \\cdot \\left\\langle \\psi \\right| IZ \\left| \\psi \\right\\rangle + 0.8 \\cdot \\left\\langle \\psi \\right| XY \\left| \\psi \\right\\rangle.$$\n",
    "\n",
    "How one can see the $\\left\\langle H \\right\\rangle$ expectation value could be computed by adding the expectation values of its parts (Pauli terms). The algorithm does exactly that. It constructs a quantum circuit for each Pauli term and computes the expectation value of the corresponding Pauli term. Then, the algorithm sums all calculated expectation values of Pauli terms and obtains the expectation value of $H$. In this algorithm, we will do this routine of estimating the expectation value of $H$ over and over again for different trial wavefunctions (ansatz states) $\\left| \\psi \\right\\rangle$. \n",
    "\n",
    "It is known that the eigenvector $\\left| \\psi_g \\right\\rangle$ that minimizes the expectation value $\\left\\langle H \\right\\rangle$ corresponds to the eigenvector of $H$ that has the smallest eigenvalue [[1]]. So, basically we can try all possible trial wavefunctions $\\left| \\psi \\right\\rangle$s to find the $\\left| \\psi_g \\right\\rangle$ that has the smallest expectation value. Here the question is how we create those trial states? In the algorithm, the trial states are created from a parametrized circuit. By changing the parameters one obtains different wavefunctions (ansatz states) [[1]]. If your circuit with its parameters is good enough you will have access to the subspace of the states [[2]] that includes the $\\left| \\psi_g \\right\\rangle$. Otherwise, if the circuit will not have a possibility to generate our desired $\\left| \\psi_g \\right\\rangle$ it will be impossible to find the right solution.\n",
    "\n",
    "The parameters of the state preparation circuit are controlled by a classical computer. At each step, the classical computer will change the parameters by using some optimization method in order to create an ansatz state that will have a smaller expectation value then previous ansatz states had. This way the classical computer and the quantum computer are working together to archive the goal of the algorithm (to find the ground state energy). That's way, VQE is a quantum-classical hybrid algorithm.\n",
    "\n",
    "<img src=\"images/vqe_parts.png\" alt=\"drawing\" width=\"600\"/>\n",
    "\n",
    "This image is taken from the [[1]] paper. Here are the main parts of the VQE algorithm. The *Classical feedback decision* is some optimization method that changes the parameters of the *Quantum state preparation*. With different *Quantum modules* the algorithm calculates expectation values of each Pauli term and then it sums them by using *Classical adder* via classical computer. Then the algorithm returns to *Classical feedback decision* to choose better parameters for the *Quantum state preparation*. VQE repeats this procedure until the optimization method is satisfied with the obtained result. Note, that we will use these names in the image as method names in our code.\n",
    "\n",
    "  [1]: https://www.nature.com/articles/ncomms5213?origin=ppub\n",
    "  [2]: https://www.mustythoughts.com/variational-quantum-eigensolver-explained\n",
    "  [3]: 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\n",
    "  [4]: https://qiskit.org/textbook/ch-applications/vqe-molecules.html\n",
    "  [5]: https://github.com/DavitKhach/quantum-algorithms-tutorials/blob/master/iterative_quantum_phase_estimation.ipynb\n",
    "  [6]: https://github.com/DavitKhach/quantum-algorithms-tutorials/blob/master/quantum_phase_estimation.ipynb\n",
    "  [7]: https://github.com/DavitKhach/quantum-algorithms-tutorials\n",
    "  [8]: https://nbviewer.jupyter.org/github/DavitKhach/quantum-algorithms-tutorials/blob/master/variational_quantum_eigensolver.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see how the algorithm 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 scipy.optimize import minimize\n",
    "\n",
    "from qiskit import *\n",
    "from qiskit.circuit.library.standard_gates import U2Gate\n",
    "from qiskit.aqua.operators import WeightedPauliOperator\n",
    "from qiskit.aqua.algorithms import NumPyEigensolver"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simplicity we will consider one qubit Hamiltonian consisting of a sum of 4 Pauli operators:\n",
    "\n",
    "$$H = a \\cdot I + b \\cdot Z + c \\cdot X + d \\cdot Y,$$\n",
    "\n",
    "here are matrix representations of each Pauli operators:\n",
    "\n",
    "$$\n",
    "I = \\begin{pmatrix}\n",
    "1 & 0\\\\\n",
    "0 & 1\n",
    "\\end{pmatrix},\n",
    "\\qquad\n",
    "Z = \\begin{pmatrix}\n",
    "1 & 0\\\\\n",
    "0 & -1\n",
    "\\end{pmatrix},\n",
    "\\qquad\n",
    "X = \\begin{pmatrix}\n",
    "0 & 1\\\\\n",
    "1 & 0\n",
    "\\end{pmatrix},\n",
    "\\qquad\n",
    "Y = \\begin{pmatrix}\n",
    "0 & -i\\\\\n",
    "i & 0\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "The following method creates $H$ for a given $a$, $b$, $c$, $d$ coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hamiltonian_operator(a, b, c, d):\n",
    "    \"\"\"\n",
    "    Creates a*I + b*Z + c*X + d*Y pauli sum \n",
    "    that will be our Hamiltonian operator.\n",
    "    \n",
    "    \"\"\"\n",
    "    pauli_dict = {\n",
    "        'paulis': [{\"coeff\": {\"imag\": 0.0, \"real\": a}, \"label\": \"I\"},\n",
    "                   {\"coeff\": {\"imag\": 0.0, \"real\": b}, \"label\": \"Z\"},\n",
    "                   {\"coeff\": {\"imag\": 0.0, \"real\": c}, \"label\": \"X\"},\n",
    "                   {\"coeff\": {\"imag\": 0.0, \"real\": d}, \"label\": \"Y\"}\n",
    "                   ]\n",
    "    }\n",
    "    return WeightedPauliOperator.from_dict(pauli_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the code the $a$, $b$, $c$, $d$ coefficients are chosen as random real numbers from the range $[0, 10]$. Thus we will have randomly generated Hamiltonian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "scale = 10\n",
    "a, b, c, d = (scale*random(), scale*random(), \n",
    "              scale*random(), scale*random())\n",
    "H = hamiltonian_operator(a, b, c, d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we are using Qiskit's ```NumPyEigensolver``` class (like was done here [[4]]) that will find the smallest eigenvalue of the given Hamiltonian with a classical algorithm. This value will be used as a reference for comparing it with the VQE result.\n",
    "\n",
    "  [1]: https://www.nature.com/articles/ncomms5213?origin=ppub\n",
    "  [2]: https://www.mustythoughts.com/variational-quantum-eigensolver-explained\n",
    "  [3]: 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\n",
    "  [4]: https://qiskit.org/textbook/ch-applications/vqe-molecules.html\n",
    "  [5]: https://github.com/DavitKhach/quantum-algorithms-tutorials/blob/master/iterative_quantum_phase_estimation.ipynb\n",
    "  [6]: https://github.com/DavitKhach/quantum-algorithms-tutorials/blob/master/quantum_phase_estimation.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The exact ground state energy is: -5.775991418998328\n"
     ]
    }
   ],
   "source": [
    "exact_result = NumPyEigensolver(H).run()\n",
    "reference_energy = min(np.real(exact_result.eigenvalues))\n",
    "print('The exact ground state energy is: {}'.format(reference_energy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to create \"good\" ansatz states. First of all, note that every possible wavefunction $\\left| \\psi \\right\\rangle$ can be presented as a vector: \n",
    "$$\n",
    "\\left| \\psi \\right\\rangle = \\begin{pmatrix}\n",
    "\\cos{\\left( \\theta/2 \\right)}\\\\\n",
    "e^{i \\varphi} \\cdot \\sin{\\left( \\theta/2 \\right)}\n",
    "\\end{pmatrix},\n",
    "$$\n",
    "\n",
    "where the numbers $\\theta$ and $\\varphi$ define a point on the unit three-dimensional sphere that is called Bloch sphere and is presented in the following picture from this book [[3]].\n",
    "\n",
    "<img src=\"images/bloch_sphere.png\" alt=\"drawing\" width=\"300\"/>\n",
    "\n",
    "For a random one qubit Hamiltonian, a \"good\" quantum state preparation circuit should be able to generate all possible states in the Bloch sphere [[2]]. Before quantum state preparation, our qubit is in $\\left| 0 \\right\\rangle = \\begin{pmatrix}\n",
    "1\\\\\n",
    "0\n",
    "\\end{pmatrix}$\n",
    "state. This corresponds to the vertical position of the vector in the Bloch sphere. In order to generate any possible $\\left| \\psi \\right\\rangle$ we will apply $R_x(t_1)$ and $R_y(t_2)$ gates on the $\\left| 0 \\right\\rangle$ initial state: $R_y(t_2) R_x(t_1) \\left| 0 \\right\\rangle = \\left| \\psi \\right\\rangle$. $R_x(t_1)$ corresponds to the rotation in the Bloch sphere around the *x* axis and $R_y(t_2)$ the rotation around the *y* axis. With these two rotations, one can have access to any point in the Bloch sphere. Here we show the matrix forms of $R_x(t_1)$ and $R_y(t_2)$ gates:\n",
    "\n",
    "$$\n",
    "R_x(t_1) = \\begin{pmatrix}\n",
    "cos(\\frac{t_1}{2}) & -i \\cdot sin(\\frac{t_1}{2})\\\\\n",
    "-i \\cdot sin(\\frac{t_1}{2}) & cos(\\frac{t_1}{2})\n",
    "\\end{pmatrix},\n",
    "\\qquad\n",
    "R_y(t_2) = \\begin{pmatrix}\n",
    "cos(\\frac{t_2}{2}) & -sin(\\frac{t_2}{2})\\\\\n",
    "sin(\\frac{t_2}{2}) & cos(\\frac{t_2}{2})\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "These two gates with there parameters ($t_1$ and $t_2$) will generate for us the trial (ansatz) wavefunctions. The two parameters will be in control of the Classical Computer and its optimization model. Here is the corresponding circuit:\n",
    "\n",
    "<img src=\"images/vqe_one_qubit.png\" alt=\"drawing\" width=\"300\"/>\n",
    "\n",
    "Note that in the general case if we know for some reason that $\\left| \\psi_g \\right\\rangle$ is in the $M$ subspace of the Hilbert space, we will not need to create a parameterized circuit that will have access to all possible states in the Hilbert space. It will be sufficient to have an ansatz circuit that will be able to create all quantum states in the $M$ subspace. \n",
    "\n",
    "The following function corresponds to the quantum state preparation circuit for given parameters ($t_1$ and $t_2$).\n",
    "\n",
    "  [1]: https://www.nature.com/articles/ncomms5213?origin=ppub\n",
    "  [2]: https://www.mustythoughts.com/variational-quantum-eigensolver-explained\n",
    "  [3]: 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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quantum_state_preparation(circuit, parameters):\n",
    "    q = circuit.qregs[0] # q is the quantum register where the info about qubits is stored\n",
    "    circuit.rx(parameters[0], q[0]) # q[0] is our one and only qubit XD\n",
    "    circuit.ry(parameters[1], q[0])\n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we come to the part of the algorithm where we should do individual/separate measurements of expectation values of the Pauli terms. For one qubit case, there are only four possible Pauli terms: 4 basic Pauli operators $I$, $X$, $Y$, $Z$.\n",
    "\n",
    "$$\\left\\langle H \\right\\rangle = \\left\\langle \\psi \\right| H \\left| \\psi \\right\\rangle = a \\cdot \\left\\langle \\psi \\right| I \\left| \\psi \\right\\rangle + b \\cdot \\left\\langle \\psi \\right| Z \\left| \\psi \\right\\rangle + c \\cdot \\left\\langle \\psi \\right| X \\left| \\psi \\right\\rangle + d \\cdot \\left\\langle \\psi \\right| Y \\left| \\psi \\right\\rangle.$$\n",
    "\n",
    "For $I$ operator the  expectation value is always unity: $\\left\\langle \\psi \\right| I \\left| \\psi \\right\\rangle = \\left\\langle \\psi \\right| \\left| \\psi \\right\\rangle = 1$. So, its contribution to the overall expectation value $\\left\\langle H \\right\\rangle$ will be equal to $a \\cdot \\left\\langle \\psi \\right| I \\left| \\psi \\right\\rangle = a$.\n",
    "\n",
    "For rest of the Pauli operators, we should make the following remark: every one qubit quantum state $\\left| \\psi \\right\\rangle$ can be represented via different sets of basis vectors:\n",
    "\n",
    "$$\\left| \\psi \\right\\rangle = c_1^z \\cdot \\left| 0 \\right\\rangle + c_2^z \\cdot \\left| 1 \\right\\rangle = c_1^x \\cdot \\left| + \\right\\rangle + c_2^x \\cdot \\left| - \\right\\rangle = c_1^y \\cdot \\left| +i \\right\\rangle + c_2^y \\cdot \\left| -i \\right\\rangle,$$\n",
    "\n",
    "where\n",
    "\n",
    "\\begin{align*}\n",
    "&\\text{Z eigenvectors} \\qquad\n",
    "\\left| 0 \\right\\rangle = \\begin{pmatrix}\n",
    "1\\\\\n",
    "0\n",
    "\\end{pmatrix},\n",
    "&&\\left| 1 \\right\\rangle = \\begin{pmatrix}\n",
    "0\\\\\n",
    "1\n",
    "\\end{pmatrix},\n",
    "\\\\\n",
    "&\\text{X eigenvectors} \\qquad\n",
    "\\left| + \\right\\rangle = \\frac{1}{\\sqrt{2}} \\begin{pmatrix}\n",
    "1\\\\\n",
    "1\n",
    "\\end{pmatrix},\n",
    "&&\\left| - \\right\\rangle = \\frac{1}{\\sqrt{2}} \\begin{pmatrix}\n",
    "1\\\\\n",
    "-1\n",
    "\\end{pmatrix},\n",
    "\\\\\n",
    "&\\text{Y eigenvectors} \\qquad\n",
    "\\left| +i \\right\\rangle = \\frac{1}{\\sqrt{2}} \\begin{pmatrix}\n",
    "1\\\\\n",
    "i\n",
    "\\end{pmatrix}, \n",
    "&&\\left| -i \\right\\rangle = \\frac{1}{\\sqrt{2}} \\begin{pmatrix}\n",
    "1\\\\\n",
    "-i\n",
    "\\end{pmatrix}.\n",
    "\\end{align*}\n",
    "\n",
    "The first presented eigenvectors for each Pauli has an eigenvalue equal to $+1$: $Z \\left| 0 \\right\\rangle = +1\\left| 0 \\right\\rangle$, $X \\left| + \\right\\rangle = +1\\left| + \\right\\rangle$, $Y \\left| +i \\right\\rangle = +1\\left| +i \\right\\rangle$. And the second presented eigenvectors for each Pauli has an eigenvalue equal to $-1$: $Z \\left| 1 \\right\\rangle = -1\\left| 1 \\right\\rangle$, $X \\left| - \\right\\rangle = -1\\left| - \\right\\rangle$, $Y \\left| -i \\right\\rangle = -1\\left| -i \\right\\rangle$. Now, let's calculate the expectation values of these Pauli operators: \n",
    "\n",
    "\\begin{align*}\n",
    "\\left\\langle \\psi \\right| Z \\left| \\psi \\right\\rangle &= \\left( {c_1^z}^* \\cdot \\left\\langle 0 \\right| + {c_2^z}^* \\cdot \\left\\langle 1 \\right| \\right) Z \\left( c_1^z \\cdot \\left| 0 \\right\\rangle + c_2^z \\cdot \\left| 1 \\right\\rangle \\right) = {\\left| c_1^z \\right|}^2 - {\\left| c_2^z \\right|}^2,\n",
    "\\\\\n",
    "\\left\\langle \\psi \\right| X \\left| \\psi \\right\\rangle &= \\left( {c_1^x}^* \\cdot \\left\\langle + \\right| + {c_2^x}^* \\cdot \\left\\langle - \\right| \\right) X \\left( c_1^x \\cdot \\left| + \\right\\rangle + c_2^x \\cdot \\left| - \\right\\rangle \\right) = {\\left| c_1^x \\right|}^2 - {\\left| c_2^x \\right|}^2,\n",
    "\\\\\n",
    "\\left\\langle \\psi \\right| Y \\left| \\psi \\right\\rangle &= \\left( {c_1^y}^* \\cdot \\left\\langle +i \\right| + {c_2^y}^* \\cdot \\left\\langle -i \\right| \\right) Y \\left( c_1^y \\cdot \\left| +i \\right\\rangle + c_2^y \\cdot \\left| -i \\right\\rangle \\right) = {\\left| c_1^y \\right|}^2 - {\\left| c_2^y \\right|}^2,\n",
    "\\end{align*}\n",
    "\n",
    "where we take into account that the inner product of orthonormal vectors is 0 (e.g. $\\left\\langle 0 \\right| \\left| 1 \\right\\rangle = 0$, $\\left\\langle + \\right| \\left| - \\right\\rangle = 0$, $\\left\\langle +i \\right| \\left| -i \\right\\rangle = 0$). But what are these $\\left| c \\right|^2$s? The ${\\left| c_1^z \\right|}^2$ and ${\\left| c_2^z \\right|}^2$ are by definition the probabilities that after Z basis measurement (measuring is it $\\left| 0 \\right\\rangle$ or is it $\\left| 1 \\right\\rangle$) the quantum state $\\left| \\psi \\right\\rangle$ will become $\\left| 0 \\right\\rangle$ or $\\left| 1 \\right\\rangle$ respectively. In order to find that value, we should run our program with our trial $\\left| \\psi \\right\\rangle$  wavefunction and do $Z$ measurement on the qubit $N$ times (its named ```shots``` in the code). The probability of finding the qubit after measurment in $\\left| 0 \\right\\rangle$ state will be equal to ${\\left| c_1^z \\right|}^2 = \\frac{n_0}{N}$, where $n_0$ is the number of the $\\left| 0 \\right\\rangle$ state measurments. Similarly, ${\\left| c_2^z \\right|}^2 = \\frac{n_1}{N}$, where $n_1$ is the number of the $\\left| 1 \\right\\rangle$ state measurments. Thus, the final expectation value will be $\\left\\langle Z \\right\\rangle = \\frac{n_0 - n_1}{N}$.\n",
    "\n",
    "For $\\left\\langle X \\right\\rangle = \\frac{n_+ - n_-}{N}$ and $\\left\\langle Y \\right\\rangle = \\frac{n_{+i} - n_{-i}}{N}$ the expectation value estimation procedure stays the same. Here $n_+$ and $n_-$ are numbers of measurements in X basis that corresponds to $\\left| + \\right\\rangle$ or $\\left| - \\right\\rangle$ outcomes respectively. And $n_{+i}$ and $n_{-i}$ are numbers of measurements in $Y$ basis that corresponds to $\\left| +i \\right\\rangle$ or $\\left| -i \\right\\rangle$ outcomes respectively. The difficulty comes from the fact that one may have the possibility to measure only in the $Z$ basis. To solve this difficulty we still do a $Z$ basis measurement, but, before that, we apply specific operators to the $\\left| \\psi \\right\\rangle$ state. We try to apply such an operator that after measuring the probability of $\\left| 0 \\right\\rangle$ outcome will be equal to the probability of $\\left| + \\right\\rangle$ $\\left( \\left| +i \\right\\rangle \\right)$ outcome. And the probability of $\\left| 1 \\right\\rangle$ outcome will be equal to the probability of $\\left| - \\right\\rangle$ $\\left( \\left| -i \\right\\rangle \\right)$ outcome. Let's define such operators:\n",
    "\n",
    "$$H_{gate} \\left| \\psi \\right\\rangle = H_{gate} \\left( c_1^x \\cdot \\left| + \\right\\rangle + c_2^x \\cdot \\left| - \\right\\rangle \\right) = c_1^x \\cdot \\left| 0 \\right\\rangle + c_2^x \\cdot \\left| 1 \\right\\rangle,$$\n",
    "\n",
    "where $H_{gate}$ is an operator such that $H_{gate} \\left| + \\right\\rangle = \\left| 0 \\right\\rangle$ and $H_{gate} \\left| - \\right\\rangle = \\left| 1 \\right\\rangle$.\n",
    "\n",
    "$$Y_{gate} \\left| \\psi \\right\\rangle = Y_{gate} \\left( c_1^y \\cdot \\left| +i \\right\\rangle + c_2^y \\cdot \\left| -i \\right\\rangle \\right) = c_1^y \\cdot \\left| 0 \\right\\rangle + c_2^y \\cdot \\left| 1 \\right\\rangle,$$\n",
    "\n",
    "where $Y_{gate}$ is an operator such that $Y_{gate} \\left| +i \\right\\rangle = \\left| 0 \\right\\rangle$ and $Y_{gate} \\left| -i \\right\\rangle = \\left| 1 \\right\\rangle$.\n",
    "\n",
    "This kind of operators can be easily found:\n",
    "\n",
    "$$\n",
    "H_{gate} = \\frac{1}{\\sqrt{2}}\\begin{pmatrix}\n",
    "1 & 1\\\\\n",
    "1 & -1\n",
    "\\end{pmatrix},\n",
    "\\qquad\n",
    "Y_{gate} = \\frac{1}{\\sqrt{2}}\\begin{pmatrix}\n",
    "1 & -i\\\\\n",
    "1 & i\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "BTW the $H_{gate}$ is the well known Hadamard gate (its common notation is $H$, but we don't use this notation because we already have used the $H$ for the Hamiltonian operator). In the following code you can see an implementation of these gates via qiskit's ```u2``` gate:\n",
    "\n",
    "$$\n",
    "u2(\\varphi, \\lambda) = \\frac{1}{\\sqrt{2}}\\begin{pmatrix}\n",
    "1 & -e^{i \\lambda}\\\\\n",
    "e^{i \\varphi} & e^{i (\\varphi + \\lambda)}\n",
    "\\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H_gate:\n",
      "[[ 1.+0.j  1.-0.j]\n",
      " [ 1.+0.j -1.+0.j]]\n",
      "Y_gate:\n",
      "[[ 1.+0.j -0.-1.j]\n",
      " [ 1.+0.j  0.+1.j]]\n"
     ]
    }
   ],
   "source": [
    "H_gate = U2Gate(0, np.pi).to_matrix()\n",
    "print(\"H_gate:\")\n",
    "print((H_gate * np.sqrt(2)).round(5))\n",
    "\n",
    "Y_gate = U2Gate(0, np.pi/2).to_matrix()\n",
    "print(\"Y_gate:\")\n",
    "print((Y_gate * np.sqrt(2)).round(5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we are ready to go! In this code, we define a method that will create a circuit with its specific measurement ($Z$ or $X$ or $Y$ measurments)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vqe_circuit(parameters, measure):\n",
    "    \"\"\"\n",
    "    Creates a device ansatz circuit for optimization.\n",
    "    :param parameters_array: list of parameters for constructing ansatz state that should be optimized.\n",
    "    :param measure: measurement type. E.g. 'Z' stands for Z measurement.\n",
    "    :return: quantum circuit.\n",
    "    \"\"\"\n",
    "    q = QuantumRegister(1)\n",
    "    c = ClassicalRegister(1)\n",
    "    circuit = QuantumCircuit(q, c)\n",
    "\n",
    "    # quantum state preparation\n",
    "    circuit = quantum_state_preparation(circuit, parameters)\n",
    "\n",
    "    # measurement\n",
    "    if measure == 'Z':\n",
    "        circuit.measure(q[0], c[0])\n",
    "    elif measure == 'X':\n",
    "        circuit.u2(0, np.pi, q[0])\n",
    "        circuit.measure(q[0], c[0])\n",
    "    elif measure == 'Y':\n",
    "        circuit.u2(0, np.pi/2, q[0])\n",
    "        circuit.measure(q[0], c[0])\n",
    "    else:\n",
    "        raise ValueError('Not valid input for measurement: input should be \"X\" or \"Y\" or \"Z\"')\n",
    "\n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the quantum module part. The ```quantum_module``` finds the expectation values of a Pauli operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quantum_module(parameters, measure):\n",
    "    # measure\n",
    "    if measure == 'I':\n",
    "        return 1\n",
    "    elif measure == 'Z':\n",
    "        circuit = vqe_circuit(parameters, 'Z')\n",
    "    elif measure == 'X':\n",
    "        circuit = vqe_circuit(parameters, 'X')\n",
    "    elif measure == 'Y':\n",
    "        circuit = vqe_circuit(parameters, 'Y')\n",
    "    else:\n",
    "        raise ValueError('Not valid input for measurement: input should be \"I\" or \"X\" or \"Z\" or \"Y\"')\n",
    "    \n",
    "    shots = 8192\n",
    "    backend = BasicAer.get_backend('qasm_simulator')\n",
    "    job = execute(circuit, backend, shots=shots)\n",
    "    result = job.result()\n",
    "    counts = result.get_counts()\n",
    "    \n",
    "    # expectation value estimation from counts\n",
    "    expectation_value = 0\n",
    "    for measure_result in counts:\n",
    "        sign = +1\n",
    "        if measure_result == '1':\n",
    "            sign = -1\n",
    "        expectation_value += sign * counts[measure_result] / shots\n",
    "        \n",
    "    return expectation_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ```pauli_operator_to_dict``` creates a  dictionary from the ```WeightedPauliOperator``` object (the Hamiltonian is encoded in this object). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pauli_operator_to_dict(pauli_operator):\n",
    "    \"\"\"\n",
    "    from WeightedPauliOperator return a dict:\n",
    "    {I: 0.7, X: 0.6, Z: 0.1, Y: 0.5}.\n",
    "    :param pauli_operator: qiskit's WeightedPauliOperator\n",
    "    :return: a dict in the desired form.\n",
    "    \"\"\"\n",
    "    d = pauli_operator.to_dict()\n",
    "    paulis = d['paulis']\n",
    "    paulis_dict = {}\n",
    "\n",
    "    for x in paulis:\n",
    "        label = x['label']\n",
    "        coeff = x['coeff']['real']\n",
    "        paulis_dict[label] = coeff\n",
    "\n",
    "    return paulis_dict\n",
    "pauli_dict = pauli_operator_to_dict(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is the main method that takes parameters for ansatz state preparation and returns the corresponding expectation value of the Hamiltonian. For each Pauli term, we create separate quantum modules that calculate the expectation value of each Pauli. Then, all expectation values of Pauli operators multiplied by there corresponding coefficients ($a$, $b$, $c$, $d$) are summed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vqe(parameters):\n",
    "        \n",
    "    # quantum_modules\n",
    "    quantum_module_I = pauli_dict['I'] * quantum_module(parameters, 'I')\n",
    "    quantum_module_Z = pauli_dict['Z'] * quantum_module(parameters, 'Z')\n",
    "    quantum_module_X = pauli_dict['X'] * quantum_module(parameters, 'X')\n",
    "    quantum_module_Y = pauli_dict['Y'] * quantum_module(parameters, 'Y')\n",
    "    \n",
    "    # summing the measurement results\n",
    "    classical_adder = quantum_module_I + quantum_module_Z + quantum_module_X + quantum_module_Y\n",
    "    \n",
    "    return classical_adder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the final part of the code where we are minimizing (optimizing) the returned value from the ```vqe``` method by changing parameters for the quantum state preparation circuit (trial/ansatz $\\left| \\psi \\right\\rangle$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The exact ground state energy is: -5.775991418998328\n",
      "The estimated ground state energy from VQE algorithm is: -5.771199193113992\n"
     ]
    }
   ],
   "source": [
    "parameters_array = np.array([np.pi, np.pi])\n",
    "tol = 1e-3 # tolerance for optimization precision.\n",
    "\n",
    "vqe_result = minimize(vqe, parameters_array, method=\"Powell\", tol=tol)\n",
    "print('The exact ground state energy is: {}'.format(reference_energy))\n",
    "print('The estimated ground state energy from VQE algorithm is: {}'.format(vqe_result.fun))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is it! VQE algorithm for one qubit. But what should be done if there is a Hamiltonian for multiple qubits? For example, the Hamiltonian for the two-qubit system could look like this:\n",
    "\n",
    "$$\\left\\langle H \\right\\rangle = \\left\\langle \\psi \\right| H \\left| \\psi \\right\\rangle = 0.4 \\cdot \\left\\langle \\psi \\right| IX \\left| \\psi \\right\\rangle + 0.6 \\cdot \\left\\langle \\psi \\right| IZ \\left| \\psi \\right\\rangle + 0.8 \\cdot \\left\\langle \\psi \\left| XY \\right| \\psi \\right\\rangle.$$\n",
    "\n",
    "This time the algorithm should calculate the expectation value of each Pauli product. The most not trivial case is the $XY$ Pauli term and we will consider only this case. Note that for tensor products of operators like $XY$ Pauli term (actually $XY = X \\otimes Y$):\n",
    "\n",
    "$$U_1 \\otimes U_2 \\left| \\psi_1 \\right\\rangle \\otimes \\left| \\psi_2 \\right\\rangle =  U_1 \\left| \\psi_1 \\right\\rangle \\otimes U_2 \\left| \\psi_2 \\right\\rangle,$$\n",
    "\n",
    "were $U_1$ and $U_2$ are some unitary operators acting on two qubits with separate $\\left| \\psi_1 \\right\\rangle$ and $\\left| \\psi_2 \\right\\rangle$ wavefunctions. This way it can be shown that tensor products of the eigenvectors of $X$ and $Y$ are the eigenvectors of $XY$ Pauli term. It is true for any Pauli product term. Also, notice that all eigenvectors of  Pauli terms have eigenvalues equal to either $+1$ or $-1$. The eigenvectors of $XY$ that have $+1$ eigenvalue are $\\left| + \\right\\rangle \\otimes \\left| +i \\right\\rangle = \\left| ++i \\right\\rangle$, $\\left| - \\right\\rangle \\otimes \\left| -i \\right\\rangle = \\left| --i \\right\\rangle$, and the eigenvectors that have $-1$ eigenvalue are $\\left| + \\right\\rangle \\otimes \\left|-i \\right\\rangle = \\left| +-i \\right\\rangle$, $\\left| - \\right\\rangle \\otimes \\left| +i \\right\\rangle = \\left| -+i \\right\\rangle$. As an example:\n",
    "\n",
    "$$XY \\left| + -i \\right\\rangle= X \\otimes Y \\left| + \\right\\rangle \\otimes \\left| -i \\right\\rangle =  X \\left| + \\right\\rangle \\otimes Y \\left| -i \\right\\rangle = \n",
    "-\\left| + \\right\\rangle \\otimes \\left| -i \\right\\rangle = - \\left| + -i \\right\\rangle.$$\n",
    "\n",
    "The measurement logic stays the same. So, for $\\left\\langle XY \\right\\rangle$ one should apply $H_{gate}$ on the first qubit and $Y_{gate}$ on the second qubit (in tensor product notation $H_{gate} \\otimes Y_{gate}$) before $Z$ basis measurement. To see that we should represent the combined quantum state $\\left| \\psi \\right\\rangle = \\left| \\psi_1 \\right\\rangle \\otimes \\left| \\psi_2 \\right\\rangle$ in the basis of the eigenvectors of $XY$ operator:\n",
    "\n",
    "$$H_{gate} \\otimes Y_{gate} \\left| \\psi \\right\\rangle = H_{gate} \\otimes Y_{gate} \\left( c_1^{XY}\\left| ++i \\right\\rangle + c_2^{XY}\\left| +-i \\right\\rangle + c_3^{XY}\\left| -+i \\right\\rangle + c_4^{XY}\\left| --i \\right\\rangle \\right) \\\\\n",
    "= c_1^{XY}\\left| 00 \\right\\rangle + c_2^{XY}\\left| 01 \\right\\rangle + c_3^{XY}\\left| 10 \\right\\rangle + c_4^{XY}\\left| 11 \\right\\rangle,$$\n",
    "\n",
    "where $\\left| c \\right|^2$ are probabilities of obtaining the corresponding eigenvectors of $XY$ after the measurement. By doing similar calculations as we did for single Paulis, it can be shown that the expectation value $\\left\\langle \\psi \\right| XY \\left| \\psi \\right\\rangle$ is equal to the sum of $\\left| c \\right|^2$s with $-1$ coefficient for those eigenvectors which eigenvalue is $-1$ and with $+1$ coefficient for those eigenvectors which eigenvalue is $+1$.\n",
    "\n",
    "This way one can scale this solution for bigger Hamiltonians."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**[[Homepage][7]]**\n",
    "\n",
    "<h3 align=\"left\">\n",
    "References\n",
    "</h3>\n",
    "\n",
    "[[1]] [A. Peruzzo et al., Nature Communications, \"A variational eigenvalue solver on a photonic quantum processor\" (2014).][1]\n",
    "\n",
    "[[2]] [Michał Stęchły, \"Variational Quantum Eigensolver explained\"][2].\n",
    "\n",
    "[[3]] [M.A. Nielsen, I.L. Chuang, Cambridge University Press New York, \"Quantum Computation and Quantum Information: 10th Anniversary Edition\n",
    "10th\" (2011)][3]\n",
    "\n",
    "[[4]] [Learn Quantum Computation using Qiskit: Simulating Molecules using VQE][4]\n",
    "\n",
    "  [1]: https://www.nature.com/articles/ncomms5213?origin=ppub\n",
    "  [2]: https://www.mustythoughts.com/variational-quantum-eigensolver-explained\n",
    "  [3]: 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\n",
    "  [4]: https://qiskit.org/textbook/ch-applications/vqe-molecules.html\n",
    "  [5]: https://github.com/DavitKhach/quantum-algorithms-tutorials/blob/master/iterative_quantum_phase_estimation.ipynb\n",
    "  [6]: https://github.com/DavitKhach/quantum-algorithms-tutorials/blob/master/quantum_phase_estimation.ipynb\n",
    "  [7]: https://github.com/DavitKhach/quantum-algorithms-tutorials\n",
    "  [8]: https://nbviewer.jupyter.org/github/DavitKhach/quantum-algorithms-tutorials/blob/master/variational_quantum_eigensolver.ipynb"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}