Quantum_physics/02_Time_evolution_and_ground_states.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Time evolution and ground states with matrix product states (MPS)\n",
    "\n",
    "In this notebook, we define a nearest-neighbour Hamiltonian on a finite, linear chain. We compute time evolution under the Hamiltonian and the ground state of the Hamiltonian.\n",
    "\n",
    "Refer e.g. to https://arxiv.org/abs/1008.3477 for an introduction and for original references describing the algorithms used below.\n",
    "\n",
    "Related keywords:\n",
    "\n",
    "- density matrix renormalization group (DMRG)\n",
    "\n",
    "- time-dependent density matrix renormalization group (t-DMRG)\n",
    "\n",
    "- time-evolving block decimation (TEBD)\n",
    "\n",
    "- Variational matrix product state (MPS) algorithms\n",
    "\n",
    "- Tensor trains (TT)\n",
    "\n",
    "The code below uses Python and `mpnum` (http://mpnum.readthedocs.io/). This notebook is part of https://github.com/milan-hl/mpnum-examples/.\n",
    "\n",
    "Below, very basic algorithms are used. Refer to the literature for better available algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definitions\n",
    "\n",
    "Open ``init.ipynb`` to see which MPSs and MPOs it defines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run init.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use qubits - local dimension 2\n",
    "ldim = 2\n",
    "# Run on 6 qubits\n",
    "n_sites = 6\n",
    "width = 2  # nearest-neighbour\n",
    "# Parameter of the Hamiltonian\n",
    "h = .3\n",
    "# Evolution time\n",
    "t = 0.9\n",
    "steps = 5\n",
    "tau = t / steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the following nearest-neighbour Hamiltonian of $n$ qubits:\n",
    "\n",
    "\\begin{align}\n",
    "H &= \\sum_{i=0}^{n-2} h_{i,i+1} \\\\\n",
    "h_{i, i+1} &= \\sigma_x^i \\sigma_x^{i+1} + \\sigma_y^i \\sigma_y^{i+1} + h \\sigma_z^i +  \\delta_{i,n-1} h \\sigma_z^{i+1}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_term = mp.chain((mpx, mpx)) + mp.chain((mpy, mpy)) + h * mp.chain((mpz, mid))\n",
    "h_terms = [h_term.copy() for _ in range(n_sites - 1)]\n",
    "h_terms[-1] += h * mp.chain((mid, mpz))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numerical compression (without modifying the Hamiltonian terms) does not reduce the bond dimension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(3,), (3,), (3,), (3,), (4,)]\n",
      "[(3,), (3,), (3,), (3,), (4,)]\n"
     ]
    }
   ],
   "source": [
    "print([t.ranks for t in h_terms])\n",
    "[t.compress(method='svd', relerr=1e-10) for t in h_terms]\n",
    "print([t.ranks for t in h_terms])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time evolution with sparse matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert local terms (MPOs on two qubits) to full matrices on two qubits\n",
    "h_terms_arr = [t.to_array_global().reshape([ldim**width] * 2) for t in h_terms]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = 0\n",
    "\n",
    "for pos, term in enumerate(h_terms_arr):\n",
    "    H += spa.kron(\n",
    "        spa.eye(ldim**pos), \n",
    "        spa.kron(\n",
    "            term,\n",
    "            spa.eye(ldim**(n_sites - width - pos))\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "psi_0 = mp.chain(([mup, mdown] * ((n_sites + 1) // 2))[:n_sites])\n",
    "psi_0_arr = psi_0.to_array().ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Norm of the evolved state: 1.0\n"
     ]
    }
   ],
   "source": [
    "psi_t_spa = spa.linalg.expm_multiply(-1j * H * t, psi_0_arr)\n",
    "print('Norm of the evolved state:', np.linalg.norm(psi_t_spa))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Algorithm used by `expm_multiply()` is referenced at https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.expm_multiply.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time evolution with MPS\n",
    "\n",
    "We do a first order trotter approximation:\n",
    "$$\n",
    "\\exp(-i H t) \\approx [ \\exp(-i H_o \\tau) \\exp(-i H_e \\tau) ]^M, \\quad \\tau = t / M\n",
    "$$\n",
    "where\n",
    "$$\n",
    "H = H_e + H_o, \\quad\n",
    "H_e = \\sum_{i=0}^{n/2} h_{2i,2i+1}, \\quad\n",
    "H_o = \\sum_{i=0}^{n/2} h_{2i+1, 2i+2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_evos = [\n",
    "    # Convert results from matrix to MPO (on 2 qubits)\n",
    "    mp.MPArray.from_array_global(\n",
    "        # Compute exp(-i h_{i,i+1} t) as full matrix (on 2 qubits)\n",
    "        sp.linalg.expm(\n",
    "            (-1j * tau)\n",
    "            # Convert MPO on 2 qubits to square matrix\n",
    "            * term.to_array_global().reshape([ldim**width] * 2)\n",
    "        )\n",
    "        .reshape([ldim] * (2 * width)),\n",
    "        ndims=2\n",
    "    )\n",
    "    for term in h_terms\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_evo_even = mp.chain(h_evos[0::2])\n",
    "h_evo_odd = mp.chain(h_evos[1::2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Norm of the time-evolved state: 1.0\n",
      "Overlap with solution from sparse matrix: 0.959732582153\n"
     ]
    }
   ],
   "source": [
    "psi = psi_0\n",
    "for i in range(steps):\n",
    "    psi = mp.partialdot(h_evo_even, psi, 0)\n",
    "    psi = mp.partialdot(h_evo_odd, psi, 1)\n",
    "    psi.compress(method='svd', relerr=1e-10)\n",
    "    \n",
    "psi_t_arr = psi.to_array().ravel()\n",
    "print('Norm of the time-evolved state:', np.linalg.norm(psi_t_arr))\n",
    "print('Overlap with solution from sparse matrix:', abs(np.vdot(psi_t_arr, psi_t_spa)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notes:\n",
    "- Only 5 Trotter steps cause small overlap\n",
    "\n",
    "- For larger numbers of sites, a larger compression error needs to be tolerated (otherwise, bond dimension and memory usage will grow)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ground state with sparse matrices\n",
    "\n",
    "Compute the eigenvalue which has the smallest real part and a corresponding eigenvector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigval_spa, eigvec_spa = spa.linalg.eigsh(H, k=1, which='SR')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Algorithm used by `eigsh()` is referenced at https://docs.scipy.org/doc/scipy/reference/generated/generated/scipy.sparse.linalg.eigsh.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ground state with matrix products states (MPS) and operators (MPO)\n",
    "\n",
    "The following code performs a variational search for eigenvectors of a matrix product operator (MPO) Hamiltonian whose eigenvalue has a small real part, with the search restricted to matrix product states (MPS) of a certain bond dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalue difference: [ 0.01046594]\n",
      "Eigenvector overlap: 0.99901663645\n",
      "Relative variance: 0.00126286092335\n"
     ]
    }
   ],
   "source": [
    "# Convert local terms to MPO\n",
    "H_mpo = mp.local_sum(h_terms)\n",
    "\n",
    "# Look for the eigenvalue with smallest real part\n",
    "import functools as ft\n",
    "eigs_fun = ft.partial(spa.linalg.eigsh, k=1, which='SR')\n",
    "\n",
    "eigval, eigvec = mp.linalg.eig(H_mpo, num_sweeps=5, startvec_rank=4, eigs=eigs_fun)\n",
    "eigvec_arr = eigvec.to_array().ravel()\n",
    "\n",
    "# Todo: Verify that discarded imaginary part is small\n",
    "Hpsi = mp.dot(H_mpo, eigvec)\n",
    "ept_Hsq = mp.norm(Hpsi)**2\n",
    "ept_H = mp.inner(eigvec, Hpsi).real\n",
    "rel_var = (ept_Hsq - ept_H**2) / ept_H**2\n",
    "\n",
    "print('Eigenvalue difference:', eigval - eigval_spa)\n",
    "print('Eigenvector overlap:', abs(np.vdot(eigvec_spa, eigvec_arr)))\n",
    "print('Relative variance:', rel_var)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notes:\n",
    "\n",
    "- **No guarantee** that the ground state has been found\n",
    "- Number of iterations (`num_sweeps`) must be chosen be the user\n",
    "- For possible convergence criteria, refer to the literature (e.g. https://arxiv.org/abs/1008.3477)\n",
    "- To verify whether the result is an eigenvector, compute e.g. $(\\langle H^2 \\rangle - \\langle H \\rangle^2) / \\langle H \\rangle^2$\n",
    "- Larger bond dimensions (`startvec_rank`) enables more precise results\n",
    "- Documentation is at http://mpnum.readthedocs.io/en/latest/mpnum.html#mpnum.linalg.eig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}