tensor_networks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "63ef3c93",
   "metadata": {},
   "source": [
    "# Efficient quantum simulation using Tensor Network states"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b737328",
   "metadata": {},
   "source": [
    "The aim of this notebook is to demonstrate the impact of the Tensor Network state representation\n",
    "on computational complexity in quantum computing simulations. During the morning lecture, the asymptotic complexity of\n",
    "TN methods was discussed in an abstract manner, and now we will see, through some concrete examples via the ITensors package\n",
    "precisely what the practical difference in runtimes is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "32dd2d09",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "using ITensors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b06c335f",
   "metadata": {},
   "source": [
    "# Exercise 1: 5 minutes "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ddef59e",
   "metadata": {},
   "source": [
    "Complete the code in the following function to compute the Kronecker product ⊗ of vectors/matrices a,b.\n",
    "Note that the Kronecker product gives the abstract tensor product expressed in a particular basis.\n",
    "https://en.wikipedia.org/wiki/Kronecker_product"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0bf4f90a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function tensor(a,b)\n",
    "    return nothing #fill this in so the test cell below passes\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16ee9e64",
   "metadata": {},
   "source": [
    "run the below cell to test your implementation of tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5679cee3",
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "AssertionError: tensor([1 2; 3 4], [5 6; 7 8]) == [5 6 10 12; 7 8 14 16; 15 18 20 24; 21 24 28 32]",
     "output_type": "error",
     "traceback": [
      "AssertionError: tensor([1 2; 3 4], [5 6; 7 8]) == [5 6 10 12; 7 8 14 16; 15 18 20 24; 21 24 28 32]",
      "",
      "Stacktrace:",
      " [1] test_tensor()",
      "   @ Main ./In[3]:2",
      " [2] top-level scope",
      "   @ In[3]:7",
      " [3] eval",
      "   @ ./boot.jl:368 [inlined]",
      " [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
      "   @ Base ./loading.jl:1428"
     ]
    }
   ],
   "source": [
    "function test_tensor()\n",
    "    @assert tensor([1 2; 3 4], [5 6; 7 8]) == [5 6 10 12; 7 8 14 16; 15 18 20 24; 21 24 28 32]\n",
    "    @assert tensor([1, 2], [3, 4]) == [3, 4, 6, 8]\n",
    "    return nothing\n",
    "end\n",
    "\n",
    "test_tensor()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebd1c307",
   "metadata": {},
   "source": [
    "hint: there is an inbuilt kron function in the language"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed16c698",
   "metadata": {},
   "source": [
    "## Key takeaway: taking the tensor product multiplies all the dimensions of the inputs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4e4fffb",
   "metadata": {},
   "source": [
    "# Exercise 2: 5 minutes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1ba7336",
   "metadata": {},
   "source": [
    "Using the represenation |0> = [1, 0] and |1> = [0, 1], compute for n qubits,\n",
    "GHZ state:      (|0> ⊗ ... ⊗ |0>) + (|1> ⊗ ... |1>)\n",
    "Hadamard state: (|0> + |1>) ⊗ ... ⊗ (|0> + |1>)>\n",
    "Note that both these states are very relevant in quantum computing:\n",
    "The GHZ state is the unmagnetized ground state of the Ising model,\n",
    "The Hadamard state is the input in many early quantum algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b25b362a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "hadamard_state (generic function with 1 method)"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function ghz_state(n)\n",
    "    return nothing\n",
    "end\n",
    "\n",
    "function hadamard_state(n)\n",
    "    return nothing\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de6be729",
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "AssertionError: ghz_state(3) == [1, 0, 0, 0, 0, 0, 0, 1]",
     "output_type": "error",
     "traceback": [
      "AssertionError: ghz_state(3) == [1, 0, 0, 0, 0, 0, 0, 1]",
      "",
      "Stacktrace:",
      " [1] test_states()",
      "   @ Main ./In[67]:2",
      " [2] top-level scope",
      "   @ In[67]:7",
      " [3] eval",
      "   @ ./boot.jl:368 [inlined]",
      " [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
      "   @ Base ./loading.jl:1428"
     ]
    }
   ],
   "source": [
    "function test_states()\n",
    "    @assert ghz_state(3) == [1, 0, 0, 0, 0, 0, 0, 1]\n",
    "    @assert hadamard_state(3) == [1, 1, 1, 1, 1, 1, 1, 1]\n",
    "    return nothing\n",
    "end\n",
    "\n",
    "test_states()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c307db04",
   "metadata": {},
   "source": [
    "hint: recursively call tensor, for example by using foldl(tensor, .)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d698325",
   "metadata": {},
   "source": [
    "## Key takeaway:  The GHZ state is very sparse in this representation, the Hadamard state is not"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0466a7c5",
   "metadata": {},
   "source": [
    "# Exercise 3: 5 minutes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35833e67",
   "metadata": {},
   "source": [
    "Print the GHZ and Hadamard states for 30 qubits.\n",
    "Compute the time it takes to do the inner product <hadamard_state(n), hadamard_state(n)>\n",
    "for n ranging up to 25-30 (depending on the speed of your machine)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4cb35918",
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "do some stuff here to print/plot the timings",
     "output_type": "error",
     "traceback": [
      "do some stuff here to print/plot the timings",
      "",
      "Stacktrace:",
      " [1] error(s::String)",
      "   @ Base ./error.jl:35",
      " [2] top-level scope",
      "   @ In[68]:6",
      " [3] eval",
      "   @ ./boot.jl:368 [inlined]",
      " [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
      "   @ Base ./loading.jl:1428"
     ]
    }
   ],
   "source": [
    "function time_inner(n)\n",
    "    t = @timed 1+1 #replace this with a computation of inner product\n",
    "    return t[:time]\n",
    "end\n",
    "\n",
    "error(\"do some stuff here to print/plot the timings\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a3737c1",
   "metadata": {},
   "source": [
    "hint: use LinearAlgebra.dot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf1fe27f",
   "metadata": {},
   "source": [
    "## Key takeaway:  \n",
    "## The length of both the GHZ and Hadamard states increases as 2^n\n",
    "## The time to compute the inner product increases exponentially in n.\n",
    "## Even on very big computers, going above 40-50 qubits is impossible."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65362811",
   "metadata": {},
   "source": [
    "To be able to simulate large numbers of qubits, a better representation of the relevant states in the Hilbert space is needed.\n",
    "For a state like the GHZ state, which is sparse in the basis above, it should be obvious that such a represenation exists.\n",
    "That there is a good representation that encapsulates all interesting states (such as the Hadamard state) is not obvious.\n",
    "\n",
    "During the morning lecture, it was explained that the Tensor Network representation seems to fit this bill, and why.\n",
    "We will now redo the above exercises in the Tensor Network representation, and see the benefits it provides."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48fb24e7",
   "metadata": {},
   "source": [
    "# Exercise 4: 15 minutes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c2c5377",
   "metadata": {},
   "source": [
    "Construct tensor network representations of the\n",
    "GHZ state:      (|0> ⊗ ... ⊗ |0>) + (|1> ⊗ ... |1>)\n",
    "Hadamard state: (|0> + |1>) ⊗ ... ⊗ (|0> + |1>)>\n",
    "In ITensors such a representation is called an MPS (matrix product state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d66cf468",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tn_hadamard_state (generic function with 1 method)"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function tn_ghz_state(n)\n",
    "    return nothing\n",
    "end\n",
    "\n",
    "function tn_hadamard_state(n)\n",
    "    return nothing\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c890fa37",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "test_tn_states (generic function with 1 method)"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function test_tn_states()\n",
    "    N = 5   \n",
    "\n",
    "    function elt(ϕ, inds)\n",
    "        s = siteinds(ϕ)\n",
    "        V = ITensor(1.)\n",
    "        for j=1:N\n",
    "            V *= (ϕ[j]*state(s[j],inds[j]))\n",
    "        end\n",
    "        return scalar(V)\n",
    "    end\n",
    "\n",
    "    ψ = tn_ghz_state(N)\n",
    "    @assert elt(ψ, [1,1,1,1,1]) ≈ 1    \n",
    "    @assert elt(ψ, [2,2,2,2,2]) ≈ 1    \n",
    "    @assert elt(ψ, [1,1,1,2,1]) ≈ 0\n",
    "    \n",
    "\n",
    "    ψ = tn_hadamard_state(N)\n",
    "    @assert elt(ψ, [1,1,1,2,1]) ≈ 1    \n",
    "    @assert elt(ψ, [1,2,1,2,2]) ≈ 1\n",
    "\n",
    "    return nothing\n",
    "end\n",
    "test_tn_states()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2f2b481",
   "metadata": {},
   "source": [
    "hint1:  You can create an MPS using the states constructed in Exercise 2, and an array of siteinds\n",
    "        https://itensor.github.io/ITensors.jl/stable/examples/MPSandMPO.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9665529",
   "metadata": {},
   "source": [
    "## Key takeaway:  The GHZ and Hadamard state have representations with low bond dimensions. This is what makes computations with them very efficient"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "054fba57",
   "metadata": {},
   "source": [
    "# Exercise 5: 15 minutes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "110ae1eb",
   "metadata": {},
   "source": [
    "Print the GHZ and Hadamard states for 30 qubits\n",
    "Compute the time it takes to do the inner product <hadamard_state(n), hadamard_state(n)>\n",
    "for n ranging up to 100\n",
    "Do this using the TN representation of the states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "657598b0",
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "do some stuff here to print/plot the timings",
     "output_type": "error",
     "traceback": [
      "do some stuff here to print/plot the timings",
      "",
      "Stacktrace:",
      " [1] error(s::String)",
      "   @ Base ./error.jl:35",
      " [2] top-level scope",
      "   @ In[71]:1",
      " [3] eval",
      "   @ ./boot.jl:368 [inlined]",
      " [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
      "   @ Base ./loading.jl:1428"
     ]
    }
   ],
   "source": [
    "error(\"do some stuff here to print/plot the timings\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4decb51",
   "metadata": {},
   "source": [
    "hint:\n",
    "Use ITensors.inner.\n",
    "If your code takes too long, construct the GHZ state using the solution code of Exercise 4."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de28eec8",
   "metadata": {},
   "source": [
    "# Exercise 6: 15 minutes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2538e35f",
   "metadata": {},
   "source": [
    "Construct a tensor network representation of the HWA quantum circuit from the lecture for 3 qubits.\n",
    "That is, the following sequence of gates (Ry, Rx, Ry, CNOTeven, CNOTodd)\n",
    "where, in the 3 qubit case, CNOTeven is a CNOT gate between qubits 1 and 2, and CNOTodd is a CNOT gate between qubits 2 and 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0dadd4d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "hwa_layer (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function hwa_layer(sites, angles::Vector{<:Number})\n",
    "  return nothing\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "700ab565",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "test_hwa_layer (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function test_hwa_layer()\n",
    "    N = 3\n",
    "    function elt(ϕ, inds)\n",
    "        s = siteinds(ϕ)\n",
    "        V = ITensor(1.)\n",
    "        for j=1:N\n",
    "            V *= (ϕ[j]*state(s[j],inds[j]))\n",
    "        end\n",
    "        return scalar(V)\n",
    "    end\n",
    "    sites = siteinds(\"S=1/2\", N)\n",
    "    in = MPS(sites, [\"Dn\", \"Up\", \"Dn\"])\n",
    "    hwa = hwa_layer(sites, [pi/2, pi/2, -pi/2, pi/2, pi/2, -pi/2, pi/2, pi/2, -pi/2])\n",
    "    out = apply(hwa, in)\n",
    "    @assert elt(out,[2,2,1]) ≈ (sqrt(2)+im*sqrt(2))/2\n",
    "    \n",
    "    return nothing\n",
    "end\n",
    "\n",
    "test_hwa_layer()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "188d759d",
   "metadata": {},
   "source": [
    "hint1: You can see how to chain elementary gates into a circuit here: https://github.com/ITensor/ITensors.jl/blob/main/examples/autodiff/ops/vqe.jl\n",
    "hint2: In the above file look at layer, variational_circuit, and lines 77-79"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "003a4139",
   "metadata": {},
   "source": [
    "## Key takeaway:  Every quantum circuit has a natural Tensor Network representation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a81455e",
   "metadata": {},
   "source": [
    "# Exercise 7: optional"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9acd467",
   "metadata": {},
   "source": [
    "Create a randomly initialized MPS, and orthogonalize it (put it in canonical form).\n",
    "Experiment with the effect on the bond dimension of the state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7b7bae7",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "Create a random state, and orthogonalize it.",
     "output_type": "error",
     "traceback": [
      "Create a random state, and orthogonalize it.",
      "",
      "Stacktrace:",
      " [1] error(s::String)",
      "   @ Base ./error.jl:35",
      " [2] top-level scope",
      "   @ In[72]:1",
      " [3] eval",
      "   @ ./boot.jl:368 [inlined]",
      " [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
      "   @ Base ./loading.jl:1428"
     ]
    }
   ],
   "source": [
    "error(\"Create a random state, and orthogonalize it.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e056788e",
   "metadata": {},
   "source": [
    "hint:\n",
    "call orthogonalize or orthogonalize!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4dd25b9",
   "metadata": {},
   "source": [
    "## key takeaway: The bond dimension of a random state can already be quite well reduced by using the canonical form, demonstrating how efficient the Tensor Network representation is."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d8ced3d",
   "metadata": {},
   "source": [
    "# Exercise 8: optional"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75188bf9",
   "metadata": {},
   "source": [
    "Verify that all the tensors to the left of the orthogonality center are left-orthogonal (contracting the tensor with itself along the site index, and the bond index pointing away from the center yields the identity matrix), and the tensors to the right of the orthogonality center are right-orthogonal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35d881b9",
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "do the contractions",
     "output_type": "error",
     "traceback": [
      "do the contractions",
      "",
      "Stacktrace:",
      " [1] error(s::String)",
      "   @ Base ./error.jl:35",
      " [2] top-level scope",
      "   @ In[73]:1",
      " [3] eval",
      "   @ ./boot.jl:368 [inlined]",
      " [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
      "   @ Base ./loading.jl:1428"
     ]
    }
   ],
   "source": [
    "error(\"do the contractions\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f481ffa2",
   "metadata": {},
   "source": [
    "hint:\n",
    "you can copy any tensor, and then prime some of its indexes. Contracting the tensor with its copy will contract along the unprimed indices. https://itensor.github.io/ITensors.jl/stable/ITensorType.html#Priming_and_tagging_ITensor"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4edeb9e",
   "metadata": {},
   "source": [
    "## Key takeaway: the orthogonality properties allow one to save computational resources, since many of the contractions occuring on measurement of local operators are guaranteed to be trivial"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.1",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}