Benchmarking the Hafnian¶

This tutorial shows how to use Hafnian, a C (masquerading as Python) library to calculate the Hafnian.

The Hafnian¶

The hafnian of an $$n$$-by-$$n$$ symmetric matrix $$A = A^T$$ is defined as

\begin{align}\label{eq:hafA} \text{haf}(A) = \sum_{M \in \text{PMP}(n)} \prod_{\scriptscriptstyle (i, j) \in M} A_{i, j} % = \sum_{\mu \in \text{PMP}(n)} \prod_{j=1}^n A_{\mu(2j-1),\mu(2j)} \end{align}

where PMP$$(n)$$ stands for the set of perfect matching permutations of $$n$$ (even) objects.

Using the library¶

Import the library in the usual way:

[1]:
from thewalrus import hafnian

To use it we need to pass square numpy arrays as arguments, thus we also must import NumPy:

[2]:
import numpy as np

The library provides functions to compute hafnians of real and complex matrices. The functions arguments must be passed as the NumPy arrays or matrices.

[3]:
size = 10
nth = 4
matrix = np.ones([size,size])
hafnian(matrix)
[3]:
945.0
[4]:
size = 10
nth = 4
matrix = 1j*np.ones([size,size])
hafnian(matrix)
[4]:
945j

Not surprisingly, the hafnian of a matrix containing only ones is given by $$(n-1)!! = \frac{n!}{(n/2)! 2^{n/2}}$$

[5]:
from math import factorial
factorial(size)/(factorial(size//2)*2**(size//2))
[5]:
945.0

Note that when doing floating point computations with large numbers, precision can be lost.

Benchmarking the performance of the code¶

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.

[6]:
a0 = 1000.
anm1 = 2.
n = 20
r = (anm1/a0)**(1./(n-1))
nreps = [(int)(a0*(r**((i)))) for i in range(n)]
[7]:
nreps
[7]:
[1000,
721,
519,
374,
270,
194,
140,
101,
73,
52,
37,
27,
19,
14,
10,
7,
5,
3,
2,
2]

The following function generates random Haar unitaries of dimensions $$n$$

[8]:
from scipy import diagonal, randn
from scipy.linalg import qr
def haar_measure(n):
'''A Random matrix distributed with Haar measure
See https://arxiv.org/abs/math-ph/0609050
How to generate random matrices from the classical compact groups
z = (randn(n,n) + 1j*randn(n,n))/np.sqrt(2.0)
q,r = qr(z)
d = diagonal(r)
ph = d/np.abs(d)
q = np.multiply(q,ph,q)
return q

Now let’s benchmark the scaling of the calculation with the matrix size

[9]:
import time
times = np.empty(n)
for ind,reps in enumerate(nreps):
start = time.time()
for i in range(reps):
size = 2*(ind+1)
nth = 1
matrix = haar_measure(size)
A = matrix @ matrix.T
A = 0.5*(A+A.T)
res = hafnian(A)
end = time.time()
times[ind] = (end - start)/reps
print(2*(ind+1), times[ind])
2 0.00031961894035339354
4 0.00020129009357933858
6 0.0010492227440394878
8 0.001197782429781827
10 0.0013135018172087494
12 0.0017131505553255376
14 0.0013083321707589286
16 0.0031093134738431117
18 0.0061694171330700185
20 0.013150962499471812
22 0.032473312841879355
24 0.06553327595746075
26 0.15734933551989103
28 0.3482480730329241
30 0.7330920457839966
32 1.6255977153778076
34 3.6591072559356688
36 8.095563332239786
38 18.329116582870483
40 41.31634557247162

We can now plot the (average) time it takes to calculate the hafnian vs. the size of the matrix:

[10]:
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_formats=['svg']
[11]:
plt.semilogy(2*np.arange(1,n+1),times,"+")
plt.xlabel(r"Matrix size $n$")
plt.ylabel(r"Time in seconds for 4 threads")
[11]:
Text(0, 0.5, 'Time in seconds for 4 threads')

The specs of the computer on which this benchmark was performed are:

[12]:
processor       : 0
vendor_id       : AuthenticAMD
cpu family      : 21
model           : 101
model name      : AMD A12-9800 RADEON R7, 12 COMPUTE CORES 4C+8G
stepping        : 1
microcode       : 0x6006118
cpu MHz         : 3992.225
cache size      : 1024 KB
physical id     : 0
siblings        : 4
core id         : 0
cpu cores       : 2
apicid          : 16
initial apicid  : 0
fpu             : yes
fpu_exception   : yes
cpuid level     : 13
wp              : yes

Note