notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benasque Spring School on Near-Term Quantum Computing\n",
    "## Quantum Circuit simulation with Tensor Networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with creating a simple quantum circuit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Quac\n",
    "\n",
    "n_qubits = 8\n",
    "circ = Circuit(n_qubits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "push!(circ, H(1))\n",
    "circ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's create a simple quantum circuits with `10` qubits and `8` layers of random gates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "n = 10\n",
    "depth = 8\n",
    "\n",
    "# Define the quantum circuit\n",
    "circ = Circuit(n)\n",
    "gates = [Swap, CX, CY, CZ, CRx, CRy, CRz]\n",
    "\n",
    "Random.seed!(1)\n",
    "\n",
    "# We want to apply some entangling gates\n",
    "for _ in 1:depth\n",
    "    for i in 1:n\n",
    "        gate = rand(gates)\n",
    "        j = rand(setdiff(1:n, i))\n",
    "        push!(circ, gate(i, j))\n",
    "    end\n",
    "end\n",
    "\n",
    "circ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exact simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have created a circuit, let's convert it to a quantum tensor network with `Qrochet`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Qrochet\n",
    "\n",
    "U = Quantum(circ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Tenet\n",
    "\n",
    "tn = TensorNetwork(U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `CairoMakie` to visualize the tensor network:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using CairoMakie\n",
    "using GraphMakie\n",
    "using NetworkLayout\n",
    "\n",
    "plot(tn, layout=Stress())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get a `TensorNetwork` object that contains the tensor network representation of the circuit. It has $2N$ open indices, where $N$ is the number of qubits in the circuit, since we need to connect the input and output  qubits of the circuit to the tensor network. We will now contract the tensor network with the initial and final states to get the expectation value of the circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ψᵢₙ = Product(fill([1, 0], n))\n",
    "ψₒᵤₜ = Product(fill([1, 0], n))\n",
    "\n",
    "qtn = merge(Quantum(ψᵢₙ), U, Quantum(ψₒᵤₜ)') # < ψᵢₙ | U | ψₒᵤₜ >\n",
    "tn = TensorNetwork(qtn)\n",
    "\n",
    "plot(tn, layout=Stress())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `EinExprs.jl` to see the contraction path of the tensor network:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using EinExprs\n",
    "\n",
    "path = einexpr(tn, optimizer=Greedy())\n",
    "@show total_flops = mapreduce(flops, +, Branches(path))\n",
    "@show max_dims = maximum(ndims, Branches(path));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Largest Tensor in the computation: $2^{14} * 16$ bytes $= 0.262144$ MB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the contraction path:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(path, layout=Stress())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay! Let's contract the circuit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contract(TensorNetwork(qtn); path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Approximate simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a more complex circuit that is made of two-qubit local gates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Random.seed!(1)\n",
    "\n",
    "n = 50\n",
    "depth = 14\n",
    "\n",
    "# Define the quantum circuit\n",
    "gates = [Rxx, Ryy, Rzz]\n",
    "circ = Circuit(n)\n",
    "for i in 1:depth\n",
    "    for j in 1:n-1\n",
    "        rand_gate = rand(gates)\n",
    "        push!(circ, rand_gate(j, j+1; θ=rand()*2π))\n",
    "    end\n",
    "end\n",
    "\n",
    "print(circ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the `ψᵢₙ ` and `ψₒᵤₜ` product states:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ψᵢₙ = Product(fill([1, 0], n))\n",
    "ψₒᵤₜ = Product(fill([1, 0], n))\n",
    "\n",
    "plot(TensorNetwork(ψᵢₙ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = Quantum(circ)\n",
    "\n",
    "qtn = merge(Quantum(ψᵢₙ), U, Quantum(ψₒᵤₜ)') # < ψᵢₙ | U | ψₒᵤₜ >\n",
    "tn = TensorNetwork(qtn)\n",
    "\n",
    "plot(tn; layout=Spring())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = einexpr(tn, optimizer=Greedy())\n",
    "\n",
    "@show total_flops = mapreduce(flops, +, Branches(path))\n",
    "@show max_dims = maximum(ndims, Branches(path));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Largest Tensor in the computation: $2^{26} * 16$ bytes $= 1.07$ GB !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(path, layout=Stress())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@time c = contract(tn; path) # Exact result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much larger circuit! We need an approximation! We will use `Matrix Product States (MPS)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ψᵢₙ = convert(Chain, Product(fill([1, 0], n)))\n",
    "ψₒᵤₜ = Product(fill([1, 0], n))\n",
    "\n",
    "@time begin\n",
    "    # Apply gates using `evolve!` function  \n",
    "    for gate in circ\n",
    "        evolve!(ψᵢₙ, gate; maxdim=2, renormalize=false)\n",
    "    end\n",
    "\n",
    "    overlap(ψᵢₙ, ψₒᵤₜ) # <ψ_evolved | ψ_out>\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "\n",
    "@show fidelity = norm(ψᵢₙ);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can limit the bond dimension of the MPS to the computational complexity of the operation. A limit in the bond dimension can be thought as a limit in the entanglement of the state. The bond dimension is a trade-off between accuracy and computational cost. The larger the bond dimension, the more accurate the approximation, but the more expensive the computation.\n",
    "\n",
    "Let's use a reduced bond dimension $\\chi$ and compare it with the exact result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's evolve for different bond dimensions χ and compare to the exact result\n",
    "results = ComplexF32[]\n",
    "χ_list = [2, 4, 8, 16, 32, 64, 128]\n",
    "ψₒᵤₜ = Product(fill([1, 0], n))\n",
    "\n",
    "for χ in χ_list\n",
    "    @show χ\n",
    "    ψᵢₙ = convert(Chain, Product(fill([1, 0], n)))\n",
    "    canonize!(ψᵢₙ)\n",
    "\n",
    "    for gate in circ\n",
    "        evolve!(ψᵢₙ, gate; maxdim=χ, iscanonical=true, renormalize=true)\n",
    "    end\n",
    "\n",
    "    push!(results, only(overlap(ψᵢₙ, ψₒᵤₜ)))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = Figure()\n",
    "ax = Axis(fig[1, 1], xlabel=\"χ\", ylabel=\"|Amplitude|\", title=\"Expectation value\", xscale=log2)\n",
    "\n",
    "scatterlines!(ax, Float32.(χ_list), abs.(results), color=:red, linestyle=:dash, marker=:circle, label=\"Approximation\")\n",
    "lines!(ax, Float32.(χ_list), abs(only(c))*ones(length(χ_list)), color=:blue, linestyle=:solid, label=\"Exact value\")\n",
    "\n",
    "# Adjust the legend manually\n",
    "legend = Legend(fig[1, 2], ax, valign=:top)\n",
    "fig[1, 2] = legend\n",
    "\n",
    "fig # Display the figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time Evolution of a Matrix Product State\n",
    "\n",
    "We will use the `Time Evolving Block Decimation (TEBD)` algorithm to simulate the time evolution of a Matrix Product State (MPS) circuit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Transverse-field Ising Model\n",
    "\n",
    "The Hamiltonian of transverse field Ising model is given by:\n",
    "$H = -\\sum_{\\braket{i,j}} \\sigma_{X}^{(i)} \\sigma_{X}^{(j)} + \\lambda \\sum_i \\sigma_Z^{(i)}$\n",
    "\n",
    "We will simulate the time evolution of a quantum state under this Hamiltonian, and we will measure some `observables` of the system.\n",
    "\n",
    "In the `TEBD` algorithm we will use the Suzuki-Trotter decomposition to the first order, which approximates $e^{-iHt}$ as $e^{-iH_1 \\Delta t} e^{-iH_2 \\Delta t} ... e^{-iH_n \\Delta t}$, where $H = H_1 + H_2 + ... + H_n$, where $\\Delta t$ is the time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observables = Dense.([Z(5)])\n",
    "\n",
    "n = 10\n",
    "δₜ = 0.1\n",
    "timesteps = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start the quench from a product state where all qubits are in the $|+\\rangle$ state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ket₊ = 1/√2 * [1, 1]\n",
    "ψ = convert(Chain, Product(fill(ket₊, n)))\n",
    "\n",
    "plot(TensorNetwork(ψ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can find the energy of the system $|\\psi\\rangle=| + + \\dots + \\rangle$ by measuring the expectation value of the Hamiltonian $\\langle \\psi | H | \\psi \\rangle$.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "\n",
    "λ = 0.3\n",
    "\n",
    "ket₊ = 1/√2 * [1, 1]\n",
    "ψ = convert(Chain, Product(fill(ket₊, n)))\n",
    "\n",
    "function energy_Ising(ψ, λ)\n",
    "    id = Matrix(1.0I, 2, 2)\n",
    "    H_ij = -1 * kron(Matrix(X()), Matrix(X())) + λ * (kron(Matrix(Z()), id) + kron(id, Matrix(Z())))\n",
    "    H_ij = reshape(H_ij,2,2,2,2)\n",
    "\n",
    "    e_i = []\n",
    "    for i in 1:n-1\n",
    "        H_gate = Dense(Qrochet.Operator(), H_ij; sites=Site[Site(i), Site(i+1), Site(i, dual=true), Site(i+1, dual=true)])\n",
    "        push!(e_i, only(expect(ψ, [H_gate])))\n",
    "    end\n",
    "    return e_i\n",
    "end\n",
    "\n",
    "energy_Ising(ψ, λ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define two gates that will be used in the TEBD. The `trotter_XX` gate is a Trotterized version of the XX gate, and the `trotter_Z` gate is a Trotterized version of the Z gate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function trotter_XX(i,j; δₜ=δₜ)\n",
    "    mat = kron(Matrix(X()), Matrix(X()))\n",
    "    mat = cis(δₜ * mat)\n",
    "    mat = reshape(mat,2,2,2,2)\n",
    "    Dense(Qrochet.Operator(), mat; sites=Site[Site(i), Site(j), Site(i, dual=true), Site(j, dual=true)])\n",
    "end\n",
    "\n",
    "function trotter_Z(i; λ, δₜ=δt)\n",
    "    mat = Matrix(Z())\n",
    "    mat = cis(- λ * δₜ * mat)\n",
    "\n",
    "    Dense(Qrochet.Operator(), mat; sites=Site[Site(i), Site(i, dual=true)])\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can start the time evolution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ket₊ = 1/√2 * [1, 1]\n",
    "ψ = Product(fill(ket₊, n))\n",
    "ψ = convert(Chain, ψ)\n",
    "\n",
    "# Time evolution with the TEBD algorithm:\n",
    "expected_vals = Float64[]\n",
    "@time for it in 1:timesteps\n",
    "    # First apply the gates on even qubits\n",
    "    for (i,j) in Iterators.filter(==(2) ∘ length, Iterators.partition(1:n,2))\n",
    "        evolve!(ψ, trotter_XX(i,j; δₜ))\n",
    "    end\n",
    "\n",
    "    # Then apply the gates on odd qubits\n",
    "    for (i,j) in Iterators.filter(==(2) ∘ length, Iterators.partition(2:n,2))\n",
    "        evolve!(ψ, trotter_XX(i,j; δₜ))\n",
    "    end\n",
    "\n",
    "    # Finally apply the one-qubit gates\n",
    "    for i in 1:n\n",
    "        evolve!(ψ, trotter_Z(i; λ=0.3, δₜ))\n",
    "    end\n",
    "\n",
    "    # Compute the expectation value of the observable\n",
    "    push!(expected_vals, abs(only(expect(ψ, observables))))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the expectation value:\n",
    "fig = Figure(fontsize=18)\n",
    "ax = Axis(fig[1, 1], xlabel=\"time\", ylabel=\"Observable\")\n",
    "lines!(ax, (1:timesteps) .* δₜ, expected_vals, color=:red, linestyle=:solid, marker=:circle, label=\"Exact result\", fontsize=12)\n",
    "\n",
    "legend = Legend(fig[1, 2], ax, valign=:top)\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But what if we had a large system? We can limit the bond dimension $\\chi$!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expectedvals = Dict()\n",
    "\n",
    "ket₊ = 1/√2 * [1, 1]\n",
    "ψ_original = convert(Chain, Product(fill(ket₊, n)))\n",
    "\n",
    "observables = Dense.([Z(5)])\n",
    "\n",
    "χ_list = [1, 2, 4, 8, 16]\n",
    "\n",
    "for χ in χ_list\n",
    "    # Reset ψ to its original state for each χ iteration\n",
    "    ψ = deepcopy(ψ_original)\n",
    "    canonize!(ψ)\n",
    "    expectedvals_χ = Float64[]\n",
    "\n",
    "    @show χ\n",
    "\n",
    "    for it in 1:timesteps\n",
    "        for (i,j) in Iterators.filter(==(2) ∘ length, Iterators.partition(1:n,2))\n",
    "            evolve!(ψ, trotter_XX(i,j; δₜ); maxdim=χ, iscanonical=true)\n",
    "        end\n",
    "        for (i,j) in Iterators.filter(==(2) ∘ length, Iterators.partition(2:n,2))\n",
    "            evolve!(ψ, trotter_XX(i,j; δₜ); maxdim=χ, iscanonical=true)\n",
    "        end\n",
    "        for i in 1:n\n",
    "            evolve!(ψ, trotter_Z(i; λ=0.3, δₜ); maxdim=χ, iscanonical=true)\n",
    "        end\n",
    "        push!(expectedvals_χ, abs(only(expect(ψ, observables))))\n",
    "    end\n",
    "\n",
    "    expectedvals[χ] = expectedvals_χ\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using CairoMakie\n",
    "\n",
    "# Assuming `expectedvals` and `maxdims` are defined as before\n",
    "\n",
    "fig = Figure()\n",
    "ax = Axis(fig[1, 1], xlabel=\"Timestep\", ylabel=\"Observable\", title=\"Observables for Different maxdim\")\n",
    "\n",
    "# Reverse the order of these arrays to match the reversed maxdims\n",
    "color_palette = reverse([:red, :green, :blue, :purple, :orange])\n",
    "line_styles = reverse([(:solid, 2), (:dash, 2), (:dot, 2), (:dashdot, 2), (:dashdotdot, 2)])\n",
    "markers = reverse([:circle, :rect, :diamond, :cross, :star5])\n",
    "\n",
    "for (index, χ) in enumerate(χ_list)\n",
    "    # Convert expectedvals[maxdim] to a Float32 array explicitly\n",
    "    y_values = Float32.(expectedvals[χ])\n",
    "    x_values = Float32.(1:length(y_values)) # Ensure x values are also Float32 for consistency\n",
    "\n",
    "    lines!(ax, x_values, y_values,\n",
    "          color=color_palette[index],\n",
    "          linestyle=line_styles[index][1],\n",
    "          linewidth=line_styles[index][2],\n",
    "          marker=markers[index],\n",
    "          markersize=10,\n",
    "          label=\"χ = $χ\")\n",
    "end\n",
    "\n",
    "legend = Legend(fig[1, 2], ax, valign=:top)\n",
    "fig[1, 2] = legend\n",
    "\n",
    "fig # Display the figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise !\n",
    "\n",
    "Use the `TEBD` algorithm to simulate the time evolution of a Matrix Product State (MPS) under the XXZ Hamiltonian: $H = -\\frac{1}{2}\\sum^{N}_{i}(\\sigma_{i}^{x}\\sigma_{i+1}^{x}+\\sigma_{i}^{y}\\sigma_{i+1}^{y}+\\Delta\\sigma_{i}^{z}\\sigma_{i+1}^{z})$ with $\\Delta=0$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Measure the magnetization of the system during the quench, which is defined as $M = \\frac{1}{N} \\sum_{i=1}^{N}\\langle\\sigma_{i}^{z}\\rangle$.\n",
    "Start with a system where all spins are in the $|\\uparrow\\rangle$ state. Use the `TEBD` algorithm with a bond dimension $\\chi$ to simulate the time evolution of the system.\n",
    "\n",
    "Tip: Start with a small system (e.g. $N=10$) and a small number of time steps to test your implementation.\n",
    "\n",
    "- Plot the average magnetization of the system during the quench using `Makie`.\n",
    "- How does the magnetization of the system change with the bond dimension $\\chi$?\n",
    "- Do a 2D plot of the magnetization per site $m_{i}$ of the system as a function of time and site. See https://docs.makie.org/stable/reference/plots/heatmap/#example_4a31556a for an example of how to plot a heatmap with `Makie`.\n",
    "- Compute the energy and verify that it is conserved during the time evolution. How does that change if we use a larger time step $\\delta t$?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.2",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}