Rational hybrid Monte Carlo with block solvers and multiple pseudofermions

The dominant cost of most lattice QCD simulations is the inversion of the Dirac operator required to calculate the force term in the rational hybrid Monte Carlo (RHMC) update. One way to improve this situation is to use multiple pseudofermions, which reduces the size and variance of this force and hence allows a larger integration step size to be used. This means fewer force term calculations are required, but at the cost of having to invert the Dirac operator for each pseudofermion field. This bottleneck can be addressed: recently there has been renewed interest in the use of block Krylov solvers, which can solve multiple right-hand-side vectors with significantly fewer iterations than are required if each vector is solved using a separate Krylov solver. We combine these two ideas, achieving a significant speed-up of RHMC lattice QCD simulations.

Figure 1

Top: Number of Dirac operator calls for the block SBCGrQ solver to converge compared to the SCG solver, for a range of fermion masses. As the mass is made lighter, the block solver improvement increases. Bottom: The square root of the ratio of the lowest eigenvalue to the $n_{\mathrm{pf}}\mathrm{th}$ eigenvalue of the Dirac operator. The bound on the convergence rate of Eq. (19) is determined by this quantity, and qualitatively it also seems to describe the actual convergence of the block solver quite well.

Figure 2

Top: Residual of shifted solution versus shift $\sigma$, for fixed $k=400$ solver iterations. We see a dramatic decrease in the residual for small shifts as $n_{\mathrm{pf}}$ is increased. Bottom: Same quantity but keeping fixed the unshifted residual $\left|r\right|/\left|r_0\right|={10}^{-7}$. We see that the reduction in iterations leads to a relative increase in the residuals of the larger shifts with $n_{\mathrm{pf}}$.

Figure 3

Measured expectation values of acceptance rate (solid lines), compared with the predicted acceptance rate from measured force variances using Eq. (13) (dotted lines). For high acceptance and small $\delta \tau$ the agreement is reasonable; it turns out the difference between the prediction and the measured values is largely due to the neglected correlation between initial and final force terms not being negligible in these data, so increasing the trajectory length would improve the agreement.

Figure 4

Gauge and fermion force norms versus $n_{\mathrm{pf}}$, with large-$n_{\mathrm{pf}}$ scaling predictions. Top: Force norms with a fit to Eq. (16). Bottom: Fourth root of variance of force norms, approximately proportional to the number of integration steps required for the OMF2 integrator, along with a fit to Eq. (17). Force (I) is measured at every integration step along the trajectory, while Force (II) is measured on the same set of 2000 thermalized configurations.

Figure 5

Top: Histogram of the rms pseudofermion force norm, $\sqrt{F^2}$, for different $n_{\mathrm{pf}}$. For $n_{\mathrm{pf}}=1$ the distribution is very non-Gaussian, with a long tail of large values. Bottom: Histogram of $\sqrt{n_{\mathrm{pf}}{F}^2}$ wich shows an approximate $n_{\mathrm{pf}}$ invariance for intermediate values of $n_{\mathrm{pf}}$, due to the $c_1$ and $c_3$ terms dominating Eqs. (16) and (17) for these values of $n_{\mathrm{pf}}$.

Figure 6

Histogram of $e^{-\mathrm{\Delta }H}$ for $n_{\mathrm{pf}}=1$ to 6 with the acceptance rate tuned to $\simeq 90\%$. The black dotted line shows the prediction for a Gaussian distribution of $\mathrm{\Delta }H$ with the same acceptance rate. For $n_{\mathrm{pf}}=1$ (top left) the distribution is very far from gaussian, with an excess of very small values, but as $n_{\mathrm{pf}}$ is increased the distribution approaches the Gaussian one.

Figure 7

Trajectory speed-up versus $n_{\mathrm{pf}}$, normalized to 1 for $n_{\mathrm{pf}}=1$. At each step the multishift CG (SCG) solver is used $n_{\mathrm{pf}}$ times. The black dashed line is the simple prediction from the variance of the pseudofermion force norm using Eq. ($14$), and the green dotted line is the simple prediction from the condition number of the Dirac operator using Eq. (15). Going from $n_{\mathrm{pf}}=1$ to $n_{\mathrm{pf}}=2$ gives a significant cost reduction, but increasing $n_{\mathrm{pf}}$ further results in a larger cost per trajectory.

Figure 8

Top: Runtime of Dirac operator acting on vectors in block form, normalized to the nonblock form. Bottom: Solver runtime per Dirac operator call versus $n_{\mathrm{pf}}$. For small $n_{\mathrm{pf}}$ one iteration of block multishift SBCGrQ is much faster than multishift SCG since the block Dirac operator is faster. For large enough $n_{\mathrm{pf}}$, however, the SBCGrQ solver overhead that grows $\propto n_{\mathrm{pf}}{N}_{\mathrm{shifts}}$ eventually dominates the cost.

Figure 9

Cost of calculating the pseudofermion force term versus $n_{\mathrm{pf}}$, using either multishift (SCG) or block multishift (SBCGrQ) solvers with stopping criterion ${10}^{-6}$ or ${10}^{-12}$. The block solver is a significant improvement, moreover it allows the use of a very tight stopping criterion without significant extra cost, which reduces possible reversibility violations.

Figure 10

Cost of calculating the pseudofermion action versus $n_{\mathrm{pf}}$, using either multishift (SCG) or block multishift (SBCGrQ) solvers with stopping criterion ${10}^{-14}$. This is done twice per trajectory: at the start for the heat bath and at the end for the accept/reject step. The block solver significantly reduces the cost of this step.

Figure 11

The cost of generating an accepted $\tau =1$ trajectory using the block multishift (SBCGrQ) inverter for different $n_{\mathrm{pf}}$ versus the number of integrator steps $n_{\mathrm{steps}}=\tau /\delta \tau$. Top: cost in Dirac operator calls; bottom: computer runtime cost.

Figure 12

Trajectory speed-up versus $n_{\mathrm{pf}}$ at $\simeq 90\%$ acceptance, normalised to 1 for $n_{\mathrm{pf}}=1$, using either multishift (SCG) or block multishift (SBCGrQ) solvers. Top: speed-up in Dirac operator calls; bottom: speed-up in computer runtime. The black dashed line is the simple prediction using the force norm variance of Eq. (14), and the green dotted line is the simple prediction from the condition number of the Dirac operator using Eq. (15). The optimal $n_{\mathrm{pf}}$ and the overall gain are both significantly increased by the use of block methods.