{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Benchmarking the Hafnian\n",
"This tutorial shows how to use Hafnian, a C (masquerading as Python) library to calculate the Hafnian."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The Hafnian\n",
"The hafnian of an $n$-by-$n$ symmetric matrix $A = A^T$ is defined as\n",
"\n",
"\\begin{align}\\label{eq:hafA}\n",
"\\text{haf}(A) = \\sum_{M \\in \\text{PMP}(n)} \\prod_{\\scriptscriptstyle (i, j) \\in M} A_{i, j}\n",
"% = \\sum_{\\mu \\in \\text{PMP}(n)} \\prod_{j=1}^n A_{\\mu(2j-1),\\mu(2j)}\n",
"\\end{align}\n",
"where PMP$(n)$ stands for the set of perfect matching permutations of $n$ (even) objects. \n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Using the library"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import the library in the usual way:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from thewalrus import hafnian"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To use it we need to pass square numpy arrays as arguments, thus we also must import NumPy:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The library provides functions to compute hafnians of real and complex matrices. The functions arguments must be passed as the NumPy arrays or matrices."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"945.0"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"size = 10\n",
"nth = 4\n",
"matrix = np.ones([size,size])\n",
"hafnian(matrix)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"945j"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"size = 10\n",
"nth = 4\n",
"matrix = 1j*np.ones([size,size])\n",
"hafnian(matrix)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Not surprisingly, the hafnian of a matrix containing only ones is given by $(n-1)!! = \\frac{n!}{(n/2)! 2^{n/2}}$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"945.0"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from math import factorial\n",
"factorial(size)/(factorial(size//2)*2**(size//2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that when doing floating point computations with large numbers, precision can be lost."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Benchmarking the performance of the code"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For sizes $n=2,30$ we will generate random symmetric matrices and measure the (average) amount of time it takes to calculate their hafnian. The number of samples for each will be geometrically distributed, with 1000 samples for size $n=2$ and 10 samples for $n=30$. The unitaries will be random Haar distributed."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"a0 = 1000.\n",
"anm1 = 2.\n",
"n = 20\n",
"r = (anm1/a0)**(1./(n-1))\n",
"nreps = [(int)(a0*(r**((i)))) for i in range(n)]"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[1000,\n",
" 721,\n",
" 519,\n",
" 374,\n",
" 270,\n",
" 194,\n",
" 140,\n",
" 101,\n",
" 73,\n",
" 52,\n",
" 37,\n",
" 27,\n",
" 19,\n",
" 14,\n",
" 10,\n",
" 7,\n",
" 5,\n",
" 3,\n",
" 2,\n",
" 2]"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nreps"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following function generates random Haar unitaries of dimensions $n$"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"from scipy import diagonal, randn\n",
"from scipy.linalg import qr\n",
"def haar_measure(n):\n",
" '''A Random matrix distributed with Haar measure\n",
" See https://arxiv.org/abs/math-ph/0609050\n",
" How to generate random matrices from the classical compact groups\n",
" by Francesco Mezzadri '''\n",
" z = (randn(n,n) + 1j*randn(n,n))/np.sqrt(2.0)\n",
" q,r = qr(z)\n",
" d = diagonal(r)\n",
" ph = d/np.abs(d)\n",
" q = np.multiply(q,ph,q)\n",
" return q"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's benchmark the scaling of the calculation with the matrix size"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2 0.00031961894035339354\n",
"4 0.00020129009357933858\n",
"6 0.0010492227440394878\n",
"8 0.001197782429781827\n",
"10 0.0013135018172087494\n",
"12 0.0017131505553255376\n",
"14 0.0013083321707589286\n",
"16 0.0031093134738431117\n",
"18 0.0061694171330700185\n",
"20 0.013150962499471812\n",
"22 0.032473312841879355\n",
"24 0.06553327595746075\n",
"26 0.15734933551989103\n",
"28 0.3482480730329241\n",
"30 0.7330920457839966\n",
"32 1.6255977153778076\n",
"34 3.6591072559356688\n",
"36 8.095563332239786\n",
"38 18.329116582870483\n",
"40 41.31634557247162\n"
]
}
],
"source": [
"import time\n",
"times = np.empty(n)\n",
"for ind,reps in enumerate(nreps):\n",
" start = time.time()\n",
" for i in range(reps):\n",
" size = 2*(ind+1)\n",
" nth = 1\n",
" matrix = haar_measure(size)\n",
" A = matrix @ matrix.T\n",
" A = 0.5*(A+A.T)\n",
" res = hafnian(A)\n",
" end = time.time()\n",
" times[ind] = (end - start)/reps\n",
" print(2*(ind+1), times[ind])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now plot the (average) time it takes to calculate the hafnian vs. the size of the matrix:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"%config InlineBackend.figure_formats=['svg']"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'Time in seconds for 4 threads')"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
"