meeting_05.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework 1 solution\n",
    "\n",
    "Let's use the following convention for numbering legs:\n",
    "```\n",
    " 1--A--3\n",
    "    |\n",
    "    2\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using TensorOperations\n",
    "using LinearMaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'rand_UMPS :: Tuple{Any,Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "rand_UMPS"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    rand_UMPS(d, D; keep_it_real=true)\n",
    "\n",
    "Return a random three-valent tensor A, that defines a uniform MPS (UMPS).\n",
    "The bond dimension of the physical leg should be d, and the bond dimension\n",
    "of the two \"virtual\" legs (the horizontal ones) should be D.\n",
    "keep_it_real is keyword argument, for whether the matrix should be real or\n",
    "complex.\n",
    "\n",
    "This means you can call\n",
    "`rand_UMPS(2, 9)`\n",
    "or\n",
    "`rand_UMPS(2, 9; keep_it_real=true)`\n",
    "and they both give a you real A, but you can also call\n",
    "`rand_UMPS(2, 9; keep_it_real=false)`\n",
    "to get a complex A.\n",
    "\"\"\"\n",
    "function rand_UMPS(d, D; keep_it_real=true)\n",
    "    shp = (D, d, D)\n",
    "    if keep_it_real\n",
    "        A = randn(shp)\n",
    "    else\n",
    "        A_real = randn(shp)\n",
    "        A_imag = randn(shp)\n",
    "        A = complex.(A_real, A_imag) / sqrt(2)\n",
    "    end\n",
    "    return A\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'tm :: Tuple{Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tm"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    tm(A)\n",
    "\n",
    "Return the transfer matrix of A:\n",
    " --A---\n",
    "   |  \n",
    " --A*--\n",
    "\"\"\"\n",
    "function tm(A)\n",
    "    @tensor T[i1,i2,j1,j2] := A[i1,p,j1]*conj(A)[i2,p,j2]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'tm_eigs_dense :: Tuple{Any,Any,Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tm_eigs_dense"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function eig_and_trunc(T, nev; by=identity, rev=false)\n",
    "    S, U = eig(T)\n",
    "    perm = sortperm(S; by=by, rev=rev)\n",
    "    S = S[perm]\n",
    "    U = U[:, perm]\n",
    "    S = S[1:nev]\n",
    "    U = U[:, 1:nev]\n",
    "    return S, U\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "    tm_eigs(A, dirn, nev)\n",
    "\n",
    "Return some of the eigenvalues and vectors of the transfer matrix of A.\n",
    "dirn should be \"L\", \"R\" or \"BOTH\", and determines which eigenvectors to return.\n",
    "nev is the number of eigenpairs to return (starting with the eigenvalues with\n",
    "largest magnitude).\n",
    "\"\"\"\n",
    "function tm_eigs_dense(A, dirn, nev)\n",
    "    T = tm(A)\n",
    "    D = size(T, 1)\n",
    "    T = reshape(T, (D^2, D^2))\n",
    "    nev = min(nev, D^2)\n",
    "    \n",
    "    result = ()\n",
    "    if dirn == \"R\" || dirn == \"BOTH\"\n",
    "        SR, UR = eig_and_trunc(T, nev; by=abs, rev=true)\n",
    "        UR = [reshape(UR[:,i], (D, D)) for i in 1:nev]\n",
    "        result = tuple(result..., SR, UR)\n",
    "    end\n",
    "    if dirn == \"L\" || dirn == \"BOTH\"\n",
    "        SL, UL = eig_and_trunc(T', nev; by=abs, rev=true)\n",
    "        UL = [reshape(UL[:,i], (D, D)) for i in 1:nev]\n",
    "        result = tuple(result..., SL, UL)\n",
    "    end\n",
    "    return result\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'tm_l :: Tuple{Any,Any}' in module 'Main'.\u001b[39m\n",
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'tm_r :: Tuple{Any,Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tm_r"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    tm_l(A, x)\n",
    "\n",
    "Return y, where\n",
    "/------   /------A--\n",
    "|       = |      |  \n",
    "\\- y* -   \\- x* -A*-\n",
    "\"\"\"\n",
    "function tm_l(A, x)\n",
    "    @tensor y[i, j] := (x[a, b] * A[b, p, j]) * conj(A[a, p, i])\n",
    "    return y\n",
    "end\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "    tm_r(A, x)\n",
    "\n",
    "Return y, where\n",
    "-- y -\\   --A-- x -\\\n",
    "      | =   |      |\n",
    "------/   --A*-----/\n",
    "\"\"\"\n",
    "function tm_r(A, x)\n",
    "    @tensor y[i, j] := A[i, p, a] * (conj(A[j, p, b]) * x[a, b])\n",
    "    return y\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tm_eigs_sparse (generic function with 1 method)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function tm_eigs_sparse(A, dirn, nev)\n",
    "    if dirn == \"BOTH\"\n",
    "        SR, UR = tm_eigs_sparse(A, \"R\", nev)\n",
    "        SL, UL = tm_eigs_sparse(A, \"L\", nev)\n",
    "        return SR, UR, SL, UL\n",
    "    else\n",
    "        D = size(A, 1)\n",
    "        x = zeros(eltype(A), (D, D))\n",
    "        if dirn == \"L\"\n",
    "            f = v -> vec(tm_l(A, copy!(x, v)))\n",
    "        else\n",
    "            f = v -> vec(tm_r(A, copy!(x, v)))\n",
    "        end\n",
    "\n",
    "        fmap = LinearMap{eltype(A)}(f, D^2)\n",
    "        S, U, nconv, niter, nmult, resid = eigs(fmap, nev=nev, which=:LM, ritzvec=true)\n",
    "        U = [reshape(U[:,i], (D, D)) for i in 1:size(U, 2)]\n",
    "\n",
    "        return S, U\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tm_eigs (generic function with 1 method)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function tm_eigs(A, dirn, nev; max_dense_D=10)\n",
    "    D = size(A, 1)\n",
    "    if D <= max_dense_D || nev >= D^2\n",
    "        return tm_eigs_dense(A, dirn, nev)\n",
    "    else\n",
    "        return tm_eigs_sparse(A, dirn, nev)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'normalize! :: Tuple{Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "normalize!"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    normalize!(A)\n",
    "\n",
    "Normalize the UMPS defined by A, and return the dominant left and right\n",
    "eigenvectors l and r of its transfer matrix, normalized so that l'*r = 1.\n",
    "\"\"\"\n",
    "function normalize!(A)\n",
    "    SR, UR, SL, UL = tm_eigs(A, \"BOTH\", 1)\n",
    "    S1 = SR[1]\n",
    "    A ./= sqrt(S1)\n",
    "    \n",
    "    l = UL[1]\n",
    "    r = UR[1]  \n",
    "    #We need this to be 1\n",
    "    n = vec(l)'*vec(r)\n",
    "    abs_n = abs(n)\n",
    "    phase_n = abs_n/n\n",
    "    sfac = 1.0/sqrt(abs_n)\n",
    "    l .*= sfac/phase_n\n",
    "    r .*= sfac\n",
    "    return l, r\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework 2 solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'tm_l_op :: Tuple{Any,Any,Any}' in module 'Main'.\u001b[39m\n",
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'tm_r_op :: Tuple{Any,Any,Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tm_r_op"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    tm_l_op(A, O, x)\n",
    "\n",
    "Return y, where\n",
    "/------   /------A--\n",
    "|         |      |  \n",
    "|       = |      O  \n",
    "|         |      |  \n",
    "\\- y* -   \\- x* -A*-\n",
    "\"\"\"\n",
    "function tm_l_op(A, O, x)\n",
    "    @tensor y[i, j] := (x[a, b] * A[b, p2, j]) * (conj(A[a, p1, i]) * conj(O[p1, p2]))\n",
    "    return y\n",
    "end\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "    tm_r_op(A, O, x)\n",
    "\n",
    "Return y, where\n",
    "-- y -\\   --A-- x -\\\n",
    "      |     |      |\n",
    "      | =   O      |\n",
    "      |     |      |\n",
    "------/   --A*-----/\n",
    "\"\"\"\n",
    "function tm_r_op(A, O, x)\n",
    "    @tensor y[i, j] := (A[i, p1, a] * O[p1, p2]) * (conj(A[j, p2, b]) * x[a, b])\n",
    "    return y\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'expect_local :: NTuple{4,Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "expect_local"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    expect_local(A, O, l, r)\n",
    "\n",
    "Return the expectation value of the one-site operator O for the UMPS state\n",
    "defined by the tensor A.\n",
    "\"\"\"\n",
    "function expect_local(A, O, l, r)\n",
    "    l = tm_l_op(A, O, l)\n",
    "    expectation = vec(l)'*vec(r)\n",
    "    return expectation\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'correlator_twopoint :: NTuple{6,Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "correlator_twopoint"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    correlator_twopoint(A, O1, O2, m, l, r)\n",
    "\n",
    "Return the (connected) two-point correlator of operators O1 and O2 for the\n",
    "state UMPS(A), when O1 and O2 are i sites apart, where i ranges from 1 to m. In\n",
    "other words, return <O1_0 O2_i> - <O1> <O2>, for all i = 1,...,m, where the\n",
    "expectation values are with respect to the state |UMPS(A)>.\n",
    "\"\"\"\n",
    "function correlator_twopoint(A, O1, O2, m, l, r)\n",
    "    local_O1 = expect_local(A, O1, l, r)\n",
    "    local_O2 = expect_local(A, O2, l, r)\n",
    "    disconnected = local_O1 * local_O2\n",
    "    \n",
    "    l = tm_l_op(A, O1, l)\n",
    "    r = tm_r_op(A, O2, r)\n",
    "    \n",
    "    result = zeros(eltype(A), m)\n",
    "    result[1] = vec(l)'*vec(r) - disconnected\n",
    "    for i in 1:m\n",
    "        r = tm_r(A, r)\n",
    "        result[i] = vec(l)'*vec(r) - disconnected\n",
    "    end\n",
    "    return result\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'correlation_length :: Tuple{Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "correlation_length"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    correlation_length(A)\n",
    "\n",
    "Return the correlation length ξ of the UMPS defined by A. ξ = - 1/ln(|lambda[2]|),\n",
    "where lambda[2] is the eigenvalue of the MPS transfer matrix with second largest\n",
    "magnitude. (We assume here that UMPS(A) is normalized.)\n",
    "\"\"\"\n",
    "function correlation_length(A)\n",
    "    S, U = tm_eigs(A, \"L\", 2)\n",
    "    s2 = S[2]\n",
    "    ξ = -1/log(abs(s2))\n",
    "    return ξ\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# But I don't care about random MPSes!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that you know how to compute correlation functions for MPS states, we should start thinking about how do we get a useful MPS state in the first place. So far all our states have been defined by tensors $A$ with random elements.\n",
    "\n",
    "Some fun little analytical examples exist, such as how to represent GHZ or W states as MPSes. You should google them up. The real question is though, given an arbitrary local Hamiltonian, can we find an MPS for the ground state. We know that for gapped Hamiltonians (i.e. most Hamiltonians) correlations in the ground state decay exponentially and entanglement is limited, so there's hope that the answer would be \"yes we can\".\n",
    "\n",
    "You can think about this problem in a few different ways. One is to ask, what's the lowest energy state we can find on the manifold of MPS states. The question is then a matter of optimizing for minimal expectation value $\\langle \\psi | H | \\psi \\rangle$ under the constraint that the state $|\\psi \\rangle$ is an MPS. Another way to approach this would be writing down some expression for the exact ground state (an expression we obviously can't just go an evaluate, otherwise there would be no problem), and then try to massage that expression into an MPS.\n",
    "\n",
    "Various algorithms for implementing these different approaches in different situations (finite-size, infinite-size, translation invariant or not, what kind of boundary conditions) exist. Some examples are the Density Matrix Renormalization Group (DMRG), Time-Dependent Variational Principle (TDVP) and Time-Evolving Block Decimation (TEBD). We'll implement an infinite-lattice version of TEBD. Why TEBD? Because it's simple, and because it you can also use it to do time-evolution of MPS states. (I maintain that the choice to pick TEBD is independent of the fact that Guifre invented it.)\n",
    "\n",
    "Implementing iTEBD will be the rest of this minicourse. However, we'll only really get to that next time. First, we need to talk a bit about gauge freedom in MPS and how to fix this freedom, since TEBD will be making use of a specific gauge."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gauge freedom and canonical forms\n",
    "In any tensor network, I can always take any contracted leg, and insert on that leg a pair of invertibale matrices $g g^{-1} = \\mathbb{1}$. This doesn't affect the value the network contracts to in anyway: I've just changed the basis in which I perform one of the sums (traces).\n",
    "\n",
    "This means that an MPS representation of a state is not unique. I can insert any $g g^{-1}$ on any of the virtual legs of the MPS, and the state stays the same. Obviously, I would like to stay within the framework of translation invariant MPSes, so I allow only for gauge transformations (as these changes of basis are called) of the form\n",
    "```\n",
    "--A-- -> --gi--A--g--\n",
    "  |            |\n",
    "```\n",
    "where by `gi` I mean the inverse of $g$.\n",
    "\n",
    "Since we have some freedom in picking $A$, we should probably try to make use of it somehow, while fixing the freedom. How? Well, we could try to impose some nice properties on $A$, if such a property can be achieved by a gauge transformation. And it turns out there's a really useful property that we can impose, which goes as follows.\n",
    "\n",
    "First, we have to complicate our setup slightly. Previously our how MPS was defined by $A$. Now, let's think of an MPS that is instead written as\n",
    "```\n",
    "... --λ--Γ--λ--Γ--λ--Γ--λ-- ...\n",
    "         |     |     |\n",
    "```\n",
    "Here $\\Gamma$ is some tensor, and $\\lambda$ is a diagonal matrix. Obviously we could just absorb $\\lambda$ into $\\Gamma$ and a get single tensor $A$. The point is, by using the gauge freedom, we can choose $\\Gamma$ and $\\lambda$ in such a way, that *$\\lambda$ always has on its diagonal the Schmidt coefficients of partitioning the state into two parts at $\\lambda$*. What are Schmidt coefficients? https://en.wikipedia.org/wiki/Schmidt_decomposition\n",
    "\n",
    "Another way of stating the same condition is to say that the left dominant eigenvector of\n",
    "```\n",
    "--λ--Γ--\n",
    "     |\n",
    "--λ--Γ--\n",
    "```\n",
    "and the right dominant eigenvector of\n",
    "```\n",
    "--Γ--λ--\n",
    "  |\n",
    "--Γ--λ--\n",
    "```\n",
    "are both the identity matrix.\n",
    "\n",
    "Choosing our gauge so that this property holds, turns out to be quite handy. You'll see next time how it's used in the TEBD algorithm, but you can probably already imagine that it's kinda nice that on each virtual legs you have these Schmidt coefficients $\\lambda$, that immediately tell you everything you would want to know about the entanglement between the two sides of the state separated by this leg. The bond dimension is also obviously the Schmidt rank, which gives us a very concrete, physical characterization of what high bond dimension and low bond dimension mean.\n",
    "\n",
    "So, given $A$, how do we find a gauge transformation such that the above property holds? I'm getting tired of typing, so just go read Section II of https://arxiv.org/pdf/0711.3960.pdf, that will explain it all. Then, it's time to define your homework."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Homework 3\n",
    "#### -1)\n",
    "Make sure that, for the expectation values and correlation functions,<br>\n",
    "a) Your expectation values for Hermitian operators are real.<br>\n",
    "b) Your correlation functions decay like they should (make a plot).<br>\n",
    "c) You can generate two-point correlators with $D=100$ for distances $1, \\dots, 100$ or so in a few seconds.<br>\n",
    "Feel free to make use of my code above, as long as you understand what's going on.\n",
    "\n",
    "#### 0)\n",
    "Read Section II of https://arxiv.org/pdf/0711.3960.pdf. Note that there they start from some $\\Gamma$ and $\\lambda$, that are not the canonical form, and try to turn them into $\\Gamma'$ and $\\lambda'$ that would be in the canonical form. We, on the other hand, start from just $A$. You can easily translate to their language by setting their $\\Gamma$ to your $A$, and their $\\lambda$ to just the identity matrix.\n",
    "\n",
    "#### 1)\n",
    "As explained in above reference, the form of MPS transfer matrix $T$, i.e. the symmetry it has because of how it consists of $A$ and $A^\\ast$, guarantees that the dominant eigenvectors $l$ and $r$ are Hermitian and positive semi-definite. However, as we discussed in the last meeting, when we numerically ask for the eigenvectors, they come with an arbitrary phase. Last time we modified `normalize!` so that it partially fixes this phase (and also the norm of $l$ and $r$), imposing `vec(l)'*vec(r) == 1`. Your first job is to further modify `normalize!` so that it fully fixes the phases, so that the `l` and `r` it returns are Hermitian and positive semi-definite, in addition to having the property that `vec(l)'*vec(r) == 1`.\n",
    "\n",
    "How do you do that? Well, you should definitely start by getting the `normalize!` function from above, and only make small modifications to that. Since we know that $r$ is actually $a \\cdot \\tilde{r}$, where $\\tilde{r}$ is Hermitian and pos. semi-def. and $a$ is a scalar, we can just take the trace $\\mathrm{Tr} r = a \\mathrm{Tr} \\tilde{r}$, and we know that $\\mathrm{Tr} \\tilde{r}$ is something real and positive. Any part of $\\mathrm{Tr} r$ that is not real and positive, is $a$, and you divide $r$ by $a$ to get $\\tilde{r}$. And similarly for $l$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'normalize! :: Tuple{Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "normalize!"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    normalize!(A)\n",
    "\n",
    "Normalize the UMPS defined by A, and return the dominant left and right\n",
    "eigenvectors l and r of its transfer matrix, normalized so that they are\n",
    "both Hermitian and positive semi-definite (when thought of as matrices),\n",
    "and l'*r = 1.\n",
    "\"\"\"\n",
    "function normalize!(A)\n",
    "    SR, UR, SL, UL = tm_eigs(A, \"BOTH\", 1)\n",
    "    S1 = SR[1]\n",
    "    A ./= sqrt(S1)\n",
    "    \n",
    "    l = UL[1]\n",
    "    r = UR[1]  \n",
    "    # We want both l and r to be Hermitian and pos. semi-def.\n",
    "    # We know they are that, up to a phase.\n",
    "    # We can find this phase, and divide it away, because it is also the\n",
    "    # phase of the trace of l (respectively r).\n",
    "    r_tr = trace(r)\n",
    "    phase_r = r_tr/abs(r_tr)\n",
    "    r ./= phase_r\n",
    "    l_tr = trace(l)\n",
    "    phase_l = l_tr/abs(l_tr)\n",
    "    l ./= phase_l\n",
    "    # Finally divide them by a real scalar that makes\n",
    "    # their inner product be 1.\n",
    "    n = vec(l)'*vec(r)\n",
    "    abs_n = abs(n)\n",
    "    phase_n = n/abs_n\n",
    "    (phase_n ≉ 1) && warn(\"In normalize! phase_n = \", phase_n, \" ≉ 1\")\n",
    "    sfac = sqrt(abs_n)\n",
    "    l ./= sfac\n",
    "    r ./= sfac\n",
    "    return l, r\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2)\n",
    "Now that we can trust our $l$ and $r$ to be Hermitian and $\\geq 0$, we can go ahead and implement a function that takes in $A$, $l$ and $r$, and outputs $\\Gamma$ and $\\lambda$, that define the same MPS, but now in canonical form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mreplacing docs for 'canonical_form :: Tuple{Any,Any,Any}' in module 'Main'.\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "canonical_form"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    canonical_form(A, l, r)\n",
    "\n",
    "Return a three-valent tensor Γ and a vector λ, that define the canonical\n",
    "of the UMPS defined by A. l and r should be the normalized dominant\n",
    "left and right eigenvectors of A.\n",
    "\"\"\"\n",
    "function canonical_form(A, l, r)\n",
    "    l_H = 0.5*(l + l')\n",
    "    r_H = 0.5*(r + r')\n",
    "    (l_H ≉ l) && warn(\"In canonical_form, l is not Hermitian: \", l)\n",
    "    (r_H ≉ r) && warn(\"In canonical_form, r is not Hermitian: \", r)\n",
    "    evl, Ul = eig(Hermitian(l_H))\n",
    "    evr, Ur = eig(Hermitian(r_H))\n",
    "    X = Ur * Diagonal(sqrt.(complex.(evr)))\n",
    "    YT = Diagonal(sqrt.(complex.(evl))) * Ul'\n",
    "    U, λ, V = svd(YT*X)\n",
    "    Xi = Diagonal(sqrt.(complex.(1./evr))) * Ur'\n",
    "    YTi = Ul * Diagonal(sqrt.(complex.(1 ./ evl)))\n",
    "    @tensor Γ[x,i,y] := (V'[x,a] * Xi[a,b]) * A[b,i,c] * (YTi[c,d] * U[d,y])\n",
    "    return Γ, λ\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3)\n",
    "Prove to yourself that<br>\n",
    "a)<br>\n",
    "the way you've chosen to transform the MPS, going form $A$ to $(\\Gamma, \\lambda)$, is actually a gauge transformation. In other words, prove that they both define the same state.\n",
    "\n",
    "b)<br>\n",
    "that the MPS defined $(\\Gamma, \\lambda)$ is actually in canonical form, i.e., that the left dominant eigenvector of\n",
    "```\n",
    "--λ--Γ--\n",
    "     |\n",
    "--λ--Γ--\n",
    "```\n",
    "and the right dominant eigenvector of\n",
    "```\n",
    "--Γ--λ--\n",
    "  |\n",
    "--Γ--λ--\n",
    "```\n",
    "are both the identity matrix. Do this both numerically and analytically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "true\n",
      "true\n",
      "true\n",
      "true\n",
      "true\n",
      "2.3936273716601304e-15\n",
      "2.404546554598619e-15\n",
      "true\n",
      "true\n"
     ]
    }
   ],
   "source": [
    "let\n",
    "    d = 2\n",
    "    D = 20\n",
    "    A = rand_UMPS(d, D; keep_it_real=false)\n",
    "    \n",
    "    # Test the new normalization function.\n",
    "    l, r = normalize!(A)\n",
    "    SR, UR, SL, UL = tm_eigs(A, \"BOTH\", 1)\n",
    "    # The MPS is normalized.\n",
    "    println(SR[1] ≈ 1.0)\n",
    "    println(SL[1] ≈ 1.0)\n",
    "    # l and r are really eigenvectors, and have their inner product be 1.\n",
    "    println(l ≈ tm_l(A, l))\n",
    "    println(r ≈ tm_r(A, r))\n",
    "    println(vec(l)'*vec(r) ≈ 1.0)\n",
    "    # l and r are really Hermitian.\n",
    "    println(vecnorm(l - l'))\n",
    "    println(vecnorm(r - r'))\n",
    "    # l and r are really pos. semi-def.\n",
    "    evl = eig(l)[1]\n",
    "    evr = eig(r)[1]\n",
    "    allposl = all(abs(imag(i)) < 1e-14 && real(i) >= 0 for i in evl)\n",
    "    allposr = all(abs(imag(i)) < 1e-14 && real(i) >= 0 for i in evr)\n",
    "    println(allposl)\n",
    "    println(allposr)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "vecnorm(should_be_unity_l - eye(D, D)) = 2.503750070484813e-14\n",
      "vecnorm(should_be_unity_r - eye(D, D)) = 2.6864148124101486e-14\n",
      "true\n",
      "true\n"
     ]
    }
   ],
   "source": [
    "let\n",
    "    d = 2\n",
    "    D = 10\n",
    "    A = rand_UMPS(d, D; keep_it_real=false)\n",
    "    \n",
    "    l, r = normalize!(A)\n",
    "    Γ, λ = canonical_form(A, l, r)\n",
    "    λm = diagm(λ)\n",
    "    \n",
    "    # Check that identity really is the dominant eigenvector as it should.\n",
    "    @tensor should_be_unity_l[x,y] := (λm[a,b] * Γ[b,i,x]) * (conj(λm)[a,c] * conj(Γ)[c,i,y])\n",
    "    @tensor should_be_unity_r[x,y] := (Γ[x,i,b] * λm[b,a]) * (conj(Γ)[y,i,c] * conj(λm)[c,a])\n",
    "    @show vecnorm(should_be_unity_l - eye(D,D))\n",
    "    @show vecnorm(should_be_unity_r - eye(D,D))\n",
    "    \n",
    "    # Check that we get the same two-point correlators for the original A,\n",
    "    # and for a couple of different ways of contracting together λ and Γ.\n",
    "    O1 = randn(d, d)\n",
    "    O2 = randn(d, d)\n",
    "    O1 = (O1 + O1')/2\n",
    "    O2 = (O2 + O2')/2\n",
    "    \n",
    "    m = 3\n",
    "    # The original A.\n",
    "    corrsA = correlator_twopoint(A, O1, O2, m, l, r)\n",
    "    # Contract λ into Γ from the right.\n",
    "    @tensor A2[x,i,y] := Γ[x,i,a] * λm[a,y]\n",
    "    l2, r2 = normalize!(A2)\n",
    "    corrsA2 = correlator_twopoint(A2, O1, O2, m, l2, r2)\n",
    "    println(all(corrsA .≈ corrsA2))\n",
    "    # Contract λ into Γ from the left.\n",
    "    @tensor A3[x,i,y] := λm[x,a] * Γ[a,i,y]\n",
    "    l3, r3 = normalize!(A3)\n",
    "    corrsA3 = correlator_twopoint(A3, O1, O2, m, l3, r3)\n",
    "    println(all(corrsA .≈ corrsA3))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## iTEBD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we have arrived at a point where we have everything we need to implement iTEBD, and use it to find the ground state for a given Hamiltonian. We'll first go through the Suzuki-Trotter decomposition of the imaginary time-development operator, see how that leads to an algorithm (iTEBD) for finding an MPS for the ground state, and then finally see why we should keep our MPS in the canonical form while applying this algorithm.\n",
    "\n",
    "Given a Hamiltonian $H$, here's one way to define the ground state:\n",
    "$$ |E_0\\rangle \\sim e^{-\\beta H} | \\psi \\rangle, \\quad \\beta >> 1. $$\n",
    "Here $\\sim$ means $=$ up to normalization. $| \\psi \\rangle$ can be any state, as long as it's not orthogonal to $|E_0 \\rangle$ (just pick one at random). If $\\beta$ is large enough, $e^{\\beta H}$ becomes proportional to the projector to $| E_0 \\rangle$, and the above equation holds.\n",
    "\n",
    "Of course, $e^{- \\beta H}$ is a matrix exponentially big in the system size. Instead of trying to evaluate it, what we'll do is decompose it into a circuit of local gates. This is called the Suzuki-Trotter decomposition.\n",
    "\n",
    "For this, we'll need to assume $H$ is a sum over local terms. Let's assume it consists of two-site terms:\n",
    "$$ H = \\sum_i h_{i,i+1}. $$\n",
    "Then we can write \n",
    "$$ H = \\sum_{i \\text{ odd}} h_{i,i+1} + \\sum_{i \\text{ even}} h_{i,i+1} = H_{\\text{odd}} + H_{\\text{even}}. $$\n",
    "The useful point about this is that all the terms in $H_{\\text{odd}}$ operate on different sites, and thus commute with each other, and similarly for $H_{\\text{even}}$. Next, introduce a small parameter $\\tau$, and use the Baker–Campbell–Hausdorff formula as follows:\n",
    "$$ e^{-\\beta H} = \\big[ e^{-\\tau H} \\big]^{\\frac{\\beta}{\\tau}}\n",
    "= \\big[ e^{-\\tau H_{\\text{odd}} \\; -\\tau H_{\\text{even}}} \\big]^{\\frac{\\beta}{\\tau}}\n",
    "= \\big[ e^{-\\tau H_{\\text{odd}}} \\; e^{-\\tau  H_{\\text{even}}} \\; e^{-\\tau^2 (\\text{commutators of $H_{\\text{odd}}$ and $H_{\\text{even}}$)}} \\big]^{\\frac{\\beta}{\\tau}}.$$\n",
    "Making $\\tau$ small enough so that anything $O(\\tau^2)$ can be ignored, makes this\n",
    "$$ \\big[ \\prod_{i \\text{ odd}} e^{-\\tau h_{i, i+1}} \\; \\prod_{i \\text{ even}} e^{-\\tau h_{i, i+1}} \\big]^{\\frac{\\beta}{\\tau}}.$$\n",
    "Diagrammatically we can write this as\n",
    "\n",
    "<img src=\"fig/trotter.svg\">\n",
    "\n",
    "Based on this, we can formulate the basic idea of the iTEBD algorithm (note that the i in iTEBD stands for infinite systems, but so far everything we've done works the same way for finite systems): Take the expression $e^{-\\beta H} | \\psi \\rangle$. Decompose $e^{-\\beta H}$ into a circuit using the above Suzuki-Trotter decomposition, with some choice of the parameter $\\tau$ that is small. Take the initial state $| \\psi \\rangle$ to be an MPS. The expression for the ground state then looks like\n",
    "\n",
    "<img src=\"fig/trotter_gs.svg\">\n",
    "\n",
    "With these steps, we've gotten rid of all the exponentially large objects in $e^{-\\beta H} | \\psi \\rangle$: The evolution operator is described by local gates, and the state is defined by local MPS tensors. We can then absorb the layers of the circuit into the MPS, layer by layer, contracting only local things at a time. If we assume translation invariance this can all be done at the thermodynamic limit, i.e., with an infinite system. By \"absorbing\" a layer I mean this:\n",
    "\n",
    "<img src=\"fig/trotter_absorb.svg\">\n",
    "\n",
    "At the last step the four-valent tensor has been split into two MPS matrices. This can be done by considering it as a matrix from the two left-most legs to the two right-most legs, and performing and a singular value decomposition (SVD).\n",
    "\n",
    "With this method we can keep absorbing layers of the circuit $e^{-\\beta H}$ into our MPS, effectively increasing $\\beta$, until the state doesn't change significantly any more, i.e., until the procedure converges. At the end, our MPS should be the ground state of the Hamiltonian, subject to the caveat that since our $\\tau$ is small, and not exactly zero, and our $\\beta$ is large, and not exactly infinite, we only find the ground state approximately (the error coming from ignoring the $\\tau^2$ terms is called the Trotter error).\n",
    "\n",
    "Seems great! So why don't we get coding? Well, there's one problem here. Every time we use an SVD to split\n",
    "\n",
    "<img src=\"fig/trotter_split.svg\">\n",
    "\n",
    "we are SVDing a $dD \\times dD$ matrix. This means there are $dD$ singular values, meaning the bond dimension of our MPS after this procedure will be $dD$ (whereas a moment before it was just $D$). So our bond dimension will grow during this procedure, and exponentially too! This is not good.\n",
    "\n",
    "However, we never claimed we wanted the exact ground state. What we wanted was a low bond dimension MPS that approximates the ground state. So somehow we should reduce the bond dimension back to what it was, while not changing the state our MPS describes too much (or as little as possible). We should probably take our MPS tensors, that now have some indices ranging over $dD$ values, and throw out some of the elements (like throwing out some columns of a matrix) so that all the dimensions are at most $D$. But which ones do we throw out?\n",
    "\n",
    "Turns out, you should always throw out the parts that correspond to the smallest Schmidt values. This is the same result as the fact that the optimal low-rank approximation to matrix is given by its truncated SVD (see Low rank approximation in https://en.wikipedia.org/wiki/Singular-value_decomposition#Applications_of_the_SVD), since the Schmidt decomposition is just an SVD in disguise. Thus we should be performing this whole iTEBD procedure in the canonical form, so that on each leg we have the Schmidt value corresponding to that bipartition of the state. Then every time we perform one of the SVDs that increases the bond dimension, we can go back to the canonical form and remove from the neighbouring tensors the parts that correspond to the smallest Schmidt values (which has the same effect as setting those values to zero).\n",
    "\n",
    "Thus a naive summary of the iTEBD algorithm would be as follows:\n",
    "1. Absorb a layer of the circuit for $e^{-\\beta H}$ into your MPS.\n",
    "2. Use SVDs to transform the network back into an MPS, but now with a larger bond dimension.\n",
    "3. Gauge transform this network into the canonical form.\n",
    "4. Throw out the smallest Schmidt values on some of the legs, to reduce the bond dimension back to $D$.\n",
    "5. Go to 1., unless your state has already converged, i.e., isn't changing any more, in which case we are done.\n",
    "\n",
    "Note that every second time we go through this process, we should absorb a layer of $e^{-\\tau H_{\\text{odd}}}$ and every second time a layer of $e^{-\\tau H_{\\text{even}}}$.\n",
    "\n",
    "I call the above summary \"naive\", because I've ignored some details. For instance, the canonical form of your MPS will now have translation symmetry only by two sites, which means you'll have to keep to tensors $\\Gamma_A$ and $\\Gamma_B$ and corresponding two vectors of Schmidt values $\\lambda_A$ and $\\lambda_B$. It also turns out it's not important to be canonicalizing your MPS at every single step: If you start with an MPS in canonical form, you can perform the absorption, SVD and truncation all without recanonicalizing in between, since your MPS will stay close enough to being canonical. You can even repeat this for a few iterations before having to canonicalize again. These kinds of details you should check up on your own from the literature. I mainly recommend<br>\n",
    "https://arxiv.org/pdf/cond-mat/0605597.pdf<br>\n",
    "for understanding the iTEBD algorithm, and Section II of<br>\n",
    "https://arxiv.org/pdf/0711.3960.pdf<br>\n",
    "for how to deal with the canonical form (including the one with two-site translation symmetry).\n",
    "Other resources exist, which you can try if you want to, especially given that the first paper above is quite terse. There's even a wikipedia article, https://en.wikipedia.org/wiki/Time-evolving_block_decimation, although I must admit I haven't read it.\n",
    "\n",
    "Your homework for next time is to go and implement this! You should start by reading https://arxiv.org/pdf/cond-mat/0605597.pdf, and then try putting all of it into code. Don't worry if you don't get it all done by the next meeting. What I want to see is a serious attempt to get as far as you can (including writing some code!). Once you hit a roadblock and don't know what to do anymore, you should try to understand what exactly it is you don't understand, and then come to the next meeting with questions. Next time we'll go through what you got done and clarify what needs to be clarified. We'll hopefully also talk a bit about things like performance, how to measure convergence, etc.\n",
    "\n",
    "I've written below some templates that you can use as a basis of your own code. They also include some hints as to how to go about implementing things."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "itebd_optimize (generic function with 1 method)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# U should be the four-valent tensor e^(-\\tau h).\n",
    "# D is the bond dimension to be used.\n",
    "# You probably want to add some other arguments for things like\n",
    "# how exactly to decide when convergence has been reached.\n",
    "function itebd_optimize(U, D)\n",
    "    d = size(U, 1)\n",
    "    ΓA, λA, ΓB, λB = itebd_random_initial(d, D)\n",
    "    # Here there should be a loop, that applies one step of iTEBD,\n",
    "    # i.e., absorbs U into ΓA, λA, ΓB, λB, and splits with an SVD\n",
    "    # to get back to an MPS defined by some other ΓA2, λA2, ΓB2, λB2,\n",
    "    # and truncates the SVD back to bond dimension D.\n",
    "    # Your ΓA2, λA2, ΓB2, λB2 will not be exactly in the canonical form\n",
    "    # any more after all this, so you should probably recanonicalize.\n",
    "    #\n",
    "    # Finally, you should check whether doing all this changed your MPS\n",
    "    # significantly: If not, the iTEBD procedure has converged, and we can\n",
    "    # quit and return. Convergence can be checked in many ways, but one\n",
    "    # good one is to compare the Schmidt values λA2, λB2 to the previous\n",
    "    # Schmidt values λA, λB, and see if they've changed a lot.\n",
    "    return ΓA, λA, ΓB, λB\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "itebd_step (generic function with 1 method)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function itebd_step(ΓA, λA, ΓB, λB, U)\n",
    "    # This is the heart of the algorithm: Absorb, split, truncate.\n",
    "    # See Figure 3 of https://arxiv.org/pdf/cond-mat/0605597.pdf\n",
    "    # of what should go in here. Remember that you should do the same\n",
    "    # thing twice, first U operating of odd pairs of sites and then on\n",
    "    # even pairs of sites.\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "itebd_random_initial (generic function with 1 method)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function itebd_random_initial(d, D)\n",
    "    # Generate a random initial guess for the MPS. It should have\n",
    "    # translation symmetry by two sites, i.e., it should consist of\n",
    "    # tensors ΓA, λA, ΓB, λB, and be in the canonical form.\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "truncate_svd"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# A handy function for truncating SVDs. You'll need this in implementing itebd_step.\n",
    "\"\"\"\n",
    "    truncate_svd(U, S, V, D)\n",
    "\n",
    "Given an SVD of some matrix M as M = U*diagm(S)*V', truncate this\n",
    "SVD, keeping only the D largest singular values.\n",
    "\"\"\"\n",
    "function truncate_svd(U, S, V, D)\n",
    "    U = U[:, 1:D]\n",
    "    S = S[1:D]\n",
    "    V = V[:, 1:D]\n",
    "    return U, S, V\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "double_canonicalize (generic function with 1 method)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function double_canonicalize(ΓA, λA, ΓB, λB)\n",
    "    # Take in an MPS defined by ΓA, λA, ΓB, λB, and do a gauge transformation to\n",
    "    # put it in a canonical form. See the end of Section II of\n",
    "    # https://arxiv.org/pdf/0711.3960.pdf for how to do this.\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's test this algorithm of ours!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "build_ising_ham (generic function with 2 methods)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function build_ising_ham(h=1.0)\n",
    "    X = [0 1; 1 0]\n",
    "    Z = [1 0; 0 -1]\n",
    "    I2 = eye(2)\n",
    "    XX = kron(X, X)\n",
    "    ZI = kron(Z, I2)\n",
    "    IZ = kron(I2, Z)\n",
    "    H = -(XX + h/2*(ZI+IZ))\n",
    "    return H\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mMethodError: no method matching start(::Void)\u001b[0m\nClosest candidates are:\n  start(\u001b[91m::SimpleVector\u001b[39m) at essentials.jl:258\n  start(\u001b[91m::Base.MethodList\u001b[39m) at reflection.jl:560\n  start(\u001b[91m::ExponentialBackOff\u001b[39m) at error.jl:107\n  ...\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mMethodError: no method matching start(::Void)\u001b[0m\nClosest candidates are:\n  start(\u001b[91m::SimpleVector\u001b[39m) at essentials.jl:258\n  start(\u001b[91m::Base.MethodList\u001b[39m) at reflection.jl:560\n  start(\u001b[91m::ExponentialBackOff\u001b[39m) at error.jl:107\n  ...\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1mitebd_optimize\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Array{Float64,4}, ::Int64\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./In[31]:7\u001b[22m\u001b[22m",
      " [2] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:515\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "let\n",
    "    # Here's a piece that you can use as a basis for testing your iTEBD.\n",
    "    # You should probably implement function that measures the per site ground\n",
    "    # state energy (very much like your previous expect_local, but now for an\n",
    "    # MPS of the form ΓA, λA, ΓB, λB), and see if you are getting it right.\n",
    "    # The exact per site ground state energy for the Ising model at criticality\n",
    "    # (h=1.0) is -4/π, so if you are getting something close to that, you're doing\n",
    "    # great.\n",
    "    τ = 1e-2\n",
    "    D = 50\n",
    "    h = build_ising_ham()\n",
    "    U = expm(-τ*h)\n",
    "    U = reshape(U, (2,2,2,2))\n",
    "    h = reshape(h, (2,2,2,2))\n",
    "    ΓA, λA, ΓB, λB = itebd_optimize(U, D)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}