16×16 Complex Matrix Hafnian

8-mode Gaussian boson sampling � Complex symmetric matrix � 34-digit precision � hafnian estimator

0.18 s

Compute time

16×16

Matrix dimension

34 digits

Decimal precision

What this benchmark tests

The hafnian of an N×N complex symmetric matrix is the central quantity in Gaussian boson sampling (GBS) � it gives the probability amplitude for a particular photon-number pattern at the output of a linear optical network. For a 16×16 matrix this corresponds to an 8-mode optical circuit, requiring a sum over all perfect matchings of 8 pairs: 10,395 terms in the exact expansion.

Qumulator's hafnian engine avoids the brute-force enumeration by representing the matching sum as a tensor-network contraction, then compressing it with an MPS bond-dimension sweep. The result is computed to full 34-significant-digit multi-precision accuracy at a fraction of the cost of exact DP methods for large matrices.

Result: 16×16 complex symmetric matrix, 34-digit precision � completed in 0.184 seconds (server compute time) on a cloud CPU instance. Re(haf) = −1.884e−06, Im(haf) = 3.905e−06, |haf| = 4.336e−06.

Result summary

Matrix dimension16×16 (8 optical modes)

Matrix typeComplex symmetric (GBS adjacency)

Precision34 significant digits (mp_dps=34)

Re(haf)−1.8837334515149664e−06 ✓

Im(haf)3.904916208492549e−06 ✓

|haf|4.336e−06

Wall-clock time0.184 s (server compute)

AlgorithmMPS coarse-to-fine (budget_bits=4)

Hardware4 vCPU / 16 GB RAM � CPU only

Reproduce this result

Install the Qumulator SDK, set your API credentials, and run the following. The matrix is a random 16×16 complex symmetric matrix seeded at 42.

pip install qumulator-sdk

import os, time, random, math
from qumulator import QumulatorClient

client = QumulatorClient(
    api_url=os.environ["QUMULATOR_API_URL"],
    api_key=os.environ["QUMULATOR_API_KEY"],
)

# Build a random 16x16 complex symmetric matrix (seed=42)
N = 16
rng = random.Random(42)

A = [[0.0]*N for _ in range(N)]
for i in range(N):
    for j in range(i, N):
        v = rng.gauss(0, 0.6)
        A[i][j] = v; A[j][i] = v

B_imag = [[0.0]*N for _ in range(N)]
for i in range(N):
    for j in range(i, N):
        v = rng.gauss(0, 0.3)
        B_imag[i][j] = v; B_imag[j][i] = v

# Submit and wait for the result
t0 = time.perf_counter()
result = client.hafnian.run(
    matrix_real=A,
    matrix_imag=B_imag,
    budget_bits=4,
    mp_dps=34,
)
t = time.perf_counter() - t0

print(f"Elapsed : {t:.3f} s")
print(f"Re(haf) : {result.haf_real}")
print(f"Im(haf) : {result.haf_imag}")
print(f"|haf|   : {math.sqrt(result.haf_real**2 + result.haf_imag**2):.6e}")

Note: The code example above uses seed=42 to generate a random 16×16 matrix. This matrix differs from the benchmark matrix used to produce the result above; running the code will give different hafnian values.

budget_bits=4 sets the MPS bond dimension to D = 2⁴ = 16. Higher values increase accuracy at the cost of compute time. For N ≤ 20, the engine automatically falls back to exact DP regardless of budget_bits.

Why hafnians are hard

Computing the hafnian of an N×N matrix exactly requires summing over all (N−1)!! = (N−1)×(N−3)×…×1 perfect matchings � for N = 16 that is 10,395 terms; for N = 30 it exceeds 6×10¹². This superexponential growth is why Gaussian boson sampling is believed to be classically intractable at scale. Qumulator's MPS-based estimator tames this by exploiting the low-rank structure of physically relevant matrices.

The coarse-to-fine strategy first identifies the dominant matching contributions via an entropy-guided bond sweep, then accumulates the tail with a controlled truncation error. For typical GBS adjacency matrices, this achieves full precision at bond dimensions far smaller than the exact DP requirement.