{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the VQE algorithm to work it is assumed that the Hamiltonian can be decomposed into a sum of Pauli constituents. Inspired by an implementation for 2 qubits on the quantum Slack channel I thought I'd try to implement a general decomposition algorithm in Julia.\n",
    "\n",
    "## The Pauli matrices\n",
    "\n",
    "From [Wikipedia](https://en.wikipedia.org/wiki/Pauli_matrices): \n",
    "\n",
    "> ...the Pauli matrices are a set of three 2 × 2 complex matrices which are Hermitian and unitary.\n",
    "\n",
    "> ...together with the identity matrix $I$ (sometimes considered as the zeroth Pauli matrix $σ_0$), the Pauli matrices form a basis for the real vector space of 2 × 2 Hermitian matrices. This means that any 2 × 2 Hermitian matrix can be written in a unique way as a linear combination of Pauli matrices, with all coefficients being real numbers.\n",
    "\n",
    "First we need to define the Pauli matrices.\n",
    "\n",
    "$$\n",
    "\\sigma_0 = I = \\begin{pmatrix} 1 & 0 \\\\ 0 & 1 \\end{pmatrix},\n",
    "\\sigma_1 = \\sigma_x = \\begin{pmatrix} 0 & 1 \\\\ 1 & 0 \\end{pmatrix},\n",
    "\\sigma_2 = \\sigma_y = \\begin{pmatrix} 0 &-i \\\\ i & 0 \\end{pmatrix},\n",
    "\\sigma_3 = \\sigma_z = \\begin{pmatrix} 1 & 0 \\\\ 0 &-1 \\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "In Julia code this becomes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "σ₀=[1 0; 0 1]\n",
    "σ₁=[0 1; 1 0]\n",
    "σ₂=[0 -im; im 0]\n",
    "σ₃=[1 0; 0 -1]\n",
    "\n",
    "σ = [σ₀, σ₁, σ₂, σ₃];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decomposing a 2x2 Hermitian Matrix\n",
    "\n",
    "Taking the wikipedia article at face value, for a 2x2 matrix $H$ the decomposition would be:\n",
    "\n",
    "$$\n",
    "H = \\sum_{i=0,x,y,z} a_{i} \\sigma_i,\n",
    "\\quad\n",
    "a_{i,j} = \\frac{1}{2} tr\\left( \\sigma_i H \\right)\n",
    "$$\n",
    "\n",
    "Let's try it with a random matrix with real entries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra # For the trace and kron operations\n",
    "\n",
    "temp = rand(2,2)\n",
    "H = 0.5(temp+temp')\n",
    "\n",
    "a = [tr(σ[i]*H) for i in 1:4]/2\n",
    "\n",
    "@assert H ≈ sum(a[i]*σ[i] for i in 1:4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So it works!\n",
    "\n",
    "What about a Hermitian matrix with complex entries?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp = rand(ComplexF64, 2,2)\n",
    "H = 0.5(temp+temp')\n",
    "\n",
    "a = [tr(σ[i]*H) for i in 1:4]/2\n",
    "\n",
    "@assert H ≈ sum(a[i]*σ[i] for i in 1:4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That works too!\n",
    "\n",
    "## Generalising to larger dimensions\n",
    "\n",
    "What about larger dimensions? From [this blog post](https://michaelgoerz.net/notes/decomposing-two-qubit-hamiltonians-into-pauli-matrices.html) I found the formula:\n",
    "\n",
    "$$\n",
    "H = \\sum_{i,j=0,x,y,z} a_{i,j} \\left( \\sigma_i \\otimes \\sigma_j \\right),\n",
    "\\quad\n",
    "a_{i,j} = \\frac{1}{4} tr\\left[\\left( \\sigma_i \\otimes \\sigma_j \\right) H \\right]\n",
    "$$\n",
    "\n",
    "This means that $a$ is a $4\\times4$ matrix containing all the combinations of two Pauli matrices. This translates very easily into Julia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "⊗ = kron\n",
    "\n",
    "temp = rand(ComplexF64, 4,4)\n",
    "H = 0.5(temp+temp')\n",
    "\n",
    "a = [tr((σ[i]⊗σ[j])*H) for i in 1:4, j in 1:4]/4\n",
    "\n",
    "@assert H ≈ sum(a[i,j]*(σ[i]⊗σ[j]) for i in 1:4, j in 1:4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far so good. Now we want to generalise it to $n$ dimensions. \n",
    "\n",
    "The pattern seems clear. To find the coefficients, $a_{1,2,3,...}$, calculate the tensor product of the $n$ combinations of Pauli matrices. Multiply by the matrix to be decomposed and then take the trace and divide by the normalisation factor, $2^n$.\n",
    "\n",
    "The result will be an $n$ dimensional array with $4^n$ entries. Julia's `CartesianIndices` type, gives us a way to iterate over all the dimensions in a linear way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" The tensor product of pauli matrices for the given index of arbitrary length [i,j,k,...]\"\"\"\n",
    "pauli_product(idx) = foldr(⊗, σ[idx[j]] for j in 1:length(idx))\n",
    "\n",
    "\"\"\" Returns multidimensional array of Pauli coefficients, depending on the size of H \"\"\"\n",
    "function pauli_decomposition(H)\n",
    "    @assert ishermitian(H)\n",
    "    \n",
    "    n = Int(log(2, size(H, 1)))\n",
    "    a = Array{Float64}(undef, fill(4, n)...) # 4^n array\n",
    "    for idx in CartesianIndices(a)\n",
    "        a[idx] = real(tr(pauli_product(idx)*H)) # Ignore numerical error in imaginary part\n",
    "    end\n",
    "    a/2^n # Normalisatin factor\n",
    "end\n",
    "\n",
    "\"\"\" Sums the pauli terms with the given coefficients \"\"\"\n",
    "function pauli_sum(a)\n",
    "    sum(a[idx]*pauli_product(idx) for idx in CartesianIndices(a))\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time to test it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.422849 seconds (1.31 M allocations: 105.146 MiB, 4.67% gc time)\n"
     ]
    }
   ],
   "source": [
    "n=5\n",
    "temp = rand(ComplexF64, 2^n,2^n)\n",
    "H = 0.5(temp+temp')\n",
    "\n",
    "@time a = pauli_decomposition(H)\n",
    "@assert H ≈ pauli_sum(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Works on my machine.. I've tested it up to n=7 and it works perfectly! \n",
    "\n",
    "(After n=7 the matrices the calculation takes quite a while and consumes the memory on my laptop.)\n",
    "\n",
    "We can also test it with a matrix with a known result: $$H = \\frac{1}{2}(I\\otimes I-\\sigma_x\\otimes \\sigma_x-\\sigma_y\\otimes \\sigma_y+\\sigma_z\\otimes \\sigma_z)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Array{Float64,2}:\n",
       " 0.5   0.0   0.0  0.0\n",
       " 0.0  -0.5   0.0  0.0\n",
       " 0.0   0.0  -0.5  0.0\n",
       " 0.0   0.0   0.0  0.5"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = [1  0  0 0;\n",
    "     0  0 -1 0;\n",
    "     0 -1  0 0;\n",
    "     0  0  0 1]\n",
    "\n",
    "a = pauli_decomposition(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which matches perfectly what we expected :)"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.5.1",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}