Improving Polynomial-filtered Hybrid Monte Carlo With Hasenbusch

TL;DR

This paper compares and combines polynomial filtering with Hasenbusch preconditioning in Hybrid Monte Carlo for lattice QCD, demonstrating a speedup and easier tuning on a small lattice with Wilson fermions.

Contribution

It introduces a novel combination of polynomial and Hasenbusch filters, improving efficiency and simplifying parameter tuning in HMC simulations.

Findings

Combined filters outperform standard Hasenbusch filters in speed.

The new method simplifies the tuning process.

Speedup observed on a small lattice with Wilson fermions.

Abstract

The predominant method for generating Lattice QCD configurations is Hybrid Monte Carlo (HMC). In order to speed up this generation, a wide range of preconditioning techniques that modify the lattice action have been devised. This work compares the performance of the well-known Hasenbusch preconditioning technique with the polynomial filtering technique on a small $16^3 \times 32$ lattice with two flavours of Wilson fermions at a pion mass $M_{\pi} \sim 400~\mathrm{MeV}$. We explore a novel method of combining polynomial and Hasenbusch filters, revealing a speedup when compared to the standard two Hasenbusch filters. This comes with the added advantage of simplified tuning.

Equations

$$S=S_G[U]+S_F[U,\psi ,\overline{\psi }],$$

$$H[U,\varphi ,{\varphi }^{†}]=\sum \mathrm{tr}[P^2]+S[U,\varphi ,{\varphi }^{†}]$$

$$\frac{\partial S_F}{\partial U}=\frac{\partial }{\partial U}[{\varphi }^{†}(M^{†}M{)}^{-1}\varphi ]=-{\varphi }^{†}(M^{†}M{)}^{-1}\frac{\partial M^{†}M}{\partial U}(M^{†}M{)}^{-1}\varphi .$$

$$S=S_{UV}+S_{IR},$$

$$S_{MP}={\varphi }_1^{†}{J}^{-1}{\varphi }_1+{\varphi }_2^{†}J{K}^{-1}{\varphi }_2.$$

$$S_{PF}={\varphi }_1^{†}P(K){\varphi }_1+{\varphi }_2^{†}[P(K)K{]}^{-1},$$

$$C=\frac{N_{mat}}{P_{acc}}$$

$$S=S_G+{\varphi }_1^{†}P(K){\varphi }_1+{\varphi }_2^{†}[P(K)K{]}^{-1}{\varphi }_2,$$

$$I_{LPF}[\tau ]={\left(\hat{S}\left[\frac{h}{2}\right]\hat{T}\left[h\right]\hat{S}\left[\frac{h}{2}\right]\right)}^n,$$

$$I[\tau ]={\hat{S}}_1\left[\frac{\tau }{4}\right]{\hat{S}}_2\left[\frac{\tau }{6}\right]\hat{T}\left[\frac{\tau }{3}\right]{\hat{S}}_2\left[\frac{\tau }{3}\right]\hat{T}\left[\frac{\tau }{6}\right]{\hat{S}}_1\left[\frac{\tau }{2}\right]\hat{T}\left[\frac{\tau }{6}\right]{\hat{S}}_2\left[\frac{\tau }{3}\right]\hat{T}\left[\frac{\tau }{3}\right]{\hat{S}}_2\left[\frac{\tau }{6}\right]{\hat{S}}_1\left[\frac{\tau }{4}\right]$$

$$S_{2PF}={\varphi }_1^{†}{P}_1(K){\varphi }_1+{\varphi }_2^{†}Q(K){\varphi }_2+{\varphi }_3^{†}[P_2(K)K{]}^{-1}{\varphi }_3$$

$$S_{2MF}={\varphi }_1^{†}{J}_1^{-1}{\varphi }_1+{\varphi }_2^{†}{J}_1{J}_2^{-1}{\varphi }_2+{\varphi }_3^{†}{J}_2{K}^{-1}{\varphi }_3.$$

$$S_{\mathrm{PF}-\mathrm{MP}}={\varphi }_1^{†}P(J){\varphi }_1+{\varphi }_2^{†}[P(J)J{]}^{-1}{\varphi }_2+{\varphi }_3^{†}J{K}^{-1}{\varphi }_3$$

$$S_{1PF-2MP}={\varphi }_1^{†}P(J_1){\varphi }_1+{\varphi }_2^{†}[P(J_1)J_1{]}^{-1}{\varphi }_2+{\varphi }_3^{†}{J}_1{J}_2^{-1}{\varphi }_3+{\varphi }_4^{†}{J}_2,K^{-1}{\varphi }_4$$

$$S_{2PF-1MP}={\varphi }_1^{†}P(J){\varphi }_1++{\varphi }_2^{†}Q(J){\varphi }_2+{\varphi }_3^{†}[P_2(J)J{]}^{-1}{\varphi }_3+{\varphi }_4^{†}J{K}^{-1}{\varphi }_4,$$

Full text

1 Introduction and Motivation

Lattice QCD is the non-perturbative method of choice when dealing with strong interactions. In order to measure observables on the lattice, we numerically evaluate path integrals via an ensemble average over a large number of configurations, each of which are described by the state of the gauge field $U$ and fermion field $\psi$. In order to generate these configurations, we most commonly use Hybrid Monte Carlo (HMC), which involves repeated inversions of the Dirac matrix $M$ that describes the strong force between fermions at any two lattice sites. The large number of inversions required and the size of the matrix involved mean that HMC is very computational intensive, making it difficult to simulate near physical quark masses.

A large variety of optimization techniques have been applied to HMC. These techniques reduce the number of matrix inversions required or improve the condition number of the fermion matrix such that the overall cost is reduced. In our work [1], we examine mass preconditioning [2] and polynomial filtering [3] on a $16^{3}\times 32$ lattice to compare their performance, before investigating a more novel technique where these two methods are combined. In section 2, we explain the formulation of both methods before the initial comparison in section 3. Then in section 4 we investigate the performance benefits of using polynomial filters on top of mass preconditioners.

2 Theory

The standard method for generating lattice configurations with dynamical fermions is Hybrid Monte Carlo (HMC). The starting point is the lattice action

[TABLE]

consisting of a gluonic part $S_{G}$ that depends purely on the gauge fields $U$ and a fermionic part $S_{F}$ that also depends on the fermion fields $\psi,\overline{\psi}$. We wish to find configurations $(U,\psi,\overline{\psi})$ distributed according to $\exp(-S[U,\psi,\overline{\psi}])$. To do this, we first use Wick’s theorem to convert the fermions field $\psi$ to bosonic pseudofermion fields $\phi$ to make them computationally friendly, modifying the fermion action in the process. We then need to sample $\phi$ from the distribution $\exp[-S_{F}(U,\phi,\phi^{\dagger})]$, which can be done in practise by relating to Gaussian noise vectors $\chi$ distributed according to $\exp[-\chi^{\dagger}\chi]$. To generate correctly distributed gauge fields $U$, we introduce a fictitious conjugate momentum field $P$ distributed according to $\exp[-\sum\mathrm{tr}[P^{2}]]$, then produce configurations which preserve the Hamiltonian

[TABLE]

by using Hamilton’s equations, leading to integration steps {IEEEeqnarray}rrcl &^T[ϵ]: (P,U) →(P, e^iϵP U) \IEEEyesnumber\IEEEyessubnumber

and ^S[ϵ]: (P,U) →(P - ϵF, U), \IEEEyessubnumber

where $F=\frac{\partial S}{\partial U}$ is the force term. We combine a sequence of these steps into a trajectory to take an existing configuration $(P,U)$ and produce a new state $(P^{\prime},U^{\prime})$, such that we get to a new candidate configuration $(U^{\prime},\phi^{\prime})$. This candidate undergoes a Metropolis acceptance step to ensure detailed balance and hence that the process will converge to the target distribution $\exp[-S]$.

The main cost arises from the fermionic force term, which in the two-flavour degenerate case $S_{F}=\phi^{\dagger}(M^{\dagger}M)^{-1}\phi$ takes the form

[TABLE]

The issue here is the inversion of the fermion matrix $K\equiv M^{\dagger}M$. It is a very large, sparse matrix, so it takes iterative solvers (e.g. conjugate gradient) many matrix multiplications to invert. It can also give rise to large force terms, which then necessitates a reduction in the step-size to keep a good Metropolis acceptance rate and thus increases the number of required inversions. Therefore, modern lattice simulations use a variety of HMC improvements to reduce the force terms (thus reducing the number of required inversions) or to make the matrix $K$ easier to invert.

A template for achieving a reduction in computational cost was proposed by [4], which considers splitting the action into two terms

[TABLE]

such that $S_{UV}$ captures the high-energy modes ($\sim$ large forces) of the action whilst $S_{IR}$ captures the low-energy modes ($\sim$ small forces), and $F_{UV}$ is relatively cheap to calculate compared to $F_{IR}$. Such terms can then be placed on different time-scales using a multiple time-scale integrator, which allows the expensive $F_{IR}$ to be evaluated less often and hence improve the cost.

One of the standard HMC improvement techniques is mass preconditioning (MP) [2], where we factorize the fermion action $S_{F}$ into two terms as follows:

[TABLE]

Here, $J$ is a fermion matrix just like $K$, but with a modified mass parameter $\kappa^{\prime}<\kappa$ giving rise to a ‘heavier’ quark mass. The first action term $S_{1}$ captures the high-energy modes but has a cheaper force term than $S_{2}$, so we can use multiple time-scales to reduce the overall cost.

We compare this technique to polynomial filtering (PF) [3], whereby the fermion action is filtered via

[TABLE]

where $P(K)$ is a low-order polynomial approximating $K^{-1}$. In our work, we use Chebyshev polynomials such that the only parameter to tune is the polynomial order $p$. The polynomial filter term $S_{1}$ has a very cheap force term due to a lack of matrix inversions, and it captures the high-energy modes of the system by virtue of $P(K)$ approximating the inverse. Hence, we can put $S_{1}$ and $S_{2}$ on separate time-scales to achieve a cost reduction.

3 Comparison of filtering methods

3.1 Setup

Our initial study compares the computational cost of polynomial filtering against mass preconditioning. We use a $n_{f}=2$, $16^{3}\times 32$ lattice with a Wilson fermion action with hopping parameter $\kappa=0.15825$, and even-odd preconditioning. This lattice has a pion mass of $\sim 400\ \mathrm{MeV}$ and a lattice spacing of $\sim 0.08\ \mathrm{fm}$.

The metric we use for the cost is

[TABLE]

where $N_{mat}$ is the number of fermion matrix $K$ multiplications per trajectory and $P_{acc}$ is the Metropolis acceptance rate. We average this quantity over at least 2000 trajectories for each displayed data point in the graphs that follow.

As noted in equation (4), we split the different action terms onto different integration time-scales. By way of example, the single polynomial filter action (1PF) takes the form

[TABLE]

and we set the step-sizes $h_{0}=h_{G}<h_{1}<h_{2}$ to reflect the relative cost and frequency scales. The corresponding number of steps $n_{i}$ at each scale are given by the relation $\tau=h_{i}n_{i}$ where $\tau$ is the trajectory length; we have fixed $\tau=1$ as is standard. We use a generalized integration scheme (see section 3.2) that gives a very flexible choice of step-sizes, and the ‘balanced forces’ method to tune these step-sizes. See the paper [1] for further details.

3.2 Generalized multi-scale integration

To create multiple time-scales in a HMC integration, one typically uses a nested integration scheme. Finer integration scales are built by substituting integration schemes for the time updates in the original integration scheme. For example, consider the space-time-space leapfrog integrator

[TABLE]

where $n=\tau/h$ and $\hat{T},\hat{S}$ are as defined in (2), (2). To add a nested integration scale to this scheme, we can replace the time updates $\hat{T}[h]$ with $m$ leapfrog steps in the second term $S_{2}$: {IEEEeqnarray}rCl I_2[h] & = ( ^S_2 [h2] ^T[h] ^S_2[h2] )^m,

s.t.  I_nested[τ] = ( ^S_1 [ h12 ] I_2[h_1] ^S_1[ h12 ] )^n.

This ensures that $S_{1}$ is integrated with step-size $h_{1}=\tau/n$ and $S_{2}$ is integrated with finer step-size $h_{2}=h_{1}/m$. However, nested schemes force the step-sizes of each scale to evenly divide those on each coarser scale.

We have devised an generalized scheme for constructing multiple integration time-scales without any relative step-size restriction. The basic idea is to note that when we integrate some Hamiltonian $H=T+S=T+S_{1}+S_{2}+\ldots$, we always have just one kind of time step $\hat{T}[h]$ that is applied in a uniform direction $h>0$. It thus makes sense to parametrize the progress of time-steps via a time parameter $t$ that increases monotonically from [math] to $\tau$. We can then combine integration schemes for Hamiltonians for each action term $H_{i}=T+S_{i}$ into a generalized scheme for the full Hamiltonian $H=T+S$ by integrating with time steps from [math] to $\tau$, inserting the action step updates $\hat{S}_{i}[h]$ at their respective ‘times’ $t$ in the composite integrators. For example, consider a two-step and a three-step leapfrog integrator for the two action terms: {IEEEeqnarray}rCl I_1[τ] & = ^S_1 [ τ4 ] ^T [ τ2 ] ^S_1 [ τ2 ] ^T [ τ2 ] ^S_1 [ τ4 ], \IEEEyesnumber\IEEEyessubnumber

and  I_2[τ] = ^S_2 [ τ6 ] ^T [ τ3 ] ^S_2 [ τ3 ] ^T [ τ3 ] ^S_2 [ τ3 ] ^T [ τ3 ] ^S_2 [ τ6 ]. \IEEEyessubnumber We can combine these two schemes via the generalized scheme by overlaying the space updates based on their position in ‘time’. This produces the integrator

[TABLE]

Further details are given in [1], where we also prove that this method produces an integration scheme that meets the requirements for detailed balance, given that the composite integrators also do so.

3.3 Results

Figure 1 shows the cost function $C$ (7) for a single polynomial filter and for a single mass preconditioner. As one can see, the mass preconditioned action performs much better: a cost of $C=43,800\pm 3,500$ at $\kappa^{\prime}=0.1545$ compared with $C=87,500\pm 7,400$ at $p=10$. This suggests that a short order polynomial of order $10$ can’t capture as much of the dynamics as a mass preconditioner that requires $80$ or more iterations to invert.

We can improve the performance of both methods by using two filters instead of one. In the case of polynomial filtering, we take two polynomials $P_{1}(K)$ and $P_{2}(K)$ with orders $p_{2}>p_{1}$ that factorize into a polynomial $Q(K)=P_{2}(K)/P_{1}(K)$, such that we can use action

[TABLE]

and thus keep the cost of evaluating the intermediate force $F_{2}$ low. In the case of mass preconditioning, we use two modified mass parameters $\kappa_{1}<\kappa_{2}<\kappa$, and action

[TABLE]

To keep the parameter space manageable, we fixed the cheapest term $S_{1}$ in both cases for our analysis – namely, $p_{1}=4$ and $\kappa_{1}=0.145$.

Figure 2 shows the cost of these two methods side by side. In both cases we see improved results over the one-filter case. Two polynomial filters reduce the cost approximately as much as a single mass preconditioner: a cost of $C=47,700\pm 3,700$ at $p_{2}=54$. Two mass preconditioners reduce the cost further, to $C=31,100\pm 2,200$ at $\kappa_{2}=0.1555$.

4 Combined filters

4.1 Combining polynomial filtering and mass-preconditioning

It is informative to compare the two techniques’ relative efficiency. Polynomial filtering works better at capturing the high-frequency (energy) modes of the system, because we can then use a low polynomial order. To capture lower frequency modes, we have to increase the polynomial order and the cost to evaluate the force terms increases dramatically. On the other hand, Hasenbusch filtering is best suited to small changes $\Delta\kappa=\kappa^{\prime}-\kappa$ in the quark mass, because this implies that $JK^{-1}\approx I$ and hence the correction term is easier to invert. However, the filter term $\phi_{1}^{\dagger}J^{-1}\phi_{1}$ would then have a more expensive force term due to the lighter mass.

To combine the benefits of both of these methods, we place a polynomial filter on top of a mass preconditioner, resulting in the action

[TABLE]

where $J=K(\kappa^{\prime}),\kappa^{\prime}>\kappa$ and $P(J)$ is a polynomial. The idea here is that the high frequency modes are captured by the cheap polynomial $P(J)$, such that the relative mass of the Hasenbusch correction term $\Delta\kappa$ can be kept small.

We tested this on our lattice with the polynomial $P(K)$ fixed to order $p=4$ and $\kappa^{\prime}$ varying. The cost function for various $\kappa^{\prime}$ is shown in Figure 3 in the second column, which can be compared to the 2MP action in the first column. As the graph shows, the cost for PF-MP is very similar to that of 2MP.

4.2 3-filter actions

Note that in the cost function for 2MP and PF-MP, the cost increases significantly depending on the choice of the intermediate term parameter $\kappa_{2}/\kappa^{\prime}$, so a degree of tuning is required to achieve a minimum. This is not a cheap procedure, as we need at least 500 trajectories to get a decent handle on the acceptance rate $P_{acc}$ for each set of parameters and hence the cost. This is complicated by the fact that the Hasenbusch parameter $\kappa_{2}/\kappa^{\prime}$ is a real parameter whose optimal value is heavily dependent on the quark mass $\kappa$, so we’d have to tune it again for different quark masses. This makes it difficult to optimize an action with three mass preconditioners effectively, because there are too many parameters $\kappa_{1},\kappa_{2},\kappa_{3}$ which require fine tuning.

On the other hand, polynomial filtering only depends on a single integral parameter $p$ once a class of polynomials is chosen, and the same polynomial order should filter out similar proportions of the action regardless of $\kappa$. Hence, a polynomial filtering term requires very little to no tuning once a good set of polynomials are found. We can thus add more polynomial filters to our action without increasing the time it takes to tune.

We consider one polynomial filter on top of two mass preconditioners (1PF-2MP)

[TABLE]

where we fix $p=4$, $\kappa_{1}=0.145$; and two polynomial filters on top of a single mass preconditioner (2PF-1MP)

[TABLE]

where we fix $p=4$, $q=20$. The cost function for these two actions are shown in columns three and four respectively in Figure 3.

As can be seen in Figure 3, the optimal cost for the two three-filter actions are very similar to those for PF-MP and 2MP. The extra polynomial filter in 1PF-2MP on top of 2MP has no effect on the performance. However, for 2PF-1MP the cost has a much lower dependence on the choice of the Hasenbusch parameter $\kappa$, as can be seen at the $\kappa^{\prime}=0.1565$ data point. Since the polynomial filter orders $(p,q)$ don’t require much tuning, this indicates that the 2PF-1MP action doesn’t require much tuning to reach close-to-optimum computational costs. In fact, it requires even less tuning than 2MP.

5 Conclusion and Outlook

We have compared the performance of polynomial filtering and mass preconditioning on a modest $16^{3}\times 32$ lattice before combining the two techniques to try to achieve lower costs. Whilst we didn’t observe any significant improvements in cost over the ‘standard’ 2MP algorithm, we did find that 2PF-1MP cost much less to tune due to the cost’s low dependence on $\kappa^{\prime}$ and the fact that the polynomial filters don’t require very much tuning. It remains to be seen if this behaviour still holds for larger lattices and for smaller quark masses.