16×16 Complex Matrix Hafnian

8-mode Gaussian boson sampling · Complex symmetric matrix · 34-digit precision · KLT coarse-to-fine estimator

1.5 s

Wall-clock 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 KLT coarse-to-fine hafnian engine avoids the brute-force enumeration by representing the matching sum as a tensor-network contraction, then compressing it with a KLT-guided 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 1.537 seconds end-to-end on a cloud CPU instance. Re(haf) = 7.4023, Im(haf) = −0.9097, |haf| = 7.458.

Result summary

Matrix dimension16×16 (8 optical modes)

Matrix typeComplex symmetric (GBS adjacency)

Precision34 significant digits (mp_dps=34)

Re(haf)7.402315403377568 ✓

Im(haf)−0.909662989687380 ✓

|haf|7.458000e+00

Wall-clock time1.537 s

AlgorithmKLT MPS coarse-to-fine (budget_bits=4)

Hardware4 vCPU / 8 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

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: 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 KLT 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.