[Paper Review] A faster hafnian formula for complex matrices and its benchmarking on a supercomputer

Andreas Björklund, Brajesh Gupt+1 more|arXiv (Cornell University)|May 31, 2018

Matrix Theory and Algorithms6 citations

TL;DR

This paper presents a new O(n³2ⁿ/²) algorithm for computing the hafnian of complex matrices using Schur decomposition, significantly accelerating exact hafnian evaluation. The method achieves a factor-of-n speedup over prior ring-based algorithms and demonstrates near-perfect weak scaling on the Titan supercomputer, enabling hafnian computation for matrices up to 56×56 with high accuracy and parallel efficiency.

ABSTRACT

We introduce new and simple algorithms for the calculation of the number of perfect matchings of complex weighted, undirected graphs with and without loops. Our compact formulas for the hafnian and loop hafnian of $n \ imes n $ complex matrices run in $O(n^3 2^{n/2})$ time, are embarrassingly parallelizable and, to the best of our knowledge, are the fastest exact algorithms to compute these quantities. Despite our highly optimized algorithm, numerical benchmarks on the Titan supercomputer with matrices up to size $56 \ imes 56$ indicate that one would require the 288000 CPUs of this machine for about a month and a half to compute the hafnian of a $100 \ imes 100$ matrix.

Motivation & Objective

• To develop a faster exact algorithm for computing the hafnian of complex matrices, a key quantity in quantum optics and combinatorics.
• To reduce the computational bottleneck in simulating Gaussian Boson Sampling, a quantum supremacy benchmark.
• To enable scalable, high-accuracy hafnian evaluation on supercomputers using hybrid MPI+OpenMP parallelism.
• To establish the practical limits of classical simulation for large-scale quantum sampling problems.

Proposed method

• The algorithm replaces the O(n⁴) dynamic programming tabulation in Cygan and Pilipczuk’s ring-based hafnian method with an O(n³) Schur decomposition of complex matrices.
• It uses compact formulas for the hafnian and loop hafnian of n×n complex matrices, derived via elementary linear algebra over the complex field.
• The method is embarrassingly parallelizable, enabling efficient distribution across thousands of CPU cores using hybrid MPI and OpenMP threading.
• Numerical benchmarks are performed on the Titan supercomputer using up to 288,000 CPUs, with weak scaling tests at matrix sizes up to 56×56.
• The algorithm maintains high numerical accuracy across distributed computations, with negligible error variation across different MPI/OpenMP process configurations.
• The approach is designed to minimize polynomial prefactors in exponential scaling, crucial for assessing quantum supremacy thresholds.

Experimental results

Research questions

• RQ1 Can the hafnian of a complex matrix be computed in O(n³2ⁿ/²) time using only standard linear algebra operations?
• RQ2 How does the proposed algorithm scale in practice on a large-scale supercomputer with hybrid parallelism?
• RQ3 What is the practical upper bound for exact hafnian computation on classical hardware, given current algorithmic and architectural constraints?
• RQ4 To what extent does the algorithm’s accuracy remain stable under distributed, non-commutative summation across multiple processors?

Key findings

• The proposed algorithm achieves O(n³2ⁿ/²) time complexity, representing a factor-of-n speedup over the previous best-known ring-based hafnian algorithm.
• On the Titan supercomputer, computing a 56×56 hafnian using 288,000 CPUs takes approximately 1.5 months, indicating that a 100×100 hafnian would require at least 1.5 months with perfect scaling.
• Weak scaling experiments show that total computation time plateaus as matrix size and processor count are doubled, confirming good weak scaling with a fixed problem-to-processor ratio.
• The algorithm exhibits near-perfect strong scaling across up to 288,000 CPUs, with computation time halving as processor count doubles.
• Percentage error in hafnian results remains negligible across different MPI and OpenMP process configurations, confirming numerical robustness under distributed summation.
• Despite optimization, exact hafnian computation for n > 60 remains intractable on classical hardware, with n=100 requiring ~2155 years on 16 processors.