Exact Local Bosonic Algorithm for Dynamical Quarks

TL;DR

This paper introduces an exact local bosonic algorithm for simulating dynamical fermions in lattice QCD, offering a potential alternative to the widely used Hybrid Monte Carlo method.

Contribution

The paper presents a novel exact local bosonic algorithm that is a feasible alternative to existing methods for lattice QCD simulations.

Findings

Algorithm is exact and local.

Viable alternative to Hybrid Monte Carlo.

Potential improvements in simulation efficiency.

Abstract

We present an exact local bosonic algorithm for the simulation of dynamical fermions in lattice QCD. We show that this algorithm is a viable alternative to the Hybrid Monte Carlo algorithm.

Equations

$$Z=\int [dU]|\det D{|}^2{e}^{-S_G[U]}$$

$$|\det D{|}^2=\det D^{†}\det D\simeq \frac{1}{|\det P(D){|}^2}$$

$$Z\simeq \int [dU][d\varphi ][d{\varphi }^{†}]e^{-S_L[U,\varphi ]}$$

$$Z=\int [dU][d\varphi ][d{\varphi }^{†}]|\det (DP(D)){|}^2{e}^{-S_L(U,\varphi )}$$

$$Z=\int [dU][d\eta ][d{\eta }^{†}][d\varphi ][d{\varphi }^{†}]e^{-S_{exact}(U,\varphi ,\eta )}$$

$$S_{exact}(U,\varphi ,\eta )=S_L(U,\varphi )+S_C(U,\eta )$$

$$P_{acc}((U,\varphi )\to (U^{\prime },{\varphi }^{\prime }))=$$

$$\displaystyle\left\{\begin{array}[]{ll}\min\left(1,e^{-|W\chi|^{2}+|\chi|^{2}}\right)&\mbox{{\footnotesize if $S_{G}(U)\geq S_{G}(U^{\prime})$}}\\ \min\left(1,e^{+|W^{-1}\chi|^{2}-|\chi|^{2}}\right)&\mbox{{\footnotesize if $S_{G}(U)<S_{G}(U^{\prime})$}}\end{array}\right.$$

$$<P_{acc}>=erfc\left(f\sqrt{96V}{\left(\frac{k}{k_c}\right)}^{n+1}\right)$$

Full text

1 Introduction

The local bosonic algorithm was proposed by Lüscher [1] in a hermitian version as a new method to simulate full QCD. Here we present several significant improvements of this algorithm: an inexpensive stochastic Metropolis accept/reject test [2, 3] to make the algorithm exact; a non-hermitian polynomial approximation [2, 3] and a simple even-odd preconditioning [4, 3]. We illustrate the effectiveness of our algorithm for two-flavor QCD and compare its performance with the Hybrid Monte Carlo (HMC) algorithm in the relevant regime of small quark mass. We successfully model its static properties (the Monte Carlo acceptance), and try to disentangle its dynamics. From our analysis it appears that this version of the local bosonic algorithm is indeed a viable alternative to the HMC.

2 Description of the algorithm

The full QCD partition function with two fermion flavors is given by

[TABLE]

where $D$ represents the fermion matrix and $S_{G}$ the gauge action. This partition function can be approximated by introducing a polynomial approximation of the inverse of the fermion matrix.

[TABLE]

The polynomial $P(z)=\prod_{k=1}^{n}(z-z_{k})$ of degree $n$ is defined in the complex plane and approximates the inverse of $z$. Since we are investigating full QCD with two flavors, the determinant $|\det P(D)|^{2}$ manifestly factorizes into positive pairs, so that the $\frac{1}{|\det P(D)|^{2}}$ term of the approximation (2) can be expressed as a Gaussian integral over a set of boson fields $\phi_{k}$ ($k=1,...,n$) with color and Dirac indices. The full QCD partition function (1) is then approximated by

[TABLE]

where $\phi$ represents the set of all boson field families, $S_{L}$ is the local action $S_{L}=S_{G}+S_{b}$ and $S_{b}=\sum_{k=1}^{n}|(D-z_{k})\phi_{k}|^{2}$. Making use of the locality of $S_{L}$ we may now simulate the partition function (3) by locally updating the boson fields and the gauge fields, using heat-bath and over-relaxation algorithms.

The simulation of full QCD can be obtained from (3) by correcting the errors due to the approximation through a Metropolis test at the end of each trajectory. Introducing the error term in the partition function we obtain

[TABLE]

The correction term $|\det(DP(D))|^{2}$ can be evaluated by expressing the determinant $|\det(DP(D))|^{2}$ as a Gaussian integral over an auxiliary field $\eta$ and incorporating it in the partition sum (4)

[TABLE]

by defining a new exact action

[TABLE]

where $S_{C}(U,\eta)=|[DP(D)]^{-1}\eta|^{2}$ is the correction action. In carrying out the simulation one generates configurations of $(U,\phi)$ and $\eta$ such that the probability of finding a particular configuration is proportional to $\exp(-S_{exact})$. The strategy is to alternatively update the $(U,\phi)$ fields and the $\eta$ field. The Metropolis acceptance probability $P_{acc}$ is defined using the exact action $S_{exact}$ so that the transition probability of the algorithm satisfies detailed balance with respect to $S_{exact}$ [3]. To be more precise:

[TABLE]

where $\chi$ is a Gaussian distributed spinor (used to generate $\eta$ with a global heat-bath by $\eta=DP(D)\chi$ or $\eta=D^{\prime}P(D^{\prime})\chi$) and $W=[D^{\prime}P(D^{\prime})]^{-1}DP(D)$. Here $D$ and $D^{\prime}$ denote the fermion matrix with the $U$ and $U^{\prime}$ gauge configurations,

respectively. The algorithm can be summarized as follows:

$\bullet$ Generate a new Gaussian spinor $\chi$ with variance one;

$\bullet$ Update locally the boson and gauge fields $m$ times (in a reversible and optimized order [4]) according to the approximate partition function (3);

$\bullet$ Accept/reject the new configuration according to the Metropolis acceptance probability (7).

In order to evaluate the Metropolis acceptance probability we have to solve a linear system involving $DP(D)$ or $D^{\prime}P(D^{\prime})$, for which we use the BiCGstab algorithm [5]. This linear system is very well conditioned because $P(D^{\prime})$ (or $P(D)$) approximates the inverse of $D^{\prime}$ (or $D$). The cost for solving it is minimal and scales like the local updating algorithms in the volume and quark mass.

3 Numerical results

Numerical simulations using the exact local bosonic algorithm have been performed for different lattice parameters. The majority of the simulations are reported in [3]. We explored the acceptance of the Metropolis correction test and the number of iterations of the BiCGStab algorithm used in that test for inverting $DP(D)$. The study was performed by varying the degree $n$ of the polynomial and the hopping parameter $k$. One observes that the acceptance increases quite rapidly with the degree of the polynomial, above some threshold [3] (see Fig.1). On the other hand, the number of iterations needed to invert $DP(D)$ remains very low for high enough acceptance. The data show clearly that the overhead due to the Metropolis test remains negligible provided that the degree of the polynomial is tuned to have sufficient acceptance. Results obtained from simulations using even-odd preconditioning confirm, as expected, that the improvement of the approximation reduces the required number of bosonic fields by at least a factor two [3, 4]. This improvement is trivially implemented by a simple redefinition of the roots $z_{k}$ [3].

We have also performed a simulation for parameters relevant for physical interesting conditions. On a $8^{4}$ lattice at $\beta=5.6$ and $k=0.1585$ (very near $k_{c}$) with Schrödinger functional boundary conditions we have compared the performance of our a

lgorithm with HMC. In Table 1 we present the plaquette and its integrated autocorrelation time in units of fermion-matrix multiplications ($[D\phi]$) The simulation with the local bosonic algorithm was performed using 48 boson fields. Its average acceptance is $<P_{acc}>=0.589(1)$. The HMC data is obtained from [6].These results show that for light quarks the peformance of our algorithm is comparable to the performance of HMC.

$<\Box>$ $\tau_{int}(\Box)$

HMC 0.5726(1) 33750(6000) $[D\phi]$

multiboson 0.57275(15) 23040(4000) $[D\phi]$

Table 1

4 Optimisation

Using some general assumptions and the known error bounds on the Chebyshev-like approximation of $D^{-1}$, we can obtain [3] an ansatz for the acceptance probability, which takes the form

[TABLE]

where $V$ is the lattice volume, $k$ is the hopping parameter of the Wilson fermion matrix $D$ and $k_{c}$ is the critical hopping parameter at a given $\beta$. We expect the fitting parameter $f$ to depend smoothly on $\beta$, but very little on $n,V$ and $k$. In [3] we show that our ansatz (11) is very satisfactory (see Fig.1). At other values of $\beta$, one can fix $f$ by a test on a small lattice, and then predict the acceptance for larger volumes and different quark masses. We can use our ansatz to analyze and optimize the cost of a simulation. From the analysis it turns out that the optimal average acceptance is around 0.7-0.8.

5 Scaling

The total complexity with the volume and the quark mass is related to the scaling of the number of bosonic fields and the autocorrelation of the observables. As the volume $V$ increases, the number $n$ of bosonic fields should grow like $\log V$. As the quark mass $m_{q}$ decreases, the number of fields $n$ must grow like $m_{q}^{-1}$. This is a consequence of the exponential convergence of the polynomial approximation. We expect the autocorrelation to behave like $\tau\sim nm_{q}^{-\alpha}$ [4, 3, 7] so that the complexity is proportional to $n^{2}m_{q}^{-\alpha}$. The exponent $\alpha$ can only be derived from MC data since there is no clear theoretical understanding of the coupled dynamics of the gauge and boson fields. From an exploratory simulation [3] on a $4^{4}$ lattice at $\beta=0$ we obtained that $\alpha$ is near $1$. This is only indicative, because the dependence of $\alpha$ on the volume and the coupling has yet to be explored.

6 Conclusion

We have presented an alternative algorithm to HMC for simulating dynamical quarks in lattice QCD. This algorithm is based on a local bosonic action. A non-hermitian polynomial approximation of the inverse of the quark matrix is used to define the local bosonic action. The addition of a global Metropolis test corrects the systematic errors. The overhead of the correction test is minimal. This algorithm is exact for any choice of the polynomial approximation. No critical tuning of the approximation parameters is needed. Only the efficiency of the algorithm, which can be monitored, will be affected by the choice of parameters. The cost of the algorithm scales similarly with that of HMC with the volume $V$ of the lattice and with the inverse of the quark mass. Finally, for light quarks the peformance of our algorithm is comparable to the performance of HMC provided that there is enough memory to store the needed boson fields.