Demo_TFI_FF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demo: Simulation of QAOA on large transverse field Ising chains via free fermions\n",
    "\n",
    "Here we show how to use the functionalities contained in tfi_ff.py\n",
    "They are based on the mapping of the TFI to free fermions as discussed in [our paper](https://arxiv.org/abs/2004.14666). This solution is not original work done by us but the implementation is. Selected prior references for the mapping - using slightly different notation/bases - are \n",
    "* [Wang, Hadfield, Jiang, Rieffel - 2017](https://arxiv.org/abs/1706.02998) - also published in PRA\n",
    "* [ E. Lieb, T. Schultz, and D. Mattis, Annals of Physics16,407  (1961)](https://www.math.ucdavis.edu/~bxn/lieb1961.pdf) - original work, not as close to the presentation here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import tfi_ff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Energy evaluation \n",
    "For applications that require only the energy expectation value for given parameters, use the _eval\\_energy_ method:\n",
    "Choose a number of qubits N and a number of QAOA blocks p, where the ansatz  \n",
    "\\begin{equation}|\\psi(\\vec{\\theta},\\vec{\\varphi})\\rangle = \\prod_{j=p}^{1} L_X(\\varphi_j)L_{ZZ}(\\theta_j)|\\bar{\\psi}\\rangle\\end{equation}\n",
    "will contain the ground state of the TFI model  \n",
    "\\begin{equation}H_{TFI} = -\\sum_{k=1}^{N} Z^{(k)}Z^{(k+1)}-t\\sum_{k=1}^{N}X^{(k)}\\end{equation}\n",
    "for $p\\geq \\left\\lfloor N/2\\right\\rfloor$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The energy of the state prepared with the angles\n",
      "[-0.597 -0.031 -0.209  0.516  0.694  0.068 -0.058 -0.57   0.174  1.396\n",
      "  0.446  1.3    1.176  0.402 -0.66  -0.142  1.101  0.537  0.021  0.792]\n",
      "has the energy per site E=0.2027 in the TFI on 20 spins with transverse field 1.2\n"
     ]
    }
   ],
   "source": [
    "# Fix number of qubits\n",
    "N = 20\n",
    "# Fix number of QAOA blocks. For $p\\geq N/2$\n",
    "p = 10\n",
    "# Choose a transverse field strength \n",
    "t = 1.2\n",
    "# Generate 2 random parameters per block (theta and phi). Periodicity of the ansatz in all parameters is pi,\n",
    "# here we choose to center the parameters around 0 but with i.i.d. distribution over the whole hypercube \n",
    "# [-pi/2, pi/2]^(2p)\n",
    "param = np.random.random(2*p)*np.pi-np.pi/2\n",
    "# Pass the parameters, the system size and the transverse field to the evaluation function\n",
    "E = tfi_ff.eval_energy(param, N, t)\n",
    "print(f'The energy of the state prepared with the angles\\n{np.round(param,3)}\\n'\n",
    "      f'has the energy per site E={E:.4f} in the TFI on {N} spins with transverse field {t}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that our QAOA ansatz always contains an even number of parameters by design. \n",
    "In case an odd number of parameters is given, the last parameter is discarded.\n",
    "\n",
    "One can skip the arguments $N$ and $t$ which are set to len(param) and 1 by default"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient and Fubini-Study matrix evaluation\n",
    "The gradient of the TFI energy can be computed via\n",
    "* the standard finite difference method\n",
    "* [the parameter shift rule](https://arxiv.org/abs/1811.11184)\n",
    "* matrix elements of the Hamiltonian between the QAOA state and the corresponding derivative state (see below)\n",
    "* [an ancilla qubit expectation value](https://arxiv.org/abs/1804.03023)\n",
    "\n",
    "Our code makes use of the third variant, which is feasible only because all computations decompose into two-dimensional matrix-vector operations and only a small bit of memory is needed to buffer several states. In particular the combination with the Fubini-Study matrix computation, which requires the derivative states as well, makes this a practical and fast choice.  \n",
    "The derivative of the energy expectation value - where we look at $\\theta_j$ exemplarily -\n",
    "\\begin{equation}E(\\vec{\\theta},\\vec{\\varphi})=\\langle\\psi(\\vec{\\theta},\\vec{\\varphi})|H_{TFI}|\\psi(\\vec{\\theta},\\vec{\\varphi})\\rangle\\end{equation}\n",
    "w.r.t $\\theta_j$ is given by\n",
    "\\begin{equation}\\partial_{\\theta_j}E(\\vec{\\theta},\\vec{\\varphi})=2\\mathfrak{Re}\\langle\\psi(\\vec{\\theta},\\vec{\\varphi})|H|\\partial_{\\theta_j}\\psi(\\vec{\\theta},\\vec{\\varphi})\\rangle\\end{equation}\n",
    "which in the free fermion picture with $r=\\lfloor N/2\\rfloor$ fermions decomposes into \n",
    "\\begin{align}\n",
    "\\partial_{\\theta_j}E(\\vec{\\theta},\\vec{\\varphi})&=\\sum_{q=1}^r 2\\mathfrak{Re}\n",
    "\\langle \\psi^{(q)}| H^{(q)} |\\partial_{\\theta_j}\\psi^{(q)}\\rangle \\\\\n",
    "|\\partial_{\\theta_j}\\psi^{(q)}\\rangle &= \\left(\\prod_{k=N}^j R_{Z}^{(q)}(\\varphi_k)R_{ax}^{(q)}(\\theta_k)\\right)W^{(q)}\\left(\\prod_{k=j-1}^1 R_{Z}^{(q)}(\\varphi_k)R_{ax}^{(q)}(\\theta_k)\\right)|\\bar{\\psi}\\rangle\n",
    "\\end{align}\n",
    "\n",
    "with $W^{(q)}=\\cos \\alpha_q Z + \\sin\\alpha_q Y$ and some fermion-specific angles $\\alpha_q$. This means that while evolving the $q$-th fermion state with the QAOA unitaries, we can split off a copy at the $2j$-th unitary in order to apply the respective $W^{(q)}$ (for $\\partial_{\\theta_j}$) or $Z$ (for $\\partial_{\\varphi_j}$) operator to it and continue evolving this state with the remaining QAOA unitaries in parallel to the original, derivative-free state.\n",
    "\n",
    "Evidently, we need these states to compute the Fubini-Study matrix ${}^{[1]}$\n",
    "\\begin{equation}\n",
    "F_{ij} = \\mathfrak{Re}\\left\\{\\langle\\partial_i\\psi|\\partial_j\\psi\\rangle\\right\\}-\\langle\\partial_i\\psi|\\psi\\rangle\\langle\\psi|\\partial_j\\psi\\rangle\n",
    "\\end{equation}\n",
    "such that the gradient uses resources that are required by $F$ in any case but computes some matrix elements of $H$ with them rather than the overlaps in F, which are independent of the Hamiltonian.\n",
    "\n",
    "This also is the reason why the gradient and Fubini-Study matrix functionality is implemented as one method, namely _eval\\_all_, whereas in scenarios requiring the gradient but not $F$, it will be more efficient to use a finite difference method or parameter shift rule of your choice.\n",
    "\n",
    "\n",
    "[1] For convenience, relabel $\\partial_j=\\partial_{\\theta_{j/2}}$ for even $j$ and $\\partial_j=\\partial_{\\varphi_{(j-1)/2}}$ for odd $j$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The energy of the state prepared with the angles\n",
      "[ 0.794  0.07   1.105  0.759  0.767  0.054  0.701 -1.394  1.263  0.976]\n",
      "has the energy per site E=-0.2743 in the TFI on 10 spins with transverse field 1.2.\n",
      "\n",
      "The gradient is [-0.065 -0.049 -0.123  0.987  0.026  0.448  0.056  0.372  1.302  0.026]\n",
      "and the Fubini-Study matrix reads\n",
      "[[ 2.5   -0.     2.482  0.07   1.408  1.689  1.368  1.016 -1.463  0.288]\n",
      " [-0.     3.145 -0.083 -0.388  2.256  0.095  2.184  0.853  0.767 -0.834]\n",
      " [ 2.482 -0.083  2.476 -0.091  1.353  1.602  1.305  1.005 -1.502  0.344]\n",
      " [ 0.07  -0.388 -0.091  2.919 -0.558  1.325 -0.38  -0.357  0.352 -0.445]\n",
      " [ 1.408  2.256  1.353 -0.558  3.171 -0.022  3.099  0.764 -0.222 -0.807]\n",
      " [ 1.689  0.095  1.602  1.325 -0.022  3.313 -0.057  2.204 -0.819  0.833]\n",
      " [ 1.368  2.184  1.305 -0.38   3.099 -0.057  3.046  0.572 -0.189 -0.881]\n",
      " [ 1.016  0.853  1.005 -0.357  0.764  2.204  0.572  3.295 -0.492  0.99 ]\n",
      " [-1.463  0.767 -1.502  0.352 -0.222 -0.819 -0.189 -0.492  4.151  0.462]\n",
      " [ 0.288 -0.834  0.344 -0.445 -0.807  0.833 -0.881  0.99   0.462  3.942]]\n"
     ]
    }
   ],
   "source": [
    "# As before specify number of qubits N, number of blocks p and transverse field\n",
    "N = 10; p=5; t=1.2;\n",
    "# Generate random parameters as before \n",
    "param = np.random.random(2*p)*np.pi-np.pi/2\n",
    "# The entire computation is obscured in one method here because the computation of the gradient and F \n",
    "# are more efficient when performed together.\n",
    "E, grad, F = tfi_ff.eval_all(param, N, t)\n",
    "print(f'The energy of the state prepared with the angles\\n{np.round(param,3)}\\n'\n",
    "      f'has the energy per site E={E:.4f} in the TFI on {N} spins with transverse field {t}.\\n\\n'\n",
    "      f'The gradient is {np.round(grad,3)}\\n'\n",
    "      f'and the Fubini-Study matrix reads\\n{np.round(F,3)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For N=10 qubits and p= 5 QAOA blocks the computation took 0.01 seconds.\n",
      "For N=20 qubits and p=10 QAOA blocks the computation took 0.07 seconds.\n",
      "For N=30 qubits and p=15 QAOA blocks the computation took 0.16 seconds.\n",
      "For N=40 qubits and p=20 QAOA blocks the computation took 0.43 seconds.\n",
      "For N=50 qubits and p=25 QAOA blocks the computation took 1.01 seconds.\n",
      "For N=60 qubits and p=30 QAOA blocks the computation took 1.78 seconds.\n"
     ]
    }
   ],
   "source": [
    "# We restricted ourselves to 10 qubits here for readability of the output, but regarding the computation time, \n",
    "# there are not strong constraints as the algorithm scales polynomially in the system size and the number of \n",
    "# parameters, with very very small prefactors:\n",
    "import time\n",
    "\n",
    "# Choose some system sizes and fix p=N/2\n",
    "for N in 10*np.arange(1,7):\n",
    "    start = time.process_time()\n",
    "    E, g, F = tfi_ff.eval_all(np.random.random(N)*np.pi-np.pi/2)\n",
    "    end = time.process_time()\n",
    "    print(f'For N={N} qubits and p={N//2:2} QAOA blocks the computation took {end-start:.2f} seconds.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Possible improvements\n",
    "The code is fast because the problem is easy. There will be room for improvement for the unlikely case that a faster performance is required. Some ideas for improvement are:\n",
    "* pre-compute the QAOA blocks as derivative generators W^q can be pulled out of these blocks such that all time evolutions only require the blocks and not the separate layers.\n",
    "* the blocks in the previous point could additionally be precomputed analytically, saving many of the matrix multiplications.\n",
    "* port the code to C/C++, which should be easily doable as it only contains matrix-vector and vector-vector multiplications.\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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}