AnsatzCapabilities.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Checking expressibility of ansatze\n",
    "\n",
    "$\\newcommand{\\ket}[1]{\\lvert #1 \\rangle}$\n",
    "$\\newcommand{\\bra}[1]{\\langle #1 \\rvert}$\n",
    "$\\newcommand{\\braket}[1]{\\langle #1 \\rangle}$\n",
    "\n",
    "[1905.10876](https://arxiv.org/abs/1905.10876) is much better on quantifying the capabilities of ansatze, but I didn't understand most of it. Instead I'll use the following much simpler quantity that appears in the paper: \n",
    "$$A = \\int_{Haar} \\ket{\\psi} \\bra{\\psi} d\\psi - \\int_\\theta U(\\theta) \\ket{s} \\bra{s} U^\\dagger(\\theta) d\\theta$$\n",
    "where $\\ket{s}$ means the starting state of our ansatz and $U$ is the ansatz. \n",
    "\n",
    "The Haar integral is a way to take the integral over a group like the unitaries. Instead, we'll approximate it with Monte-Carlo Numerical Integration (guide on [Jarrod McClean's blog](https://jarrodmcclean.com/integrating-over-the-unitary-group/)). \n",
    "\n",
    "Finally, we'll define the expressibility of an ansatz as the Frobenius norm of $A$. The lower it is, the more expressible the ansatz is. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "\n",
    "import cirq\n",
    "import openfermioncirq\n",
    "\n",
    "from tools import utils"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define Hubbard and get ansatz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normal stuff for Hubbard..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermion.utils import HubbardSquareLattice\n",
    "# HubbardSquareLattice parameters\n",
    "x_n = 2\n",
    "y_n = 2\n",
    "n_dofs = 1 # 1 degree of freedom for spin, this might be wrong. Having only one dof means ordered=False. \n",
    "periodic = 0 # Not sure what this is, periodic boundary conditions?\n",
    "spinless = 0 # Has spin\n",
    "\n",
    "lattice = HubbardSquareLattice(x_n, y_n, n_dofs=n_dofs, periodic=periodic, spinless=spinless)\n",
    "\n",
    "from openfermion.hamiltonians import FermiHubbardModel\n",
    "from openfermion.utils import SpinPairs\n",
    "tunneling = [('neighbor', (0, 0), 1.)] # Not sure if this is right\n",
    "interaction = [('onsite', (0, 0), 2., SpinPairs.DIFF)] # Not sure if this is right\n",
    "potential = [(0, 1.)]\n",
    "mag_field = 0. \n",
    "particle_hole_sym = False # Not sure if this is right\n",
    "\n",
    "hubbard = FermiHubbardModel(lattice , tunneling_parameters=tunneling, interaction_parameters=interaction, \n",
    "                            potential_parameters=potential, magnetic_field=mag_field, \n",
    "                            particle_hole_symmetry=particle_hole_sym)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(256, 256)\n",
      "Ground state energy:  -6.828427124746191\n"
     ]
    }
   ],
   "source": [
    "from openfermion import get_sparse_operator, get_ground_state\n",
    "hub_sparse = get_sparse_operator(hubbard.hamiltonian())\n",
    "print(hub_sparse.shape)\n",
    "genergy, gstate = get_ground_state(hub_sparse)\n",
    "print(\"Ground state energy: \", genergy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermioncirq import SwapNetworkTrotterHubbardAnsatz\n",
    "\n",
    "ansatz = SwapNetworkTrotterHubbardAnsatz(x_n, y_n, 1., 2., periodic=False, iterations=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "starting_state = gstate # CHANGE THIS\n",
    "# starting_state = np.array([1] + [0 for i in range(255)])\n",
    "M = 500 \n",
    "N = len(starting_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate the first integral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_unitary():\n",
    "    \"\"\"Generate uniformly random unitary from U(N). \n",
    "       Credit to Jarrod McClean\"\"\"\n",
    "    Z = np.random.randn(N, N) + 1.0j * np.random.randn(N, N)\n",
    "    [Q, R] = sp.linalg.qr(Z)\n",
    "    D = np.diag(np.diagonal(R) / np.abs(np.diagonal(R)))\n",
    "    return np.dot(Q, D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_unitary_density = 0\n",
    "\n",
    "for i in range(M):\n",
    "    psi = random_unitary() @ starting_state \n",
    "#     random_unitary_density += utils.adjoint([psi]) @ [psi] # Sorry for weird formatting\n",
    "    random_unitary_density += psi\n",
    "random_unitary_density /= M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.06371926278274356"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Testing to make sure this works as expected\n",
    "total_2 = 0\n",
    "\n",
    "for i in range(M):\n",
    "    psi = random_unitary() @ starting_state \n",
    "#     total_2 += utils.adjoint([psi]) @ [psi]\n",
    "    total_2 += psi\n",
    "total_2 /= M\n",
    "\n",
    "# Should be low\n",
    "np.linalg.norm(random_unitary_density - total_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uniformly choose parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "ansatz_density = 0\n",
    "for i in range(M): \n",
    "    params = np.random.uniform(-2, 2, 15)\n",
    "    circuit = cirq.resolve_parameters(ansatz.circuit, ansatz.param_resolver(params))\n",
    "    final_state = circuit.final_wavefunction(starting_state)\n",
    "#     density = cirq.density_matrix_from_state_vector(final_state)\n",
    "#     ansatz_density += density \n",
    "    ansatz_density += final_state\n",
    "ansatz_density /= M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.190856809406871"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(random_unitary_density - ansatz_density)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Identity ansatz\n",
    "\n",
    "It's hard to contextualize this result, so I'll try an ansatz which is just the identity $I$. This is a terrible ansatz and the score should reflect that. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999073322440026"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i_density = 0\n",
    "\n",
    "for i in range(M):\n",
    "#     i_density += utils.adjoint([starting_state]) @ [starting_state]\n",
    "    i_density += starting_state\n",
    "i_density /= M\n",
    "\n",
    "np.linalg.norm(random_unitary_density - i_density)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "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
}