VQE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variational Quantum Eigensolver (VQE) Example\n",
    "\n",
    "### Table of Contents\n",
    "1. [Introduction](#Introduction)\n",
    "2. [Qubits, gates, and all that](#Qubits,-gates,-and-all-that)\n",
    "3. [Make that Hamiltonian!](#Make-that-Hamiltonian!)\n",
    "4. [A first attempt at a quantum circuit](#A-first-attempt-at-a-quantum-circuit)\n",
    "5. [A \"real\" quantum circuit](#A-\"real\"-quantum-circuit)\n",
    "6. [A \"real\" measurement of the energy](#A-\"real\"-measurement-of-the-energy)\n",
    "7. [Appendix: All together now](#Appendix:-All-together-now)\n",
    "\n",
    "\n",
    "### Introduction\n",
    "\n",
    "The [variational quantum eigensolver (VQE)](https://en.wikipedia.org/wiki/Quantum_algorithm#Variational_quantum_eigensolver) is a hybrid classical-quantum algorithm that variationally determines the ground state energy of a Hamiltonian.\n",
    "\n",
    "It's quantum in the sense that the expectation value of the energy is computed via a quantum algorithm, but it is classical in the sense that the energy is minimized with a classical optimization algorithm.\n",
    "\n",
    "From a molecular electronic structure perspective, it is equivalent to computing the Full Configuration Interaction (FCI) for a given basis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quantum computing can be a little unintuitive, so it helps to see a working example. Although we won't actually do any \"real\" quantum computing, the methods can be understood from a linear algebra perspective. So that's what we are going to try and do for the VQE.\n",
    "\n",
    "The problem we want to tackle is to implement the VQE to compute the energy of molecular hydrogen (H$_2$) in a minimal basis. We will base this implementation off the really neat paper\n",
    "\n",
    "> O’Malley, Peter JJ, et al. \"[Scalable quantum simulation of molecular energies.](https://link.aps.org/doi/10.1103/PhysRevX.6.031007)\" Physical Review X 6.3 (2016): 031007.\n",
    "    \n",
    "In the O'Malley paper, they implement the VQE on a real quantum computer to compute the potential energy surface of H$_2$. The schematic we are going to follow can be seen below, and we are going to implement the \"software\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"circuit_for_h2.png\" width=\"800\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the big, overarching plan:\n",
    "\n",
    "1. Put the Hamiltonian in the computational (qubit) basis.\n",
    "2. Obtain a variational ansatz to parameterize the wave function.\n",
    "3. Represent this ansatz as a quantum circuit.\n",
    "4. Given this circuit, measure the expectation value of Hamiltonian (energy).\n",
    "5. Vary the circuit parameters until the energy is minimized.\n",
    "\n",
    "We'll look at these in turn, but first, let's look at single qubit states and some common matrices used to operate on them. This will set the groundwork for building up our VQE procedure.\n",
    "\n",
    "### Qubits, gates, and all that\n",
    "\n",
    "First, let's define some of the common quantum operators, gates, and states we want to work with. For starters, here are the identity and Pauli spin matrices:\n",
    "\n",
    "$$\n",
    "\\mathbf{I} = \\begin{pmatrix} \n",
    "1 & 0 \\\\\n",
    "0 & 1 \n",
    "\\end{pmatrix}\\quad\n",
    "{X} = \\begin{pmatrix} \n",
    "0 & 1 \\\\\n",
    "1 & 0 \n",
    "\\end{pmatrix}\\quad\n",
    "{Y} = \\begin{pmatrix} \n",
    "0 & -i \\\\\n",
    "i & 0 \n",
    "\\end{pmatrix}\\quad\n",
    "{Z} = \\begin{pmatrix} \n",
    "1 & 0 \\\\\n",
    "0 & -1 \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Two other important gates are the Hadamard matrix and phase matrix:\n",
    "\n",
    "$$\n",
    "{H} = \\frac{1}{\\sqrt{2}}\\begin{pmatrix} \n",
    "1 & 1 \\\\\n",
    "1 & -1 \n",
    "\\end{pmatrix}\\quad\n",
    "{S} = \\begin{pmatrix} \n",
    "1 & 0 \\\\\n",
    "0 & i \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Since for a single qubit the two basis states are\n",
    "$$\n",
    "|\\psi\\rangle = \\begin{pmatrix} \n",
    "|0\\rangle\\\\\n",
    "|1\\rangle \n",
    "\\end{pmatrix}, \\quad \\mathrm{meaning} \\qquad\n",
    "|0\\rangle = \\begin{pmatrix} \n",
    "1\\\\\n",
    "0 \n",
    "\\end{pmatrix}\\quad\n",
    "|1\\rangle = \\begin{pmatrix} \n",
    "0\\\\\n",
    "1 \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "we can define projection matrices, which are useful for defining, among other things, controlled gates (like CNOT, which we will see in a bit).\n",
    "\n",
    "$$\n",
    "|0\\rangle\\langle 0| = \\begin{pmatrix} \n",
    "1 & 0 \\\\\n",
    "0 & 0\n",
    "\\end{pmatrix}\\quad\n",
    "|1\\rangle\\langle 1| = \\begin{pmatrix} \n",
    "0 & 0\\\\\n",
    "0 & 1 \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Also useful for single qubits, are the rotation matrices:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "{R_x(\\theta)} &= \\begin{pmatrix} \n",
    "\\mathrm{cos}(\\theta/2) & -i\\cdot\\mathrm{sin}(\\theta/2) \\\\\n",
    "-i\\cdot\\mathrm{sin}(\\theta/2) & \\mathrm{cos}(\\theta/2) \n",
    "\\end{pmatrix}\\\\\n",
    "{R_y(\\theta)} &= \\begin{pmatrix} \n",
    "\\mathrm{cos}(\\theta/2) & -\\mathrm{sin}(\\theta/2) \\\\\n",
    "\\mathrm{sin}(\\theta/2) & \\mathrm{cos}(\\theta/2) \n",
    "\\end{pmatrix}\\\\\n",
    "{R_z(\\theta)} &= \\begin{pmatrix} \n",
    "\\mathrm{exp}(-i\\theta/2) & 0 \\\\\n",
    "0 & \\mathrm{exp}(i\\theta/2) \n",
    "\\end{pmatrix}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "\n",
    "So putting all these together for future use, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.set_printoptions(precision=4,suppress=True)\n",
    "\n",
    "\n",
    "# Pauli matrices\n",
    "I  = np.array([[ 1, 0],\n",
    "               [ 0, 1]])\n",
    "Sx = np.array([[ 0, 1],\n",
    "               [ 1, 0]])\n",
    "Sy = np.array([[ 0,-1j],\n",
    "               [1j, 0]])\n",
    "Sz = np.array([[ 1, 0],\n",
    "               [ 0,-1]])\n",
    "\n",
    "# Hadamard matrix\n",
    "H = (1/np.sqrt(2))*np.array([[ 1, 1],\n",
    "                             [ 1,-1]])\n",
    "\n",
    "# Phase matrix\n",
    "S = np.array([[ 1, 0],\n",
    "              [ 0,1j]])\n",
    "\n",
    "# single qubit basis states |0> and |1>\n",
    "q0 = np.array([[1],\n",
    "               [0]])\n",
    "q1 = np.array([[0],\n",
    "               [1]])\n",
    "\n",
    "# Projection matrices |0><0| and |1><1|\n",
    "P0  = np.dot(q0,q0.conj().T)\n",
    "P1  = np.dot(q1,q1.conj().T)\n",
    "\n",
    "\n",
    "# Rotation matrices as a function of theta, e.g. Rx(theta), etc.\n",
    "Rx = lambda theta : np.array([[    np.cos(theta/2),-1j*np.sin(theta/2)],\n",
    "                              [-1j*np.sin(theta/2),    np.cos(theta/2)]])\n",
    "Ry = lambda theta : np.array([[    np.cos(theta/2),   -np.sin(theta/2)],\n",
    "                              [    np.sin(theta/2),    np.cos(theta/2)]])\n",
    "Rz = lambda theta : np.array([[np.exp(-1j*theta/2),                0.0],\n",
    "                              [                0.0, np.exp(1j*theta/2)]])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, using these single qubit operations, we can build up some common two-qubit states and operations. By and large, moving from single to multiple qubits just involves tensor products of single qubit states and gates. Note that we number our qubits from bottom to top in this example, **so qubit-1 is on top and qubit-0 is on bottom**; e.g. in the figure:\n",
    "\n",
    "<img src=\"qubit_order.png\" width=\"100\">\n",
    "\n",
    "So for our two qubits the basis looks like\n",
    "\n",
    "$$\n",
    "|\\psi\\rangle = |\\hbox{qubit-1}\\rangle \\otimes |\\hbox{qubit-0}\\rangle = \n",
    "\\begin{pmatrix} \n",
    "|0\\rangle\\otimes|0\\rangle\\\\\n",
    "|0\\rangle\\otimes|1\\rangle\\\\\n",
    "|1\\rangle\\otimes|0\\rangle\\\\\n",
    "|1\\rangle\\otimes|1\\rangle\\\\\n",
    "\\end{pmatrix} = \n",
    "\\begin{pmatrix} \n",
    "|00\\rangle\\\\\n",
    "|01\\rangle\\\\\n",
    "|10\\rangle\\\\\n",
    "|11\\rangle\\\\\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "so that, e.g.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "|00\\rangle &= \\begin{pmatrix} \n",
    "1\\\\\n",
    "0\\\\\n",
    "\\end{pmatrix} \\otimes \\begin{pmatrix} \n",
    "1\\\\\n",
    "0\\\\\n",
    "\\end{pmatrix} = \\begin{pmatrix} \n",
    "1\\\\\n",
    "0\\\\\n",
    "0\\\\\n",
    "0\\\\\n",
    "\\end{pmatrix},\\\\\n",
    "|01\\rangle &= \n",
    "\\begin{pmatrix} \n",
    "1\\\\\n",
    "0\\\\\n",
    "\\end{pmatrix} \\otimes \\begin{pmatrix} \n",
    "0\\\\\n",
    "1\\\\\n",
    "\\end{pmatrix} = \\begin{pmatrix} \n",
    "0\\\\\n",
    "1\\\\\n",
    "0\\\\\n",
    "0\\\\\n",
    "\\end{pmatrix}, \\mathrm{~etc.}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "For two qubits, useful are the CNOT$_{01}$ and CNOT$_{10}$, where the first digit is the \"control\" qubit and the second is the \"target\" qubit. (Again, recall that in the O'Malley paper, we define qubit-1 as the \"top\" qubit, and qubit-0 as the \"bottom\" qubit, so if you change the qubit numbering the definitions will be swapped).\n",
    "\n",
    "They can be defined as the sum of two tensor products:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathrm{CNOT}_{10} &= (|0\\rangle\\langle 0| \\otimes \\mathbf{I}) + (|1\\rangle\\langle 1| \\otimes {X}) = \\begin{pmatrix} \n",
    "1 & 0 & 0 & 0\\\\\n",
    "0 & 1 & 0 & 0\\\\\n",
    "0 & 0 & 0 & 1\\\\\n",
    "0 & 0 & 1 & 0\n",
    "\\end{pmatrix}\\\\\n",
    "\\mathrm{CNOT}_{01} &= (\\mathbf{I} \\otimes |0\\rangle\\langle 0|) + (X \\otimes |1\\rangle\\langle 1|) = \\begin{pmatrix} \n",
    "1 & 0 & 0 & 0\\\\\n",
    "0 & 0 & 0 & 1\\\\\n",
    "0 & 0 & 1 & 0\\\\\n",
    "0 & 1 & 0 & 0\n",
    "\\end{pmatrix}\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The SWAP gate does what you'd expect, swapping qubit 1 and qubit 0:\n",
    "\n",
    "$$\n",
    "\\mathrm{SWAP} = \\begin{pmatrix} \n",
    "1 & 0 & 0 & 0\\\\\n",
    "0 & 0 & 1 & 0\\\\\n",
    "0 & 1 & 0 & 0\\\\\n",
    "0 & 0 & 0 & 1\n",
    "\\end{pmatrix}\\\\\n",
    "$$\n",
    "\n",
    "Most other operations are just simple tensor products of the single qubit operations. Like if you want to apply a Pauli X gate to qubit 1 and a Pauli Y gate to qubit 0, it's just $X_1 \\otimes Y_0$:\n",
    "\n",
    "$$\n",
    "{X_1 \\otimes Y_0} = \n",
    "\\begin{pmatrix} \n",
    "0 & 1 \\\\\n",
    "1 & 0 \\\\\n",
    "\\end{pmatrix}\\otimes\n",
    "\\begin{pmatrix} \n",
    "0 & -i \\\\\n",
    "i & 0 \\\\\n",
    "\\end{pmatrix} =\n",
    "\\begin{pmatrix} \n",
    "0 & 0 & 0 & -i\\\\\n",
    "0 & 0 & i & 0\\\\\n",
    "0 & -i & 0 & 0\\\\\n",
    "i & 0 & 0 & 0\\\\\n",
    "\\end{pmatrix}\\\\\n",
    "$$\n",
    "\n",
    "Okay, let's get some of these coded up for later, and if we need any products, we know we can just `np.kron` them later as needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import block_diag\n",
    "\n",
    "# CNOTij, where i is control qubit and j is target qubit\n",
    "CNOT10 = np.kron(P0,I) + np.kron(P1,Sx) # control -> q1, target -> q0\n",
    "CNOT01 = np.kron(I,P0) + np.kron(Sx,P1) # control -> q0, target -> q1\n",
    "\n",
    "SWAP   = block_diag(1,Sx,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make that Hamiltonian!\n",
    "\n",
    "Now we have our building blocks, we need to think about how to represent the Hamiltonian, which is usually in the fermion basis, in the basis of qubits. There are a few ways to do this, but the most common are the [Jordan-Wigner (JW)](https://en.wikipedia.org/wiki/Jordan%E2%80%93Wigner_transformation) transformation and the [Bravyi-Kitaev (BK)](https://arxiv.org/abs/1208.5986) transformation.\n",
    "\n",
    "According to the O'Malley paper, using the BK-transformed Hamiltonian and exploiting some of the symmetry in the H$_2$ molecule, the Hamiltonian can be represented by only two qubits. (Nice.) The Hamiltonian has five components:\n",
    "\n",
    "$$\n",
    "\\hat{H}_{\\mathrm{BK}} = g_0 \\mathbf{I} + g_1 Z_0 + g_2 Z_1 + g_3 Z_0Z_1 + g_4 Y_0Y_1 + g_5 X_0 X_1\n",
    "$$\n",
    "\n",
    "It general, the coefficients $g_i$ can be obtained from a cheap [Hartree-Fock](https://en.wikipedia.org/wiki/Hartree%E2%80%93Fock_method) calculation.\n",
    "\n",
    "So let's build a Hamiltonian for the H$_2$ molecule with a bond length of 0.75 A. We will also also add the nuclear repulsion energy `nuclear_repulsion`. The parameters we use can be found in Table 1 of the Appendix in the O'Malley [paper](https://journals.aps.org/prx/pdf/10.1103/PhysRevX.6.031007), though for the nuclear repulsion energy you have to calculate yourself (I used Gaussian). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# See DOI: 10.1103/PhysRevX.6.031007\n",
    "# Here, we use parameters given for H2 at R=0.75A\n",
    "g0 = -0.4804\n",
    "g1 = +0.3435\n",
    "g2 = -0.4347\n",
    "g3 = +0.5716\n",
    "g4 = +0.0910\n",
    "g5 = +0.0910\n",
    "\n",
    "nuclear_repulsion = 0.7055696146"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With all this, we can build the Hamiltonian, `Hmol`, explicitly in matrix form by taking tensor products, also known as Kronecker products (using `np.kron`), of the single-qubit matrix operators we built previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Hmol = (g0 * np.kron( I, I) + # g0 * I\n",
    "        g1 * np.kron( I,Sz) + # g1 * Z0\n",
    "        g2 * np.kron(Sz, I) + # g2 * Z1\n",
    "        g3 * np.kron(Sz,Sz) + # g3 * Z0Z1\n",
    "        g4 * np.kron(Sy,Sy) + # g4 * Y0Y1\n",
    "        g5 * np.kron(Sx,Sx))  # g5 * X0X1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's take a look at the Hamiltonian matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.0000+0.j  0.0000+0.j  0.0000+0.j  0.0000+0.j]\n",
      " [ 0.0000+0.j -1.8302+0.j  0.1820+0.j  0.0000+0.j]\n",
      " [ 0.0000+0.j  0.1820+0.j -0.2738+0.j  0.0000+0.j]\n",
      " [ 0.0000+0.j  0.0000+0.j  0.0000+0.j  0.1824+0.j]]\n"
     ]
    }
   ],
   "source": [
    "print(Hmol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have the Hamiltonian in the computational basis, let's just diagonalize it to get the energy (lowest eigenvalue). By adding the nuclear repulsion energy to the result we should get the same result as a Full Configuration Interaction (FCI) calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Classical diagonalization: -1.1456295 Eh\n",
      "Exact (from G16):          -1.1457417 Eh\n"
     ]
    }
   ],
   "source": [
    "electronic_energy = np.linalg.eigvalsh(Hmol)[0] # take the lowest value\n",
    "print(\"Classical diagonalization: {:+2.8} Eh\".format(electronic_energy + nuclear_repulsion))\n",
    "print(\"Exact (from G16):          {:+2.8} Eh\".format(-1.1457416808))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considering that the Hamiltonian elements had a precision of 1E-04, this is very good agreement. However, this approach utilizes a classical algorithm to obtain the eigenvalues. We want to see if we can obtain the eigenvalues using a quantum circuit.\n",
    "\n",
    "\n",
    "### A first attempt at a quantum circuit\n",
    "\n",
    "The usual input for quantum algorithms is to start in the $|00\\cdots\\rangle$ state. This is represented by a zero vector, with the first element set to 1. Because our Hamiltonian for H$_2$ only requires two qubits, we will start with the state $|01\\rangle$. To obtain this from $|00\\rangle$, we just need to act on the zeroth qubit with the Pauli X operator. This is the first step in the quantum circuit in the figure I showed above from the O'Malley paper. (They apply X$_{\\pi}$ -- same thing.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.]\n",
      " [ 1.]\n",
      " [ 0.]\n",
      " [ 0.]]\n"
     ]
    }
   ],
   "source": [
    "# initial basis, put in |01> state with Sx operator on q0\n",
    "psi0 = np.zeros((4,1))\n",
    "psi0[0] = 1\n",
    "psi0 = np.dot(np.kron(I,Sx),psi0)\n",
    "print(psi0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We haven't defined our VQE ansatz yet, but before we do, let's write a function to return the expected value of the Hamiltonian `Hmol` given an `ansatz`, its parameter `theta`, and the initial state `psi0`. This ansatz will eventually be encoded by the quantum circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def expected(theta,ansatz,Hmol,psi0):\n",
    "    circuit = ansatz(theta[0])\n",
    "    psi = np.dot(circuit,psi0)\n",
    "    return np.real(np.dot(psi.conj().T,np.dot(Hmol,psi)))[0,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the expectation value in had, we now define an ansatz. In the O'Malley paper, they utilize the Unitary Coupled Cluster (UCC) ansatz, which in this case depends only on a single parameter $\\theta$:\n",
    "\n",
    "$U(\\theta) = \\mathrm{exp}\\left(-i\\theta X_0Y_1\\right)$\n",
    "\n",
    "so that a parameterized wave function $|\\psi(\\theta)\\rangle$ for the ground state of H$_2$ is given as\n",
    "\n",
    "$|\\psi(\\theta)\\rangle = \\mathrm{exp}\\left(-i\\theta X_0Y_1\\right)|01\\rangle$\n",
    "\n",
    "and $X_0Y_1$ is the tensor product of the Pauli-X on qubit 0 and Pauli-Y on qubit 1.\n",
    "\n",
    "Before thinking about how we might represent $U(\\theta)$ as a series of quantum gates, let's plug it into the expression\n",
    "\n",
    "$E(\\theta) = \\frac{\\langle \\psi | U^{\\dagger}(\\theta)\\hat{H}_{\\mathrm{mol}}U(\\theta)|\\psi\\rangle}{\\langle \\psi | U^{\\dagger}(\\theta)U(\\theta)|\\psi\\rangle} = \\langle \\psi | U^{\\dagger}(\\theta)\\hat{H}_{\\mathrm{mol}}U(\\theta)|\\psi\\rangle$\n",
    "\n",
    "(Note that as long as $\\psi$ is normalized and $U$ is unitary, we can ignore the normalization $\\langle \\psi | U^{\\dagger}(\\theta)U(\\theta)|\\psi\\rangle$, since it always equals 1.)\n",
    "\n",
    "Given the `ansatz` and the initial state `psi0`, we can minimize `expected()` using the classical optimizers in the `scipy` package. So straightforwardly plugging in and minimizing yields the lazy result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lazy VQE: \n",
      "  [+] theta:  -0.11487186 deg\n",
      "  [+] energy: -1.1456295 Eh\n"
     ]
    }
   ],
   "source": [
    "from scipy.linalg import expm\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "# our UCC ansatz, not yet represented in terms of quantum gates\n",
    "ansatz = lambda theta: expm(-1j*np.array([theta])*np.kron(Sy,Sx))\n",
    "\n",
    "# initial guess for theta\n",
    "theta  = [0.0]\n",
    "result = minimize(expected,theta,args=(ansatz,Hmol,psi0))\n",
    "theta  = result.x[0]\n",
    "val    = result.fun\n",
    "\n",
    "print(\"Lazy VQE: \")\n",
    "print(\"  [+] theta:  {:+2.8} deg\".format(theta))\n",
    "print(\"  [+] energy: {:+2.8} Eh\".format(val + nuclear_repulsion))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which equals the result we got from diagonalization of the Hamiltonian. So we know that the UCC ansatz works! But we were lazy and didn't bother to think about how the quantum computer would compute the exponential. So, how do we represent the UCC ansatz in terms of quantum gates acting on the initial qubit state?\n",
    "\n",
    "### A \"real\" quantum circuit\n",
    "\n",
    "So before we saw the UCC ansatz works, but we cheated by keeping it as a matrix exponential. This is not a suitable form for a quantum computer. Let's do better.\n",
    "\n",
    "According to the O'Malley paper , we can represent $U(\\theta)$ for this problem as:\n",
    "\n",
    "<img src=\"circuit_close.png\" width=\"500\">\n",
    "\n",
    "This means that we should change our ansatz to read"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read right-to-left (bottom-to-top?)\n",
    "\n",
    "ansatz = lambda theta: (np.dot(np.dot(np.kron(-Ry(np.pi/2),Rx(np.pi/2)),\n",
    "                        np.dot(CNOT10, \n",
    "                        np.dot(np.kron(I,Rz(theta)),\n",
    "                               CNOT10))),\n",
    "                               np.kron(Ry(np.pi/2),-Rx(np.pi/2))))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that while you read the the circuit diagram left-to-right, when you read the matrix expression above it is better read right-to-left. The right-most matrices are applied to the state first, so that the first gates we apply are -R$_x(\\pi/2)$ to qubit-0 and R$_y(\\pi/2)$ to qubit-1. Also note that when we apply these two gates, it is simultaneous so the \"total\" gate is really -R$_y(\\pi/2) \\otimes $R$_x(\\pi/2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "VQE: \n",
      "  [+] theta:  +2.9118489 deg\n",
      "  [+] energy: -1.1456295 Eh\n"
     ]
    }
   ],
   "source": [
    "theta  = [0.0]\n",
    "result = minimize(expected,theta,args=(ansatz,Hmol,psi0))\n",
    "theta  = result.x[0]\n",
    "val    = result.fun\n",
    "\n",
    "print(\"VQE: \")\n",
    "print(\"  [+] theta:  {:+2.8} deg\".format(theta))\n",
    "print(\"  [+] energy: {:+2.8} Eh\".format(val + nuclear_repulsion))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which is the correct answer! Since we can now compute the expectation value of our Hamiltonian using quantum gates, we can pass the computed energy to a classical optimizer, which gives new parameters for the quantum gates. When this process is repeated until convergence, we obtain the FCI ground state energy. Also, once we have the optimized wave function parameters, the ground state can be easily reconstructed for additional simulations, etc.\n",
    "\n",
    "You might have noticed, though, that the above is still not sufficient for a quantum computer. The reason is that although we have represented our wave function with quantum gates, the *measurement* of the expectation value is still poorly defined as a physical operation. Even if you have prepared your qubits to represent a molecular wave function, measuring the expectation value of the Hamiltonian is not simply accomplished physically by applying a \"Hamiltonian operation\". \n",
    "\n",
    "An analogy: similar to classical computation, you might want a string, or float, or whatever as the \"true\" output of your function, but for the computer to compute it for you -- it needs to ultimately be in binary. Same thing for the quantum computer. Our function should ultimately return the energy, but it needs to process this in terms of quantum bits. \n",
    "\n",
    "\n",
    "### A \"real\" measurement of the energy\n",
    "\n",
    "All that is to say that we were cheating again. Experimentally, the \"only\" measurements we can make are those which probe the final quantum state of the qubits. What we need a way to connect measurements of qubits to an expression for the expectation value of the molecular electronic Hamiltonian.\n",
    "\n",
    "Put another way, the problem stems back to our defintion of the expected value:\n",
    "\n",
    "```python\n",
    "def expected(theta,ansatz,Hmol,psi0):\n",
    "    circuit = ansatz(theta[0])\n",
    "    psi = np.dot(circuit,psi0)\n",
    "    return np.real(np.dot(psi.conj().T,np.dot(Hmol,psi)))[0,0]\n",
    "```\n",
    "\n",
    "Simply dotting in `Hmol` with the wave function will not work, because physically we don't have a measuring apparatus for \"energy\". We can, however, measure the state of each qubit by measuring the spin ($\\hat{\\mathrm{S}}_z$) of each qubit. We need to reduce the Hamiltonian's expected value into these types of \"easy\" projective measurements that can be done in the computational basis. These are sums of Pauli measurements.\n",
    "\n",
    "Now in some respects, we are already halfway there. For a normalized wave function $|\\psi'\\rangle$:\n",
    "\n",
    "$$E = \\langle \\psi'|\\hat{H}_\\mathrm{mol}|\\psi'\\rangle$$\n",
    "\n",
    "and using the definition of our H$_2$ Hamiltonian in the computational basis we have:\n",
    "\n",
    "$$\n",
    "\\begin{align} E &= \\langle \\psi'|g_0 \\mathbf{I} + g_1 Z_0 + g_2 Z_1 + g_3 Z_0Z_1 + g_4 Y_0Y_1 + g_5 X_0 X_1|\\psi'\\rangle \\\\\n",
    "&= g_0\\langle \\mathbf{I} \\rangle + g_1\\langle Z_0 \\rangle + g_2 \\langle Z_1 \\rangle + g_3 \\langle Z_0Z_1 \\rangle + g_4 \\langle Y_0Y_1 \\rangle + g_5\\langle X_0 X_1 \\rangle \\\\\n",
    "&= \\sum_i g_i \\langle \\hat{O}_i \\rangle\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "meaning that, given our wave function in the computational basis, if we can compute the expected value of the (products of) Pauli operators, we can relate this to the expected value of the Hamiltonian through the sum given above. This is given in that figure:\n",
    "\n",
    "<img src=\"expectation.png\" width=\"150\">\n",
    "\n",
    "Let's go a step further, though. It would be even better if we could relate the above expression to a *single type* of Pauli measurement, that is, measuring the spin of just *one* qubit. Then we don't need to have multiple measurement apparatus.\n",
    "\n",
    "Thankfully, there is a way to do this. The trick is to apply an additional unitary transformation at the end of the circuit so that, by measuring the spin of the top qubit, we can obtain any Pauli measurement. In our case, that means relating each of the the $\\langle \\hat{O}_i\\rangle$ quantities to the expected value of $\\langle Z_1 \\otimes \\mathbf{I}\\rangle$ by some appropriate unitary operator. This is what is meant by the R$_t$ gates in this part of the figure. The R$_t$ are the unitaries we are talking about. \n",
    "\n",
    "<img src=\"measurement.png\" width=\"150\">\n",
    "\n",
    "The little measuring gauge means we finally apply the measurement. Because we apply a particular form of the unitaries, we only need to measure the state of qubit-1 like we talked about above, but that's not necessarily the only way to go about it.\n",
    "\n",
    "You can find a table of some of these transformations [here](https://docs.microsoft.com/en-us/quantum/concepts/pauli-measurements), but here are a few examples:\n",
    "\n",
    "For example, the simplest case is if you want to measure $Z_1 \\otimes \\mathbf{I}$. Then you don't have to do anything:\n",
    "\n",
    " $\\qquad Z_1 \\otimes \\mathbf{I} = (\\mathbf{I} \\otimes \\mathbf{I})^{\\dagger} \\otimes (Z_1 \\otimes \\mathbf{I}) \\otimes (\\mathbf{I} \\otimes \\mathbf{I})$\n",
    "\n",
    "But if you want to measure $Y_1 \\otimes \\mathbf{I}$, then\n",
    "\n",
    "$\\qquad Y_1 \\otimes \\mathbf{I} = ({HS^{\\dagger}} \\otimes \\mathbf{I})^{\\dagger} \\otimes (Z_1 \\otimes \\mathbf{I}) \\otimes ({HS^{\\dagger}} \\otimes \\mathbf{I})$\n",
    "\n",
    "For $\\mathbf{I} \\otimes Z_0$, we have to apply the SWAP gate\n",
    "\n",
    "$\\qquad \\mathbf{I} \\otimes Z_0 = (\\mathrm{SWAP})^{\\dagger} \\otimes (Z_1 \\otimes \\mathbf{I}) \\otimes (\\mathrm{SWAP})$\n",
    "\n",
    "And as a final example, for $X_1Z_0$, (*caveat: we have different qubit ordering compared to the Microsoft documentation, so our CNOT$_{10}$ vs CNOT$_{01}$ definitions are swapped*)\n",
    "\n",
    "$\\qquad X_1 \\otimes Z_0 = (\\mathrm{CNOT}_{01}(H \\otimes \\mathbf{I}))^{\\dagger} \\otimes (Z_1 \\otimes \\mathbf{I}) \\otimes (\\mathrm{CNOT}_{01}(H \\otimes \\mathbf{I}))$\n",
    "\n",
    "Since you can think of these unitary transformations either acting on the operator, *or* acting on the state, the end result is that by applying the particular transformation and *then* measuring $Z_1$ you can get any Pauli measurement you want.\n",
    "\n",
    "It might be easier to see by example. Let's see how this plays out for our Hamiltonian. \n",
    "\n",
    "We will update the function `expected()` to `projective_expected()`, and remove the `Hmol` argument. In our case, this function will include the hard-coded representation of the Hamiltonian, in terms of measuring $Z_1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def projective_expected(theta,ansatz,psi0):\n",
    "    # this will depend on the hard-coded Hamiltonian + coefficients\n",
    "    circuit = ansatz(theta[0])\n",
    "    psi = np.dot(circuit,psi0)\n",
    "    \n",
    "    # for 2 qubits, assume we can only take Pauli Sz measurements (Sz \\otimes I)\n",
    "    # we just apply the right unitary for the desired Pauli measurement\n",
    "    measureZ = lambda U: np.dot(np.conj(U).T,np.dot(np.kron(Sz,I),U))\n",
    "    \n",
    "    energy = 0.0\n",
    "    \n",
    "    # although the paper indexes the hamiltonian left-to-right (0-to-1) \n",
    "    # qubit-1 is always the top qubit for us, so the tensor pdt changes\n",
    "    # e.g. compare with the \"exact Hamiltonian\" we explicitly diagonalized\n",
    "    \n",
    "    # <I1 I0>\n",
    "    energy += g0 # it is a constant\n",
    "\n",
    "    # <I1 Sz0>\n",
    "    U = SWAP\n",
    "    energy += g1*np.dot(psi.conj().T,np.dot(measureZ(U),psi))\n",
    "\n",
    "    # <Sz1 I0>\n",
    "    U = np.kron(I,I)\n",
    "    energy += g2*np.dot(psi.conj().T,np.dot(measureZ(U),psi))\n",
    "\n",
    "    # <Sz1 Sz0>\n",
    "    U = CNOT01\n",
    "    energy += g3*np.dot(psi.conj().T,np.dot(measureZ(U),psi))\n",
    "\n",
    "    # <Sx1 Sx0>\n",
    "    U = np.dot(CNOT01,np.kron(H,H))\n",
    "    energy += g4*np.dot(psi.conj().T,np.dot(measureZ(U),psi))\n",
    "\n",
    "    # <Sy1 Sy0>\n",
    "    U = np.dot(CNOT01,np.kron(np.dot(H,S.conj().T),np.dot(H,S.conj().T)))\n",
    "    energy += g5*np.dot(psi.conj().T,np.dot(measureZ(U),psi))\n",
    "\n",
    "    return np.real(energy)[0,0]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the expectation value now defined in terms of measuring the spin of the zero-th qubit, let's carry out the VQE procedure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "VQE: \n",
      "  [+] theta:  +2.9118489 deg\n",
      "  [+] energy: -1.1456295 Eh\n"
     ]
    }
   ],
   "source": [
    "theta  = [0.0]\n",
    "result = minimize(projective_expected,theta,args=(ansatz,psi0))\n",
    "theta  = result.x[0]\n",
    "val    = result.fun\n",
    "\n",
    "print(\"VQE: \")\n",
    "print(\"  [+] theta:  {:+2.8} deg\".format(theta))\n",
    "print(\"  [+] energy: {:+2.8} Eh\".format(val + nuclear_repulsion))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Success! We get the same energy (and theta) as the previous VQE circuit, but now all measurements are related back to the result of measuring the spin of the qubit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Appendix: All together now\n",
    "\n",
    "Here's all of the pieces together in one place"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Classical diagonalization: -1.1456295 Eh\n",
      "Exact (from G16):          -1.1457417 Eh\n",
      "VQE: \n",
      "  [+] theta:  +2.9118489 deg\n",
      "  [+] energy: -1.1456295 Eh\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import block_diag\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "np.set_printoptions(precision=4,suppress=True)\n",
    "\n",
    "\n",
    "# Pauli matrices\n",
    "I  = np.array([[ 1, 0],\n",
    "               [ 0, 1]])\n",
    "Sx = np.array([[ 0, 1],\n",
    "               [ 1, 0]])\n",
    "Sy = np.array([[ 0,-1j],\n",
    "               [1j, 0]])\n",
    "Sz = np.array([[ 1, 0],\n",
    "               [ 0,-1]])\n",
    "\n",
    "# Hadamard matrix\n",
    "H = (1/np.sqrt(2))*np.array([[ 1, 1],\n",
    "                             [ 1,-1]])\n",
    "\n",
    "# Phase matrix\n",
    "S = np.array([[ 1, 0],\n",
    "              [ 0,1j]])\n",
    "\n",
    "# single qubit basis states |0> and |1>\n",
    "q0 = np.array([[1],\n",
    "               [0]])\n",
    "q1 = np.array([[0],\n",
    "               [1]])\n",
    "\n",
    "# Projection matrices |0><0| and |1><1|\n",
    "P0  = np.dot(q0,q0.conj().T)\n",
    "P1  = np.dot(q1,q1.conj().T)\n",
    "\n",
    "\n",
    "# Rotation matrices as a function of theta, e.g. Rx(theta), etc.\n",
    "Rx = lambda theta : np.array([[    np.cos(theta/2),-1j*np.sin(theta/2)],\n",
    "                              [-1j*np.sin(theta/2),    np.cos(theta/2)]])\n",
    "Ry = lambda theta : np.array([[    np.cos(theta/2),   -np.sin(theta/2)],\n",
    "                              [    np.sin(theta/2),    np.cos(theta/2)]])\n",
    "Rz = lambda theta : np.array([[np.exp(-1j*theta/2),                0.0],\n",
    "                              [                0.0, np.exp(1j*theta/2)]])\n",
    "\n",
    "# CNOTij, where i is control qubit and j is target qubit\n",
    "CNOT10 = np.kron(P0,I) + np.kron(P1,Sx) # control -> q1, target -> q0\n",
    "CNOT01 = np.kron(I,P0) + np.kron(Sx,P1) # control -> q0, target -> q1\n",
    "\n",
    "SWAP   = block_diag(1,Sx,1)\n",
    "\n",
    "# See DOI: 10.1103/PhysRevX.6.031007\n",
    "# Here, we use parameters given for H2 at R=0.75A\n",
    "g0 = -0.4804\n",
    "g1 = +0.3435\n",
    "g2 = -0.4347\n",
    "g3 = +0.5716\n",
    "g4 = +0.0910\n",
    "g5 = +0.0910\n",
    "\n",
    "nuclear_repulsion = 0.7055696146\n",
    "\n",
    "Hmol = (g0 * np.kron( I, I) + # g0 * I\n",
    "        g1 * np.kron( I,Sz) + # g1 * Z0\n",
    "        g2 * np.kron(Sz, I) + # g2 * Z1\n",
    "        g3 * np.kron(Sz,Sz) + # g3 * Z0Z1\n",
    "        g4 * np.kron(Sy,Sy) + # g4 * Y0Y1\n",
    "        g5 * np.kron(Sx,Sx))  # g5 * X0X1\n",
    "\n",
    "electronic_energy = np.linalg.eigvalsh(Hmol)[0] # take the lowest value\n",
    "print(\"Classical diagonalization: {:+2.8} Eh\".format(electronic_energy + nuclear_repulsion))\n",
    "print(\"Exact (from G16):          {:+2.8} Eh\".format(-1.1457416808))\n",
    "\n",
    "# initial basis, put in |01> state with Sx operator on q0\n",
    "psi0 = np.zeros((4,1))\n",
    "psi0[0] = 1\n",
    "psi0 = np.dot(np.kron(I,Sx),psi0)\n",
    "\n",
    "\n",
    "# read right-to-left (bottom-to-top?)\n",
    "ansatz = lambda theta: (np.dot(np.dot(np.kron(-Ry(np.pi/2),Rx(np.pi/2)),\n",
    "                        np.dot(CNOT10, \n",
    "                        np.dot(np.kron(I,Rz(theta)),\n",
    "                               CNOT10))),\n",
    "                               np.kron(Ry(np.pi/2),-Rx(np.pi/2))))\n",
    "\n",
    "def projective_expected(theta,ansatz,psi0):\n",
    "    # this will depend on the hard-coded Hamiltonian + coefficients\n",
    "    circuit = ansatz(theta[0])\n",
    "    psi = np.dot(circuit,psi0)\n",
    "    \n",
    "    # for 2 qubits, assume we can only take Pauli Sz measurements (Sz \\otimes I)\n",
    "    # we just apply the right unitary for the desired Pauli measurement\n",
    "    measureZ = lambda U: np.dot(np.conj(U).T,np.dot(np.kron(Sz,I),U))\n",
    "    \n",
    "    energy = 0.0\n",
    "    \n",
    "    # although the paper indexes the hamiltonian left-to-right (0-to-1) \n",
    "    # qubit-1 is always the top qubit for us, so the tensor pdt changes\n",
    "    # e.g. compare with the \"exact Hamiltonian\" we explicitly diagonalized\n",
    "    \n",
    "    # <I1 I0> \n",
    "    energy += g0 # it is a constant\n",
    "\n",
    "    # <I1 Sz0>\n",
    "    U = SWAP\n",
    "    energy += g1*np.dot(psi.conj().T,np.dot(measureZ(U),psi))\n",
    "\n",
    "    # <Sz1 I0>\n",
    "    U = np.kron(I,I)\n",
    "    energy += g2*np.dot(psi.conj().T,np.dot(measureZ(U),psi))\n",
    "\n",
    "    # <Sz1 Sz0>\n",
    "    U = CNOT01\n",
    "    energy += g3*np.dot(psi.conj().T,np.dot(measureZ(U),psi))\n",
    "\n",
    "    # <Sx1 Sx0>\n",
    "    U = np.dot(CNOT01,np.kron(H,H))\n",
    "    energy += g4*np.dot(psi.conj().T,np.dot(measureZ(U),psi))\n",
    "\n",
    "    # <Sy1 Sy0>\n",
    "    U = np.dot(CNOT01,np.kron(np.dot(H,S.conj().T),np.dot(H,S.conj().T)))\n",
    "    energy += g5*np.dot(psi.conj().T,np.dot(measureZ(U),psi))\n",
    "\n",
    "    return np.real(energy)[0,0]\n",
    "\n",
    "theta  = [0.0]\n",
    "result = minimize(projective_expected,theta,args=(ansatz,psi0))\n",
    "theta  = result.x[0]\n",
    "val    = result.fun\n",
    "\n",
    "# check it works...\n",
    "assert np.allclose(val + nuclear_repulsion,-1.1456295)\n",
    "\n",
    "print(\"VQE: \")\n",
    "print(\"  [+] theta:  {:+2.8} deg\".format(theta))\n",
    "print(\"  [+] energy: {:+2.8} Eh\".format(val + nuclear_repulsion))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}