ApplicationChallenge.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from qiskit import *\n",
    "%config InlineBackend.figure_format = 'svg'\n",
    "from scipy.optimize import minimize, minimize_scalar\n",
    "\n",
    "from qiskit.visualization import *\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "### Some global variables \n",
    "shots = 5000\n",
    "visualization = True "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Disclaimer: If you view this on the Github website, my \\ket & \\bra Latex macros don't seem to be rendering. To view properly, please run the Markdown cells locally.**\n",
    "\n",
    "<!--- Defining Latex macros -->\n",
    "$\\newcommand{\\ket}[1]{\\lvert #1 \\rangle}$\n",
    "$\\newcommand{\\bra}[1]{\\langle #1 \\rvert}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decomposing the Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "hamiltonian = [[0, 0, 0, 0], \n",
    "             [0, -1, 1, 0], \n",
    "             [0, 1, -1, 0], \n",
    "             [0, 0, 0, 0]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our Hamiltonian is a 2-qubit operator, so let's compute 2-tensored Pauli matrices to see what we can decompose it into. Using the Kronecker product, this isn't so bad. \n",
    "\n",
    "$$I \\otimes I = \\begin{bmatrix} 1 & 0 & 0 & 0\\\\ 0 & 1 & 0 & 0\\\\ 0 & 0 & 1 & 0\\\\ 0 & 0 & 0 & 1 \\end{bmatrix}$$\n",
    "\n",
    "$$ X \\otimes X = \\begin{bmatrix} 0 & 0 & 0 & 1\\\\ 0 & 0 & 1 & 0\\\\ 0 & 1 & 0 & 0\\\\ 1 & 0 & 0 & 0 \\end{bmatrix}$$\n",
    "\n",
    "$$ Y \\otimes Y = \\begin{bmatrix} 0 & 0 & 0 & -1\\\\ 0 & 0 & 1 & 0\\\\ 0 & 1 & 0 & 0\\\\ -1 & 0 & 0 & 0 \\end{bmatrix}$$\n",
    "\n",
    "$$ Z \\otimes Z = \\begin{bmatrix} 1 & 0 & 0 & 0\\\\ 0 & -1 & 0 & 0\\\\ 0 & 0 & -1 & 0\\\\ 0 & 0 & 0 & 1 \\end{bmatrix}$$\n",
    "\n",
    "It is clear from this that the 4 most basic 2-qubit operators all have 0's everywhere except for a main diagonal. This is useful because our Hamiltonian has only non-zero values in the 4 center spots. \n",
    "\n",
    "Notice the Hamiltonian is $$H = \\frac{1}{2} ( X \\otimes X + Y \\otimes Y + Z \\otimes Z - I \\otimes I )$$\n",
    "\n",
    "## Creating the ansatz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_ansatz(theta):\n",
    "    \"\"\"Returns QuantumCircuit, QuantumRegister, and ClassicalRegister components of ansatz \"\"\"\n",
    "    # Recommended ansatz: (IX) CX (RZ I) (HI) |00> \n",
    "    qc = QuantumCircuit()\n",
    "    qr = QuantumRegister(2,'qreg')\n",
    "    qc.add_register( qr )\n",
    "    qc.h( qr[0] )\n",
    "    qc.rz(theta, qr[0] )\n",
    "    qc.cx( qr[0], qr[1] )\n",
    "    qc.x( qr[1] )\n",
    "    cr = ClassicalRegister(2,'creg')\n",
    "    qc.add_register(cr)\n",
    "\n",
    "    #return qc.draw(output='mpl')\n",
    "    return qc, qr, cr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       "  \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Created with matplotlib (https://matplotlib.org/) -->\n",
       "<svg height=\"168.9334pt\" version=\"1.1\" viewBox=\"0 0 355.2825 168.9334\" width=\"355.2825pt\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       " <defs>\n",
       "  <style type=\"text/css\">\n",
       "*{stroke-linecap:butt;stroke-linejoin:round;}\n",
       "  </style>\n",
       " </defs>\n",
       " <g id=\"figure_1\">\n",
       "  <g id=\"patch_1\">\n",
       "   <path d=\"M 0 168.9334 \n",
       "L 355.2825 168.9334 \n",
       "L 355.2825 0 \n",
       "L 0 0 \n",
       "z\n",
       "\" style=\"fill:#ffffff;\"/>\n",
       "  </g>\n",
       "  <g id=\"axes_1\">\n",
       "   <g id=\"line2d_1\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 139.0079 43.5608 \n",
       "L 343.5374 43.5608 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-linecap:square;stroke-width:2;\"/>\n",
       "   </g>\n",
       "   <g id=\"line2d_2\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 139.0079 89.0118 \n",
       "L 343.5374 89.0118 \n",
       "\" style=\"fill:none;stroke:#000000;stroke-linecap:square;stroke-width:2;\"/>\n",
       "   </g>\n",
       "   <g id=\"line2d_3\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 139.0079 132.985642 \n",
       "L 343.5374 132.985642 \n",
       "\" style=\"fill:none;stroke:#778899;stroke-linecap:square;stroke-width:2;\"/>\n",
       "   </g>\n",
       "   <g id=\"line2d_4\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 139.0079 135.939957 \n",
       "L 343.5374 135.939957 \n",
       "\" style=\"fill:none;stroke:#778899;stroke-linecap:square;stroke-width:2;\"/>\n",
       "   </g>\n",
       "   <g id=\"line2d_5\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 252.6354 89.0118 \n",
       "L 252.6354 43.5608 \n",
       "\" style=\"fill:none;stroke:#6fa4ff;stroke-linecap:square;stroke-width:2;\"/>\n",
       "   </g>\n",
       "   <g id=\"line2d_6\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 166.2785 139.0079 \n",
       "L 170.8236 129.9177 \n",
       "\" style=\"fill:none;stroke:#778899;stroke-linecap:square;stroke-width:1.5;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_2\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 146.961825 58.332375 \n",
       "L 176.504975 58.332375 \n",
       "L 176.504975 28.789225 \n",
       "L 146.961825 28.789225 \n",
       "z\n",
       "\" style=\"fill:#6fa4ff;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_3\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 192.412825 58.332375 \n",
       "L 221.955975 58.332375 \n",
       "L 221.955975 28.789225 \n",
       "L 192.412825 28.789225 \n",
       "z\n",
       "\" style=\"fill:#bb8bff;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_4\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 252.6354 47.992272 \n",
       "C 253.81064 47.992272 254.937904 47.525345 255.768924 46.694324 \n",
       "C 256.599945 45.863304 257.066873 44.73604 257.066873 43.5608 \n",
       "C 257.066873 42.38556 256.599945 41.258296 255.768924 40.427276 \n",
       "C 254.937904 39.596255 253.81064 39.129327 252.6354 39.129327 \n",
       "C 251.46016 39.129327 250.332896 39.596255 249.501876 40.427276 \n",
       "C 248.670855 41.258296 248.203928 42.38556 248.203928 43.5608 \n",
       "C 248.203928 44.73604 248.670855 45.863304 249.501876 46.694324 \n",
       "C 250.332896 47.525345 251.46016 47.992272 252.6354 47.992272 \n",
       "z\n",
       "\" style=\"fill:#6fa4ff;stroke:#6fa4ff;stroke-linejoin:miter;stroke-width:1.5;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_5\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 252.6354 99.351902 \n",
       "C 255.377627 99.351902 258.007909 98.262404 259.946957 96.323357 \n",
       "C 261.886004 94.384309 262.975503 91.754027 262.975503 89.0118 \n",
       "C 262.975503 86.269573 261.886004 83.639291 259.946957 81.700243 \n",
       "C 258.007909 79.761196 255.377627 78.671697 252.6354 78.671697 \n",
       "C 249.893173 78.671697 247.262891 79.761196 245.323843 81.700243 \n",
       "C 243.384796 83.639291 242.295298 86.269573 242.295298 89.0118 \n",
       "C 242.295298 91.754027 243.384796 94.384309 245.323843 96.323357 \n",
       "C 247.262891 98.262404 249.893173 99.351902 252.6354 99.351902 \n",
       "z\n",
       "\" style=\"fill:#6fa4ff;stroke:#6fa4ff;stroke-linejoin:miter;stroke-width:2;\"/>\n",
       "   </g>\n",
       "   <g id=\"patch_6\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 283.314825 103.783375 \n",
       "L 312.857975 103.783375 \n",
       "L 312.857975 74.240225 \n",
       "L 283.314825 74.240225 \n",
       "z\n",
       "\" style=\"fill:#05bab6;\"/>\n",
       "   </g>\n",
       "   <g id=\"line2d_7\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 252.6354 94.92043 \n",
       "L 252.6354 83.10317 \n",
       "\" style=\"fill:none;stroke:#ffffff;stroke-linecap:square;stroke-width:2;\"/>\n",
       "   </g>\n",
       "   <g id=\"line2d_8\">\n",
       "    <path clip-path=\"url(#p29205ff115)\" d=\"M 246.72677 89.0118 \n",
       "L 258.54403 89.0118 \n",
       "\" style=\"fill:none;stroke:#ffffff;stroke-linecap:square;stroke-width:2;\"/>\n",
       "   </g>\n",
       "   <g id=\"text_1\">\n",
       "    <g clip-path=\"url(#p29205ff115)\">\n",
       "     <!-- $H$ -->\n",
       "     <defs>\n",
       "      <path d=\"M 16.890625 72.90625 \n",
       "L 26.8125 72.90625 \n",
       "L 21 43.015625 \n",
       "L 56.78125 43.015625 \n",
       "L 62.59375 72.90625 \n",
       "L 72.515625 72.90625 \n",
       "L 58.296875 0 \n",
       "L 48.390625 0 \n",
       "L 55.171875 34.71875 \n",
       "L 19.390625 34.71875 \n",
       "L 12.59375 0 \n",
       "L 2.6875 0 \n",
       "z\n",
       "\" id=\"DejaVuSans-Oblique-72\"/>\n",
       "     </defs>\n",
       "     <g transform=\"translate(156.7934 47.147987)scale(0.13 -0.13)\">\n",
       "      <use transform=\"translate(0 0.09375)\" xlink:href=\"#DejaVuSans-Oblique-72\"/>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "   <g id=\"text_2\">\n",
       "    <g clip-path=\"url(#p29205ff115)\">\n",
       "     <!-- $R_z$ -->\n",
       "     <defs>\n",
       "      <path d=\"M 25.203125 64.796875 \n",
       "L 20.21875 38.921875 \n",
       "L 32.90625 38.921875 \n",
       "Q 40.375 38.921875 44.984375 43.046875 \n",
       "Q 49.609375 47.171875 49.609375 53.8125 \n",
       "Q 49.609375 59.125 46.5 61.953125 \n",
       "Q 43.40625 64.796875 37.59375 64.796875 \n",
       "z\n",
       "M 43.3125 35.015625 \n",
       "Q 46.4375 34.28125 48.515625 31.390625 \n",
       "Q 50.59375 28.515625 53.328125 19.921875 \n",
       "L 59.515625 0 \n",
       "L 49.125 0 \n",
       "L 43.40625 18.703125 \n",
       "Q 41.21875 25.921875 38.328125 28.359375 \n",
       "Q 35.453125 30.8125 29.5 30.8125 \n",
       "L 18.609375 30.8125 \n",
       "L 12.59375 0 \n",
       "L 2.6875 0 \n",
       "L 16.890625 72.90625 \n",
       "L 39.109375 72.90625 \n",
       "Q 49.21875 72.90625 54.609375 68.328125 \n",
       "Q 60.015625 63.765625 60.015625 55.171875 \n",
       "Q 60.015625 47.5625 55.421875 41.984375 \n",
       "Q 50.828125 36.421875 43.3125 35.015625 \n",
       "z\n",
       "\" id=\"DejaVuSans-Oblique-82\"/>\n",
       "      <path d=\"M 11.625 54.6875 \n",
       "L 54.296875 54.6875 \n",
       "L 52.6875 46.484375 \n",
       "L 11.53125 7.171875 \n",
       "L 45.515625 7.171875 \n",
       "L 44.09375 0 \n",
       "L -0.296875 0 \n",
       "L 1.3125 8.203125 \n",
       "L 42.484375 47.515625 \n",
       "L 10.203125 47.515625 \n",
       "z\n",
       "\" id=\"DejaVuSans-Oblique-122\"/>\n",
       "     </defs>\n",
       "     <g transform=\"translate(200.0994 42.716515)scale(0.13 -0.13)\">\n",
       "      <use transform=\"translate(0 0.09375)\" xlink:href=\"#DejaVuSans-Oblique-82\"/>\n",
       "      <use transform=\"translate(69.482422 -16.3125)scale(0.7)\" xlink:href=\"#DejaVuSans-Oblique-122\"/>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "   <g id=\"text_3\">\n",
       "    <g clip-path=\"url(#p29205ff115)\">\n",
       "     <!-- $\\pi$/2 -->\n",
       "     <defs>\n",
       "      <path d=\"M 9.125 54.6875 \n",
       "L 61.53125 54.6875 \n",
       "L 59.765625 45.703125 \n",
       "L 52.875 45.703125 \n",
       "L 46.34375 12.109375 \n",
       "Q 45.65625 8.59375 46.578125 7.03125 \n",
       "Q 47.46875 5.515625 50.140625 5.515625 \n",
       "Q 50.875 5.515625 51.953125 5.671875 \n",
       "Q 53.078125 5.765625 53.421875 5.8125 \n",
       "L 52.15625 -0.6875 \n",
       "Q 50.34375 -1.3125 48.484375 -1.609375 \n",
       "Q 46.578125 -1.90625 44.78125 -1.90625 \n",
       "Q 38.921875 -1.90625 37.3125 1.265625 \n",
       "Q 35.6875 4.5 37.359375 13.09375 \n",
       "L 43.703125 45.703125 \n",
       "L 23.53125 45.703125 \n",
       "L 14.65625 0 \n",
       "L 5.46875 0 \n",
       "L 14.359375 45.703125 \n",
       "L 7.375 45.703125 \n",
       "z\n",
       "\" id=\"DejaVuSans-Oblique-960\"/>\n",
       "      <path d=\"M 25.390625 72.90625 \n",
       "L 33.6875 72.90625 \n",
       "L 8.296875 -9.28125 \n",
       "L 0 -9.28125 \n",
       "z\n",
       "\" id=\"DejaVuSans-47\"/>\n",
       "      <path d=\"M 19.1875 8.296875 \n",
       "L 53.609375 8.296875 \n",
       "L 53.609375 0 \n",
       "L 7.328125 0 \n",
       "L 7.328125 8.296875 \n",
       "Q 12.9375 14.109375 22.625 23.890625 \n",
       "Q 32.328125 33.6875 34.8125 36.53125 \n",
       "Q 39.546875 41.84375 41.421875 45.53125 \n",
       "Q 43.3125 49.21875 43.3125 52.78125 \n",
       "Q 43.3125 58.59375 39.234375 62.25 \n",
       "Q 35.15625 65.921875 28.609375 65.921875 \n",
       "Q 23.96875 65.921875 18.8125 64.3125 \n",
       "Q 13.671875 62.703125 7.8125 59.421875 \n",
       "L 7.8125 69.390625 \n",
       "Q 13.765625 71.78125 18.9375 73 \n",
       "Q 24.125 74.21875 28.421875 74.21875 \n",
       "Q 39.75 74.21875 46.484375 68.546875 \n",
       "Q 53.21875 62.890625 53.21875 53.421875 \n",
       "Q 53.21875 48.921875 51.53125 44.890625 \n",
       "Q 49.859375 40.875 45.40625 35.40625 \n",
       "Q 44.1875 33.984375 37.640625 27.21875 \n",
       "Q 31.109375 20.453125 19.1875 8.296875 \n",
       "z\n",
       "\" id=\"DejaVuSans-50\"/>\n",
       "     </defs>\n",
       "     <g transform=\"translate(201.0244 54.631245)scale(0.08 -0.08)\">\n",
       "      <use transform=\"translate(0 0.78125)\" xlink:href=\"#DejaVuSans-Oblique-960\"/>\n",
       "      <use transform=\"translate(60.205078 0.78125)\" xlink:href=\"#DejaVuSans-47\"/>\n",
       "      <use transform=\"translate(90.271484 0.78125)\" xlink:href=\"#DejaVuSans-50\"/>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "   <g id=\"text_4\">\n",
       "    <g clip-path=\"url(#p29205ff115)\">\n",
       "     <!-- $X$ -->\n",
       "     <defs>\n",
       "      <path d=\"M 13.71875 72.90625 \n",
       "L 23.6875 72.90625 \n",
       "L 36.1875 46.578125 \n",
       "L 58.796875 72.90625 \n",
       "L 70.3125 72.90625 \n",
       "L 40.28125 38.09375 \n",
       "L 58.40625 0 \n",
       "L 48.484375 0 \n",
       "L 33.796875 30.71875 \n",
       "L 7.171875 0 \n",
       "L -4.296875 0 \n",
       "L 29.78125 39.203125 \n",
       "z\n",
       "\" id=\"DejaVuSans-Oblique-88\"/>\n",
       "     </defs>\n",
       "     <g transform=\"translate(293.6014 92.598987)scale(0.13 -0.13)\">\n",
       "      <use transform=\"translate(0 0.09375)\" xlink:href=\"#DejaVuSans-Oblique-88\"/>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "   <g id=\"text_5\">\n",
       "    <g clip-path=\"url(#p29205ff115)\">\n",
       "     <!-- ${qreg}_{0}$ -->\n",
       "     <defs>\n",
       "      <path d=\"M 41.703125 8.203125 \n",
       "Q 38.09375 3.46875 33.171875 1.015625 \n",
       "Q 28.265625 -1.421875 22.3125 -1.421875 \n",
       "Q 14.015625 -1.421875 9.296875 4.171875 \n",
       "Q 4.59375 9.765625 4.59375 19.578125 \n",
       "Q 4.59375 27.484375 7.5 34.859375 \n",
       "Q 10.40625 42.234375 15.828125 48.09375 \n",
       "Q 19.34375 51.90625 23.90625 53.953125 \n",
       "Q 28.46875 56 33.5 56 \n",
       "Q 39.546875 56 43.453125 53.609375 \n",
       "Q 47.359375 51.21875 49.125 46.390625 \n",
       "L 50.6875 54.59375 \n",
       "L 59.71875 54.59375 \n",
       "L 45.125 -20.609375 \n",
       "L 36.078125 -20.609375 \n",
       "z\n",
       "M 13.921875 20.90625 \n",
       "Q 13.921875 13.671875 16.9375 9.890625 \n",
       "Q 19.96875 6.109375 25.6875 6.109375 \n",
       "Q 34.1875 6.109375 40.1875 14.234375 \n",
       "Q 46.1875 22.359375 46.1875 33.984375 \n",
       "Q 46.1875 41.015625 43.078125 44.75 \n",
       "Q 39.984375 48.484375 34.1875 48.484375 \n",
       "Q 29.9375 48.484375 26.3125 46.5 \n",
       "Q 22.703125 44.53125 20.015625 40.71875 \n",
       "Q 17.1875 36.71875 15.546875 31.34375 \n",
       "Q 13.921875 25.984375 13.921875 20.90625 \n",
       "z\n",
       "\" id=\"DejaVuSans-Oblique-113\"/>\n",
       "      <path d=\"M 44.578125 46.390625 \n",
       "Q 43.21875 47.125 41.453125 47.515625 \n",
       "Q 39.703125 47.90625 37.703125 47.90625 \n",
       "Q 30.515625 47.90625 25.140625 42.453125 \n",
       "Q 19.78125 37.015625 18.015625 27.875 \n",
       "L 12.5 0 \n",
       "L 3.515625 0 \n",
       "L 14.203125 54.6875 \n",
       "L 23.1875 54.6875 \n",
       "L 21.484375 46.1875 \n",
       "Q 25.046875 50.921875 30 53.453125 \n",
       "Q 34.96875 56 40.578125 56 \n",
       "Q 42.046875 56 43.453125 55.828125 \n",
       "Q 44.875 55.671875 46.296875 55.28125 \n",
       "z\n",
       "\" id=\"DejaVuSans-Oblique-114\"/>\n",
       "      <path d=\"M 48.09375 32.234375 \n",
       "Q 48.25 33.015625 48.3125 33.84375 \n",
       "Q 48.390625 34.671875 48.390625 35.5 \n",
       "Q 48.390625 41.453125 44.890625 44.921875 \n",
       "Q 41.40625 48.390625 35.40625 48.390625 \n",
       "Q 28.71875 48.390625 23.578125 44.15625 \n",
       "Q 18.453125 39.9375 15.828125 32.171875 \n",
       "z\n",
       "M 55.90625 25.203125 \n",
       "L 14.109375 25.203125 \n",
       "Q 13.8125 23.34375 13.71875 22.265625 \n",
       "Q 13.625 21.1875 13.625 20.40625 \n",
       "Q 13.625 13.625 17.796875 9.90625 \n",
       "Q 21.96875 6.203125 29.59375 6.203125 \n",
       "Q 35.453125 6.203125 40.671875 7.515625 \n",
       "Q 45.90625 8.84375 50.390625 11.375 \n",
       "L 48.6875 2.484375 \n",
       "Q 43.84375 0.53125 38.6875 -0.4375 \n",
       "Q 33.546875 -1.421875 28.21875 -1.421875 \n",
       "Q 16.84375 -1.421875 10.71875 4.015625 \n",
       "Q 4.59375 9.46875 4.59375 19.484375 \n",
       "Q 4.59375 28.03125 7.640625 35.375 \n",
       "Q 10.6875 42.71875 16.609375 48.484375 \n",
       "Q 20.40625 52.09375 25.65625 54.046875 \n",
       "Q 30.90625 56 36.8125 56 \n",
       "Q 46.09375 56 51.578125 50.4375 \n",
       "Q 57.078125 44.875 57.078125 35.5 \n",
       "Q 57.078125 33.25 56.78125 30.6875 \n",
       "Q 56.5 28.125 55.90625 25.203125 \n",
       "z\n",
       "\" id=\"DejaVuSans-Oblique-101\"/>\n",
       "      <path d=\"M 59.625 54.6875 \n",
       "L 50.296875 6.78125 \n",
       "Q 47.609375 -7.125 40.015625 -13.953125 \n",
       "Q 32.421875 -20.796875 19.578125 -20.796875 \n",
       "Q 14.84375 -20.796875 10.78125 -20.09375 \n",
       "Q 6.734375 -19.390625 3.21875 -17.921875 \n",
       "L 4.890625 -9.1875 \n",
       "Q 8.203125 -11.328125 11.90625 -12.34375 \n",
       "Q 15.625 -13.375 19.828125 -13.375 \n",
       "Q 28.375 -13.375 33.859375 -8.703125 \n",
       "Q 39.359375 -4.046875 41.109375 4.6875 \n",
       "L 41.890625 8.796875 \n",
       "Q 38.140625 4.5 33.15625 2.25 \n",
       "Q 28.171875 0 22.40625 0 \n",
       "Q 14.109375 0 9.34375 5.484375 \n",
       "Q 4.59375 10.984375 4.59375 20.609375 \n",
       "Q 4.59375 28.171875 7.46875 35.421875 \n",
       "Q 10.359375 42.671875 15.578125 48.296875 \n",
       "Q 19.046875 52 23.65625 54 \n",
       "Q 28.265625 56 33.296875 56 \n",
       "Q 38.8125 56 42.90625 53.4375 \n",
       "Q 47.015625 50.875 49.03125 46.1875 \n",
       "L 50.59375 54.6875 \n",
       "z\n",
       "M 46.09375 34.625 \n",
       "Q 46.09375 41.265625 42.96875 44.875 \n",
       "Q 39.84375 48.484375 34.078125 48.484375 \n",
       "Q 30.515625 48.484375 27.296875 47.0625 \n",
       "Q 24.078125 45.65625 21.78125 43.109375 \n",
       "Q 18.0625 38.921875 15.984375 33.234375 \n",
       "Q 13.921875 27.546875 13.921875 21.484375 \n",
       "Q 13.921875 14.75 17.0625 11.125 \n",
       "Q 20.21875 7.515625 26.125 7.515625 \n",
       "Q 34.671875 7.515625 40.375 15.25 \n",
       "Q 46.09375 23 46.09375 34.625 \n",
       "z\n",
       "\" id=\"DejaVuSans-Oblique-103\"/>\n",
       "      <path d=\"M 31.78125 66.40625 \n",
       "Q 24.171875 66.40625 20.328125 58.90625 \n",
       "Q 16.5 51.421875 16.5 36.375 \n",
       "Q 16.5 21.390625 20.328125 13.890625 \n",
       "Q 24.171875 6.390625 31.78125 6.390625 \n",
       "Q 39.453125 6.390625 43.28125 13.890625 \n",
       "Q 47.125 21.390625 47.125 36.375 \n",
       "Q 47.125 51.421875 43.28125 58.90625 \n",
       "Q 39.453125 66.40625 31.78125 66.40625 \n",
       "z\n",
       "M 31.78125 74.21875 \n",
       "Q 44.046875 74.21875 50.515625 64.515625 \n",
       "Q 56.984375 54.828125 56.984375 36.375 \n",
       "Q 56.984375 17.96875 50.515625 8.265625 \n",
       "Q 44.046875 -1.421875 31.78125 -1.421875 \n",
       "Q 19.53125 -1.421875 13.0625 8.265625 \n",
       "Q 6.59375 17.96875 6.59375 36.375 \n",
       "Q 6.59375 54.828125 13.0625 64.515625 \n",
       "Q 19.53125 74.21875 31.78125 74.21875 \n",
       "z\n",
       "\" id=\"DejaVuSans-48\"/>\n",
       "     </defs>\n",
       "     <g transform=\"translate(71.2699 48.011777)scale(0.1625 -0.1625)\">\n",
       "      <use xlink:href=\"#DejaVuSans-Oblique-113\"/>\n",
       "      <use transform=\"translate(63.476562 0)\" xlink:href=\"#DejaVuSans-Oblique-114\"/>\n",
       "      <use transform=\"translate(104.589844 0)\" xlink:href=\"#DejaVuSans-Oblique-101\"/>\n",
       "      <use transform=\"translate(166.113281 0)\" xlink:href=\"#DejaVuSans-Oblique-103\"/>\n",
       "      <use transform=\"translate(229.589844 -16.40625)scale(0.7)\" xlink:href=\"#DejaVuSans-48\"/>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "   <g id=\"text_6\">\n",
       "    <g clip-path=\"url(#p29205ff115)\">\n",
       "     <!-- ${qreg}_{1}$ -->\n",
       "     <defs>\n",
       "      <path d=\"M 12.40625 8.296875 \n",
       "L 28.515625 8.296875 \n",
       "L 28.515625 63.921875 \n",
       "L 10.984375 60.40625 \n",
       "L 10.984375 69.390625 \n",
       "L 28.421875 72.90625 \n",
       "L 38.28125 72.90625 \n",
       "L 38.28125 8.296875 \n",
       "L 54.390625 8.296875 \n",
       "L 54.390625 0 \n",
       "L 12.40625 0 \n",
       "z\n",
       "\" id=\"DejaVuSans-49\"/>\n",
       "     </defs>\n",
       "     <g transform=\"translate(71.2699 93.462777)scale(0.1625 -0.1625)\">\n",
       "      <use xlink:href=\"#DejaVuSans-Oblique-113\"/>\n",
       "      <use transform=\"translate(63.476562 0)\" xlink:href=\"#DejaVuSans-Oblique-114\"/>\n",
       "      <use transform=\"translate(104.589844 0)\" xlink:href=\"#DejaVuSans-Oblique-101\"/>\n",
       "      <use transform=\"translate(166.113281 0)\" xlink:href=\"#DejaVuSans-Oblique-103\"/>\n",
       "      <use transform=\"translate(229.589844 -16.40625)scale(0.7)\" xlink:href=\"#DejaVuSans-49\"/>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "   <g id=\"text_7\">\n",
       "    <g clip-path=\"url(#p29205ff115)\">\n",
       "     <!-- 2 -->\n",
       "     <g transform=\"translate(161.7334 127.754825)scale(0.104 -0.104)\">\n",
       "      <use xlink:href=\"#DejaVuSans-50\"/>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "   <g id=\"text_8\">\n",
       "    <g clip-path=\"url(#p29205ff115)\">\n",
       "     <!-- $creg$ -->\n",
       "     <defs>\n",
       "      <path d=\"M 53.609375 52.59375 \n",
       "L 51.8125 43.703125 \n",
       "Q 48.578125 46.046875 44.9375 47.21875 \n",
       "Q 41.3125 48.390625 37.40625 48.390625 \n",
       "Q 33.109375 48.390625 29.21875 46.875 \n",
       "Q 25.34375 45.359375 22.703125 42.578125 \n",
       "Q 18.5 38.328125 16.203125 32.609375 \n",
       "Q 13.921875 26.90625 13.921875 20.796875 \n",
       "Q 13.921875 13.421875 17.609375 9.8125 \n",
       "Q 21.296875 6.203125 28.8125 6.203125 \n",
       "Q 32.515625 6.203125 36.6875 7.328125 \n",
       "Q 40.875 8.453125 45.40625 10.6875 \n",
       "L 43.703125 1.8125 \n",
       "Q 39.796875 0.203125 35.671875 -0.609375 \n",
       "Q 31.546875 -1.421875 27.203125 -1.421875 \n",
       "Q 16.3125 -1.421875 10.453125 4.015625 \n",
       "Q 4.59375 9.46875 4.59375 19.578125 \n",
       "Q 4.59375 28.078125 7.640625 35.234375 \n",
       "Q 10.6875 42.390625 16.703125 48.09375 \n",
       "Q 20.796875 52 26.3125 54 \n",
       "Q 31.84375 56 38.375 56 \n",
       "Q 42.1875 56 45.9375 55.140625 \n",
       "Q 49.703125 54.296875 53.609375 52.59375 \n",
       "z\n",
       "\" id=\"DejaVuSans-Oblique-99\"/>\n",
       "     </defs>\n",
       "     <g transform=\"translate(72.9924 139.803972)scale(0.195 -0.195)\">\n",
       "      <use xlink:href=\"#DejaVuSans-Oblique-99\"/>\n",
       "      <use transform=\"translate(54.980469 0)\" xlink:href=\"#DejaVuSans-Oblique-114\"/>\n",
       "      <use transform=\"translate(96.09375 0)\" xlink:href=\"#DejaVuSans-Oblique-101\"/>\n",
       "      <use transform=\"translate(157.617188 0)\" xlink:href=\"#DejaVuSans-Oblique-103\"/>\n",
       "     </g>\n",
       "    </g>\n",
       "   </g>\n",
       "  </g>\n",
       " </g>\n",
       " <defs>\n",
       "  <clipPath id=\"p29205ff115\">\n",
       "   <rect height=\"154.5334\" width=\"340.8825\" x=\"7.2\" y=\"7.2\"/>\n",
       "  </clipPath>\n",
       " </defs>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<Figure size 451.5x204.68 with 1 Axes>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gen_ansatz(np.pi / 2)[0].draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating expectation\n",
    "\n",
    "Using [Davit Khachatryan's explanation of the 1-qubit example](https://github.com/DavitKhach/quantum-algorithms-tutorials/blob/master/variational_quantum_eigensolver.ipynb), we can do similar calculations to figure out how to get the expectation of 2-qubit observables: \n",
    "\n",
    "Since we know that all the Pauli operators have the same eigenvalues, once we calculate what we need to measure for one Pauli operator, we know it'll be the same thing for the others in their respective bases. For simplicity, I'll do the calculations in the $Z \\otimes Z$ basis because it's the computational basis. \n",
    "\n",
    "For some qubit $\\ket{\\psi}$ we have \n",
    "$$\\ket{\\psi} = c_{00}^z \\ket{00} + c_{01}^z \\ket{01} + c_{10}^z \\ket{10} + c_{11}^z \\ket{11}$$\n",
    "\n",
    "Calculating expectation: \n",
    "$$ \\bra{\\psi} Z \\otimes Z \\ket{\\psi} = \\bra{\\psi} \\big( c_{00}^z \\ket{00} - c_{01}^z \\ket{01} - c_{10}^z \\ket{10} + c_{11}^z \\ket{11} \\big) = \\mid c_{00}^z \\mid^2 - \\mid c_{01}^z \\mid^2 - \\mid c_{10}^z \\mid^2 + \\mid c_{11}^z \\mid^2$$\n",
    "\n",
    "So this means we need to measure the ansatz in the computational basis and then combine the probabilities like so: $$prob(00) + prob(11) - prob(01) - prob(10)$$. \n",
    "\n",
    "Since the Pauli $X$ and $Y$ also have +1, -1 eigenvalues for their respective bases, we know we'll get the same outcome for expectation. This means our above equation for calculating expectation from probabilities works for them too. \n",
    "\n",
    "Notice we don't have to calculate the expectation for observable $I \\otimes I$ because \n",
    "$$ \\langle I \\otimes I \\rangle = \\bra{\\psi} I \\otimes I \\ket{\\psi} = \\langle \\psi \\rvert \\psi \\rangle = 1 $$\n",
    "so instead we can just subtract $\\frac{1}{2}$ from the remaining expectation.\n",
    "\n",
    "## Changing basis\n",
    "\n",
    "Finally we have to measure the Pauli observables in their respective eigenbases. \n",
    "\n",
    "$Z \\otimes Z$ has the computational basis as its eigenbasis, so we don't need to do anything. \n",
    "\n",
    "$X \\otimes X$ has the Hadamard basis as its eigenbasis, since $X \\ket{+} = \\ket{+}$ and $X \\ket{-} = - \\ket{-}$. So we need to apply a basis transformation that maps $\\ket{+}, \\ket{-} \\rightarrow \\ket{0}, \\ket{1}$. Of course, the Hadamard gate does this, so we need to just apply it to both qubits in the ansatz. \n",
    "\n",
    "$Y \\otimes Y$ is a little less conventional. $Y$ has eigenvectors \n",
    "$$ \\ket{+i} = \\frac{1}{\\sqrt{2}} \\begin{bmatrix} 1\\\\ i \\end{bmatrix} \\text{ and } \\ket{-i} = \\frac{1}{\\sqrt{2}} \\begin{bmatrix} 1\\\\ -i \\end{bmatrix} $$\n",
    "So we need a unitary matrix to map $\\ket{+i}, \\ket{-i} \\rightarrow \\ket{0}, \\ket{1}$. This is simpler than it seems with the following lemma: \n",
    "\n",
    "**Lemma 1:** For orthonormal bases $\\ket{v_i}, \\ket{w_i}$ such that $\\ket{w_i} = U \\ket{v_i}$, $U = \\sum_i  \\ket{w_i} \\bra{v_i}$. \n",
    "\n",
    "*Proof:* (from Nielsen and Chuang p71) $$ \\sum_i \\ket{w_i} \\bra{v_i} = \\sum_i U \\ket{v_i} \\bra{v_i} = U I = U $$ where the second equality follows from the Completeness Relation. \n",
    "\n",
    "Using this, the unitary we need to measure in $Y$ basis is simply: $$ U = \\ket{0} \\bra{+i} + \\ket{1} \\bra{-i} = \\frac{1}{\\sqrt{2}} \\begin{bmatrix} 1 & -i\\\\ 1 & i \\end{bmatrix}  $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def measure(basis, theta):\n",
    "    \"\"\"We change to appropriate basis and measure.\n",
    "        basis: 'XX', 'YY', or 'ZZ'\n",
    "    \"\"\"\n",
    "\n",
    "    qc, qr, cr = gen_ansatz(theta)\n",
    "\n",
    "    if basis == 'XX':\n",
    "        qc.h( qr[0] )\n",
    "        qc.h( qr[1] )\n",
    "    elif basis == 'YY': \n",
    "        qc.u3( np.pi/2, 0, np.pi/2, qr[0] ) # Qiskit's U3 gate can describe any unitary - it defines our Y gate\n",
    "        qc.u3( np.pi/2, 0, np.pi/2, qr[1] )\n",
    "    elif basis == 'ZZ': \n",
    "        pass\n",
    "    else: \n",
    "        raise ValueError(\"Basis must be one of 'XX', 'YY', 'ZZ'.\")\n",
    "    \n",
    "    qc.measure( qr[0],cr[0] )\n",
    "    qc.measure( qr[1],cr[1] )\n",
    "    \n",
    "    backend = Aer.get_backend('qasm_simulator')\n",
    "    job = execute( qc, backend, shots=shots )\n",
    "    \n",
    "    return job.result().get_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "thetas = []\n",
    "\n",
    "# This is our function to be minimized so it needs to output a float (the upper bound on ground state energy)\n",
    "def vqe(theta): \n",
    "    if visualization:\n",
    "        thetas.append(theta)\n",
    "    \n",
    "    # Measure in XX basis \n",
    "    xx = measure('XX', theta)\n",
    "    \n",
    "    # Measure in YY basis \n",
    "    yy = measure('YY', theta)\n",
    "    \n",
    "    # Measure in ZZ basis \n",
    "    zz = measure('ZZ', theta)\n",
    "    \n",
    "    # Calculate expectation by combining probabilities like so: 00 + 11 - 01 - 10 \n",
    "    exp = 0\n",
    "    for dic in [xx, yy, zz]:\n",
    "        for ket in ['00', '11']:\n",
    "            exp += dic.get(ket, 0) / shots\n",
    "        for ket in ['01', '10']: \n",
    "            exp -= dic.get(ket, 0) / shots \n",
    "        \n",
    "    # Subtract 1 for II, apply coefficient of 1/2, then return \n",
    "    return 0.5 * (exp - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/warren/anaconda3/lib/python3.6/site-packages/scipy/optimize/_minimize.py:522: RuntimeWarning: Method Powell cannot handle constraints nor bounds.\n",
      "  RuntimeWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "VQE came up with an upper bound of -2.0\n"
     ]
    }
   ],
   "source": [
    "sol = minimize_scalar(vqe, bounds=(0, 2*np.pi), method='bounded')\n",
    "\n",
    "print(\"VQE came up with an upper bound of \" + str(sol.fun))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(3.12975511)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the eigenvalues using NumPy so we can check our answers... The lowest eigenvalue is -2 which is what we got too!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 0., -2.,  0.,  0.]),\n",
       " array([[ 0.        ,  0.        ,  1.        ,  0.        ],\n",
       "        [ 0.70710678,  0.70710678,  0.        ,  0.        ],\n",
       "        [ 0.70710678, -0.70710678,  0.        ,  0.        ],\n",
       "        [ 0.        ,  0.        ,  0.        ,  1.        ]]))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.eig(hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Investigating our arrival at the solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: -2.0\n",
       " message: 'Solution found.'\n",
       "    nfev: 14\n",
       "  status: 0\n",
       " success: True\n",
       "       x: 3.1440226061471668"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Here's the full OptimizeResult object SciPy gives us after minimizing - notice nfev is the number of inputs the minimizer had to try to reach the the optimum\n",
    "sol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since there weren't too many iterations with the given ansatz for us to reach the minimum, let's look at how the theta value for the ansatz progressed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.399963229728653\n",
      "3.8832220774509327\n",
      "4.799926459457306\n",
      "3.144608484320108\n",
      "3.1440163034967115\n",
      "3.1443123939084097\n",
      "3.144199297434884\n",
      "3.1441293999702373\n",
      "3.1440862009613584\n",
      "3.144059502505591\n",
      "3.1440430019524794\n",
      "3.1440328040498233\n",
      "3.1440265013993676\n",
      "3.1440226061471668\n"
     ]
    }
   ],
   "source": [
    "for theta in thetas: \n",
    "    print(theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's an ansatz from [this paper](https://arxiv.org/pdf/quant-ph/0308033.pdf) that I came across while on the [Qiskit explanation of VQE](https://qiskit.org/textbook/ch-applications/vqe-molecules.html) that's supposed to be able to get to every single 2-qubit state. \n",
    "\n",
    "Notice the crazy number of parameters which is why having good heuristics is recommended. Here's what our variational circuit looks like: \n",
    "![From Qiskit](https://qiskit.org/textbook/ch-applications/images/two_qubit_var_form.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_ansatz(theta):\n",
    "    # Overkill ansatz; notice we're redefining gen_ansatz so we don't have to rewrite vqe function\n",
    "    t0, p0, l0, t1, p1, l1, t2, p2, l2, t3, p3, l3, t4, p4, l4, t5, p5, l5, t6, p6, l6, t7, p7, l7 = theta\n",
    "    qc = QuantumCircuit()\n",
    "    qr = QuantumRegister(2, 'qreg')\n",
    "    qc.add_register( qr )\n",
    "    qc.u3(t0, p0, l0, qr[0]) \n",
    "    qc.u3(t1, p1, l1, qr[1]) \n",
    "    qc.cnot(qr[1], qr[0])\n",
    "    qc.u3(t2, p2, l2, qr[0])\n",
    "    qc.u3(t3, p3, l3, qr[1])\n",
    "    qc.cnot(qr[0], qr[1])\n",
    "    qc.u3(t4, p4, l4, qr[0])\n",
    "    qc.u3(t5, p5, l5, qr[1])\n",
    "    qc.cnot(qr[1], qr[0])\n",
    "    qc.u3(t6, p6, l6, qr[0])\n",
    "    qc.u3(t7, p7, l7, qr[1])\n",
    "    cr = ClassicalRegister(2, 'creg')\n",
    "    qc.add_register(cr)\n",
    "    return qc, qr, cr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This takes a while to run\n",
    "sol = minimize(vqe, [0 for i in range(24)], method='Powell')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(-1.9972)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol.fun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "My computer got pretty hot, but the minimizer actually still arrived at the correct solution! When I left the 'method' parameter blank, I kept getting stuck at local minima around -0.02."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol.success"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1217"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The amount of iterations it took to find the minimum \n",
    "sol.nfev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discovering excited states of the Hamiltonian\n",
    "\n",
    "While reading *Quantum Computing: An Applied Approach* by Jack Hidary, I saw a reference for [Variational Quantum Computation of Excited States](https://arxiv.org/abs/1805.08138). It's a method the authors call Variational Quantum Deflation, which allows us to find other eigenvalues of a Hamiltonian using VQE, instead of just the lowest one. The reasoning behind it is simple: we iteratively run VQE but after each iteration we only look for eigenvectors that are orthogonal to all previously-found eigenvectors. Here's the cost function: \n",
    "$$Cost(\\theta_k) = \\bra{\\psi(\\theta_k)} H \\ket{\\psi(\\theta_k)} + \\sum_{i=0}^{k-1} \\beta_i \\mid \\bra{\\psi(\\theta_k)} \\psi(\\theta_i) \\rangle \\mid^2  $$\n",
    "where $\\beta_i$ is some chosen coefficient. \n",
    "\n",
    "Notice that this is the same cost function as normal VQE except now we want the newest eigenvector to be orthogonal to all previous eigenvectors. \n",
    "\n",
    "### How to calculate inner product on a quantum computer \n",
    "\n",
    "What makes working with quantum computers interesting and challenging is that our only method of learning information is through measurement, which collapses the state. So we need to be really clever with how we do measurements to learn exactly the information we want. \n",
    "\n",
    "This paper is a good example of this because they provide a way to find out the value of $\\mid \\bra{\\psi(\\theta_k)} \\psi(\\theta_i) \\rangle \\mid^2$. \n",
    "\n",
    "One way of figuring this out is to decompose $\\psi(\\theta) = A(\\theta) \\ket{00}$, where $A$ is our ansatz circuit. Now we have \n",
    "$$ \\mid \\bra{\\psi(\\theta_k)} \\psi(\\theta_i) \\rangle \\mid^2 = \\mid \\bra{00} A(\\theta_k)^\\dagger A(\\theta_i) \\ket{00}$$\n",
    "\n",
    "Well this is easy to find out with a quantum computer: it's the probability of getting $\\ket{00}$ when measuring the state $A(\\theta_k)^\\dagger A(\\theta_i) \\ket{00}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's try it with our original ansatz\n",
    "def gen_ansatz(theta):\n",
    "    \"\"\"Returns QuantumCircuit, QuantumRegister, and ClassicalRegister components of ansatz \"\"\"\n",
    "    # Recommended ansatz: (IX) CX (RZ I) (HI) |00> \n",
    "    qc = QuantumCircuit()\n",
    "    qr = QuantumRegister(2,'qreg')\n",
    "    qc.add_register( qr )\n",
    "    qc.h( qr[0] )\n",
    "    qc.rz(theta, qr[0] )\n",
    "    qc.cx( qr[0], qr[1] )\n",
    "    qc.x( qr[1] )\n",
    "    cr = ClassicalRegister(2,'creg')\n",
    "    qc.add_register(cr)\n",
    "\n",
    "    #return qc.draw(output='mpl')\n",
    "    return qc, qr, cr\n",
    "\n",
    "# We'll need the inverse ansatz in order to do the algorithm\n",
    "def gen_combined_ansatz(theta_i, theta_k):\n",
    "    ansatz = gen_ansatz(theta_i)\n",
    "    # Append adjoint ansatz circuit to qc of ansatz\n",
    "    ansatz[0].extend(gen_ansatz(theta_k)[0].inverse())\n",
    "    return ansatz\n",
    "\n",
    "# Variational quantum deflator\n",
    "def vqd(hamiltonian, initial_guess):\n",
    "    thetas_so_far = [] \n",
    "    eigenvalues = []\n",
    "    beta = 1\n",
    "    num_eigenvectors = np.shape(hamiltonian)[0] # Assume non-degenerate \n",
    "    \n",
    "    # Cost function for VQD; we compute VQE then force orthogonality between eigenvectors\n",
    "    def cost(theta):\n",
    "        exp = vqe(theta) \n",
    "        \n",
    "        for theta_i in thetas_so_far:\n",
    "            qc, qr, cr = gen_combined_ansatz(theta_i, theta)\n",
    "            qc.measure( qr[0], cr[0] )\n",
    "            qc.measure( qr[1], cr[1] )\n",
    "            backend = Aer.get_backend('qasm_simulator')\n",
    "            job = execute( qc, backend, shots=shots )\n",
    "            exp += beta * job.result().get_counts().get('00', 0) / shots\n",
    "        return exp\n",
    "    \n",
    "    for i in range(num_eigenvectors):\n",
    "        sol = minimize(cost, initial_guess, method='Powell')\n",
    "        if sol.success:\n",
    "            thetas_so_far.append(sol.x)\n",
    "            eigenvalues.append(sol.fun)\n",
    "        else: \n",
    "            raise ValueError('Minimization failed.')\n",
    "    return eigenvalues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array(-2.), array(-1.0008), array(-0.0202), array(0.0046)]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vqd(hamiltonian, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 0., -2.,  0.,  0.]),\n",
       " array([[ 0.        ,  0.        ,  1.        ,  0.        ],\n",
       "        [ 0.70710678,  0.70710678,  0.        ,  0.        ],\n",
       "        [ 0.70710678, -0.70710678,  0.        ,  0.        ],\n",
       "        [ 0.        ,  0.        ,  0.        ,  1.        ]]))"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# For this hamiltonian, we have eigenvalues of -2, 0, 0, 0 and an orthonormal basis of eigenvectors\n",
    "np.linalg.eig(hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we weren't able to get the exact eigenvalues - instead we got -2, -1, 0, 0. \n",
    "\n",
    "I can think of 2 reasons for why our algorithms got the first eigenvalue correct but not the others: \n",
    "\n",
    "1. The beta value needs to be optimized. Varying how much we care about orthogonality compared to VQE results will dictate our results. \n",
    "\n",
    "2. The ansatz isn't able to reach all the eigenvectors of the Hamiltonian. This means VQD never gets to the global minimum. This explains why we always get the first eigenvalue of -2 correct: we know the ansatz can reach that specific eigenvector (because the ansatz was given to us as a good ansatz for VQE). \n",
    "\n",
    "When I made beta really high, I got useless upper bounds like 100, which means that enforcing orthogonality doesn't work with this ansatz. This makes me believe problem 2 is why VQD is failing. "
   ]
  }
 ],
 "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
}