2_Two/4_linear_algebra/1_linalg.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Algebra ([docs](https://docs.julialang.org/en/v1.0.0/stdlib/LinearAlgebra/))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since **quantum mechanics is linear algebra in disguise**, most quantum algorithms involve matrix operations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After `using LinearAlgebra`, **Julia speaks linear algebra fluently**.\n",
    "\n",
    "Performing linear algebra operations on a computer is, of course, an old problem. Lots of amazing libraries have been written - mostly in Fortran - which have been optimized over decades.\n",
    "\n",
    "Basically all high-level programming languages use these libraries, including R, Python, and Julia.\n",
    "\n",
    "Linear algebra in Julia is largely implemented by calling [BLAS](http://www.netlib.org/blas/)/[LAPACK](http://www.netlib.org/lapack/) functions. Sparse operations utilize functionality in [SuiteSparse](http://faculty.cse.tamu.edu/davis/suitesparse.html).\n",
    "\n",
    "As per default, Julia uses the [OpenBLAS](https://github.com/xianyi/OpenBLAS) implementation (BLAS, LAPACK, LIBM), which can be replaced by [Intel's MKL](https://software.intel.com/en-us/mkl) (BLAS, LAPACK) and [Intel's Math Library](https://software.intel.com/en-us/node/522653) (LIBM).\n",
    "\n",
    "**What is all this stuff?!?**\n",
    "\n",
    "* **BLAS**: a collection of low-level matrix and vector arithmetic operations (\"multiply two matrices\", \"multiply a matrix by vector\").\n",
    "* **LAPACK**:  a collection of higher-level linear algebra operations. Things like matrix factorizations (LU, LLt, QR, SVD, Schur, etc) that are used to do things like “find the eigenvalues of a matrix”, or “find the singular values of a matrix”, or “solve a linear system”.\n",
    "* **LIBM**: basic math functions like `sin`, `cos`, `sinh`, etcetera\n",
    "\n",
    "Sparse matrices are more difficult and there exist different collections of routines, one of which is **SuiteSparse**.\n",
    "\n",
    "**Why do I have to care?**\n",
    "\n",
    "* Switching from OpenBLAS to MKL can give you large speedups!\n",
    "* Since you might be leaving the world of Julia code, you loose easy inspectability and type genericity. The latter can be an issue for machine learning, as we'll discuss later in more detail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Taking linear algebra seriously"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia is [taking linear algebra seriously](https://www.youtube.com/watch?v=C2RO34b_oPM)! (see [here](https://github.com/JuliaLang/julia/issues/4774), and [here](https://github.com/JuliaLang/julia/issues/20978))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = rand(4,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "typeof(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Array{Float64, 2} === Matrix{Float64} # equivalent not just equal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "det(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inv(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rank(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get a vector as well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = rand(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "typeof(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Array{Float64,1} === Vector{Float64}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v^2 # can't square a vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v.^2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some things might be suprising"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1/v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But if it works, there is typically meaning to it. In this case it is calculating the [Moore-Penrose-Pseudoinverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse#Vectors) (`transpose(v)/sum(abs2,v)`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Identity matrix: `UniformScaling` operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A + 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A .+ 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `UniformScaling` operator **represents an identity matrix of any size** and is another great example of **duck typing**. It automatically gets loaded into scope when you do `using LinearAlgebra` and has the name `I`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although it never actually materializes a full identity matrix it behaves like one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A + 3I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I * A == A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, we can calculate things like, say, `A-b*I` without ever allocating a dense identity matrix, which would take up $\\mathcal{O}(n^2)$ memory.\n",
    "\n",
    "Let's benchmark the performance difference!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fullI = Matrix{Float64}(I, 4,4) # alternatively but slower, diagm(ones(4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fast(A) = A + 3*I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slow(A, fullI) = A + 3*fullI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function slower(A)\n",
    "    fullI = Matrix(1.0I, size(A)...)\n",
    "    A + 3*fullI\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using BenchmarkTools\n",
    "@btime fast($A);\n",
    "@btime slow($A, $fullI);\n",
    "@btime slower($A);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Utilizing matrix factorizations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Matrix factorizations (a.k.a. matrix decompositions)](https://en.wikipedia.org/wiki/Matrix_decomposition)\n",
    "are factorization of a matrix into a product of matrices, and are one of the central concepts\n",
    "in linear algebra.\n",
    "\n",
    "Making good use of matrix factorizations is crucial for efficient linear algebra operations.\n",
    "\n",
    "Example: Solving the linear system `Ax = b`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = rand(1:10, 5, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = rand(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve explicitly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inv(A)*b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve by left division `\\`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A\\b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using BenchmarkTools\n",
    "@btime inv($A)*$b; # it is (almost) never necessary to calculate the dense inverse\n",
    "@btime $A\\$b;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does Julia do to make this so much faster?\n",
    "\n",
    "It knows that it can perform the division much faster if it first [LU decomposes](https://en.wikipedia.org/wiki/LU_decomposition) `A`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lu(A)\\b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@btime lu($A)\\$b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's inspect the output of `lu(A)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lu(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "typeof(lu(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "supertype(LU)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### List of factorizations\n",
    "\n",
    "The following table summarizes the types of matrix factorizations that have been implemented in\n",
    "Julia. Details of their associated methods can be found in the [Standard Functions](https://docs.julialang.org/en/latest/stdlib/LinearAlgebra/#Standard-Functions-1) section\n",
    "of the Linear Algebra documentation.\n",
    "\n",
    "| Type               | Description                                                                                                    |\n",
    "|:------------------ |:-------------------------------------------------------------------------------------------------------------- |\n",
    "| `BunchKaufman`     | Bunch-Kaufman factorization                                                                                    |\n",
    "| `Cholesky`         | [Cholesky factorization](https://en.wikipedia.org/wiki/Cholesky_decomposition)                                 |\n",
    "| `CholeskyPivoted`  | [Pivoted](https://en.wikipedia.org/wiki/Pivot_element) Cholesky factorization                                  |\n",
    "| `LDLt`             | [LDL(T) factorization](https://en.wikipedia.org/wiki/Cholesky_decomposition#LDL_decomposition)                 |\n",
    "| `LU`               | [LU factorization](https://en.wikipedia.org/wiki/LU_decomposition)                                             |\n",
    "| `QR`               | [QR factorization](https://en.wikipedia.org/wiki/QR_decomposition)                                             |\n",
    "| `QRCompactWY`      | Compact WY form of the QR factorization                                                                        |\n",
    "| `QRPivoted`        | Pivoted [QR factorization](https://en.wikipedia.org/wiki/QR_decomposition)                                     |\n",
    "| `LQ`               | [QR factorization](https://en.wikipedia.org/wiki/QR_decomposition) of `transpose(A)`                           |\n",
    "| `Hessenberg`       | [Hessenberg decomposition](http://mathworld.wolfram.com/HessenbergDecomposition.html)                          |\n",
    "| `Eigen`            | [Spectral decomposition](https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix)                         |\n",
    "| `GeneralizedEigen` | [Generalized spectral decomposition](https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix#Generalized_eigenvalue_problem)                            |\n",
    "| `SVD`              | [Singular value decomposition](https://en.wikipedia.org/wiki/Singular_value_decomposition)                     |\n",
    "| `GeneralizedSVD`   | [Generalized SVD](https://en.wikipedia.org/wiki/Generalized_singular_value_decomposition#Higher_order_version) |\n",
    "| `Schur`            | [Schur decomposition](https://en.wikipedia.org/wiki/Schur_decomposition)                                       |\n",
    "| `GeneralizedSchur` | [Generalized Schur decomposition](https://en.wikipedia.org/wiki/Schur_decomposition#Generalized_Schur_decomposition) |\n",
    "\n",
    "(Taken from the Julia docs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@btime lu($A)\\$b\n",
    "@btime qr($A)\\$b\n",
    "@btime svd($A)\\$b;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the documentation (`?\\`) of the left division operator:\n",
    "\n",
    ">Matrix division using a polyalgorithm. For input matrices A and B, the result X is such that A*X == B when A is square. The solver that is used depends upon the structure of A. If A is upper or lower triangular (or diagonal), no factorization of A is required and the system is solved with either forward or backward substitution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is what the actual heuristic looks like (`@which`/`@edit` are your friends!)\n",
    "\n",
    "```julia\n",
    "function (\\)(A::AbstractMatrix, B::AbstractVecOrMat)\n",
    "    @assert !has_offset_axes(A, B)\n",
    "    m, n = size(A)\n",
    "    if m == n\n",
    "        if istril(A)\n",
    "            if istriu(A)\n",
    "                return Diagonal(A) \\ B\n",
    "            else\n",
    "                return LowerTriangular(A) \\ B\n",
    "            end\n",
    "        end\n",
    "        if istriu(A)\n",
    "            return UpperTriangular(A) \\ B\n",
    "        end\n",
    "        return lu(A) \\ B\n",
    "    end\n",
    "    return qr(A,Val(true)) \\ B\n",
    "end\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generically, a heuristic is implemented in `factorize`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "typeof(factorize(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "typeof(factorize(A+A'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fast linear algebra with multiple dispatch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, we've seen that Julia analyses the input matrix using some heuristic, factorizes it appropriately to then perform the calculation efficiently. \n",
    "\n",
    "But we can (and probably should) also be more explicit about our input to avoid this heuristic. We can encode the special structure of our matrix in a type such that we directly dispatch to the efficient method. Remember, the types decide which method is actually being run!\n",
    "\n",
    "There are many reasons to indicate what kind of matrix we have.\n",
    "\n",
    "* Don't rely on a heuristic. Not all methods have one!\n",
    "* The heurisitc comes with a small performance penalty.\n",
    "* The heurisitc isn't perfect and might fail to notice our matrix's special structure. Maybe because it's not known to base Julia. As we'll see later on, many external packages define additional special matrix types and efficient procedures for them.\n",
    "\n",
    "The following special matrix types are available out-of-the-box."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standard matrix types\n",
    "\n",
    ">[Matrices with special symmetries and structures](http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=3274)\n",
    "arise often in linear algebra and are frequently associated with various matrix factorizations.\n",
    "Julia features a rich collection of special matrix types, which allow for fast computation with\n",
    "specialized routines that are specially developed for particular matrix types.\n",
    ">\n",
    ">The following tables summarize the types of special matrices that have been implemented in Julia,\n",
    "as well as whether hooks to various optimized methods for them in LAPACK are available.\n",
    "\n",
    "| Type                          | Description                                                                                   |\n",
    "|:----------------------------- |:--------------------------------------------------------------------------------------------- |\n",
    "| [`Symmetric`](@ref)           | [Symmetric matrix](https://en.wikipedia.org/wiki/Symmetric_matrix)                            |\n",
    "| [`Hermitian`](@ref)           | [Hermitian matrix](https://en.wikipedia.org/wiki/Hermitian_matrix)                            |\n",
    "| [`UpperTriangular`](@ref)     | Upper [triangular matrix](https://en.wikipedia.org/wiki/Triangular_matrix)                    |\n",
    "| [`UnitUpperTriangular`](@ref) | Upper [triangular matrix](https://en.wikipedia.org/wiki/Triangular_matrix) with unit diagonal |\n",
    "| [`LowerTriangular`](@ref)     | Lower [triangular matrix](https://en.wikipedia.org/wiki/Triangular_matrix)                    |     |\n",
    "| [`UnitLowerTriangular`](@ref) | Lower [triangular matrix](https://en.wikipedia.org/wiki/Triangular_matrix) with unit diagonal |\n",
    "| [`UpperHessenberg`](@ref)     | Upper [Hessenberg matrix](https://en.wikipedia.org/wiki/Hessenberg_matrix)\n",
    "| [`Tridiagonal`](@ref)         | [Tridiagonal matrix](https://en.wikipedia.org/wiki/Tridiagonal_matrix)                        |\n",
    "| [`SymTridiagonal`](@ref)      | Symmetric tridiagonal matrix                                                                  |\n",
    "| [`Bidiagonal`](@ref)          | Upper/lower [bidiagonal matrix](https://en.wikipedia.org/wiki/Bidiagonal_matrix)              |\n",
    "| [`Diagonal`](@ref)            | [Diagonal matrix](https://en.wikipedia.org/wiki/Diagonal_matrix)                              |\n",
    "| [`UniformScaling`](@ref)      | [Uniform scaling operator](https://en.wikipedia.org/wiki/Uniform_scaling)                     |\n",
    "\n",
    "(Taken from the Julia docs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = Diagonal(1:5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ddense = Matrix(D) # same matrix but type doesn't indicate diagonal structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@btime $D*$b\n",
    "@btime $Ddense*$b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What method does it dispatch to?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@which D*b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@which Ddense*b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Dense Diagonal** (`Ddense*b`)\n",
    "\n",
    "```julia\n",
    "function (*)(A::AbstractMatrix{T}, x::AbstractVector{S}) where {T,S}\n",
    "    TS = promote_op(matprod, T, S)\n",
    "    mul!(similar(x,TS,axes(A,1)),A,x)\n",
    "end\n",
    "```\n",
    "\n",
    "**Diagonal** (`D*b`)\n",
    "```julia\n",
    "(*)(D::Diagonal, V::AbstractVector) = D.diag .* V\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fermions hopping on a chain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\mathcal{H} = -t\\sum_{\\langle i,j \\rangle} c_i^\\dagger c_j + \\mu \\sum_i n_i$$\n",
    "\n",
    "Here, $t$ is the hopping amplitude, $\\mu$ is the chemical potential, and $c, c^\\dagger$ are creation and annihilation operators.\n",
    "\n",
    "For simplicity, we'll consider **open boundary conditions** (not periodic), in which case the Hamiltonian is tridiagonal.\n",
    "\n",
    "Since the fermions are *not* interacting, we can work in the *single particle basis* and do not have to worry about how to construct a basis for the many-body Fock space.\n",
    "\n",
    "We use the canonical cartesian basis in which one uses $0$s to indicate empty sites and a $1$ for the particle's site, i.e. $|00100\\rangle$ represents the basis state which has the particle exclusively on the 3rd site.\n",
    "\n",
    "If you aren't familiar with second quantization just think of $\\mathcal{H}$ as any quantum mechanical operator that can be represented as a matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100 # number of sites\n",
    "t = 1\n",
    "μ = -0.5\n",
    "\n",
    "H = diagm(0 => fill(μ, N), 1 => fill(-t, N-1), -1 => fill(-t, N-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ψ = normalize(rand(N)); # some state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev(H, ψ) = ψ'*H*ψ # <φ|H|φ>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev(H, ψ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@btime ev($H, $ψ);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "typeof(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As long as the code is generic (respects the informal `AbstractArray` interface), we can use the same piece of code for completely different array types.\n",
    "\n",
    "Let's utilize the sparsity of `H` by indicating it through a type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using SparseArrays\n",
    "Hsparse = sparse(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@btime ev($Hsparse, $ψ);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a solid **30x speedup**!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our `H` isn't just sparse, but actually tridiagonal. Let's try to exploit that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Htri = Tridiagonal(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@btime ev($Htri, $ψ);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choosing the best type (and therewith an algorithm) can be tricky and one has to play around a bit. The good thing is that it's very easy to try out different types!\n",
    "\n",
    "Note that there are also great matrix types available in the ecosystem, see [JuliaMatrices](https://github.com/JuliaMatrices), for example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exact diagonalisation a.k.a Eigendecomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To diagonalize our dense \"Hamiltonian\", we simply call the built-in function `eigen`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vals, vecs = eigen(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ψ0 = vecs[:,1] # single-particle groundstate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev(H, ψ0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ev(H, ψ0) <= ev(H, ψ) # groundstate has the lowest energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "\n",
    "show_n_states = 3\n",
    "\n",
    "p = plot()\n",
    "for i in 1:show_n_states\n",
    "    plot!(p, abs2.(vecs[:,i]), xlab=\"site\", ylab=\"probability\", lab=\"n = $(i-1)\")\n",
    "end\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Periodic boundary conditions\n",
    "# T = copy(H)\n",
    "# T[1,end] = -t\n",
    "# T[end,1] = -t\n",
    "# vals, vecs = eigen(T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since Julia is using eigenproblem solvers from LAPACK (written in a low-level language) the code is, of course, **not generic**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best Julia can do, without implementing new functionality, is manually dispatch to the best LAPACK routine available.\n",
    "\n",
    "Hence, it won't work with most of our special matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigen(Htri);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we're lucky, someone has implemented a generic solver in Julia that works for a wider range of types. Example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Hbig = big.(H)\n",
    "eigen(Hermitian(Hbig));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using GenericLinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigen(Hermitian(Hbig));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arguably the most important matrix type in physics applications is a sparse matrix, i.e. `SparseMatrixCSC`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigen(Hsparse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's follow Julia's advice and take a look at [ARPACK.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Diagonalizing sparse matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[ARPACK.jl]() -  Wrapper to Fortran library [ARPACK](https://www.caam.rice.edu/software/ARPACK/) which implements **iterative** eigenvalue and singular value solvers. By far the most established sparse eigensolver.\n",
    "\n",
    "Julia implementations:\n",
    "\n",
    "* [ArnoldiMethod.jl](https://github.com/haampie/ArnoldiMethod.jl)\n",
    "* [KrylovKit.jl](https://github.com/Jutho/KrylovKit.jl)\n",
    "* [IterativeSolvers.jl](https://github.com/JuliaMath/IterativeSolvers.jl)\n",
    "* and more\n",
    "\n",
    "\n",
    "A key thing to remember is that while `eigen` is - up to numerical errors - exact, the methods in the packages above are iterative and approximative."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arpack uses a different name for the eigenvalue decomposition. They called it `eigs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Arpack\n",
    "λ, evs = eigs(Hsparse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "λ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For ArnoldiMethod, one has to go through a two-step process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using ArnoldiMethod\n",
    "decomp, history = partialschur(Hsparse)\n",
    "λ, evs = partialeigen(decomp);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "λ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In KrylovKit, they call the function `eigsolve`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using KrylovKit\n",
    "λ, evs = eigsolve(Hsparse);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "λ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Core messages of this Notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The standard libraries `LinearAlgebra` and `SparseArrays` make Julia speak linear algebra.\n",
    "* **Indicate properties and structure of a matrix**, like hermiticity or sparsity, through types. Fallback to generic types only if you run into method errors.\n",
    "* For **sparse matrix exact diagonalization**, ARPACK.jl is sort of a standard but there are great alternatives like ArnoldiMethods.jl."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# If time permits: Dude, I have a GPU!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make another case for *generic programming*, if you want to move the calculation to a GPU, chances are you only have to change the type of your matrix!\n",
    "\n",
    "Of course, this is only true if you have one :)\n",
    "\n",
    "```julia\n",
    "using CuArrays, BenchmarkTools\n",
    "\n",
    "A, B = rand(1000,1000), rand(1000,1000);\n",
    "Agpu, Bgpu = CuArray(A), CuArray(B);\n",
    "println(\"A*B (cpu)\")\n",
    "@btime A*B;\n",
    "# 9.702 ms (2 allocations: 7.63 MiB)\n",
    "\n",
    "println(\"A*B (gpu)\")\n",
    "@btime Agpu*Bgpu\n",
    "# 12.747 μs (19 allocations: 592 bytes)\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}