{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparing the ground state energies obtained by density matrix renormalization group, exact diagonalization, and an SDP hierarchy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would like to compare the ground state energy of the following spinless fermionic system [1]:\n",
    "\n",
    "$H_{\\mathrm{free}}=\\sum_{<rs>}\\left[c_{r}^{\\dagger} c_{s}+c_{s}^{\\dagger} c_{r}-\\gamma(c_{r}^{\\dagger} c_{s}^{\\dagger}+c_{s}c_{r} )\\right]-2\\lambda\\sum_{r}c_{r}^{\\dagger}c_{r},$\n",
    "\n",
    "where $<rs>$ goes through nearest neighbour pairs in a two-dimensional lattice. The fermionic operators are subject to the following constraints:\n",
    "\n",
    "$\\{c_{r}, c_{s}^{\\dagger}\\}=\\delta_{rs}I_{r}$\n",
    "\n",
    "$\\{c_r^\\dagger, c_s^\\dagger\\}=0,$\n",
    "\n",
    "$\\{c_{r}, c_{s}\\}=0.$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our primary goal is to benchmark the SDP hierarchy of Reference [2]. The baseline methods are density matrix renormalization group (DMRG) and exact diagonalization (ED), both of which are included in Algorithms and Libraries for Physics Simulations (ALPS, [3]). The range of predefined Hamiltonians is limited, so we simplify the equation by setting $\\gamma=0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prerequisites"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To run this notebook, [ALPS](http://alps.comp-phys.org/mediawiki/index.php/Main_Page), [Sympy](http://sympy.org/), [Scipy](http://scipy.org/), and [SDPA](http://sdpa.sourceforge.net/) must be installed. A recent version of [Ncpol2sdpa](https://pypi.python.org/pypi/ncpol2sdpa) is also necessary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating the ground state energy with DMRG and ED"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DMRG and ED are included in ALPS. To start the calculations, we need to import the Python interface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pyalps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For now, we are only interested in relatively small systems, we will try lattice sizes between $2\\times 2$ and $5\\times 5$. With this, we set the parameters for DMRG and ED:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lattice_range = [2, 3, 4, 5]\n",
    "parms = [{ \n",
    "  'LATTICE'        : \"open square lattice\",  # Set up the lattice\n",
    "  'MODEL'          : \"spinless fermions\",    # Select the model \n",
    "  'L'              : L,                      # Lattice dimension\n",
    "  't'              : -1 ,                    # This and the following\n",
    "  'mu'             : 2,                      # are parameters to the\n",
    "  'U'              : 0 ,                     # Hamiltonian.\n",
    "  'V'              : 0,\n",
    "  'Nmax'           : 2 ,                     # These parameters are\n",
    "  'SWEEPS'         : 20,                      # specific to the DMRG\n",
    "  'MAXSTATES'      : 300,                    # solver.\n",
    "  'NUMBER_EIGENVALUES' : 1,          \n",
    "  'MEASURE_ENERGY' : 1\n",
    "} for L in lattice_range ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will need a helper function to extract the ground state energy from the solutions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def extract_ground_state_energies(data):\n",
    "    E0 = []\n",
    "    for Lsets in data:\n",
    "        allE = []\n",
    "        for q in pyalps.flatten(Lsets):\n",
    "            allE.append(q.y[0])\n",
    "        E0.append(allE[0])\n",
    "    return sorted(E0, reverse=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We invoke the solvers and extract the ground state energies from the solutions. First we use exact diagonalization, which, unfortunately does not scale beyond a lattice size of $4\\times 4$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sparsediag comparison_sparse.in.xml\n"
     ]
    }
   ],
   "source": [
    "prefix_sparse = 'comparison_sparse'\n",
    "input_file_sparse = pyalps.writeInputFiles(prefix_sparse, parms[:-1])\n",
    "\n",
    "res = pyalps.runApplication('sparsediag', input_file_sparse)\n",
    "sparsediag_data = pyalps.loadEigenstateMeasurements(\n",
    "                     pyalps.getResultFiles(prefix=prefix_sparse)) \n",
    "\n",
    "sparsediag_ground_state_energy = extract_ground_state_energies(sparsediag_data)\n",
    "sparsediag_ground_state_energy.append(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DMRG scales to all the lattice sizes we want:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dmrg comparison_dmrg.in.xml\n"
     ]
    }
   ],
   "source": [
    "prefix_dmrg = 'comparison_dmrg'\n",
    "input_file_dmrg = pyalps.writeInputFiles(prefix_dmrg, parms)\n",
    "res = pyalps.runApplication('dmrg',input_file_dmrg)\n",
    "dmrg_data = pyalps.loadEigenstateMeasurements(\n",
    "                  pyalps.getResultFiles(prefix=prefix_dmrg)) \n",
    "dmrg_ground_state_energy = extract_ground_state_energies(dmrg_data)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating the ground state energy with SDP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ground state energy problem can be rephrased as a polynomial optimiziation problem of noncommuting variables. We use Ncpol2sdpa to translate this optimization problem to a sparse SDP relaxation [4]. The relaxation is solved with SDPA, a high-performance SDP solver that deals with sparse problems efficiently [5]. First we need to import a few more functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sympy.physics.quantum.dagger import Dagger\n",
    "from ncpol2sdpa import SdpRelaxation, generate_operators, \\\n",
    "                       fermionic_constraints, get_neighbors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set the additional parameters for this formulation, including the order of the relaxation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "level = 1\n",
    "gam, lam = 0, 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we iterate over the lattice range, defining a new Hamiltonian and new constraints in each step:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sdp_ground_state_energy = []\n",
    "for lattice_dimension in lattice_range:\n",
    "    n_vars = lattice_dimension * lattice_dimension\n",
    "    C = generate_operators('C%s' % (lattice_dimension), n_vars)\n",
    "    \n",
    "    hamiltonian = 0\n",
    "    for r in range(n_vars):\n",
    "        hamiltonian -= 2*lam*Dagger(C[r])*C[r]\n",
    "        for s in get_neighbors(r, lattice_dimension):\n",
    "            hamiltonian += Dagger(C[r])*C[s] + Dagger(C[s])*C[r]\n",
    "            hamiltonian -= gam*(Dagger(C[r])*Dagger(C[s]) + C[s]*C[r])\n",
    "    \n",
    "    substitutions = fermionic_constraints(C)\n",
    "        \n",
    "    sdpRelaxation = SdpRelaxation(C)\n",
    "    sdpRelaxation.get_relaxation(level, objective=hamiltonian, substitutions=substitutions)\n",
    "    sdpRelaxation.solve()\n",
    "    sdp_ground_state_energy.append(sdpRelaxation.primal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The level-one relaxation matches the ground state energy given by DMRG and ED."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             2          3          4          5\n",
      "DMRG -8.000000 -18.828427 -33.708204 -52.928203\n",
      "  ED -8.000000 -18.828427 -33.708204 0.000000\n",
      " SDP -8.000000 -18.828427 -33.708204 -52.928207\n"
     ]
    }
   ],
   "source": [
    "data = [dmrg_ground_state_energy,\\\n",
    "        sparsediag_ground_state_energy,\\\n",
    "        sdp_ground_state_energy]\n",
    "labels = [\"DMRG\", \"ED\", \"SDP\"]\n",
    "print (\"{:>4} {:>9} {:>10} {:>10} {:>10}\").format(\"\", *lattice_range)\n",
    "for label, row in zip(labels, data):\n",
    "    print (\"{:>4} {:>7.6f} {:>7.6f} {:>7.6f} {:>7.6f}\").format(label, *row)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] Corboz, P.; Evenbly, G.; Verstraete, F. & Vidal, G. [Simulation of interacting fermions with entanglement renormalization](http://arxiv.org/abs/0904.4151). _Physics Review A_, 2010, 81, pp. 010303.\n",
    "\n",
    "[2] Pironio, S.; Navascués, M. & Acín, A. [Convergent relaxations of polynomial optimization problems with noncommuting variables](http://arxiv.org/abs/0903.4368). _SIAM Journal on Optimization_, 2010, 20, pp. 2157-2180.\n",
    "\n",
    "[3] Bauer, B.; Carr, L.; Evertz, H.; Feiguin, A.; Freire, J.; Fuchs, S.; Gamper, L.; Gukelberger, J.; Gull, E.; Guertler, S. & others. [The ALPS project release 2.0: Open source software for strongly correlated systems](http://arxiv.org/abs/1101.2646). _Journal of Statistical Mechanics: Theory and Experiment_, IOP Publishing, 2011, 2011, P05001.\n",
    "\n",
    "[4] Wittek, P. [Ncpol2sdpa -- Sparse Semidefinite Programming Relaxations for Polynomial Optimization Problems of Noncommuting Variables](http://arxiv.org/abs/1308.6029). _arXiv:1308.6029_, 2013.\n",
    "\n",
    "[5] Yamashita, M.; Fujisawa, K. & Kojima, M. [Implementation and evaluation of SDPA 6.0 (semidefinite programming algorithm 6.0)](http://dx.doi.org/10.1080/1055678031000118482). _Optimization Methods and Software_, 2003, 18, 491-505."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python2",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}