examples/mpnum_intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introductory Notebook to mpnum\n",
    "\n",
    "mpnum implements *matrix product arrays* (MPA), which are efficient parameterizations of certain multi-partite arrays. Special cases of the MPA structure, which are omnipresent in many-body quantum physics, are *matrix product states* (MPS) and *matrix product operators* (MPO) with one and two array indices per site, respectively. In the applied math community, matrix product states are also known as *tensor trains* (TT). \n",
    "\n",
    "The main class implementing an MPA with arbitrary number of array indices (or \"physical legs\") is ``mpnum.MPArray``. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "\n",
    "import mpnum as mp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MPA and MPS basics\n",
    "\n",
    "A convenient example to deal with is a random MPA. First, we create a fixed seed, then a random MPA:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rng = np.random.RandomState(seed=42)\n",
    "mpa = mp.random_mpa(sites=4, ldim=2, rank=3, randstate=rng, normalized=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The MPA is an instance of the MPArray class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<mpnum.mparray.MPArray at 0x7f65672eaa90>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of sites:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(mpa)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Number of physical legs at each site (=number of array indices at each site):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 1, 1, 1)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpa.ndims"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the MPA has one physical leg per site, we have created a matrix product *state* (i.e. a tensor train). In the [graphical notation](http://mpnum.readthedocs.io/en/latest/intro.html#graphical-notation), this MPS looks like this\n",
    "\n",
    "![](http://mpnum.readthedocs.io/en/latest/_images/tensors_mps_no_names.png)\n",
    "\n",
    "Note that `mpnum` internally stores the local tensors of the matrix product representation on the right hand side. We see below how to obtain the \"dense\" tensor from an `MPArray`.\n",
    "Dimension of each physical leg:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((2,), (2,), (2,), (2,))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpa.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the number and dimension of the physical legs at each site can differ (altough this is rarely used in practice). \n",
    "\n",
    "Representation ranks (aka compression ranks) between each pair of sites:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 3, 2)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpa.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In physics, the representation ranks are usually called the *bond dimensions* of the representation.\n",
    "\n",
    "Dummy bonds before and after the chain are omitted in ``mpa.ranks``. (Currently, mpnum only implements open boundary conditions.)\n",
    "\n",
    "Above, we have specified ``normalized=True``.  Therefore, we have created an MPA with $\\ell_2$-norm 1.  In case the MPA does not represent a vector but has more physical legs, it is nonetheless treated as a vector. Hence, for operators <code>mp.norm</code> implements the Frobenius norm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0000000000000002"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mp.norm(mpa)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert to a dense array, which should be used with care due because the memory used increases exponentially with the number of sites:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 2, 2, 2)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr = mpa.to_array()\n",
    "arr.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting full array has one index for each physical leg.\n",
    "\n",
    "Now convert the full array back to an MPA:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpa2 = mp.MPArray.from_array(arr)\n",
    "len(mpa2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have obtained an MPA with length 1.  This is not what we expected.  The reason is that by default, all legs are placed on a single site (also notice the difference between ``mpa2.shape`` here and ``mpa.shape`` from above): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((2, 2, 2, 2),)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpa2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((2,), (2,), (2,), (2,))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpa.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We obtain the desired result by specifying the number of legs per site we want:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpa2 = mp.MPArray.from_array(arr, ndims=1)\n",
    "len(mpa2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can compute the norm distance between the two MPAs.  (Again, the Frobenius norm is used.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.2998268912398721e-16"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mp.norm(mpa - mpa2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this is an often used operation and allows for additional optimization (not implemented currently), it is advisable to use the specific <code>mp.normdist</code> for this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.2998268912398721e-16"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mp.normdist(mpa, mpa2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sums, differences and scalar multiplication of MPAs is done with the normal operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.0000000000000009"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mp.norm(3 * mpa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5000000000000011"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mp.norm(mpa + 0.5 * mpa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.50000000000000133"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mp.norm(mpa - 1.5 * mpa)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiplication with a scalar leaves the bond dimension unchanged:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 3, 2)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpa.ranks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 3, 2)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(3 * mpa).ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bond dimensions of a sum (or difference) are given by the sums of the bond dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 2, 2)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpa2 = mp.random_mpa(sites=4, ldim=2, rank=2, randstate=rng)\n",
    "mpa2.ranks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4, 5, 4)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(mpa + mpa2).ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MPO basics\n",
    "\n",
    "First, we create a random MPA with two physical legs per site:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpo = mp.random_mpa(sites=4, ldim=(3, 2), rank=3, randstate=rng, normalized=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In graphical notation, `mpo` looks like this\n",
    "\n",
    "![](http://mpnum.readthedocs.io/en/latest/_images/notebook_mpo.png)\n",
    "\n",
    "It's basic properties are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[4, (2, 2, 2, 2), (3, 3, 3)]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[len(mpo), mpo.ndims, mpo.ranks]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each site has two physical legs, one with dimension 3 and one with dimension 2.  This corresponds to a non-square full array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((3, 2), (3, 2), (3, 2), (3, 2))"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpo.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now convert the mpo to a full array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 2, 3, 2, 3, 2, 3, 2)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpo_arr = mpo.to_array()\n",
    "mpo_arr.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We refer to this arangement of axes as _local form_, since indices which correspond to the same site are neighboring. This is a natural form for the MPO representation. However, for some operations it is necessary to have row and column indices grouped together -- we refer to this as _global form_:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 3, 3, 3, 2, 2, 2, 2)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from mpnum.utils.array_transforms import local_to_global\n",
    "\n",
    "mpo_arr = mpo.to_array()\n",
    "mpo_arr = local_to_global(mpo_arr, sites=len(mpo))\n",
    "mpo_arr.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives the expected result.  Note that it is crucial to specify the correct number of sites, otherwise we do not get what we want:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 3, 2, 2, 3, 3, 2, 2)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpo_arr = mpo.to_array()\n",
    "mpo_arr = local_to_global(mpo_arr, sites=2)\n",
    "mpo_arr.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an alternative, there is the following shorthand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 3, 3, 3, 2, 2, 2, 2)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpo_arr = mpo.to_array_global()\n",
    "mpo_arr.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An array in global form can be converted into matrix-product form with the following API:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0881840590136613e-15"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpo2 = mp.MPArray.from_array_global(mpo_arr, ndims=2)\n",
    "mp.normdist(mpo, mpo2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MPO-MPS product and arbitrary MPA-MPA products\n",
    "\n",
    "We can now compute the matrix-vector product of ``mpa`` from above (which is an MPS) and ``mpo``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((2,), (2,), (2,), (2,))"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpa.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((3, 2), (3, 2), (3, 2), (3, 2))"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpo.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((3,), (3,), (3,), (3,))"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prod = mp.dot(mpo, mpa, axes=(-1, 0))\n",
    "prod.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a new MPS, with local dimension changed by ``mpo`` and looks like this:\n",
    "\n",
    "![](http://mpnum.readthedocs.io/en/latest/_images/notebook_mpo_mps.png)\n",
    "\n",
    "The ``axes`` argument is optional and defaults to ``axes=(-1, 0)`` -- i.e. contracting, at each site, the last pyhsical index of the first factor with the first physical index of the second factor. More specifically, the ``axes`` argument  specifies which _physical legs_ should be contracted:  ``axes[0]`` specifies the physical in the first argument, and ``axes[1]`` specifies the physical leg in the second argument.  This means that the same product can be achieved with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.7794893594008944e-16"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prod2 = mp.dot(mpa, mpo, axes=(0, 1))\n",
    "mp.normdist(prod, prod2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in any case, the ranks of the output of ``mp.dot`` are the *products* of the original ranks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((3, 3, 3), (2, 3, 2), (6, 9, 6))"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mpo.ranks, mpa.ranks, prod.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we compute the same product using the full arrays ``arr`` and ``mpo_arr``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(81,)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr_vec = arr.ravel()\n",
    "mpo_arr = mpo.to_array_global()\n",
    "mpo_arr_matrix = mpo_arr.reshape((81, 16))\n",
    "prod3_vec = np.dot(mpo_arr_matrix, arr_vec)\n",
    "prod3_vec.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, we need to reshape the result ``prod3_vec`` before we can convert it back to an MPA:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((3,), (3,), (3,), (3,))"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prod3_arr = prod3_vec.reshape((3, 3, 3, 3))\n",
    "prod3 = mp.MPArray.from_array(prod3_arr, ndims=1)\n",
    "prod3.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compare the two results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0433926816958574e-16"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mp.normdist(prod, prod3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compare by converting ``prod`` to a full array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0434960119970279e-16"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prod_arr = prod.to_array()\n",
    "la.norm((prod3_arr - prod_arr).reshape(81))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Converting full operators to MPOs\n",
    "\n",
    "While MPO algorithms avoid using full operators in general, we will need to convert a term acting on only two sites to an MPO in order to continue with MPO operations;  i.e. we will need to convert a full array to an MPO. \n",
    "\n",
    "First, we define a full operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "CZ = np.array([[ 1.,  0.,  0.,  0.],\n",
    "               [ 0.,  1.,  0.,  0.],\n",
    "               [ 0.,  0.,  1.,  0.],\n",
    "               [ 0.,  0.,  0., -1.]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This operator is the so-called *controlled Z* gate: Apply ``Z`` on the second qubit if the first qubit is in state ``e2``.\n",
    "\n",
    "To convert it to an MPO, we have to reshape:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "CZ_arr = CZ.reshape((2, 2, 2, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can create an MPO, being careful to specify the correct number of legs per site:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "CZ_mpo = mp.MPArray.from_array_global(CZ_arr, ndims=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test it, we apply the operator to the state which has both qubits in state ``e2``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 0, 1])"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vec = np.kron([0, 1], [0, 1])\n",
    "vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reshape and convert to an MPS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "vec_arr = vec.reshape([2, 2])\n",
    "mps = mp.MPArray.from_array(vec_arr, ndims=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compute the matrix-vector product:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.,  0.,  0., -1.])"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out = mp.dot(CZ_mpo, mps)\n",
    "out.to_array().ravel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output is as expected:  We have acquired a minus sign. \n",
    "\n",
    "We have to be careful to use ``from_array_global`` and not ``from_array`` for `CZ_mpo`, because the ``CZ_arr`` is in *global form*. Here, all physical legs have the same dimension, so we can use ``from_array`` without error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "CZ_mpo2 = mp.MPArray.from_array(CZ_arr, ndims=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, the result is not what we want:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.,  0.,  0., -1.])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out2 = mp.dot(CZ_mpo2, mps)\n",
    "out2.to_array().ravel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reason is easy to see: We have applied the following matrix to our state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.,  0.,  1.],\n",
       "       [ 0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.],\n",
       "       [ 1.,  0.,  0., -1.]])"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CZ_mpo2.to_array_global().reshape(4, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keep in mind that we have to use ``to_array_global`` before the reshape.  Using ``to_array`` would not provide us the matrix which we have applied to the state with ``mp.dot``.  Instead, it will exactly return the input:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.,  0.,  0.],\n",
       "       [ 0.,  1.,  0.,  0.],\n",
       "       [ 0.,  0.,  1.,  0.],\n",
       "       [ 0.,  0.,  0., -1.]])"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CZ_mpo2.to_array().reshape(4, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, ``from_array_global`` is just the shorthand for the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5700924586837752e-16"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from mpnum.utils.array_transforms import global_to_local\n",
    "\n",
    "CZ_mpo3 = mp.MPArray.from_array(global_to_local(CZ_arr, sites=2), ndims=2)\n",
    "\n",
    "mp.normdist(CZ_mpo, CZ_mpo3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, in the explicit version you must submit both the correct number of sites and the correct number of physical legs per site.\n",
    "Therefore, the function `MPArray.from_array_global` simplifies the conversion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating MPAs from Kronecker products"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is a frequent task to create an MPS which represents the product state of $\\vert 0 \\rangle$ on each qubit.  If the chain is very long, we cannot create the full array with ``np.kron`` and use ``MPArray.from_array`` afterwards because the array would be too large.\n",
    "\n",
    "In the following, we describe how to efficiently construct an MPA representation of a Kronecker product of vectors. The same methods can be used to efficiently construct MPA representations of Kronecker products of operators or tensors with three or more indices.\n",
    "\n",
    "First, we need the state on a single site:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 0])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e1 = np.array([1, 0])\n",
    "e1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can use ``from_kron`` to directly create an MPS representation of the Kronecker product:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 0, 0, 0, 0, 0, 0, 0])"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mps = mp.MPArray.from_kron([e1, e1, e1])\n",
    "mps.to_array().ravel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This works well for large numbers of sites because the needed memory scales linearly with the number of sites:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2000"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mps = mp.MPArray.from_kron([e1] * 2000)\n",
    "len(mps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An even more pythonic solution is the use of iterators in this example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2000"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from itertools import repeat\n",
    "\n",
    "mps = mp.MPArray.from_kron(repeat(e1, 2000))\n",
    "len(mps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do not call ``.to_array()`` on this state!\n",
    "\n",
    "The bond dimension of the state is 1, because it is a product state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 1, ..., 1, 1, 1])"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(mps.ranks)  # Convert to an array for nicer display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also create a single-site MPS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mps1 = mp.MPArray.from_array(e1, ndims=1)\n",
    "len(mps1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After that, we can use `mp.chain` to create Kronecker products of the MPS directly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mps = mp.chain([mps1, mps1, mps1])\n",
    "len(mps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It returns the same result as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 0, 0, 0, 0, 0, 0, 0])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mps.to_array().ravel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use ``mp.chain`` on the three-site MPS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "300"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mps = mp.chain([mps] * 100)\n",
    "len(mps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Note that `mp.chain` interprets the factors in the tensor product as distinct sites. Hence, the factors do not need to be of the same length or even have the same  number of indices. In contrast, there is also `mp.localouter`, which computes the tensor product of MPArrays with the same number of sites:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4, ((2,), (2,), (2,), (2,)))"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mps = mp.chain([mps1] * 4)\n",
    "len(mps), mps.shape, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4, ((2, 2), (2, 2), (2, 2), (2, 2)))"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho = mp.localouter(mps.conj(), mps)\n",
    "len(rho), rho.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Compression\n",
    "\n",
    "A typical matrix product based numerical algorithm performs many additions or multiplications of MPAs.  As mentioned above, both operations increase the rank.  If we let the bond dimension grow, the amount of memory we need grows with the number of operations we perform.  To avoid this problem, we have to find an MPA with a smaller rank which is a good approximation to the original MPA.\n",
    "\n",
    "We start by creating an MPO representation of the identity matrix on 6 sites with local dimension 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "op = mp.eye(sites=6, ldim=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((3, 3), (3, 3), (3, 3), (3, 3), (3, 3), (3, 3))"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "op.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As it is a tensor product operator, it has rank 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 1, 1, 1, 1)"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "op.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, addition increases the rank:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 3, 3, 3, 3)"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "op2 = op + op + op\n",
    "op2.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Matrix multiplication multiplies the individual ranks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9, 9, 9, 9, 9)"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "op3 = mp.dot(op2, op2)\n",
    "op3.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(NB: `compress` or `compression` below can call `canonicalize` on the MPA, which in turn could already reduce the rank to 1 in case the rank can be compressed without error.  Keep that in mind.)\n",
    "\n",
    "Keep in mind that the operator represented by ``op3`` is still the identity operator, i.e. a tensor product operator.  This means that we expect to find a good approximation with low rank easily.  Finding such an approximation is called *compression* and is achieved as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 1, 1, 1, 1)"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "op3 /= mp.norm(op3.copy())  # normalize to make overlap meaningful\n",
    "copy = op3.copy()\n",
    "overlap = copy.compress(method='svd', rank=1)\n",
    "copy.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calling ``compress`` on an MPA replaces the MPA in place with a version with smaller bond dimension.  Overlap gives the absolute value of the (Hilbert-Schmidt) inner product between the original state and the output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.99999999999999911"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "overlap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of in-place compression, we can also obtain a compressed copy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.99999999999999911, (2, 2, 2, 2, 2), (9, 9, 9, 9, 9))"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compr, overlap = op3.compression(method='svd', rank=2)\n",
    "overlap, compr.ranks, op3.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SVD compression can also be told to meet a certain truncation error (see the documentation of ``mp.MPArray.compress`` for details)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.99999999999999911, (1, 1, 1, 1, 1), (9, 9, 9, 9, 9))"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compr, overlap = op3.compression(method='svd', relerr=1e-6)\n",
    "overlap, compr.ranks, op3.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use variational compression instead of SVD compression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1+0j), (2, 2, 2, 2, 2), (9, 9, 9, 9, 9))"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compr, overlap = op3.compression(method='var', rank=2, num_sweeps=10, var_sites=2)\n",
    "# Convert overlap from numpy array with shape () to float for nicer display:\n",
    "overlap = overlap.flat[0]\n",
    "complex(overlap), compr.ranks, op3.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a reminder, it is always advisable to check whether the overlap between the input state and the compression is large enough.  In an involved algorithm, it can be useful to store the compression error at each invocation of compression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## MPO sum of local terms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A frequent task is to compute the MPO representation of a local Hamiltonian, i.e. of an operator of the form\n",
    "$$ H = \\sum_{i=1}^{n-1} h_{i, i+1} $$\n",
    "where $h_{i, i+1}$ acts only on sites $i$ and $i + 1$.  This means that $h_{i, i+1} = \\mathbb 1_{i - 1} \\otimes h'_{i, i+1} \\otimes \\mathbb 1_{n - w + 1}$ where $\\mathbb 1_k$ is the identity matrix on $k$ sites and $w = 2$ is the width of $h'_{i, i+1}$.\n",
    "\n",
    "We show how to obtain an MPO representation of such a Hamiltonian.  First of all, we need to define the local terms.  For simplicity, we choose $h'_{i, i+1} = \\sigma_Z \\otimes \\sigma_Z$ independently of $i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  0.],\n",
       "       [ 0.,  0.]])"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zeros = np.zeros((2, 2))\n",
    "zeros"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.],\n",
       "       [ 0.,  1.]])"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idm = np.eye(2)\n",
    "idm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.],\n",
       "       [ 0., -1.]])"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create a float array instead of an int array to avoid problems later\n",
    "Z = np.diag([1., -1])\n",
    "Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.,  0.,  0.],\n",
       "       [ 0., -1.,  0., -0.],\n",
       "       [ 0.,  0., -1., -0.],\n",
       "       [ 0., -0., -0.,  1.]])"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h = np.kron(Z, Z)\n",
    "h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we have to convert the local term ``h`` to an MPO:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4,)"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h_arr = h.reshape((2, 2, 2, 2))\n",
    "h_mpo = mp.MPArray.from_array_global(h_arr, ndims=2)\n",
    "h_mpo.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "h_mpo has rank 4 even though h is a tensor product.  This is far from optimal.  We improve things as follows:  (We could also compress ``h_mpo``.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1,)"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h_mpo = mp.MPArray.from_kron([Z, Z])\n",
    "h_mpo.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most simple way is to implement the formula from above with MPOs:  First we compute the $h_{i, i+1}$ from the $h'_{i, i+1}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<mpnum.mparray.MPArray at 0x7f6563b08588>,\n",
       " <mpnum.mparray.MPArray at 0x7f6563b084e0>,\n",
       " <mpnum.mparray.MPArray at 0x7f6563b086a0>,\n",
       " <mpnum.mparray.MPArray at 0x7f6563b08710>,\n",
       " <mpnum.mparray.MPArray at 0x7f6563b08630>]"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "width = 2\n",
    "sites = 6\n",
    "local_terms = []\n",
    "\n",
    "for startpos in range(sites - width + 1):\n",
    "    left = [mp.MPArray.from_kron([idm] * startpos)] if startpos > 0 else []\n",
    "    right = [mp.MPArray.from_kron([idm] * (sites - width - startpos))] \\\n",
    "        if sites - width - startpos > 0 else []\n",
    "    h_at_startpos = mp.chain(left + [h_mpo] + right)\n",
    "    local_terms.append(h_at_startpos)\n",
    "\n",
    "local_terms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we compute the sum of all the local terms and check the bond dimension of the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 5, 5, 5, 5)"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = local_terms[0]\n",
    "\n",
    "for local_term in local_terms[1:]:\n",
    "    H += local_term\n",
    "    \n",
    "H.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ranks are explained by the ranks of the local terms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1, 1, 1, 1, 1),\n",
       " (1, 1, 1, 1, 1),\n",
       " (1, 1, 1, 1, 1),\n",
       " (1, 1, 1, 1, 1),\n",
       " (1, 1, 1, 1, 1)]"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[local_term.ranks for local_term in local_terms]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We just have to add the ranks at each position.\n",
    "\n",
    "``mpnum`` provides a function which constructs ``H`` from ``h_mpo``, with an output MPO with smaller rank by taking into account the trivial action on some sites:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 3, 3, 3, 2)"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H2 = mp.local_sum([h_mpo] * (sites - width + 1))\n",
    "H2.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Without additional arguments, ``mp.local_sum()`` just adds the local terms with the first term starting on site 0, the second on site 1 and so on.  In addition, the length of the chain is chosen such that the last site of the chain coincides with the last site of the last local term.  Other constructions can be obtained by prodividing additional arguments.\n",
    "\n",
    "We can check that the two Hamiltonians are equal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.4354640488548389e-15"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mp.normdist(H, H2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, this means that we could just compress ``H``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.99999999999999889"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_comp, overlap = H.compression(method='svd', rank=3)\n",
    "overlap / mp.norm(H)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 3, 3, 3, 3)"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_comp.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also check the minimal bond dimension which can be achieved with SVD compression with small error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.99999999999999933"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_comp, overlap = H.compression(method='svd', relerr=1e-6)\n",
    "overlap / mp.norm(H)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 3, 3, 3, 2)"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_comp.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MPS, MPOs and PMPS\n",
    "\n",
    "We can represent vectors (e.g. pure quantum states) as MPS, we can represent arbitrary matrices as MPO and we can represent positive semidefinite matrices as purifying matrix product states (PMPS).  For mixed quantum states, we can thus choose between the MPO and PMPS representations.\n",
    "\n",
    "As mentioned in the introduction, MPS and MPOs are handled as MPAs with one and two physical legs per site.  In addition, PMPS are handled as MPAs with two physical legs per site, where the first leg is the \"system\" site and the second leg is the corresponding \"ancilla\" site.\n",
    "\n",
    "From MPS and PMPS representations, we can easily obtain MPO representations.  ``mpnum`` provides routines for this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4, 9, 9, 4)"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mps = mp.random_mpa(sites=5, ldim=2, rank=3, normalized=True)\n",
    "mps_mpo = mp.mps_to_mpo(mps)\n",
    "mps_mpo.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the rank of ``mps_mpo`` is the square of the rank of ``mps``.\n",
    "\n",
    "Now we create a PMPS with system site dimension 2 and ancilla site dimension 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((2, 3), (2, 3), (2, 3), (2, 3), (2, 3))"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pmps = mp.random_mpa(sites=5, ldim=(2, 3), rank=3, normalized=True)\n",
    "pmps.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9, 9, 9, 9)"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pmps_mpo = mp.pmps_to_mpo(pmps)\n",
    "pmps_mpo.ranks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the rank is squared, as expected.  We can verify that the first physical leg of each site of ``pmps`` is indeed the system site by checking the shape of ``pmps_mpo``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((2, 2), (2, 2), (2, 2), (2, 2), (2, 2))"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pmps_mpo.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Local reduced states\n",
    "\n",
    "For state tomography applications, we frequently need the local reduced states of an MPS, MPO or PMPS. We provide the following functions for this task:\n",
    "\n",
    "* ``mp.reductions_mps_as_pmps()``: Input: MPS, output: local reductions as PMPS\n",
    "* ``mp.reductions_mps_as_mpo()``: Input: MPS, output: local reductions as MPO\n",
    "* ``mp.reductions_pmps()``: Input: PMPS, output: Local reductions as PMPS\n",
    "* ``mp.reductions_mpo()``: Input: MPO, output: Local reductions as MPO\n",
    "\n",
    "The arguments of all functions are similar, e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reduction starting on site 0\n",
      "bdims: (3, 3)\n",
      "trace: 1.0\n",
      "\n",
      "Reduction starting on site 1\n",
      "bdims: (3, 3)\n",
      "trace: 1.0\n",
      "\n",
      "Reduction starting on site 2\n",
      "bdims: (3, 3)\n",
      "trace: 1.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "width = 3\n",
    "startsites = range(len(pmps) - width + 1)\n",
    "for startsite, red in zip(startsites, mp.reductions_pmps(pmps, width, startsites)):\n",
    "    print('Reduction starting on site', startsite)\n",
    "    print('bdims:', red.ranks)\n",
    "    red_mpo = mp.pmps_to_mpo(red)\n",
    "    print('trace:', mp.trace(red_mpo))\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because ``pmps`` was a normalized state, the trace of the reduced states is close to 1.\n",
    "\n",
    "You can omit the ``startsites`` argument:  The default behaviour is the first reductions starting on site 0, the second on site 1, and so on (which is just what we have requested).  The functions for reduced states can also compute different constructions by providing different arguments not described here. "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python [default]",
   "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.3"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}