vqe_random_spin_exercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational quantum eigensolving\n",
    "\n",
    "In this exercise, we exemplify how the VQE method can be used to find the approximate ground state (GS) of a given Hamiltonian. Here, we take a random Hamiltonian for concreteness.\n",
    "\n",
    "## Defining the Hamiltonian\n",
    "\n",
    "### The Heisenberg Hamiltonian\n",
    "\n",
    "We start by defining the Hamiltonian we want to find an approximate ground state of. Let us start with the Heisenberg Hamiltonian:\n",
    "\n",
    "$$H=J\\left(\\hat{S}_{1}^{x}\\hat{S}_{2}^{x}+\\hat{S}_{1}^{y}\\hat{S}_{2}^{y}+\\hat{S}_{1}^{z}\\hat{S}_{2}^{z}\\right).$$\n",
    "\n",
    "In myQLM, Hamiltonians are implemented with ``Observable`` objects (https://myqlm.github.io/qat-core.html#module-qat.core.wrappers.observable).\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from qat.core import Observable, Term\n",
    "\n",
    "nqbits = 2\n",
    "pauli_terms = #...\n",
    "\n",
    "hamiltonian = Observable(nqbits, pauli_terms=terms)\n",
    "print(\"H:\", hamiltonian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because $H$ is defined on a small number of qubits, we can afford to diagonalize it exactly and thus compute the exact GS energy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils_tuto import make_matrix\n",
    "H_mat = make_matrix(hamiltonian)\n",
    "eigvals = np.linalg.eigvalsh(H_mat)\n",
    "E0 = min(eigvals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the size of the matrix representation of $H$? For a number of qubits $n$, what is the size in general?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A slighly more complicated Hamiltonian\n",
    "\n",
    "For later purposes, we define a random Hamiltonian with 5 terms acting on 3 qubits. (Do not execute the cell if you want to work with the Heisenberg Hamiltonian first)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_random_hamiltonian(nqbits, nterms):\n",
    "    terms = []\n",
    "    for _ in range(nterms):\n",
    "        coeff = np.random.random()\n",
    "        ops = \"\".join(np.random.choice([\"X\", \"Z\"], size=nqbits))\n",
    "        qbits = np.random.choice(nqbits, size=nqbits, replace=False)\n",
    "        terms.append(Term(coefficient=coeff, pauli_op=ops, qbits=qbits))\n",
    "    hamiltonian = Observable(nqbits, pauli_terms=terms)\n",
    "    return hamiltonian\n",
    "\n",
    "\n",
    "nqbits = 3\n",
    "nterms = 5\n",
    "np.random.seed(1423543) #fixing seed to have reproducible results\n",
    "hamiltonian = make_random_hamiltonian(nqbits, nterms)\n",
    "print(\"H:\", hamiltonian)\n",
    "\n",
    "H_mat = make_matrix(hamiltonian)\n",
    "eigvals_random = np.linalg.eigvalsh(H_mat)\n",
    "E0 = min(eigvals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constructing a variational circuit\n",
    "\n",
    "The main task in VQE consists in designing a 'good' ansatz, i.e a family of variational states that can come close enough to the exact ground state. These states are constructed using quantum circuits with parameters (like rotation angles) that will be optimized to minimize the energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qat.lang.AQASM import Program, QRoutine, RY, CNOT, RX, Z, H, RZ\n",
    "\n",
    "def make_ansatz():\n",
    "    \"\"\"Function that prepares an ansatz circuit\n",
    "    \n",
    "    Returns:\n",
    "        Circuit: a parameterized circuit (i.e with some variables that are not set)\n",
    "    \"\"\"\n",
    "    \n",
    "    nparams = ... # number of parameters \n",
    "    \n",
    "    prog = Program()\n",
    "    reg = prog.qalloc(nqbits)\n",
    "    # define variables using 'new_var' \n",
    "    theta = [prog.new_var(float, '\\\\theta_%s'%i)\n",
    "             for i in range(nparams)]\n",
    "    \n",
    "    # for instance...: put a rotation on each qubit\n",
    "    for ind in range(nqbits):\n",
    "        RY(theta[ind])(reg[ind])\n",
    "        \n",
    "    # and so on\n",
    "    # ...\n",
    "    \n",
    "    circ = prog.to_circ()\n",
    "    return circ\n",
    "\n",
    "circ = make_ansatz()\n",
    "#%qatdisplay circ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Creating a variational job and a variational stack\n",
    "\n",
    "We now create a variational job from this circuit and observable.\n",
    "It is then submitted to a \"variational stack\" composed of a perfect QPU, ``LinAlg``, and a variational plugin, ``ScipyMinimizePlugin``. The latter handles parametric jobs. These are jobs whose circuit contains a parameter that is then going to be optimized, using classical optimizers, so as to minimize the value of the observable over the final state.\n",
    "\n",
    "Below, we are going to test three different classical optimizers: COBYLA, Nelder-Mead, and BFGS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qat.qpus import get_default_qpu\n",
    "from qat.plugins import ScipyMinimizePlugin #,SPSAMinimizePlugin\n",
    "\n",
    "circ = make_ansatz()\n",
    "\n",
    "job = circ.to_job(job_type=\"OBS\",\n",
    "                  observable=hamiltonian,\n",
    "                  nbshots=100)\n",
    "# note that this job contains a circuit whose variables are not set!\n",
    "# it therefore requires a special plugin so that the QPU can handle it!\n",
    "\n",
    "theta_0 = ... # initial value of the parameters\n",
    "\n",
    "linalg_qpu = get_default_qpu()\n",
    "method = \"COBYLA\"  # look at the doc of scipy.optimize.minimize for other options\n",
    "optimizer_scipy = ScipyMinimizePlugin(method=method,\n",
    "                                      tol=1e-6,\n",
    "                                      options={\"maxiter\": 200},\n",
    "                                      x0=theta_0)\n",
    "# you could also use SPSAMinimizePlugin instead\n",
    "\n",
    "# we now build a stack that can handle variational jobs\n",
    "qpu = optimizer_scipy | linalg_qpu\n",
    "\n",
    "# we submit the job\n",
    "result = qpu.submit(job)\n",
    "\n",
    "print(\"Minimum VQE energy (%s) = %s\"%(method, result.value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What do you observe? It the VQE energy close to the exact energy ``E0``? If not, where does it come from: a poor ansatz? a poor optimizer? a poor QPU? not enough ``nbshots``? Try playing with these various knobs. Hint: entanglement helps!\n",
    "\n",
    "\n",
    "## Plotting the results\n",
    "\n",
    "We can also plot the value of the variational energy over the course of the classical optimization. For this, we can retrieve information about the variational job execution in the ``meta_data`` field of the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(eval(result.meta_data['optimization_trace']),\n",
    "         label=method)\n",
    "    \n",
    "plt.plot([E0 for _ in range(300)], '--k', lw=3, label=\"exact energy\")\n",
    "    \n",
    "plt.grid()\n",
    "plt.legend(loc=\"best\")\n",
    "plt.xlabel(\"Steps\")\n",
    "plt.ylabel(\"Energy\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What can you observe in terms of convergence speed?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you have found an ansatz and an optimizer that works with a perfect QPU, can we now test what happens with a noisy QPU?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## With a noisy QPU\n",
    "\n",
    "We can now execute the circuit on a noisy QPU simulator. To achieve this, we are going to turn our perfect QPU ``qpu`` into a noisy one by simply combining it with a so-called plugin that adds \"depolarizing noise\" on the circuit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from depolarizing_plugin import DepolarizingPlugin\n",
    "depol_plugin = DepolarizingPlugin(prob_1qb=0.005, prob_2qb=0.05)\n",
    "qpu = get_default_qpu()\n",
    "noisy_qpu = depol_plugin | qpu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Depolarizing noise, as any noise, is defined by its action on the density matrix $\\rho$:\n",
    "\n",
    "$$\\mathcal{E}(\\rho) = (1-p) \\rho + \\frac{p}{3} \\left( X \\rho X + Y \\rho Y + Z \\rho Z \\right)$$\n",
    "\n",
    "with $X$, $Y$, and $Z$ the three Pauli matrices, and $p$ the probability that a depolarization event happens. (This is for one-qubit depolarizing noise; a similar expression holds for two-qubit depolarizing noise, with $X$, $Y$ and $Z$ replaced by tensor products of Pauli matrices).\n",
    "\n",
    "In other words, depolarizing noise leaves the qubit unchanged with probability $1-p$, and otherwise applies a $X$, $Y$ or $Z$ gate with probability $p/3$. $p$ is what is called ``prob_1qb`` above. \n",
    "\n",
    "This type of noise is not exactly the type of noise that is observed in actual quantum computers, but it is close enough and it allows to understand the main effects of noise.\n",
    "\n",
    "Here, the plugin adds depoalizing noise after each gate to mimic an imperfect gate.\n",
    "\n",
    "- Check that this noise preserves the trace of the density matrix\n",
    "- What are the corresponding Kraus operators?\n",
    "\n",
    "We now execute the circuit on this noisy QPU.\n",
    "\n",
    "What happens when you change the ``linalg_qpu`` object above to ``noisy_qpu``? What happens when you change the depolarizing probabilities?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bonus\n",
    "How could you test the robustness of VQE against systematic errors, without explicitly modifying the ``DepolarizingPlugin``?"
   ]
  },
  {
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}