# Numerical Challenges in Lattice QCD 2022

Europe/Zurich
Landhaus Nordhelle

#### Landhaus Nordhelle

Zum Koppenkopf 3, 58540 Meinerzhagen, Germany
Description

The purpose of this workshop is to bring together mathematicians from numerical mathematics and physicists from lattice QCD theory to exchange ideas and approaches for the great computational challenges in lattice QCD simulations. This includes linear system solvers for discretizations of the Dirac matrix, integrators for the Hybrid Monte-Carlo algorithm of lattice QCD and the computation and sampling of disconnected contributions. The workshop format puts an emphasis on stimulating interactions between mathematicians and physicists. The talks are all plenary, 30' long separated by 10' for discussions and with a coffee break of 20' every two talks for more discussions. An extended lunch break allows for longer discussions and for inspiring open-air activities.

The main topics are

1. Numerical solvers

2. Numerical integration (HMC + machine learning (based on "flow"))

3. Stochastic trace estimation

The workshop is organized as a satellite meeting to Lattice 2022, Bonn, August 08-13, and will take place from August 15-17 in Meinerzhagen, 75km from the lattice venue (1.5h drive by car or 2h by train), at the Landhaus Nordhelle, in the beautiful Sauerland, a rural, hilly area spreading across most of the south-eastern part of North Rhine-Westphalia. A contingent of single rooms has been blocked for the duration of the workshop (arrival on Sunday, departure on Wednesday).

This workshop is an event of the Research Unit 5269, working on future methods for studying confined gluons in QCD, funded by the German Research Foundation (DFG) under FOR5269, and is supported by the LatticeHadrons network of the STRONG-2020 Integrating Activity for Advanced Communities and the STIMULATE - European Joint Doctorates Program.

We are looking forward to welcome you in Meinerzhagen!

Registration
Numerical Challenges in Lattice QCD 2022 Registration Form
Participants
• Andreas Frommer
• Anthony Kennedy
• Artur Strebel
• Constantia Alexandrou
• Daniel Kressner
• Daniela Ebeling
• Eike Mueller
• Francesco Giacomo Knechtli
• Giannis Koutsou
• Gunnar Bali
• Gustavo Ramirez
• Henning Leemhuis
• Jacob Finkenrath
• Jakob Simeth
• James Brannick
• Jose Miguel Jimenez Merchan
• Juan Andres Urrea Nino
• Karsten Kahl
• Kevin Schäfers
• Kirk Soodhalter
• Liam Burke
• Lorenzo Barca
• Michael Günther
• Michael Peardon
• Michele Della Morte
• Michèle Wandelt
• Phiala Shanahan
• Piotr Bialas
• Roman Höllwieser
• Sara Collins
• Simone Bacchio
• Stefan Schaefer
• Tom Asmussen
• Tomasz Korzec
• Travis Whyte
• Urs Wenger
• Walter Wilcox
Contact: Daniela Ebeling
• Monday, August 15
• 1
Opening
Speaker: Francesco Giacomo Knechtli (Bergische Universitaet Wuppertal (DE))
• Multi-Level
Convener: Andreas Frommer
• 2
Temporal factorization of the Wilson fermion determinant and multi-level integration

When lattice QCD is formulated in sectors of fixed quark numbers, the canonical fermion determinants can be expressed explicitely in terms of transfer matrices defined at fixed time. This in turn provides a complete factorization of the fermion determinants in temporal direction. In this talk I describe this factorization for Wilson-type fermions and present explicit constructions of the transfer matrices. Possible applications of the factorization include the construction of improved estimators for generic n-point correlation functions and multi-level integration schemes. The latter is of particular interest for the calculation of disconnected correlation functions.

Speaker: Urs Wenger (Universität Bern (CH))
• 3
AMG+: Extended Multilevel Principles

By employing new extended multilevel hierarchy construction principles, AMG can be applied to many new types of problems. The extended principles include general rules for choosing relaxation, constructing the coarse-level variables, the coarse-to-fine interpolation, and coarse-level equations, and a quantitative performance predictor of the multi-level cycle convergence rate, called the mock cycle. As an illustrative example, a fast solver for the highly indefinite 1D Helmholtz equation is developed.

Speaker: Prof. James Brannick (Penn State)
• 10:20 AM
Coffee Break
• Numerical Solvers
Convener: Andreas Frommer
• 4
Disconnected Diagrams and High-Degree GMRES Polynomials

We build upon our lattice QCD POLY and HFPOLY methods by using high-degree polynomials in the context of noisy disconnected diagram evaluations. Using a new, stable form of the GMRES polynomial, we obtained a subtracted error reduction in the Wilson-Dirac scalar operator on the order of a tenth the error of the non-subtracted measurement. The associated trace of subtracted high-degree polynomials can be evaluated stochastically through a multilevel polynomial cascade method using deflation as well as a new and efficient double-polynomial construction. For demonstration purposes, our methods were tested at kappacrit in the quenched approximation.

Speaker: Walter Wilcox
• 5
Analysis of block GMRES using a $^\ast$-algebra-based approach

We discuss the challenges of extending convergence results of classical Krylov subspace methods to their block counterparts and propose a new approach to this analysis. Block KSMs such as block GMRES are generalizations of classical KSMs, and are meant to iteratively solve linear systems with multiple right-hand sides (a.k.a. a block right-hand side) all-at-once rather than individually. However, this all-at-once approach has made analysis of these methods more difficult than for classical KSMs because of the interaction of the different right-hand sides. We have proposed an approach built on interpreting the coefficient matrix and block right-hand side as being a matrix and vector over a $^\ast$-algebra of square matrices. This allows us to sequester the interactions between the right-hand sides into the elements of the $^\ast$--algebra and (in the case of GMRES) extend some classical GMRES convergence results to the block setting. We then discuss some challenges which remain and some ideas for how to proceed.

This is joint work with Marie Kubiínová from Czech Academy of Sciences, Institute of Geonics, Ostrava, Czech Republic (formerly)

[1] Marie Kubínová and Kirk M. Soodhalter. Admissible and attainable convergence behavior of block Arnoldi and GMRES. SIAM Journal on Matrix Analysis and Applications, 41 (2), pp. 464-486, 2020.

Speaker: Kirk Soodhalter (Trinity College Dublin)
• Stochastic Trace Estimation
Convener: Stefan Schaefer
• 6
Estimating the trace of a matrix function via randomized low-rank approximation

The development and analysis of randomized techniques for estimating the trace of a large matrix is a topic with a long history that has recently seen increased attention and new exciting developments. This talk is concerned with estimating the trace of a matrix function $f(A)$ for some function $f$ and a symmetric matrix $A$. Such problems feature prominently in a variety of applications in, e.g., scientific computing, statistical learning, and computational physics. The recently developed estimator Hutch++ uses randomized low-rank approximation to speed up classical stochastic trace estimation, sometimes dramatically. We present novel error bounds that are turned into an improved and adaptive variant of Hutch++. For an important subclass of functions, we derive and analyze a new randomized approach that can be significantly faster than existing methods.

Speaker: Daniel Kressner
• 7
Interpolation as a means of shift selection in multilevel Monte Carlo with lattice displacements

The calculation of disconnected diagram contributions to physical signals is
a computationally expensive task in Lattice QCD. To extract the physical
signal, the trace of the inverse Lattice Dirac operator, a large sparse matrix,
must be stochastically estimated. Because the variance of the stochastic es-
timator is typically large, variance reduction techniques must be employed.
Multilevel Monte Carlo (MLMC) methods reduce the variance of the trace
estimator by utilizing a telescoping sequence of estimators. Frequency Split-
ting is one such method that uses a sequence of inverses of shifted operators
to estimate the trace of the inverse lattice Dirac operator, however there is no
a priori way to select the shifts that minimize the cost of the multilevel trace
estimation. We present a sampling and interpolation scheme
that is able to predict the variances associated with Frequency Splitting un-
der displacements of the underlying space time lattice. The interpolation
scheme is able to predict the variances to high accuracy and therefore choose
shifts that correspond to an approximate minimum of the cost for the trace
estimation. We show that Frequency Splitting with the chosen shifts dis-
plays significant speedups over multigrid deflation, and that these shifts can
be used for multiple configurations within the same ensemble with no penalty
to performance.

Speaker: Dr Travis Whyte (College of William & Mary)
• 3:50 PM
Coffee Break
• Multi-Level
Convener: Stefan Schaefer
• 8
Multilevel Monte Carlo for quantum mechanics on a lattice

Monte Carlo simulations of quantum field theories on a lattice become increasingly expensive as the continuum limit is approached since the cost per independent sample grows with a high power of the inverse lattice spacing. Simulations on fine lattices suffer from critical slowdown, the rapid growth of autocorrelations in the Markov chain with decreasing lattice spacing a. This causes a strong increase in the number of lattice configurations that have to be generated to obtain statistically significant results. We discuss hierarchical sampling methods to tame this growth in autocorrelations; combined with multilevel variance reduction techniques, this significantly reduces the computational cost of simulations. We recently demonstrated the efficiency of this approach for two non-trivial model systems in quantum mechanics in https://arxiv.org/abs/2008.03090. This includes a topological oscillator, which is badly affected by critical slowdown due to freezing of the topological charge. On fine lattices our methods are several orders of magnitude faster than standard, single level sampling based on Hybrid Monte Carlo. For very high resolutions, multilevel Monte Carlo can be used to accelerate even the cluster algorithm, which is known to be highly efficient for the topological oscillator. Performance is further improved through perturbative matching. This guarantees efficient coupling of theories on the multilevel lattice hierarchy, which have a natural interpretation in terms of effective theories obtained by renormalisation group transformations.

Speaker: Eike Mueller (University of Bath)
• Stochastic Trace Estimation
Convener: Stefan Schaefer
• 9
Properties of the $\eta$ and $\eta^\prime$ mesons

We summarize our results for the $\eta$ and $\eta^\prime$ masses and
their four independent decay constants at the physical point as well as
their anomalous gluonic matrix elements $a_{\eta^{(\prime)}}$​.

The computation employs twenty-one $N_f=2+1$ Coordinated Lattice
Simulations (CLS) ensembles with non-perturbatively improved Wilson
fermions at four different lattice spacings and along two trajectories
in the quark mass plane, including one ensemble very close to physical
quark masses. We give details on the evaluation of the disconnected
contributions to the pseudoscalar and axialvector matrix elements. The
direct determination of both allows us to investigate the singlet PCAC
relation and to study the QCD scale dependence of the singlet currents.

Speaker: Jakob Simeth (University of Regensburg)
• 7:00 PM
Bowling Tournament
• Tuesday, August 16
• Numerical Integration
Convener: Prof. Michael Günther (University of Wuppertal)
• 10
HMC on Riemannian Manifolds

I shall give a pedagogical introduction to the Hamiltonian/Hybrid Monte Carlo algorithm (HMC) on Riemannian manifolds. I will explain how Hamiltonian systems are most naturally formulated in terms of Hamiltonian vector fields $\hat X$ on symplectic manifolds: the relationship between commutators of Hamiltonian vector fields and Poisson brackets, $[\hat X,\hat Y] = \widehat{\{X,Y\}}$; why a conserved shadow Hamiltonian exists; and why the natural Riemannian volume element $\sqrt{\det g}$ appears automatically. I will show that the symplectic integrator step $e^{\hat T\delta\tau}$ is a geodesic on a Riemannian manifold, and how this corresponds to matrix exponentiation on Lie groups. Finally, time permitting, I will explain how this leads to a practical HMC algorithm on symmetric spaces (and slightly more generally on homogeneous reductive manifolds) such as ${\mathbb CP}^n$ or $\mathbb S^2$ using Hamiltonian reduction.

Speaker: Anthony Kennedy
• 11
The potential of Padé approximations for molecular dynamics simulation
Speaker: Kevin Schaefers
• 10:20 AM
Coffee Break
• Machine Learning
Convener: Prof. Michael Günther (University of Wuppertal)
• 12
Machine learning for ensemble generation

I will discuss the use of machine learning methods to accelerate algorithms for gauge field generation, in particular via flow models.

Speaker: Phiala Shanahan
• 13
Continuous Normalizing Flow for Lattice QCD based on Trivializing Maps

In this presentation we will show the connection between Continuous Normalizing Flows (CNF) and Trivializing Maps by Luescher. Based on the latter, we will construct a CNF that can be trained to simulate lattice field theories. We discuss strategies to train the CNF for 2D and 4D SU(3) pure-gauge theories.

Speaker: Simone Bacchio
• Machine Learning
Convener: Francesco Giacomo Knechtli (Bergische Universitaet Wuppertal (DE))
• 14
Tackling critical slowing down using global correction steps with equivariant flows within the 2D Schwinger model

Discretisation of gauge theories are an elegant and successful way to solve them via supercomputers. To obtain results at the continuum, the discretised model is simulated via Monte Carlo simulations at fixed physics at different lattice spacings and then extrapolated to the continuum. In many cases the major systematic effect of the obtained result is given by the extrapolation error. To minimise this error, simulations at finer lattice spacing are necessary, which are often prevented by increasing autocorrelation times, caused by a freezing of extensive quantities, such as the topological charge.

We will discuss a potential method to overcome topological freezing in gauge models with fermions. The method combines flow-based generative models for local gauge field updates and hierarchical updates of the factorised fermion determinant. The flow-based generative models are restricted to proposing updates to gauge-fields within subdomains, thus keeping training times moderate while increasing the global volume. We apply our method to the 2-dimensional (2D) Schwinger model and show that sampling of topological sectors can be achieved also at fine lattices.
Moreover, we will discuss the potential of combining the correction steps with the Hybrid Monte Carlo method.

Speaker: Jacob Finkenrath
• 15

ecently a machine learning approach to Monte-Carlo simulations called Neural Markov Chain Monte-Carlo (NMCMC) is gaining traction. In its most popular form it uses neural networks to construct normalizing flows which are then trained to approximate the desired target distribution. The training is done using some form of gradient descent so gradient estimation is necessery. In my talk I will review several gradient estimators that can be used for this purpose and discuss their pros and cons in terms of required computational time and quality of the training.

Speaker: Prof. Piotr Białas (Jagiellonian University)
• 3:50 PM
Coffee Break
• Podium Discussion
Convener: Francesco Giacomo Knechtli (Bergische Universitaet Wuppertal (DE))
• 16
Podium discussion contribution

Some selected questions:Critical slowing down with $a\to 0$.
- Rounding issues on large volumes.
- Multiscale approaches to increase signal over noise
- Benefits of (approximate eigenvectors) and how to obtain these
- Approaches to excited state contributions: multi-state fits, GEVP, smearing.
- efficient stochastic estimation of traces, all-to-all propagators, perambulators

These are based on an analysis of challenges: * Large lattice volumes (due to small $a$ and $M_\pi$)
- Computational cost does not nexessarily translate into better signal/noise.
- critical slowing down (in general with small $a$)
- precision issues
- how to maintain good acceptance in HMC?
- algorithms that scale like $V^2$, e.g., number of eigenvectors in estimators/deflation, algorithms that scale like $V/a^2$ (smearing).
- How to beat noise of disconnected contributions on large volumes? * $n$-point functions with $n>3$:
- QED corrections to QCD, $K\to\pi\pi$ etc.
- Quasi-, Pseudo-, etc. PDFs
- transitions to scattering states/resonances * Bigger excited state problems at small quark masses
- go to larger times? How?
- improve smearing? Examples: momentum smearing, perambulators.
- several interpolators, including multi-quark states? * More mixing with lower dimensional operators than in the continuum if chiral symmetry is broken
- grandient flow?

Speaker: Gunnar Bali (Universität Regensburg)
• 17
Contribution to Podium Discussion

The following questions emerged from an e-mail discussion with Gustavo Ramirez:

1) Deflation (with approximate projection) as a multigrid method seems tricky to be ported to GPU architectures in an efficient way.

2) Can one understand why ? Is that solely due to the poor scalability of the 'little Dirac operator' ?

3) Isn't that then a general problem for multigrid methods on GPU ? The same scalability issue should be present for the near-kernel Dirac operator.

4) Are there more ideal solvers for pure GPU architectures (not hybrid CPU-GPU) ?

Speaker: Prof. Michele Della Morte
• Wednesday, August 17
• Numerical Solvers
Convener: Mike Peardon
• 18
Computational strategies for nucleon matrix elements

Nucleon matrix elements are a major lattice QCD input to the search of new physics. Their determination can be challenging, in particular at near physical quark masses. On the one hand large Euclidean time separations are necessary to determine ground state matrix elements. On the other hand the signal over noise ratio deteriorates exponentially with the time. I will detail some of the techniques our group use for the efficient calculation of quark line disconnected as well as of connected contributions to various matrix elements.

Speaker: Sara Collins
• 19
New ideas and developments for solvers for distillation?
Speaker: Juan Andres Urrea-Nino
• 10:20 AM
Coffee Break
• Stochastic Trace Estimation
Convener: Mike Peardon
• 20
Nucleon axial formfactors from lattice QCD

A precise knowledge of nucleon axial formfactors is needed for the new generation of terrestrial neutrino experiments. This is particularly challenging due to increased contamination, in this sector, from excitations, in particular from $N\pi$ scattering states. Transitions from a $N$ to a $N\pi$, mediated by an axial current are also interesting themselves, as these can be related to neutrino induced pion production. We explain how we compute the relevant Wick contractions, combining stochastic and sequentual all-to-all propagator techniques. We then employ the generalized eigenvalue approach to extract the energy levels and matrix elements of interest.

Speaker: Gunnar Bali (Universität Regensburg)
• 21
Fermion loops at the physical point

We present results obtained for disconnected fermion loops contributing to nucleon observables using simulations of twisted mass lattice QCD with physical quark masses and three lattice spacings. We focus on the methods employed, including multi-grid, hierarchical probing, low-mode deflation, and the so-called one-end trick.

Speaker: Giannis Koutsou