Stable and Efficient Algorithms for the Fermion Determinant

Abstract: Some algorithms for the numerically exact treatment of fermion determinants are summarised. This is not supposed to be a review, rather a concise handbook. The audience is expected to have a basic understanding of how to put fermions on a computer. We primarily discuss different ways to work with the fermion matrix in the "sausage" (Green's function) formulation for quantum Monte Carlo (QMC). We emphasise the need for varied approaches in different space-time volume regimes. In particular, for small spatial volumes we describe a numerically stable method based on dense matrix operations. It is designed specifically to deal with very low temperature regimes. On the other hand, for (relatively) large volumes we describe a highly efficient and scalable sparse matrix approach.

• The paper presents a compendium of algorithms that balance accuracy and efficiency for computing fermion determinants in QMC simulations.
• The analysis rigorously determines runtime and memory bounds for different parameter regimes, guiding the choice between dense and sparse matrix methods.
• The work emphasizes stabilization strategies and recursive schemes to manage numerical instabilities and optimize computational performance.

Stable and Efficient Algorithms for the Fermion Determinant: An Expert Analysis

Overview and Motivation

This work systematically presents contemporary algorithmic strategies for the stable and efficient evaluation of fermion determinants, focusing on their application in the "sausage" (Green's function) formalism commonly used in quantum Monte Carlo (QMC) simulations. The distinctive structure of the paper is that of a technical handbook rather than a conventional review, offering explicit algorithmic choices and strategic advice tailored to different parameter regimes characterized by system size and temperature.

A central challenge in the simulation of fermionic systems is the computational cost and numerical instability associated with the determinant of large, time-ordered products of non-commuting matrices. The algorithms discussed aim to navigate the high-dimensional, ill-conditioned nature of these determinants, balancing speed and stability as required by the simulation parameters.

Problem Formulation: The “Sausage” Formalism

The “sausage” formalism introduces a reduction in the computational burden of QMC approaches by compressing the time evolution (usually discretized into $N_t$ slices) into a single $V\times V$ matrix:

$\hat M = 1 + \prod_t M_t$

Here, $M_t$ are the time-slice matrices encoding potential and kinetic terms, and $\hat M$ encapsulates the full time-ordered propagation. Observables and statistical weights in the QMC ensemble depend solely on the determinant (and, for measurements, on the inverse) of $\hat M$.

This formalism is utilized throughout the work, and algorithmic strategies are tailored to regimes dictated by spatial volume $V$ and effective inverse temperature $\beta E$ (with $E$ a characteristic energy scale).

Algorithmic Regimes and Complexity Bounds

The computational landscape is strictly mapped out in terms of lower and upper bounds for both runtime and memory, associated with (i) updates (such as in HMC or BSS Monte Carlo) and (ii) measurements (i.e., Green's function calculations). The analysis rigorously identifies:

• Update Complexity: The minimal cost is $\Omega(VN_t)$ (field updates), with full dense inversion at $V\times V$0. The HMC algorithm achieves this minimal scaling for global updates, while the BSS algorithm incurs higher cost due to its local update structure.
• Measurement Complexity: Calculating all single-fermion observables (i.e., the full Green's function) has a lower bound of $V\times V$1 in both runtime and memory; naive methods scale as $V\times V$2.

This explicit complexity analysis is essential for algorithm selection and highlights that runtime scaling is tightly coupled to the volume and number of time slices.

Algorithms by Parameter Regime

Small Spatial Volume

• $V\times V$3, Any $V\times V$4: Dense linear algebra (diagonalization, matrix products) can be used for both $V\times V$5 and its inverse. The stabilised "product-of-three-matrices" scheme (embedding $V\times V$6 with careful sorting and QR decomposition [Bauer:2020stable]) is adopted for high-$V\times V$7 regimes, guaranteeing numerical stability against the exponential ill-conditioning that emerges at low temperatures. This approach is especially valuable when eigenvalues cover a wide dynamic range.
• Numerical Results: The runtime is bounded by $V\times V$8 for a single update or observable, matching the cost of the dense operations; further increases to $V\times V$9 if many temporally-resolved measurements are required.
• Claims: The method prioritizes stability and explicit control over roundoff error, albeit at the cost of increased computational expense compared to less stable alternatives.

Intermediate Volume

• $\hat M = 1 + \prod_t M_t$0, Moderate/High Temperature: Sparse matrix representations of $\hat M = 1 + \prod_t M_t$1 (with $\hat M = 1 + \prod_t M_t$2 nonzero entries per row) are exploited. A product representation for the time-evolution operator ($\hat M = 1 + \prod_t M_t$3 approximated as a product of sparse matrices) leads to a runtime $\hat M = 1 + \prod_t M_t$4, which is optimal for this regime. LU or QR-decomposition is recommended for inversion and determinant evaluation, with periodic stabilization for large $\hat M = 1 + \prod_t M_t$5.
• For Lower Temperatures ($\hat M = 1 + \prod_t M_t$6): Stabilisation intervals are introduced (accumulating products up to an optimal $\hat M = 1 + \prod_t M_t$7 before reorthogonalizing), trading off between stability and computational cost.
• Claims: This regime demonstrates how algorithmic flexibility (e.g., switching the stabilization interval) permits stable simulation without exhaustive dense linear algebra.

Large Volume

• $\hat M = 1 + \prod_t M_t$8: Dense manipulation is intractable. The text recommends switching to pseudo-fermion methods, especially in the absence of a sign problem, recasting $\hat M = 1 + \prod_t M_t$9 as a Gaussian integral over auxiliary fields. Matrix-free iterative solvers are used exclusively, with complexity scaling linearly with $M_t$0 provided the matrix-vector multiplications remain sparse.
• Sign Problems: The methods become ineffective as the sign problem worsens, as the computational requirements then become exponential in system size ($M_t$1-hard [Troyer & Wiese, 2005]).

Product Accumulation and Observables

Evaluation of observables and force terms (as in HMC) requires fast accumulation and contraction of long sequences of $M_t$2 matrices. The work presents several "pre- and suffix" storage schemes and recursive strategies, optimized for either dense or sparse algebra, to reduce the computational overhead associated with repeated products. Notably:

• Tracing observables can be reduced to $M_t$3 via sparse contractions or as low as $M_t$4 for single-time observables.
• For highly resolved measurements, recursive doubling strategies reduce redundant computation, and higher-order $M_t$5-corrections are provided for mitigating trotterization errors.

Very Low Filling Regime

In the regime of low fermion density, the paper highlights a distinct simplification: interactions become perturbative, and reweighting approaches can be applied (e.g., single-particle rescaled simulation with reweighting, as in Eq.~(84)). This can significantly improve signal-to-noise for relevant observables and is especially valuable for systems at vanishing chemical potential.

Practical and Theoretical Implications

• The presented algorithms constitute a compendium of best-in-class strategies for efficient and stable computation of fermion determinants in QMC simulations.
• The theoretically precise complexity bounds (for both runtime and memory) provide a reference for future method development and algorithmic benchmarking.
• The identification of stability bottlenecks at low temperature and large system size directly motivates further work on stable matrix-free iterative solvers, possibly integrating block decomposition, stochastic estimators, or deflation techniques.
• For physically realistic problems with severe sign problems, these algorithms remain subject to fundamental intractability; thus, alternative approaches (e.g., contour deformations, Lefschetz thimble methods, or quantum simulation) warrant further research.

Future Developments in AI and Simulation

• Advances in AI-driven symbolic regression or automated differentiation could facilitate on-the-fly selection of optimal accumulation and stabilization strategies, adjusting dynamically to system parameters.
• Learning-augmented linear algebra (e.g., surrogate solvers based on deep learning for inverse estimation) is a potential avenue to further reduce computational bottlenecks.
• Integration of these techniques into widely used QMC and lattice field theory codes will drive adoption and enable direct comparison with alternative paradigms (e.g., tensor networks, quantum algorithms).

Conclusion

This compendium rigorously formalizes the trade-offs and strategies for calculating fermion determinants in QMC, tailored to distinct simulation regimes. Through the systematic collection and presentation of algorithms, optimal computational complexity, and stability requirements, it provides essential guidance for practitioners looking to maximize both accuracy and efficiency in fermionic lattice simulations. Future developments in AI and algorithmic design that address the sign problem and further scale these techniques will be critical for simulations of correlated quantum matter and lattice field theories.