estimates.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimating the depth for classical simulability for QAOA under depolarizing noise\n",
    "\n",
    "The purpose of this notebook is to clarify how the estimats in Table X of [\"Limitations of optimization algorithms on noisy quantum devices\"](https://arxiv.org/pdf/2009.05532.pdf) were produced. In that table, we estimated the maximal depth before a noisy quantum computer performing the quantum approximate optimization algorithm (QAOA) has a performance that is similar to simple Monte Carlo algorithms when performing optimization of Ising models. The choice of models we analysed was inspired by the recent experiments of [\"Quantum Approximate Optimization of Non-Planar Graph Problems on a Planar Superconducting Processor\"](http://arxiv.org/abs/2004.04197).\n",
    "We refer to our paper [\"Limitations of optimization algorithms on noisy quantum devices\"](https://arxiv.org/pdf/2009.05532.pdf) for a technical discussion of the concepts introduced here. The notebook will solely show how to obtain a maximal depth estimate given the relevant parameters of the problem and the device. These are the noise rates of the device for layers of 1-qubit and 2-qubit gates, how many physical layers of 2-qubit gates are required for each layer of 2-qubit QAOA and...\n",
    "After that, for each example, the notebook will output a table with the estimates given by our methods. If you do not change the values presented here, it should yield the values shown in Table.\n",
    "## 3-regular MAX-CUT\n",
    "This is arguably the simplest class of examples considered in [\"Quantum Approximate Optimization of Non-Planar Graph Problems on a Planar Superconducting Processor\"](http://arxiv.org/abs/2004.04197).\n",
    "Given a graph $G=(V,E)$ that is 3-regular (each vertex has three neighbours), the goal is to minimize the Hamiltonian\n",
    "\\begin{equation}\n",
    "H_I=-\\sum\\limits_{i\\sim j}Z_iZ_j.\n",
    "\\end{equation}\n",
    "The first numbers we need to perform the estimate are how many layers of 1-qubit gates and 2-qubit gates are required to perform one layer of the QAOA unitary. In case the topology of the problem does not match the device's, it is usually necessary to perform a number of layers that scales with the number of qubits, which we will denote by $n$. We will call the number of 1-qubit layers by $L_1(n)$ and 2-qubit $L_2(n)$. For the experiment, they were given by XX, as we set below.\n",
    "Define the function of 1-qubit layers L1 and 2-qubit L2-below in terms of n below:      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def L1(n): return n \n",
    "def L2(n): return 7*n\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to input the depolarizing probability for the 1-qubit and 2-qubit gates, which we will denote by p_1 and p_2, respectively. It is also possible to define the measurement error, p_m."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_1=0.001 \n",
    "p_2=0.01\n",
    "p_m=0.01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to input what is the maximal inverse temperature $\\beta_C$ we can simulate classically efficiently for the family of Gibbs states induced by the Hamiltonian. For instance, for MAX-CUT in 3-regular graphs, the results discussed in the section \"SOLVING MAX-CUT ON NOISY DEVICES\" yield that for $\\beta\\leq1/3$ we can simulate the Gibbs state in polynomial time. Define the maximal inverse temperature beta below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beta=1/3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to define the error we tolerate. That is, the $\\epsilon$ in the statement of Propostion 1. You may want to see the result for several different values of $\\epsilon$, so we will actually define a list epsilon_list for the different values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "epsilon_list=[0.1,0.001]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to do some number crunching! Denote by $Q_k$ the quantum channel which describes the $k$ layers of noisy QAOA.\n",
    "The total relative entropy contraction for $n$ qubits and $k$ layers of the QAOA unitary will be given by:\n",
    "\\begin{equation}\n",
    "D(Q_k(\\rho)\\|\\sigma_0)\\leq (1-p_1)^{2kL_1(n)}(1-p_2)^{2kL_2(n)}(1-p_m)n.\n",
    "\\end{equation}\n",
    "Now note that we have $\\|H_I\\|\\geq \\frac{3n}{2}$. To see this, note that if we assign the value $1$ to all th spins of the graph, we achieve the value $-\\sum\\limits_{i\\sim j}1=\\frac{3n}{2}$, as the graph is 3-regular and, thus, has \\frac{3n}{2} vertices. On the other hand, a triangle inequality shows that $\\|H_I\\|\\leq \\frac{3n}{2}$ and we conclude $\\|H_I\\|=\\frac{3n}{2}$. Eq. (?) of the paper then asserts that:\n",
    "\\begin{equation}\n",
    "(1-p_1)^{2kL_1(n)}(1-p_2)^{2kL_2(n)}(1-p_m)n\\leq \\epsilon \\beta\\|H_I\\|\n",
    "\\end{equation}\n",
    "is a sufficient condition for classical simulability.\n",
    "Inserting the values discussed before for $\\beta,\\|H_I\\|$ and simplifying we obtain the condition:\n",
    "\\begin{equation}\n",
    "(1-p_1)^{2kL_1(n)}(1-p_2)^{2kL_2(n)}(1-p_m)\\leq \\frac{3}{2}\\beta\\epsilon \n",
    "\\end{equation}\n",
    "to obtain an estimate of the energy up to $\\frac{3n}{2}\\epsilon$.\n",
    "We can now define a function rel_ent which computes the L.H.S. of the equation above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rel_ent(k,n):\n",
    "    #contribution from 1-qbuit gates\n",
    "    contrib1=(1-p_1)**(2*L1(n)*k)\n",
    "    #contribution from 2-qbuit gates\n",
    "    contrib2=(1-p_2)**(2*L2(n)*k)\n",
    "    #contribution from measurement\n",
    "    contribm=(1-p_m)\n",
    "    return contrib1*contrib2*contribm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can start the number crunching! For instance, given a prescribed number of QAOA rounds, we can estimate how up to how mmany qubits we can go before we violate. Of course, for this not to return $+\\infty$, we need that the number of layers grows with the number of qubits. Let us now write a function that does exactly that:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def max_depth(QAOA_rounds,epsilon):\n",
    "    n_qubits=1\n",
    "    while rel_ent(QAOA_rounds,n_qubits)>1.5*beta*epsilon:\n",
    "        n_qubits+=1\n",
    "    return n_qubits\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now define a list of how many QAOA rounds we wish to implement and generate a table for the different values of epsilon. The values generated here correspond to row X of the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "QAOA_rounds=[7,10,11]\n",
    "Table=np.zeros([len(QAOA_rounds)+1,len(epsilon_list)+1])\n",
    "Table[1:,0]=QAOA_rounds\n",
    "Table[0,1:]=epsilon_list\n",
    "for rounds in range(1,len(QAOA_rounds)+1):\n",
    "    for eps in range(1,len(epsilon_list)+1):\n",
    "        \n",
    "        Table[rounds,eps]=max_depth(Table[rounds,0],Table[0,eps])\n",
    "df=pd.DataFrame(data=Table[1:,1:],index=Table[1:,0],columns=Table[0,1:])\n",
    "df.index.name = 'Rounds'\n",
    "df.columns.name = 'Error $\\epsilon$'\n",
    "print(\"Number of qubits before advantage is lost for error and number of QAOA rounds\")\n",
    "df.head()\n"
   ]
  },
  {
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}