# UV-filtered fermionic Monte Carlo

## Abstract

The short-range modes of the fermionic determinant can be absorbed in the gauge action using the loop expansion. The coefficients of this expansion and the zeroes of the polynomial approximating the remainder can be optimized by a simple, practical method. When the multiboson approach is used, this optimization results in a faster simulation with fewer auxiliary fields.

Dynamical fermion simulations are orders of magnitude slower than quenched ones. With Hybrid Monte Carlo (HMC), this remains true also for heavy dynamical quarks: even in the quenched regime, HMC is $`𝒪(100)`$ less efficient than a local Monte Carlo update scheme . This poor behaviour may explain why HMC computing requirements grow rather slowly as the quark mass is decreased. It clearly indicates room for improvement, which hopefully will persist down to light quarks.

The fermion determinant can be expanded in loops:

$$\mathrm{det}(\mathrm{𝟏}\kappa M)=e^{\mathrm{Tr}\mathrm{Log}(\mathrm{𝟏}\kappa M)}=e^{_l\frac{\kappa ^l}{l}\mathrm{Tr}M^l}$$
(1)   
using Wilson fermion notation. The first nonzero term ($`l=4`$) gives a shift $`\mathrm{\Delta }\beta `$ of the gauge coupling. It is remarkable that the bulk of the fermionic effects can be reabsorbed into this simple shift, which can be computed perturbatively down to rather small quark masses . Unfortunately, neither HMC nor the alternative MultiBoson (MB) method makes use of this observation. Effective loop actions () $`S_{\mathrm{eff}}`$ cannot be used to guide an exact algorithm, because $`(S_{\mathrm{exact}}S_{\mathrm{eff}})^2Volume`$, which causes an exponentially small $`e^{Volume}`$ Metropolis acceptance. Instead, I propose to make use of the identity $`e^{\mathrm{Tr}A}\times \mathrm{det}e^A=1`$ to rewrite the fermion determinant as

$`\mathrm{det}(\mathrm{𝟏}\kappa M)`$   
$`e^{_ja_j\mathrm{Tr}M^j}\times \mathrm{det}\left((\mathrm{𝟏}\kappa M)e^{+_ja_jM^j}\right)`$   
The number of nonzero coefficients $`a_j`$ and their values are subject to optimization. Note that the first four terms cause no overhead. The new operator whose determinant is to be taken will hopefully be easier to sample than the original one. In particular, for very heavy quarks the coefficients $`a_j`$ should tend to the Taylor values $`\kappa ^j/j`$, leaving only tiny higher-order terms in the determinant.

Both HMC and MB sample the determinant $`\mathrm{det}A`$ by building a polynomial $`P_n(A)`$ which approximates $`A^1`$. Gains in efficiency come from a better approximation (smaller degree or smaller error). To illustrate the advantage of using Eq.(2), in Fig.1 I show the relative error of the polynomial approximation for a real variable $`x[0.01,1]`$. The dotted line corresponds to the Chebyshev approximation used in the MB method; the solid line is the error achievable by including two terms $`a_1`$ and $`a_2`$ in the “preconditioner” $`e^{_ja_jx^j}`$, with a polynomial of the same degree $`n=20`$. Fig.2 shows the zeroes of the two polynomials in the complex plane. While the zeroes of the Chebyshev polynomial are distributed uniformly over an ellipse covering the entire approximation interval, the other zeroes are much more concentrated near the origin. This is the essence of the benefits of the “filter” $`e^{_ja_jx^j}`$: large eigenvalues (corresponding to UV modes of the Dirac operator) are effectively removed from the determinant, which frees the zeroes of the polynomial to concentrate near small eigenvalues corresponding to IR modes and genuine long-range fermionic interactions.

Choosing the coefficients $`\{a_j\}`$ and the zeroes $`\{z_k\}`$ of the polynomial approximation seems like a stiff nonlinear optimization problem, with the $`a_j`$ appearing in the exponent. Even setting $`a_j=0j`$, the choice of the zeroes $`\{z_k\}`$ has been the object of several studies . The problem is how to fix the parameters of $`W_k^n(\mathrm{𝟏}\kappa Mz_k\mathrm{𝟏})(\mathrm{𝟏}\kappa M)e^{_{j=0}^{m1}a_jM^j}`$ such that $`\mathrm{det}W1`$. Rewriting the determinant as an average over Gaussian random vectors $`\eta `$

$$\mathrm{det}^2W=\frac{𝑑\eta ^{}𝑑\eta e^{\eta ^{}W^{}W\eta }}{𝑑\eta ^{}𝑑\eta e^{\eta ^{}\eta }}=e^{|W\eta |^2+|\eta |^2}_\eta $$
(3)   
we see that a sufficient condition for $`\mathrm{det}W1`$ is $`W\eta \eta ^20`$. This optimization problem is easy to solve, using the following steps:   
1. Draw one (or more) Gaussian vectors $`\eta `$;   
2. Assign values to the $`\{a_j\}`$; compute $`\psi (\mathrm{𝟏}\kappa M)e^{_{j=0}^{m1}a_jM^j}\eta `$;   
3. Construct the polynomial $`Q_n(M)_k^n(\mathrm{𝟏}\kappa Mz_k\mathrm{𝟏})`$ which minimizes $`cQ_n(M)\psi \eta ^2`$;   
4. Compute $`\stackrel{}{}_a\left\{\frac{c}{a_j}\right\}`$ as a by-product; if $`\stackrel{}{}_a>ϵ`$, return to (2).   
Step 3 is a straightforward quadratic minimization similar to GMRES, while the loop 2–4 uses Newton’s method to solve the nonlinear minimization in the $`a_j`$. In principle, some subtle averaging over gauge fields should also be performed; in practice, different equilibrium gauge fields yield similar results.

This optimization method is simple and practical. It requires no a priori knowledge of the Dirac spectrum. It is applicable to all cases where a determinant is involved: Wilson or staggered fermions, SUSY, condensed matter. When applied to the standard MB method ($`\{a_j\}=0`$, Hermitian or non-Hermitian Wilson Dirac operator), it reveals that the usual Chebyshev approximation is not optimal.

In an actual MC simulation $`\{U_{\mathrm{old}}\}\{U_{\mathrm{new}}\}`$, the Metropolis acceptance probability is $`min(1,\frac{e^{|W_{\mathrm{old}}\eta |^2+|\eta |^2}}{e^{|W_{\mathrm{new}}\eta |^2+|\eta |^2}}_\eta )`$. This can be estimated as $`\mathrm{erfc}(cW_{\mathrm{old}}\eta \eta )`$ with $`c𝒪(1)`$ as a by-product of the above optimization. In this way, the benefits of adding more terms to the preconditioner (especially 6-link and larger loops) can be assessed without having to program the actual Monte Carlo update.

Tests of the dynamical behaviour of this UV-filtered MB algorithm have been performed for moderate quark masses: $`\beta =5.3,\kappa =0.158`$, with two flavors of Wilson fermions, on an $`8^4`$ lattice. Table I compares the efficiency of HMC, the non-Hermitian MB , and the present method. All the programs used even-odd preconditioning. The HMC program incorporates multiple step-size integration  and low-accuracy solution during the trajectory , and uses the BiCG$`\gamma _5`$ solver . The number $`n`$ of fields used by the non-Hermitian MB is consistent with the number of iterations of the HMC solver. In contrast, the UV-filtered version uses $`3`$ times fewer fields. This solves the memory bottleneck of the MB approach and dramatically reduces the work per independent configuration, which grows like $`n^2`$ . Fig.3 shows the autocorrelation function of the plaquette for HMC and the UV-filtered MB. For heavier quarks the superiority of the latter would be even greater. For lighter quarks the advantage is less pronounced, and details of the implementation and the tuning (eg., over-relaxation of the auxiliary fields) become more relevant.