Domain Wall Fermion QCD with the Exact One Flavor Algorithm

TL;DR

This paper introduces an efficient implementation of the exact one flavor algorithm (EOFA) for domain wall fermion lattice QCD, significantly reducing computational costs in simulations involving single quark flavors.

Contribution

The authors developed a new preconditioner for EOFA and optimized simulation parameters, decreasing the time per accepted trajectory by over four times.

Findings

Reduced simulation time by a factor of 4.2.

Enabled more cost-effective single-flavor QCD calculations.

Improved algorithmic efficiency for G-parity boundary condition ensembles.

Abstract

Lattice QCD calculations including the effects of one or more non-degenerate sea quark flavors are conventionally performed using the Rational Hybrid Monte Carlo (RHMC) algorithm, which computes the square root of the determinant of $\mathscr{D}^{\dagger} \mathscr{D}$, where $\mathscr{D}$ is the Dirac operator. The special case of two degenerate quark flavors with the same mass is described directly by the determinant of $\mathscr{D}^{\dagger} \mathscr{D}$ --- in particular, no square root is necessary --- enabling a variety of algorithmic developments, which have driven down the cost of simulating the light (up and down) quarks in the isospin-symmetric limit of equal masses. As a result, the relative cost of single quark flavors --- such as the strange or charm --- computed with RHMC has become more expensive. This problem is even more severe in the context of our measurements of the $\Delta I = 1/2$ $K \rightarrow \pi \pi$ matrix elements on lattice ensembles with $G$-parity boundary conditions, since $G$-parity is associated with a doubling of the number of quark flavors described by $\mathscr{D}$, and thus RHMC is needed for the isospin-symmetric light quarks as well. In this paper we report on our implementation of the exact one flavor algorithm (EOFA) introduced by the TWQCD collaboration for simulations including single flavors of domain wall quarks. We have developed a new preconditioner for the EOFA Dirac equation, which both reduces the cost of solving the Dirac equation and allows us to re-use the bulk of our existing high-performance code. Coupling these improvements with careful tuning of our integrator, the time per accepted trajectory in the production of our 2+1 flavor $G$-parity ensembles with physical pion and kaon masses has been decreased by a factor of 4.2.

Figures

Tables

Table 1. Table 1: Timings for one HMC trajectory of RBC/UKQCD’s 80 2 × 96 × 192 × 32 superscript 80 2 96 192 32 80^{2}\times 96\times 192\times 32 N f = 2 + 1 + 1 subscript 𝑁 𝑓 2 1 1 N_{f}=2+1+1 ensemble with physical quark masses and a − 1 ≈ 3 superscript 𝑎 1 3 a^{-1}\approx 3 GeV on a 12,288-node Blue Gene/Q partition McGlynn ( 2016 ) .

Action Component	Timings
Gauge	5970 s	12.0%
Light Quarks	19600 s	39.4%
Strange and Charm Quarks	24200 s	48.6%
Total	49770 s	—

Table 2. Table 2: Summary of ensembles and input parameters used in this work. Here β 𝛽 \beta is the gauge coupling, L 3 × T × L s superscript 𝐿 3 𝑇 subscript 𝐿 𝑠 L^{3}\times T\times L_{s} is the lattice volume decomposed into the length of the spatial ( L 𝐿 L ), temporal ( T 𝑇 T ), and fifth ( L s (L_{s} ) dimensions, and a ​ m l 𝑎 subscript 𝑚 𝑙 am_{l} and a ​ m h 𝑎 subscript 𝑚 ℎ am_{h} are the bare, input light and heavy quark masses. On the 16I-G, 16ID-G, and 32ID-G ensembles G 𝐺 G -parity boundary conditions are applied to the fermion fields at one or more of the spatial boundaries of the lattice; otherwise periodic boundary conditions are applied, and in all cases antiperiodic boundary conditions are used along the temporal direction.

Ensemble	Action	$\beta$	$L^3\times T\times L_s$	Möbius Scale	$G$-Parity B.C.	$a{m}_l$	$a{m}_h$
16I	DWF + I	2.13	${16}^3\times 32\times 16$	—	—	0.01	0.032
16I-G	DWF + I	2.13	${16}^3\times 32\times 16$	—	$x$	0.01	0.032
16ID-G	MDWF + ID	1.75	${16}^3\times 32\times 8$	2.00	$x$,$y$,$z$	0.01	0.045
24ID	MDWF + ID	1.633	${24}^3\times 64\times 24$	4.00	—	0.00789	0.085
32ID-G	MDWF + ID	1.75	${32}^3\times 64\times 12$	2.67	$x$,$y$,$z$	0.0001	0.045

Table 3. Table 3: Summary of spatial volumes, lattice cutoffs, and pion masses in physical units for the ensembles used in this work. All values for the 16I and 32ID-G ensembles are from Refs.  Allton  et al. ( 2007 ) and Bai  et al. ( 2015a ) , respectively. On the 16I-G (16ID-G) ensemble we assume the lattice cutoff is the same as the 16I (32ID-G) ensemble since the same action and value of β 𝛽 \beta has been used. The pion masses on the 16I-G and 16ID-G ensembles have been extracted using the fitted value of the lowest energy pion states from Table 8 and the continuum dispersion relation. Finally, the determination of the lattice scale for the 24ID ensemble was performed in Ref.  Boyle  et al. ( 2016 ) , and the determination of the pion mass in Ref.  Mawhinney ( 2016 ) .

Ensemble	$L$ (fm)	$a^{-1}$ (GeV)	$m_{\pi }$ (MeV)
16I	1.95(5)	1.62(4)	400(11)
16I-G	1.95(5)	1.62(4)	388(14)
16ID-G	2.29(1)	1.378(7)	575(11)
24ID	4.82(19)	0.981(39)	137.1(5.5)
32ID-G	4.57(2)	1.378(7)	143.1(2.0)

Table 4. Table 4: Basic integrator and HMC details for the generation of the 16I, 16I-G, and 16ID-G ensembles. We use nested Sexton-Weingarten integration schemes, detailed in Table 5 , with δ ​ τ 𝛿 𝜏 \delta\tau the coarsest time step used to evolve the outermost level. We denote the CG stopping tolerances used for the force gradient forecasting, molecular dynamics, and Monte Carlo steps by r FG subscript 𝑟 FG r_{\rm FG} , r MD subscript 𝑟 MD r_{\rm MD} , and r MC subscript 𝑟 MC r_{\rm MC} , respectively.

Ensemble	Integrator	$\delta \tau$	$r_{\mathrm{FG}}$	$r_{\mathrm{MD}}$	$r_{\mathrm{MC}}$
16I	Force Gradient QPQPQ	0.1000	${10}^{-7}$	${10}^{-8}$	${10}^{-10}$
16I-G	Omelyan ($\lambda =0.2$)	0.2000	—	${10}^{-8}$	${10}^{-10}$
16ID-G	Force Gradient QPQPQ	0.1667	${10}^{-6}$	${10}^{-7}$	${10}^{-10}$

Table 5. Table 5: Integrator layouts for the original RHMC runs. Here “Quo” is an abbreviation for the quotient action (Eqn. ( 38 )) and “ RatQuo 1 / n subscript RatQuo 1 𝑛 \mathrm{RatQuo}_{1/n} ” is an abbreviation for the rational quotient action (Eqn. ( 39 )), with a rational function approximation used to apply ( 𝒟 † ​ 𝒟 ) 1 / n superscript superscript 𝒟 † 𝒟 1 𝑛 (\mathscr{D}^{\dagger}\mathscr{D})^{1/n} and its inverse. For the EOFA reproduction runs each instance of RatQuo 1 / 2 subscript RatQuo 1 2 \mathrm{RatQuo}_{1/2} is replaced by an EOFA determinant with the same masses (Eqn. ( 10 )), while all other ensemble and integrator details are left fixed. The notation A:B for the update scheme denotes the number of steps of the next innermost integrator level (A) per step of the current level (B).

Ensemble	Level	Action	Update
16I	1	Quo(0.01,0.2) + Quo(0.2,1.0) + ${\mathrm{RatQuo}}_{1/2}$(0.032,1.0)	4:1
	2	Gauge	1:1
16I-G	1	${\mathrm{RatQuo}}_{1/2}$(0.01,0.032)	1:1
	2	${\mathrm{RatQuo}}_{1/4}$(0.032,1.0) + ${\mathrm{RatQuo}}_{1/4}$(0.032,1.0) + ${\mathrm{RatQuo}}_{1/4}$(0.032,1.0)	8:1
	3	Gauge	1:1
16ID-G	1	${\mathrm{RatQuo}}_{1/2}$(0.01,0.05) + ${\mathrm{RatQuo}}_{1/2}$(0.05,1.0) + ${\mathrm{RatQuo}}_{1/4}$(0.045,1.0)	1:1
	2	DSDR	8:1
	3	Gauge	1:1

Table 6. Table 6: Average plaquettes, quark condensates, and topological susceptibilities ( χ t subscript 𝜒 𝑡 \chi_{t} ) computed on the 16I, 16I-G and 16ID-G lattices and their corresponding EOFA reproduction ensembles. The ensemble averages on the 16I (16I-G) lattices were computed using MD trajectories 500-1500 (500-2500) after binning over 50 successive MD time units. The ensemble averages on the 16ID-G lattices were computed using MD trajectories 500:900 for the RHMC ensemble, and MD trajectories 960:1360 for the EOFA ensemble, after binning over 25 successive MD time units. We do not compute χ t subscript 𝜒 𝑡 \chi_{t} on the 16ID-G ensemble since the short 400 MD time unit measurement runs are insufficient to adequately sample the topological charge, as evidenced by the time evolutions plotted in Appendix C .

	16I		16I-G		16ID-G
Observable	RHMC	EOFA	RHMC	EOFA	RHMC	EOFA
$⟨P⟩$	0.588087(22)	0.588106(26)	0.588033(24)	0.588039(16)	0.514251(43)	0.514200(48)
$⟨{\overline{\psi }}_l{\psi }_l⟩$	0.001697(5)	0.001698(11)	0.0017151(72)	0.0017130(52)	0.005543(11)	0.005563(8)
$⟨{\overline{\psi }}_s{\psi }_s⟩$	0.0037450(31)	0.0037435(74)	0.0037541(51)	0.0037529(34)	0.0085729(82)	0.0085895(69)
$⟨{\overline{\psi }}_l{\gamma }_5{\psi }_l⟩$	-0.000015(14)	-0.000012(19)	-0.000003(15)	-0.000006(12)	0.000033(13)	-0.000001(11)
$⟨{\overline{\psi }}_s{\gamma }_5{\psi }_s⟩$	-0.000001(8)	-0.000007(12)	-0.0000004(92)	-0.0000034(81)	0.000017(10)	-0.000002(8)
${\chi }_t$	$1.03\left(19\right)\times {10}^{-5}$	$1.81\left(42\right)\times {10}^{-5}$	$2.16\left(47\right)\times {10}^{-5}$	$1.53\left(27\right)\times {10}^{-5}$	—	—

Table 7. Table 7: Low energy spectrum on the 16I ensemble computed from 100 independent measurements beginning with MD trajectory 500 and separated by 10 MD time units. Prior to fitting the correlation functions were binned over groups of 5 measurements. Corresponding effective mass plots can be found in Appendix C .

	16I
Observable	RHMC	EOFA
$a{m}_{\pi }$	0.2424(11)	0.2425(8)
$a{m}_K$	0.3252(11)	0.3253(7)
$a{m}_{\Omega }$	1.003(15)	0.994(11)
$a{m}_{\mathrm{res}}^{\prime }\left(m_l\right)$	0.0030558(80)	0.0030523(78)

Table 8. Table 8: Low energy spectra on the 16I-G and 16ID-G ensembles computed from 51 and 21 measurements, respectively. On the 16I-G ensemble we also predict the ground state pion energy using the fitted a ​ m π 𝑎 subscript 𝑚 𝜋 am_{\pi} on the 16I ensemble and the continuum dispersion relation a ​ E π pred = ( a ​ m π ) 2 + ( a ​ π / L ) 2 𝑎 superscript subscript 𝐸 𝜋 pred superscript 𝑎 subscript 𝑚 𝜋 2 superscript 𝑎 𝜋 𝐿 2 aE_{\pi}^{\rm pred}=\sqrt{\left(am_{\pi}\right)^{2}+\left(a\pi/L\right)^{2}} . Corresponding effective mass plots can be found in Appendix C .

	16I-G		16ID-G
Observable	RHMC	EOFA	RHMC	EOFA
$a{E}_{\pi }$	0.3175(43)	0.3097(48)	0.4457(101)	0.4614(72)
$a{E}_{\pi }^{\mathrm{pred}}$	0.31197(4)	0.31207(4)	—	—
$a{m}_K$	0.3271(22)	0.3272(28)	0.4343(34)	0.4382(24)
$a{m}_{\mathrm{res}}^{\prime }\left(m_l\right)$	0.003140(90)	0.003054(86)	0.00919(14)	0.00952(13)

Table 9. Table 9: The relative error ( ε 𝜀 \varepsilon ) and total running time for the EOFA heatbath on the 24ID ensemble before and after applying the tuning algorithm discussed in the text.

	$\epsilon$	Total Heatbath Time
Untuned	$1.52\times {10}^{-11}$	129.5 s
Tuned	$7.79\times {10}^{-11}$	68.9 s

Table 10. Table 10: Basic integrator and HMC details for the generation of the 24ID and 32ID-G ensembles. We denote the coarsest time step used to evolve the outermost level by δ ​ τ 𝛿 𝜏 \delta\tau , and the CG stopping tolerances used for the force gradient forecasting, molecular dynamics, and Monte Carlo steps by r FG subscript 𝑟 FG r_{\rm FG} , r MD subscript 𝑟 MD r_{\rm MD} , and r MC subscript 𝑟 MC r_{\rm MC} , respectively. We elaborate on the details of the integrator nesting in Table 11 .

Ensemble	Integrator	$\delta \tau$	$r_{\mathrm{FG}}$	$r_{\mathrm{MD}}$	$r_{\mathrm{MC}}$
24ID (RHMC)	Force Gradient QPQPQ	0.0833	${10}^{-5}$	${10}^{-7}$	${10}^{-10}$
24ID (EOFA)	Force Gradient QPQPQ	0.0833	${10}^{-5}$	${10}^{-7}$	${10}^{-10}$
32ID-G (RHMC)	Omelyan ($\lambda =0.22$)	0.0625	—	${10}^{-7}$	${10}^{-10}$
32ID-G (EOFA)	Force Gradient QPQPQ	0.1667	${10}^{-5}$	${10}^{-7}$	${10}^{-10}$

Table 11. Table 11: Integrator layouts for the 24ID and 32ID-G ensembles. The notation A:B for the update scheme denotes the number of steps of the next innermost integrator level (A) per step of the current level (B).

Ensemble	Level	Action	Update
24ID (RHMC)	1	${\mathrm{RatQuo}}_{1/2}(0.085,1.0)$	1:1
	2	$\mathrm{Quo}(0.00107,0.00789)$ + $\mathrm{Quo}(0.00789,0.0291)$ + $\mathrm{Quo}(0.0291,0.095)$ +	1:1
		$\mathrm{Quo}(0.095,0.3)$ + $\mathrm{Quo}(0.3,0.548)$ + $\mathrm{Quo}(0.548,1.0)$
	3	Gauge + DSDR	1:1
24ID (EOFA)	1	$\mathrm{EOFA}(0.085,1.0)$	1:1
	2	$\mathrm{Quo}(0.00107,0.00789)$ + $\mathrm{Quo}(0.00789,0.0291)$ + $\mathrm{Quo}(0.0291,0.095)$ +	1:1
		$\mathrm{Quo}(0.095,0.3)$ + $\mathrm{Quo}(0.3,0.548)$ + $\mathrm{Quo}(0.548,1.0)$
	3	Gauge + DSDR	1:1
32ID-G (RHMC)	1	${\mathrm{RatQuo}}_{1/2}(0.0001,0.007)$	1:1
	2	${\mathrm{RatQuo}}_{1/2}(0.007,1.0)$ + ${\mathrm{RatQuo}}_{1/4}(0.045,1.0)$	1:2
	3	DSDR	1:2
	4	Gauge	1:1
32ID-G (EOFA)	1	$\mathrm{EOFA}(0.0001,0.0058)$ + $\mathrm{EOFA}(0.0058,0.0149)$ + $\mathrm{EOFA}(0.0149,0.059)$ +	5:1
		$\mathrm{EOFA}(0.059,0.177)$ + $\mathrm{EOFA}(0.177,0.45)$ +
		$\mathrm{EOFA}(0.45,1.0)$ + ${\mathrm{RatQuo}}_{1/4}(0.045,1.0)$
	2	DSDR	1:2
	3	Gauge	1:1

Table 12. Table 12: Strange quark timings for a single MD trajectory of the 24ID ensemble on a 256-node Blue Gene/Q partition. We compare RHMC to EOFA with (“preconditioned”) and without (“dense”) Cayley-form preconditioning.

	RHMC		EOFA (Dense)		EOFA (Preconditioned)
Step	Time (s)	%	Time (s)	%	Time (s)	%
Heatbath	42.6	2.7	340.6	15.1	68.9	15.5
Force gradient integration (total)	1485.6	94.8	1840.6	81.8	355.9	80.1
Final Hamiltonian evaluation	39.4	2.5	68.8	3.1	19.8	4.4
Total	1567.6	—	2250.0	—	444.6	—
(Total RHMC) / Total	1.0	—	0.7	—	3.5	—

Table 13. Table 13: HMC details for the production ensemble generation run (1) of Ref.  Bai  et al. ( 2015a ) , as well as 13 tuning runs after switching to EOFA light quarks (2-14). We use the following notation: “O” denotes the Omelyan integrator, “FG” denotes the force gradient integrator, “ N traj subscript 𝑁 traj N_{\rm traj} ” is the number of trajectories generated for the timing run, “acceptance” is the fraction of gauge field configurations which were accepted in the final Monte Carlo step, and “efficiency” is the ratio of the total job time per trajectory for the specified integration scheme to the total job time per trajectory of the scheme used in run 1. Entries in bold correspond to the original RHMC scheme (1) and the final, fully tuned EOFA scheme (12).

Run	Integrator Type	Light Hasenbusch Masses	$\Delta \tau$	$r_{\mathrm{MD}}$	$N_{\mathrm{traj}}$	Acceptance	Efficiency
1	O	0.007	0.0625	${\mathrm{𝟏𝟎}}^{-\mathrm{𝟖}}$	850	88%	—
2	O	—	0.0625	${10}^{-8}$	10	40%	1.2
3	FG	0.043	0.0625	${10}^{-8}$	10	100%	2.0
4	FG	0.018, 0.12	0.0625	${10}^{-8}$	10	100%	1.8
5	FG	0.0118, 0.0412, 0.23	0.0625	${10}^{-8}$	10	100%	1.7
6	FG	0.0075, 0.023, 0.11, 0.4	0.0625	${10}^{-8}$	10	100%	1.7
7	FG	0.0058, 0.0149, 0.059, 0.177, 0.45	0.0625	${10}^{-8}$	10	100%	1.5
8	FG	0.0103, 0.029, 0.12, 0.41	0.1000	${10}^{-6}$	15	67%	4.0
9	FG	0.0103, 0.029, 0.12, 0.41	0.1000	${10}^{-7}$	20	95%	3.0
10	FG	0.0103, 0.029, 0.12, 0.41	0.1667	${10}^{-7}$	20	75%	4.5
11	FG	0.0058, 0.0149, 0.059, 0.177, 0.45	0.1000	${10}^{-6}$	40	80%	3.0
12	FG	0.0058, 0.0149, 0.059, 0.177, 0.45	0.1667	${\mathrm{𝟏𝟎}}^{-\mathrm{𝟕}}$	850	93%	4.2
13	FG	0.0058, 0.0149, 0.059, 0.177, 0.45	0.2000	${10}^{-7}$	60	65%	4.5
14	FG	0.0058, 0.0149, 0.059, 0.177, 0.45	0.2000	${10}^{-8}$	25	72%	3.9

Table 14. Table 14: Measured spectral range of ℳ EOFA subscript ℳ EOFA \mathcal{M}_{\rm EOFA} , heatbath relative error ( ε 𝜀 \varepsilon ), and total time for the heatbath step ( Δ ​ t HB Δ subscript 𝑡 HB \Delta t_{\rm HB} ), using N LHSB subscript 𝑁 LHSB N_{\rm LHSB} intermediate mass preconditioning steps and an order N poles subscript 𝑁 poles N_{\rm poles} rational approximation to x − 1 / 2 superscript 𝑥 1 2 x^{-1/2} , with all CG stopping tolerances set to r MC = 10 − 10 subscript 𝑟 MC superscript 10 10 r_{\rm MC}=10^{-10} . Timings are reported for a 512-node Blue Gene/Q partition.

$N_{\mathrm{LHSB}}$	Mass Ratio	${\lambda }_{\min }$	${\lambda }_{\max }$	$N_{\mathrm{poles}}$	$\epsilon$	$\Delta t_{\mathrm{HB}}$ (s)
0	0.0001/1.0	1.0	1150	11	$6.91\times {10}^{-11}$	5263.2
1	0.0001/0.043	1.0	33.3	7	$3.50\times {10}^{-11}$	2226.6
	0.043/1.0	1.0	22.8	7	$6.82\times {10}^{-12}$
2	0.0001/0.018	1.0	13.5	6	$2.13\times {10}^{-11}$	2043.8
	0.018/0.12	1.0	6.4	5	$1.11\times {10}^{-11}$
	0.12/1.0	1.0	8.3	6	$6.18\times {10}^{-12}$
3	0.0001/0.0118	1.0	8.9	6	$6.08\times {10}^{-12}$	2307.8
	0.0118/0.0412	1.0	3.3	4	$4.09\times {10}^{-11}$
	0.0412/0.23	1.0	5.5	5	$1.29\times {10}^{-11}$
	0.23/1.0	1.0	4.3	5	$1.11\times {10}^{-11}$
4	0.0001/0.0075	1.0	5.9	5	$9.63\times {10}^{-12}$	2080.7
	0.0075/0.023	1.0	2.8	4	$3.55\times {10}^{-11}$
	0.023/0.11	1.0	4.6	5	$1.00\times {10}^{-11}$
	0.11/0.4	1.0	3.6	4	$1.98\times {10}^{-11}$
	0.4/1.0	1.0	2.5	4	$2.34\times {10}^{-11}$
5	0.0001/0.0058	1.0	4.7	5	$1.11\times {10}^{-11}$	2289.0
	0.0058/0.0149	1.0	2.3	4	$1.64\times {10}^{-11}$
	0.0149/0.059	1.0	3.7	4	$9.65\times {10}^{-11}$
	0.059/0.177	1.0	3.0	4	$4.14\times {10}^{-11}$
	0.177/0.45	1.0	2.5	4	$2.71\times {10}^{-11}$
	0.45/1.0	1.0	2.2	4	$1.64\times {10}^{-11}$

Equations

$$p_{\pi }^i=\frac{\left(2n_i+1\right)\pi }{L},n_i\in \mathbb{Z},$$

$$\det \left(\frac{\mathcal{D}(m_1)}{\mathcal{D}(m_2)}\right)=\frac{1}{\det \left({\mathcal{M}}_L\right)}\cdot \frac{1}{\det \left({\mathcal{M}}_R\right)},$$

$${\left({\mathcal{D}}_{\mathrm{DWF}}\right)}_{x{x}^{\prime },s{s}^{\prime }}$$

$${\left(D_W\right)}_{x{x}^{\prime }}$$

$$L_{s{s}^{\prime }}$$

$$\left(L_{+}\right)_{ss^{\prime}}=\left(L_{-}\right)_{s^{\prime}s}=\begin{dcases}\begin{array}[]{cc}-m\delta_{L_{s}-1,s^{\prime}},&s=0\\ \delta_{s-1,s^{\prime}},&1\leq s\leq L_{s}-1\end{array}\end{dcases}.$$

$$\begin{array}[]{cc}q_{R}=P_{+}\psi_{L_{s}-1}&q_{L}=P_{-}\psi_{0}\\ \overline{q}_{R}=\overline{\psi}_{L_{s}-1}P_{-}&\overline{q}_{L}=\overline{\psi}_{0}P_{+}\\ \end{array}.$$

$${\mathcal{D}}_{\mathrm{DWF}}={\mathcal{D}}_{\mathrm{EOFA}}\cdot \tilde{\mathcal{D}},$$

$$\begin{aligned} \big{(}\mathscr{D}_{\rm EOFA}\big{)}_{xx^{\prime},ss^{\prime}}&\equiv\left(D_{W}\right)_{xx^{\prime}}\delta_{ss^{\prime}}+\delta_{xx^{\prime}}\left(M_{+}\right)_{ss^{\prime}}P_{+}+\delta_{xx^{\prime}}\left(M_{-}\right)_{ss^{\prime}}P_{-}\\ \big{(}\widetilde{\mathscr{D}}\big{)}_{ss^{\prime}}&\equiv d\left(\delta_{ss^{\prime}}-L_{ss^{\prime}}\right)+c\left(\delta_{ss^{\prime}}+L_{ss^{\prime}}\right)\end{aligned}.$$

$$\det\big{(}\widetilde{\mathscr{D}}\big{)}=\left(\left(c+d\right)^{L_{s}}+m\left(c-d\right)^{L_{s}}\right)^{12V},$$

$$\det\left[\left(\begin{array}[]{cc}A&B\\ C&D\end{array}\right)\right]=\det\left(A\right)\det\left(D-CA^{-1}B\right)=\det\left(D\right)\det\left(A-BD^{-1}C\right)$$

$${\mathcal{M}}_{\mathrm{EOFA}}\equiv 1-k{P}_{-}{\Omega }_{-}^{†}{\left[H(m_1)\right]}^{-1}{\Omega }_{-}{P}_{-}+k{P}_{+}{\Omega }_{+}^{†}{\left[H(m_2)-{\Delta }_{+}\left(m_1,m_2\right)P_{+}\right]}^{-1}{\Omega }_{+}{P}_{+}.$$

$$⟨{\mathcal{O}}_1\cdots {\mathcal{O}}_n⟩=\frac{1}{\mathcal{Z}}\int \mathcal{D}U\left(\prod _f\mathcal{D}{\psi }_f\mathcal{D}\overline{\psi }{}_f\right)\left({\mathcal{O}}_1[U]\cdots {\mathcal{O}}_n[U]\right)e^{-S[U,\overline{\psi }{}_f,{\psi }_f]}.$$

$$\frac{1}{\mathcal{Z}}\int \mathcal{D}\psi \mathcal{D}\overline{\psi }{e}^{-\overline{\psi }M\psi }=\det \left(M\right)=\frac{1}{\det \left(M^{-1}\right)}=\frac{1}{\mathcal{Z}}\int \mathcal{D}\varphi \mathcal{D}{\varphi }^{†}{e}^{-{\varphi }^{†}{M}^{-1}\varphi },$$

$$⟨{\mathcal{O}}_1\cdots {\mathcal{O}}_n⟩\approx \frac{1}{N}\sum _{i=1}^N{\mathcal{O}}_1(U_i)\cdots {\mathcal{O}}_n(U_i).$$

$$H=\frac{1}{2}{\pi }^2+S(U).$$

$$\{\begin{array}{l}{\partial }_{\tau }{U}_{\mu }(x)={\pi }_{\mu }(x)U_{\mu }(x)\\ {\partial }_{\tau }{\pi }_{\mu }(x)=-T^a{\partial }_{x,\mu }^{a}S(U)\end{array}$$

$$P_{\mathrm{accept}}=\min \left(1,e^{-\Delta H}\right).$$

$$\det \left(\mathcal{D}\right)={\left[\det \left({\mathcal{D}}^{†}\mathcal{D}\right)\right]}^{1/2}.$$

$${\left[\det \left({\mathcal{D}}^{†}\mathcal{D}\right)\right]}^{1/2}=\frac{1}{\mathcal{Z}}\int \mathcal{D}\varphi \mathcal{D}{\varphi }^{†}{e}^{-{\varphi }^{†}({\mathcal{D}}^{†}\mathcal{D}{)}^{-1/2}\varphi }\simeq \frac{1}{\mathcal{Z}}\int \mathcal{D}\varphi \mathcal{D}{\varphi }^{†}{e}^{-{\varphi }^{†}f({\mathcal{D}}^{†}\mathcal{D})\varphi }$$

$$f(x)={\alpha }_0+\sum _{k=1}^N\frac{{\alpha }_k}{{\beta }_k+x}$$

$$\det \left(\frac{{\mathcal{D}}_{\mathrm{EOFA}}(m_1)}{{\mathcal{D}}_{\mathrm{EOFA}}(m_2)}\right)=\frac{1}{\mathcal{Z}}\int \mathcal{D}\varphi \mathcal{D}{\varphi }^{†}{e}^{-{\varphi }^{†}{\mathcal{M}}_{\mathrm{EOFA}}\varphi },$$

$$\det \left(\frac{{\mathcal{D}}_{\mathrm{DWF}}(m_1)}{{\mathcal{D}}_{\mathrm{DWF}}(m_2)}\right)={\left(\frac{{\left(c+d\right)}^{L_s}+m_1{\left(c-d\right)}^{L_s}}{{\left(c+d\right)}^{L_s}+m_2{\left(c-d\right)}^{L_s}}\right)}^{12V}\det \left(\frac{{\mathcal{D}}_{\mathrm{EOFA}}(m_1)}{{\mathcal{D}}_{\mathrm{EOFA}}(m_2)}\right).$$

$${\mathcal{M}}_{\mathrm{RHMC}}={\left[{\mathcal{D}}_{\mathrm{DWF}}^{†}{\mathcal{D}}_{\mathrm{DWF}}(m_2)\right]}^{1/4}{\left[{\mathcal{D}}_{\mathrm{DWF}}^{†}{\mathcal{D}}_{\mathrm{DWF}}(m_1)\right]}^{-1/2}{\left[{\mathcal{D}}_{\mathrm{DWF}}^{†}{\mathcal{D}}_{\mathrm{DWF}}(m_2)\right]}^{1/4}$$

$$\det \left({\mathcal{M}}^{-1}\right)=\frac{1}{\mathcal{Z}}\int \mathcal{D}\varphi \mathcal{D}{\varphi }^{†}{e}^{-{\varphi }^{†}\mathcal{M}\varphi }=\frac{1}{\mathcal{Z}}\int \mathcal{D}\varphi \mathcal{D}{\varphi }^{†}{e}^{-\frac{1}{2}{\varphi }^{†}{\Sigma }^{-1}\varphi }{e}^{{\varphi }^{†}\left(\frac{1}{2}{\Sigma }^{-1}-\mathcal{M}\right)\varphi }$$

$$-\mathrm{logdet}\left({\mathcal{M}}^{-1}\right)\approx {⟨{\varphi }_i^{†}\left(\mathcal{M}-\frac{1}{2}{\Sigma }^{-1}\right){\varphi }_i⟩}_i,$$

$$\det \left(\frac{\mathcal{D}(m_1)}{\mathcal{D}(m_2)}\right)=\det \left(\frac{\mathcal{D}(m_1)}{\mathcal{D}(m_1^{\prime })}\right)\left[\prod _{i=1}^{N_m}\det \left(\frac{\mathcal{D}(m_i^{\prime })}{\mathcal{D}(m_{i+1}^{\prime })}\right)\right]\det \left(\frac{\mathcal{D}(m_{N_m}^{\prime })}{\mathcal{D}(m_2)}\right)$$

$$m_i^{\prime }=m_1+\frac{m_2-m_1}{N_m+1}i,i=1,\dots ,N_m,$$

$$\begin{gathered} {\mathcal{M}}_{\mathrm{EOFA}}^{-1/2}\simeq {\alpha }_0+\sum _{k=1}^{N_p}{\alpha }_l[\frac{1}{{\gamma }_l}-k{P}_{-}{\Omega }_{-}^{†}{\left[H(m_1)\right]}^{-1}{\Omega }_{-}{P}_{-} \\ {+k{P}_{+}{\Omega }_{+}^{†}{\left[H(m_2)-{\Delta }_{+}(m_1,m_2)P_{+}\right]}^{-1}{\Omega }_{+}{P}_{+}]}^{-1}, \end{gathered}$$

$${\left(A+BCD\right)}^{-1}=A^{-1}-A^{-1}B{\left(C^{-1}+D{A}^{-1}B\right)}^{-1}D{A}^{-1}$$

$$\mathcal{M}_{\rm EOFA}^{-1/2}\simeq\alpha_{0}+\sum_{k=1}^{N_{p}}\alpha_{l}\gamma_{l}\Big{\{}1+k\gamma_{l}P_{-}\Omega_{-}^{\dagger}\left[H(m_{1})-\gamma_{l}\Delta_{-}(m_{1},m_{2})P_{-}\right]^{-1}\Omega_{-}P_{-}\\ -k\gamma_{l}P_{+}\Omega_{+}^{\dagger}\left[H(m_{2})-\beta_{l}\gamma_{l}\Delta_{+}(m_{1},m_{2})P_{+}\right]^{-1}\Omega_{+}P_{+}\Big{\}}.$$

Full text

I Introduction

Lattice QCD simulations are typically performed using variants of the hybrid Monte Carlo (HMC) algorithm, which includes the effects of dynamical sea quarks through the determinant of a fermion matrix evaluated by stochastically sampling a discretized QCD path integral. Conventional simulations choose the Hermitian fermion matrix $\mathcal{M}=\mathscr{D}^{\dagger}\mathscr{D}$ rather than the lattice Dirac operator $\mathcal{M}=\mathscr{D}$, since the latter, in general, has a complex spectrum, and is thus less amenable to standard numerical algorithms. While $\mathscr{D}$ describes a single quark flavor, $\mathscr{D}^{\dagger}\mathscr{D}$ describes two degenerate quark flavors with the same mass. As a result the standard HMC algorithm naturally describes the light (up and down) quarks in the isospin-symmetric limit $m_{u}=m_{d}$ considered in most lattice calculations. Simulations including single quark flavors (such as the strange or charm) are typically performed by taking an overall square root of the determinant of $\mathcal{M}=\mathscr{D}^{\dagger}\mathscr{D}$, leading to the rational hybrid Monte Carlo (RHMC) algorithm. While RHMC has found widespread usage in the lattice QCD community, RHMC calculations are typically more expensive than HMC calculations for the same input quark mass, in part because many of the techniques which have been developed to accelerate HMC simulations of degenerate quark flavor pairs are not applicable to RHMC.

A number of recent developments in the HMC algorithm used by the RBC/UKQCD collaboration have driven down the cost of simulating degenerate pairs of isospin-symmetric quarks with the same mass. These developments include: extensive force tuning via Hasenbush mass preconditioning Hasenbusch and Jansen (2003), the zMöbius domain wall fermion action McGlynn (2016), reduced $L_{s}$ approximations to the light quark determinant McGlynn (2016), and the use of implicitly restarted, mixed-precision defect correction methods in the conjugate gradient algorithm 111We elaborate on the details of our defect correction solver in Section VI.. In Table 1 we list timings for a recent large-scale calculation which utilizes these techniques.

We now find that the single-flavor strange and charm quark determinants, which we simulate using the RHMC algorithm, are collectively the most expensive part of the calculation. To address this, we have turned to exploring TWQCD’s recently proposed exact one flavor algorithm (EOFA), which allows for simulating single quark flavors without the need for RHMC Chen and Chiu (2014a). Preliminary results have suggested that EOFA simulations can outperform RHMC simulations, both in terms of computer time and a reduced memory footprint, while producing exactly the same physics Chen and Chiu (2014b); Murphy (2016).

The RBC/UKQCD collaboration’s ongoing efforts to probe direct $CP$-violation in $K\rightarrow\pi\pi$ decays provide a second motivation for exploring EOFA. The collaboration has recently reported the first calculation of the $\Delta I=1/2$ $K\rightarrow\pi\pi$ decay amplitude with physical kinematics in Ref. Bai et al. (2015a), which, when combined with previous results for the $\Delta I=3/2$ amplitude Blum et al. (2015) determines the Standard Model $CP$-violating parameters $\epsilon$ and $\epsilon^{\prime}$ entirely from first principles. An important ingredient in this calculation was the introduction of $G$-parity boundary conditions for the quark fields Kelly (2013); Bai et al. (2015a): since the pion is $G$-parity odd, the pion momenta are quantized along $G$-parity directions as

[TABLE]

allowing the ensemble parameters to be tuned such that the $K\rightarrow\pi\pi$ decay has both physical kinematics and the final pions in the ground state. Since the $G$-parity transformation $G=Ce^{i\pi I_{y}}$ is the product of charge conjugation and a $180^{\circ}$ isospin rotation about the $y$-axis — at the lattice boundary the light quark doublet transforms as $(u,d)\mapsto(\overline{d},-\overline{u})$ — the $G$-parity Dirac operator inherently describes two quark flavors. The standard lattice technique for obtaining a Hermitian, positive-definite fermion matrix — by taking the square of the Dirac operator, $\mathcal{M}=\mathscr{D}^{\dagger}\mathscr{D}$ — results in a theory with four degenerate quark flavors on a $G$-parity ensemble, and a square root is required to reduce to a two-flavor simulation. Describing the light quark pair on a $G$-parity ensemble is a particularly attractive target for EOFA, since many of the techniques we use to accelerate the calculation of the light quark determinant for ensembles with periodic boundary conditions — including defect correction solvers, the forecasted force gradient integrator Yin and Mawhinney (2011), and Hasenbusch mass preconditioning — are not applicable or of limited utility for RHMC simulations, but are expected to perform well in the context of EOFA. More generally, since there is no straightforward way to start the multishift conjugate gradient solver used for RHMC with a nonzero initial guess, techniques which rely on forecasting or restarting the solver are not applicable.

In this work we discuss the RBC/UKQCD collaboration’s implementation and tests of the exact one flavor algorithm, as well as the use of EOFA in generating gauge field configurations for our ongoing first-principles calculation of the ratio of Standard Model parameters $\epsilon^{\prime}/\epsilon$ from $\Delta I=1/2$ $K\rightarrow\pi\pi$ decays with $G$-parity boundary conditions. We have independently implemented EOFA in the Columbia Physics System (CPS), BAGEL fermion sparse matrix library (BFM), and the Grid data parallel C++ QCD library (Grid), for Shamir and Möbius domain wall fermions, with periodic, anti-periodic, and $G$-parity boundary conditions. We will demonstrate in the following sections that a significant improvement over the RHMC algorithm in terms of wall clock time is indeed possible with EOFA after introducing a variety of preconditioning and tuning techniques. Early work in this direction was presented at the 34th International Symposium on Lattice Field Theory Murphy (2016); here we will elaborate on the details and discuss our first large-scale EOFA calculation.

II The Exact One Flavor Algorithm

The exact one flavor algorithm was developed by the TWQCD collaboration and used to enable efficient simulations of single quark flavors on GPU clusters, where memory usage is a significant constraint. In Ref. Ogawa et al. (2009) the authors discuss their construction of a positive-definite pseudofermion action describing a single flavor of Wilson or domain wall quark, and elaborate on the details of this construction in Ref. Chen and Chiu (2014a). The key is their observation that a ratio of determinants of domain wall fermion (DWF) Dirac operators can be factorized as

[TABLE]

with $\mathcal{M}_{L}$ and $\mathcal{M}_{R}$ Hermitian and positive-definite. In a subsequent paper the authors benchmark EOFA against RHMC for $N_{f}=1$ and $N_{f}=2+1$ lattice QCD simulations, and demonstrate a number of advantages of the EOFA formalism Chen and Chiu (2014b). These include substantial reductions in the pseudofermion force and in the memory footprint of the algorithm, since, in the context of EOFA, inversions of the Dirac operator can be performed using the ordinary conjugate gradient (CG) algorithm rather than the multishift CG used for RHMC. They ultimately find that they are able to generate HMC trajectories 15-20% faster using EOFA rather than RHMC after retuning their integration scheme to take advantage of these properties. More recently, TWQCD has used EOFA to generate $N_{f}=2+1+1$ domain wall fermion ensembles with dynamical strange and charm quarks Chen and Chiu (2017).

We note that the construction of the exact one flavor pseudofermion action has been detailed by TWQCD in Ref. Ogawa et al. (2009); Chen and Chiu (2014a) and summarized in our own formalism in Ref. Murphy (2016). We will not repeat this discussion here, other than to give a brief overview and to introduce the notation used in this work. We write the 5D Möbius domain wall fermion (MDWF) operator $\mathscr{D}_{\rm DWF}$ in terms of the 4D Wilson Dirac operator $D_{W}$ and 5D hopping matrix $L_{ss^{\prime}}$ as

[TABLE]

with

[TABLE]

Here $x$ and $s$ are spacetime indices in the 4D bulk and along the fifth dimension, respectively, with $L_{s}$ denoting the total number of $s$ sites, $P_{\pm}=(1\pm\gamma_{5})/2$ denoting the chiral projection operators, and $(R_{5})_{ss^{\prime}}\equiv\delta_{s,Ls-1-s^{\prime}}$ denoting the operator which performs a reflection in the fifth dimension. We recover four-dimensional quark fields $q$ and $\overline{q}$ with definite chiralities from the five-dimensional quark fields $\psi$ and $\overline{\psi}$ described by $\mathscr{D}_{\rm DWF}$ at the boundaries of the fifth dimension

[TABLE]

Green’s functions constructed from $q$ and $\overline{q}$ approximate continuum QCD arbitrarily well in the limit of vanishing lattice spacing and infinite 5D spacetime volume. The tunable parameters in Eqn. (3) are the domain wall parameter $M_{5}$, the bare quark mass $m$, and the Möbius scale $\alpha=2c$; the parameter $d$ is fixed at $d=1/2$. DWF with the Shamir kernel is recovered from the more general Möbius operator in the limit $\alpha\rightarrow 1$. For more detail regarding our MDWF formalism we direct the reader to Ref. Blum et al. (2016).

The construction of the exact one flavor action for domain wall fermions begins by factorizing the MDWF Dirac operator as 222We emphasize that the operator we call $\mathscr{D}_{\rm EOFA}$ is merely an alternative notation for TWQCD’s rescaled Dirac operator $D_{T}$ Chen and Chiu (2014a).

[TABLE]

with

[TABLE]

The operator $\widetilde{\mathscr{D}}$ relating $\mathscr{D}_{\rm DWF}$ and $\mathscr{D}_{\rm EOFA}$ has no dependence on the gauge field, so we are free to replace $\mathscr{D}_{\rm DWF}$ with $\mathscr{D}_{\rm EOFA}$ in Eqn. (2) without modifying physical observables described by a properly normalized path integral. In fact, it can be shown analytically using the explicit form of $\widetilde{\mathscr{D}}$ listed in Appendix A that

[TABLE]

where $V=L^{3}T$ is the 4D spacetime volume. This substitution facilitates the construction of a proper action since the operator $\gamma_{5}R_{5}\mathscr{D}_{\rm EOFA}$ is manifestly Hermitian for any choice of the Möbius scale $\alpha$, whereas $\mathscr{D}_{\rm DWF}$ satisfies a less trivial $\gamma_{5}$-Hermiticity condition when $\alpha\neq 1$ Brower et al. (2005a). However, this comes at the cost of substantially more expensive inversions, since $\mathscr{D}_{\rm EOFA}$ is dense in $ss^{\prime}$ whereas $\mathscr{D}_{\rm DWF}$ has a well-known tridiagonal block structure.

After introducing $\mathscr{D}_{\rm EOFA}$, TWQCD’s construction proceeds by applying the Schur identity

[TABLE]

to $\mathscr{D}_{\rm EOFA}$, treated as a $2\times 2$ block matrix in its spinor indices, and rearranging terms to arrive at the right-hand side of Eqn. (2). Crucially, factors of $\gamma_{5}R_{5}$ can be freely inserted under the determinant to replace $\mathscr{D}_{\rm EOFA}$ with the Hermitian operator $H\equiv\gamma_{5}R_{5}\mathscr{D}_{\rm EOFA}$, since $\det(\gamma_{5})=\det(R_{5})=1$. The final form of the exact one flavor pseudofermion action is $S_{\rm EOFA}=\phi^{\dagger}\mathcal{M}_{\rm EOFA}\phi$, with

[TABLE]

In Appendix A we collect explicit expressions for $k$, $\Omega_{\pm}$, $\Delta_{\pm}$, $\mathscr{D}_{\rm EOFA}$, and $\widetilde{\mathscr{D}}$ for Shamir and Möbius DWF, since, to the authors’ knowledge, these expressions have not previously appeared in the literature. In Ref. Chen and Chiu (2014a) these operators are constructed recursively for the more general case of DWF with weights $\rho_{s}=c\omega_{s}+d$ and $\sigma_{s}=c\omega_{s}-d$ that are allowed to vary along the fifth dimension, subject to the constraint that $\omega_{s}$ is reflection-symmetric in $s$. Shamir and Möbius DWF are simpler, special cases with $\omega_{s}=1$.

III Summary of Ensembles Used in This Work

The properties of the lattices used in this work are summarized in Tables 2 and 3. In all cases we use the Iwasaki gauge action (I) Iwasaki and Yoshié (1984), and on some ensembles supplement this with the dislocation suppressing determinant ratio (DSDR) Vranas (2006); Renfrew et al. (2008); we abbreviate the combined action including both terms as “ID”. The additional DSDR term is designed to suppress the dislocations of the gauge field associated with tunneling between topological sectors, thereby reducing the degree of residual chiral symmetry breaking. For strong coupling simulations, where these dislocations occur frequently, the DSDR term reduces the costs associated with light quark masses while still maintaining good topological sampling. We simulate $N_{f}=2+1$ quark flavors using domain wall fermions, with either the Shamir (DWF) Kaplan (1992); Shamir (1993) or Möbius (MDWF) Brower et al. (2005b, 2006, 2012) kernel. Finally, on ensembles marked “-G” we use $G$-parity boundary conditions in one or more of the spatial directions.

The 16I ensemble was first generated and used to study light meson spectroscopy with domain wall fermions in Ref. Allton et al. (2007). The 16I-G ensemble is identical to the 16I ensemble except for the boundary conditions along the $x$-direction, which have been changed from periodic to $G$-parity. Likewise, the parameters of the 16ID-G ensemble have been chosen based on a series of $\beta=1.75$ DSDR ensembles generated in Ref. Arthur et al. (2013), but have $G$-parity boundary conditions in all three spatial directions. Collectively, these three lattices are used as inexpensive, small-volume test ensembles with unphysical, heavy pion masses to perform cross-checks of the EOFA algorithm and its implementation in the BFM and CPS code libraries. The larger 24ID Mawhinney (2016) and 32ID-G Bai et al. (2015b) ensembles have physical pion masses and are currently being generated as part of production RBC/UKQCD calculations.

IV Hybrid Monte Carlo with EOFA

In lattice QCD correlation functions are computed in terms of a discretized Euclidean path integral

[TABLE]

Here $U$ is the gauge field, $\psi_{f}$ is the quark field associated with flavor $f$, and $S[U,\overline{\psi}\,_{f},\psi_{f}]$ is the action, which decomposes into a sum of contributions from the gauge field, fermions, and possibly other terms (e.g. the dislocation suppressing determinant ratio). To avoid having to deal with anticommuting Grassman variables in a computer, dynamical fermion flavors are integrated out and then reintroduced in terms of bosonic “pseudofermion” fields $\phi$ as

[TABLE]

provided $M$ is positive-definite. While pseudofermions can be represented straightforwardly in a computer, they come at the cost of applications of $M^{-1}$ rather than $M$, which is not typically available in an explicit form. Even after discretization the integration in Eqn. (11) is far too expensive to perform directly due to the enormous number of degrees of freedom on a typical lattice. Instead, Monte Carlo techniques are used to ergodically sample a sequence of representative configurations of the gauge field $\{U_{i}\}$, for which

[TABLE]

The standard Monte Carlo technique used in modern lattice QCD calculations is known as the Hybrid Monte Carlo (HMC) algorithm.

HMC generates a Markov chain of gauge field configurations $\{U_{i}\}$ by evolving a Hamiltonian system in unphysical Molecular Dynamics (MD) “time”. This Hamiltonian system is constructed by treating $U_{\mu}(x)$ as a generalized coordinate, introducing an $\mathfrak{su}(3)$-valued conjugate momentum $\pi_{\mu}(x)$, and forming the standard Hamiltonian

[TABLE]

The associated equations of motion

[TABLE]

can then be integrated using numerical integration techniques. The integration is performed over intervals of length $\Delta\tau$ — referred to as a single MD trajectory — as a sequence of $N$ small steps $\delta\tau$, with $N=\Delta\tau/\delta\tau$. Finite precision integration errors are corrected stochastically with a Metropolis accept/reject step: after every $N$ integration steps by $\delta\tau$ the total change in the Hamiltonian $\Delta H$ is computed, and the current gauge field $U_{\mu}^{\prime}(x)$ is accepted as the next configuration in the Markov chain with probability

[TABLE]

One can show that the resulting algorithm satisfies detailed balance provided the scheme used to numerically integrate Eqn. (15) is reversible Rothe (2012). Ergodicity is achieved by performing a heatbath step each time the integration is restarted to pick a new conjugate momentum $\pi_{\mu}(x)$, and thus a new trajectory in the phase space $\{(U,\pi)\}$. HMC generates a sequence of gauge field configurations whose statistical independence is governed by the length of each MD trajectory, $\Delta\tau$. The number of MD trajectories separating statistically independent gauge field configurations is typically determined ex post facto by examining the integrated autocorrelation times of representative physical observables.

The fermionic contribution to the Hamiltonian in Eqn. (14) introduces a technical obstacle for the HMC algorithm since, generically, the lattice Dirac operator $\mathscr{D}$ has a complex spectrum. Replacing $\mathscr{D}$ with the Hermitian fermion matrix $M=\mathscr{D}^{\dagger}\mathscr{D}$ in Eqn. (12) has a number of advantages. Most importantly, it allows $M^{-1}$ to be applied to pseudofermion vectors using the conjugate gradient algorithm, and it allows for a straightforward pseudofermion heatbath step: at the beginning of each MD trajectory a random Gaussian vector $\eta$ is drawn according to $P(\eta)\propto\exp(-\eta^{\dagger}\eta/2)$ and the initial pseudofermion field is seeded as $\phi=\mathscr{D}\eta$, ensuring that $\phi$ is correctly sampled as $P(\phi)\propto\exp(-\phi^{\dagger}M^{-1}\phi/2)$. However, the fermion matrix $M=\mathscr{D}^{\dagger}\mathscr{D}$ describes two degenerate quark flavors with the same mass. Single flavor simulations are typically performed by taking an overall square root of the fermion determinant,

[TABLE]

In the pseudofermion formalism applications of the operator $(\mathscr{D}^{\dagger}\mathscr{D})^{-1/2}$ are approximated by a matrix-valued function $f(\mathscr{D}^{\dagger}\mathscr{D})$,

[TABLE]

where $f(x)$ is a suitably constructed approximation to the inverse square root, valid over the spectral range of $\mathscr{D}^{\dagger}\mathscr{D}$. Variants of the HMC algorithm which construct $f$ from different classes of functions have been proposed and used in the literature; the most common is the rational HMC (RHMC) algorithm Clark (2006), where

[TABLE]

is a rational function. While rational functions are in many ways a good choice — they are economical in the sense that the inverse square root can usually be well-approximated by a modest number of terms, and the multishift CG algorithm can be used to efficiently invert $(\mathscr{D}^{\dagger}\mathscr{D}+\beta_{k})$ for all $k$ simultaneously — the additional complexity of evaluating $f(\mathscr{D}^{\dagger}\mathscr{D})$ and the associated molecular dynamics pseudofermion force makes single flavor RHMC simulations significantly more costly than degenerate two flavor HMC simulations at the same bare quark mass. This additional cost can be largely attributed to the significant linear algebra overhead associated with multishift CG.

EOFA provides an alternative construction of a single-flavor pseudofermion action through Eqn. (2): a ratio of fermion determinants can be factorized as a product of two determinants, each of which involves an operator which is Hermitian and positive-definite. This product can then be represented as a path integral over a bosonic pseudofermion field with a two-term action (Eqn. (10))

[TABLE]

leading to an algorithm which is “exact” in the sense that it avoids the numerical approximations required to implement the square root in RHMC (Eqn. (18)) and related HMC variants. EOFA is also expected to be somewhat faster than RHMC, since there is no rational approximation entering into evaluations of the Hamiltonian or the pseudofermion force, eliminating the overhead associated with multishift CG. In the remainder of this section we elaborate on the details of the action, heatbath step, and pseudofermion force entering into the Hamiltonian equations of motion (Eqn. (15)) for HMC with EOFA.

IV.1 Action

The EOFA action (Eqn. (10)) computes a ratio of determinants of $\mathscr{D}_{\rm EOFA}$ upon integrating out the pseudofermion fields (Eqn. (20)). This ratio can be related to the conventional determinant ratio computed by the RHMC algorithm through Eqns. (6) and (8), leading to the relationship

[TABLE]

We use this relationship as a test of the equivalence of RHMC and EOFA, as well as our implementation of the EOFA action, by stochastically computing the left side of Eqn. (21) with the RHMC action

[TABLE]

and the right side with the EOFA action (Eqn. (10)) on the same gauge field configuration. Observing that we can, in general, rewrite a determinant as

[TABLE]

suggests the following simple Monte Carlo integration scheme: we draw random pseudofermion vectors by independently sampling the real and imaginary parts of each component from the standard normal distribution $\mathcal{N}(\mu=0,\sigma=1)$, and compute the expectation value

[TABLE]

where the average is computed using the jackknife resampling technique. This will accurately approximate the true log determinant for finite, realistically calculable values of $N$ provided $m_{1}$ and $m_{2}$ are sufficiently close that the integrand is well-approximated by a Gaussian with unit variance. To address this latter systematic, we consider splitting Eqn. (21) as a product of determinants

[TABLE]

with equally-spaced intermediate masses

[TABLE]

and study the dependence of the result on $N_{m}$ (this procedure is identical to the method introduced in Ref. Aoki et al. (2011) for computing quark mass reweighting factors). In the upper panel of Figure 1 we plot the log determinants of $\mathcal{M}_{\rm RHMC}^{-1}$ and $\mathcal{M}_{\rm EOFA}^{-1}$ as a function of $N_{m}$, with $N=10$ stochastic evaluations, computed using a single thermalized trajectory of the 16I, 16I-G, and 16ID-G ensembles. For the case of the 16ID-G ensemble, which uses the Möbius DWF fermion action, we also include the overall constant multiplying the right side of Eqn. (21) so that in all cases we are computing the same determinant ratio of $\mathscr{D}_{\rm DWF}$ using either action.

We observe, as expected, that both formalisms agree for sufficiently large $N_{m}$. Likewise, we observe that at sufficiently small $N_{m}$ the calculation generally becomes unreliable since we do not attempt to account for the systematic error associated with approximating the integrand of the path integral by a Gaussian with unit variance (i.e. setting $\Sigma=1$ in Eqns. (23) and (24).). In both cases “sufficiently” small or large $N_{m}$ is controlled by the size of the splitting between $m_{1}$ and $m_{2}$. We also observe that, for a given choice of $N_{m}$ and $N$, both the statistical and systematic errors of the determinant ratio computed via EOFA are suppressed relative to the errors of the determinant ratio computed via RHMC. We argue that the observed error suppression can be explained by comparing the spectrum of $\mathcal{M}_{\rm RHMC}$ to the spectrum of $\mathcal{M}_{\rm EOFA}$, which we plot in the lower panels of Figure 1 for a very small lattice volume ($4^{5}$) where the complete spectrum can be computed directly. While both operators have similar condition numbers, we find that most of the spectrum of the EOFA action is concentrated into a small interval $[1,1+\Delta]$ with $\Delta\sim\mathcal{O}(0.1)$, leading to an action which is easier to estimate stochastically.

We propose that TWQCD’s EOFA construction can be thought of as a kind of preconditioning which computes the same determinant ratio as RHMC but modifies the operator inside the determinant ($\mathcal{M}_{\rm RHMC}$), mapping its spectrum onto a more compact interval. This suggests an additional application of the EOFA formalism: quark mass reweighting factors can be computed substantially more cheaply using the EOFA action than using the RHMC action, especially at light quark masses, even if the ensemble was generated using RHMC. This could be useful, for example, to include the dynamical effects of isospin breaking in ensembles generated with isospin-symmetric up and down quarks.

IV.2 Heatbath

At the beginning of each HMC trajectory we wish to draw a random pseudofermion field $\phi$ according to the distribution $P(\phi)\propto\exp(-\phi^{\dagger}\mathcal{M}_{\rm EOFA}\phi)$. To do this, we first draw a random vector $\eta$ by independently sampling the real and imaginary parts of each component from the normal distribution with $\mu=0$ and $\sigma^{2}=1/2$, and then compute $\phi=\mathcal{M}_{\rm EOFA}^{-1/2}\eta$. As with the RHMC algorithm we approximate the inverse square root by an appropriately constructed rational function, but we stress that in the context of EOFA this rational approximation enters only into the heatbath and is not necessary to compute the EOFA action itself or the associated pseudofermion force. Naively applying a rational approximation with the form of Eqn. (19) to the operator $\mathcal{M}_{\rm EOFA}$ results in

[TABLE]

where we have defined $\gamma_{l}\equiv(1+\beta_{l})^{-1}$. In this form, the nested inversions required to seed the heatbath would make EOFA prohibitively expensive. However, the Woodbury matrix identity

[TABLE]

and the cancellation between cross-terms involving products of the chiral projection operators can be used to manipulate this expression into the equivalent form

[TABLE]

With this expression the EOFA heatbath step can be performed at the cost of $2N_{p}$ CG inversions using a rational approximation with $N_{p}$ poles. Unlike the case of RHMC, multishift CG algorithms are not applicable to the EOFA heatbath since each of the $2N_{p}$ operators in Eqn. (29) generates a different Krylov space. Furthermore, since the operators $\Delta_{\pm}P_{\pm}$ have a large number of zero modes and are therefore not invertible, there is no simple transformation by which this system can be recast into a form amenable to multishift CG.

In the left panel of Figure 2 we test Eqn. (29) on a single thermalized configuration of the 16I ensemble by computing the quantity

[TABLE]

after seeding the pseudofermion field $\phi$ with a random Gaussian vector $\eta$. In exact arithmetic $\varepsilon=0$; in practice it measures the relative error in the heatbath step arising from the choice of CG stopping conditions and rational approximation to the inverse square root. We repeat this calculation, varying the number of poles in the rational approximation but keeping the stopping conditions fixed, and observe that $\varepsilon$ reaches the limits of double-precision arithmetic even with a relatively modest number of poles compared to what is typically required to compute non-integer powers of $\mathscr{D}^{\dagger}\mathscr{D}$ accurately in the context of RHMC. In the right panel of Figure 2 we demonstrate this explicitly by computing the condition numbers $\kappa=\lambda_{\rm max}/\lambda_{\min}$ of both operators as a function of the bare input quark mass. In Section VI we show how aggressive tuning of the rational approximation and stopping conditions, together with forecasting techniques for the initial CG guesses, can be combined to ameliorate the cost of the $2N_{p}$ inversions required to apply $\mathcal{M}_{\rm EOFA}^{-1/2}$.

IV.3 Pseudofermion Force

The pseudofermion force

[TABLE]

measures the back-reaction of the pseudofermions on the HMC system (Eqn. (15)) under an infinitesimal variation of the gauge field. In our notation $\{T^{a}\}$ is a basis for the Lie algebra $\mathfrak{su}(3)$, with the CPS normalization convention

[TABLE]

The EOFA pseudofermion force can be worked out explicitly by differentiating the EOFA action (Eqn. (10)) and applying the matrix identity

[TABLE]

resulting in

[TABLE]

with

[TABLE]

and

[TABLE]

Standard manipulations can be used to write a more general Dirac bilinear as

[TABLE]

allowing Eqn. (34) to be efficiently computed locally in terms of a trace over spinor indices, at the cost of the two inversions required to form $\chi_{1}$ and $\chi_{2}$. Since Dirac bilinears of the form $a^{\dagger}\left(\partial_{x,\mu}^{a}D_{W}\right)b$ enter into the pseudofermion forces associated with many of the standard pseudofermion actions for Wilson and domain wall fermions — including the RHMC action — implementing the EOFA pseudofermion force requires little new code beyond what is required to implement the EOFA Hamiltonian.

V Small Volume Reproduction Tests

To further test our implementation of EOFA we have reproduced the 16I (16I-G, 16ID-G) ensemble using EOFA for the strange quark (light quarks) in place of RHMC. For these tests we have made no serious effort to tune EOFA for performance; we have simply checked that replacing RHMC with EOFA, but leaving all other details of the simulation fixed, has no discernible impact on physical observables such as the average plaquette, topological susceptibility, and low energy spectrum.

V.1 Ensemble Generation

The details of the integrator parameters and nesting are summarized in Tables 4 and 5, respectively. We use the abbreviations

[TABLE]

and

[TABLE]

to denote the quotient and rational quotient actions, and on the EOFA reproduction ensembles replace each instance of $\mathrm{RatQuo}_{1/2}(m_{1},m_{2})$ with the EOFA action (Eqn. (10)) using the same mass parameters. The 16I and 16I-G EOFA reproduction runs were seeded with an ordered start — i.e. all gauge links were initially set to the unit matrix — and evolved for 1500 and 2500 MD time units, respectively. For the 16ID-G ensemble the last RHMC configuration (MD trajectory 908) was used to seed the start of the EOFA reproduction run, and then evolved for an additional 500 MD time units.

V.2 Basic Observables

In Table 6 we summarize results for the average plaquette $\langle P\rangle$, light quark and strange quark chiral $\langle\overline{\psi}\psi\rangle$ and pseudoscalar $\langle\overline{\psi}\gamma_{5}\psi\rangle$ condensates, and the topological susceptibility $\chi_{t}\equiv\langle Q^{2}\rangle/V$, computed on each ensemble; we observe statistically consistent results between the RHMC and EOFA ensembles in each case. Accompanying plots of the time evolution can be found in Section C. The topological charge $Q$ has been measured using the 5Li discretization introduced in Ref. de Forcrand et al. (1997) after cooling the gauge fields with 20 steps of APE smearing Albanese et al. (1987) using a smearing coefficient of 0.45. The ensemble averages were computed after binning over 50 (25) successive MD time units on the 16I and 16I-G (16ID-G) ensembles, where the bin size has been conservatively chosen based on the integrated autocorrelation times measured in Ref. Allton et al. (2007) for the 16I ensemble and Ref. Arthur et al. (2013) for a series of $\beta=1.75$ DSDR ensembles. We expect that the runs produced for this study are too short to reliably compute integrated autocorrelation times directly, but note that there is no evidence of a difference in an integrated autocorrelation time between the EOFA and RHMC ensembles in the time evolution plots of Section C.

V.3 Low Energy Spectra

In Table 7 we list results for the pion, kaon, Omega baryon, and residual masses, computed on the 16I ensemble. These calculations were performed using a measurement package previously introduced in Ref. Blum et al. (2016), and based on the all-mode averaging (AMA) technique of Ref. Blum et al. (2013). Five exact light quark propagators were computed per trajectory using a deflated, mixed-precision CG solver with 600 low-mode deflation vectors and a tight stopping residual $r=10^{-8}$, while sloppy propagators were computed for all time slices using a reduced stopping residual $r=10^{-4}$. Strange quark propagators were computed with the tight residual $r=10^{-8}$ for all time slices using ordinary CG with no deflation. AMA correlation functions were then computed by time-translational averaging of the sloppy propagators, using the available exact propagators to correct for bias. The light quark propagators were computed using Coulomb gauge-fixed wall (W) sources, with either local (L) or wall sinks; the strange quark propagators were computed using Coulomb gauge-fixed wall or $Z_{3}$ box ($Z_{3}B$) sources, and in both cases local sinks.

The pion and kaon masses were extracted by fitting to the asymptotic, large Euclidean time limit of the respective two-point correlation function,

[TABLE]

where $\mathscr{O}$ denotes the choice of interpolating operator, $X\in\{\pi,K\}$ is the ground state to which $\mathscr{O}$ couples, and $V$ and $T$ are the spatial volume and temporal extent of the lattice, respectively. In particular, we performed simultaneous fits to the $\langle PP^{LW}\rangle$, $\langle PP^{WW}\rangle$, and $\langle AP^{LW}\rangle$ correlators, with $P(x)=\overline{\psi}(x)\gamma_{5}\psi(x)$ and $A(x)=\overline{\psi}(x)\gamma_{5}\gamma_{4}\psi(x)$, and the first (second) superscript denotes the sink (source) type. The Omega baryon mass was extracted from the two-point correlation function

[TABLE]

with the interpolating operator

[TABLE]

$s_{1}=L$ and $s_{2}\in\{W,Z_{3}B\}$. The correlators were then projected onto the positive parity component

[TABLE]

and simultaneously fit to double exponential ansätze with common mass terms

[TABLE]

Finally, the residual mass was determined by fitting the ratio

[TABLE]

to a constant, where $j_{5q}^{a}$ is the five-dimensional pseudoscalar density evaluated at the midpoint of the fifth dimension, and $j_{a}^{5}$ is the physical pseudoscalar density constructed from the surface fields.

In addition, we have also measured the ground state pion energy, kaon mass, and residual mass on the 16I-G and 16ID-G ensembles. While the ground state of the kaon is at rest, the ground state of the pion has nonzero momentum $\vec{p}_{100}=(\pm\pi/L,0,0)$ on the 16I-G ensemble and $\vec{p}_{111}=(\pm\pi/L,\pm\pi/L,\pm\pi/L)$ on the 16ID-G ensemble due to the boundary conditions. These calculations make use of an extension of the AMA measurement package described above to $G$-parity ensembles; as discussed in Ref. Kelly (2013), this requires the inclusion of additional diagrams that are generated by the mixing of quark flavors at the lattice boundaries through the $G$-parity operation. We measure on 51 configurations of the 16I-G ensemble, beginning with trajectory 500 and with a separation of 40 MD time units, and use sloppy and exact CG stopping tolerances of $10^{-4}$ and $10^{-10}$, respectively, with a single exact solve per trajectory. We likewise measure on 21 configurations of the 16ID-G ensemble, beginning with trajectory 500 (960) for the RHMC (EOFA) ensemble and separated by 20 MD time units, and use the same AMA setup. We perform no additional binning for either ensemble since the separations between consecutive measurements are already comparable to the bin sizes used to compute the plaquette and quark condensates.

V.4 Pseudofermion Forces on the 16I Ensemble

TWQCD has observed that the average EOFA pseudofermion force is roughly half the size of the corresponding average RHMC pseudofermion force for a particular dynamical $N_{f}=1$ QCD simulation with domain wall quarks performed in Ref. Chen and Chiu (2014a). Following this observation, we examine the forces on the RHMC and EOFA variants of the 16I ensemble. We define a norm on the space of $\mathfrak{su}(3)$-valued pseudofermion force matrices $F_{\mu}^{a}(x)\equiv\partial_{x,\mu}^{a}S(U)$ by

[TABLE]

and consider two measures of the force associated with a given configuration of the gauge field: the first is the RMS force

[TABLE]

and the second is the maximum force

[TABLE]

in both cases taken over all lattice sites and link directions. While we expect Equation (48) to be a more pertinent definition in the context of HMC simulations — we have empirically found that acceptance is controlled by the size of $F_{\rm max}$ — both $F_{\rm RMS}$ and $F_{\rm max}$ are, a priori, reasonable measures of the pseudofermion force.

In Figure 3 we compare histograms of $F_{\rm RMS}$ and $F_{\rm max}$ between the RHMC and EOFA 16I ensembles. Each data point corresponds to a single evaluation of the pseudofermion force falling between MD trajectories 500 and 1500. We find that comparing the relative sizes of the RHMC and EOFA forces is highly dependent on whether one chooses $F_{\rm RMS}$ or $F_{\rm max}$; the mean EOFA $F_{\rm RMS}$ is roughly 30% smaller than the mean RHMC $F_{\rm RMS}$, but the distributions of $F_{\rm max}$ are nearly indistinguishable. This observation suggests that while the EOFA force distribution may have a smaller mean than the RHMC force distribution, the EOFA distribution also likely has longer tails, such that the largest forces have similar magnitudes. Since we expect the magnitude of the largest forces to correlate more strongly with the efficiency of the integrator than the magnitude of the average forces, as we have argued above, we interpret these results as suggesting that the optimal step size for an EOFA evolution should be similar to that of an RHMC simulation with the same mass parameters, even if the average force is somewhat smaller.

TWQCD has also observed a large hierarchy of scales in the pseudofermion forces associated with each of the two terms in Eqn. (34); in Ref. Chen and Chiu (2014b) they find that the average force associated with the first term — involving the left-handed component of the pseudofermion field — is more than an order of magnitude smaller than the average force associated with the second term — involving the right-handed component — for two different dynamical QCD simulations. They exploit this observation with a Sexton-Weingarten integration scheme, integrating the first term with a larger time step than the second, and find increased efficiency in their simulations. In Figure 4 we compare histograms of the RMS and maximum left-handed and right-handed forces from 1000 thermalized configurations of the 16I EOFA ensemble. Our conclusions are analogous to the comparison between the EOFA and RHMC forces: if one considers $F_{\rm RMS}$ the left-handed force contribution is indeed substantially smaller than the right-handed force contribution, but if one instead considers $F_{\rm max}$ the force distributions are very similar in both magnitude and shape. Based on the latter observation we leave both terms in Equation (34) on the same time step in our large-scale EOFA simulations.

We also note that in these small volume test runs we have not considered applying the Hasenbusch mass preconditioning technique Hasenbusch and Jansen (2003) to the EOFA formalism. Introducing a set of Hasenbusch masses $\{m^{\prime}_{i}\}_{i=1}^{N}$, with $m_{1}<m^{\prime}_{i}<m_{2}$, we can write the fermion determinant as

[TABLE]

While the left-hand side can be simulated using a single pseudofermion field, the associated forces can be large if $m_{1}\ll m_{2}$, requiring a small step size to maintain reasonable acceptance. The right-hand side, in contrast, involves $N+1$ independent pseudofermion fields, but with possibly substantially reduced forces, allowing larger step sizes to be used. For light $m_{1}$ one typically observes that the gain from increasing the step size offsets the cost of simulating extra heavy flavors, leading to a more efficient simulation. In Section VII we demonstrate that Hasenbusch preconditioning allows for a substantial speed-up in the context of the 32ID-G ensemble. We also note that in addition to reducing the size of the pseudofermion forces, the Hasenbusch technique preconditions the EOFA force in the sense that the size hierarchy between the left-handed and right-handed force contributions to a single determinant disappears in the limit $m_{i}^{\prime}\to m_{i+1}^{\prime}$. In practice, we find that the mass preconditioned simulation has comparable left-handed and right-handed force contributions even in the RMS sense.

VI Optimization and Tuning

In this section we discuss preconditioning and algorithmic techniques which reduce the cost of EOFA simulations. In some cases these are extensions of well-known lattice techniques to the EOFA formalism, while in other cases they are specific to EOFA. We illustrate these techniques using bechmark tests computed with the physical quark mass, Möbius DWF 24ID ensemble, and report timing results for code written in the Columbia Physics System (CPS) and running on 256-node or 512-node Blue Gene/Q partitions.

VI.1 Inversions of $\mathscr{D}_{\rm EOFA}$

Since the majority of the computational effort in an HMC simulation is associated with repeatedly inverting the Dirac operator, techniques to more efficiently apply the Dirac operator or to otherwise accelerate these inversions can have a dramatic impact on the overall efficiency of the integrator. To address the former, we make use of the BAGEL assembler generation library Boyle (2009) to produce highly optimized kernels and fermion solvers for the Blue Gene/Q hardware. To address the latter, we make use of multiple preconditioning techniques, as well as a mixed precision defect correction CG solver.

The first preconditioning technique we apply — “even-odd” or “red-black” preconditioning — is well-known in the lattice QCD community. Lattice sites are labeled as even if $\left(x+y+z+t\right)\equiv 0\pmod{2}$, or odd if $\left(x+y+z+t\right)\equiv 1\pmod{2}$, inducing a $2\times 2$ block structure on fermion operators

[TABLE]

Standard tricks can then be used to relate the linear system $M\psi=\phi$ to a better conditioned linear system involving only the odd sub-lattice; this preconditioned system is substantially cheaper to invert since the size of the problem has been halved. After inverting on the odd sub-lattice, the even component of $\psi$ can also be recovered at modest cost, without ever needing to explicitly invert on the even sub-lattice. The details of this construction, and its extension to EOFA, are described in Appendix B.

The second preconditioning technique we apply — Cayley-form preconditioning — is unique to EOFA, and was introduced in Ref. Murphy (2016). The generic linear system one needs to solve in the context of EOFA has the form

[TABLE]

where $H=\gamma_{5}R_{5}\mathscr{D}_{\rm EOFA}$. For Möbius domain wall fermions $\mathscr{D}_{\rm EOFA}$ is dense in $ss^{\prime}$, and thus considerably more expensive to invert than $\mathscr{D}_{\rm DWF}$, which has a tridiagonal $ss^{\prime}$ stencil, in terms of wall clock time. However, Eqn. (6) suggests that Eqn. (51) can be related to an equivalent system in terms of $\mathscr{D}_{\rm DWF}$ by using $\widetilde{\mathscr{D}}^{-1}$ as a preconditioner. We elaborate on the mathematical details in Appendix B.2, and, in particular, demonstrate that $\Delta_{\pm}\widetilde{\mathscr{D}}$ has a relatively simple, rank-one form, allowing for substantially more efficient EOFA inversions — even when $\beta\neq 0$ — by working with the preconditioned system. This technique also has the advantage that it allows for EOFA simulations which re-use existing high-performance code for applying $\mathscr{D}_{\rm DWF}$ with little modification.

Finally, we use a restarted, mixed precision defect correction solver to perform the conjugate gradient inversions of the fully preconditioned EOFA system. For memory bandwidth-limited calculations — such as applying the Dirac operator — single precision computations can be performed at approximately half the cost of full double precision computations. In the defect correction approach to mixed precision CG, the following algorithm is used:

Solve the Dirac equation in single precision arithmetic using a reduced stopping tolerance (typically $10^{-4}$ or $10^{-5}$). 2. 2.

Compute the current residual using the (single precision) solution in full double precision arithmetic. 3. 3.

If the desired final tolerance (typically $10^{-8}$ or smaller) has been reached, stop. Otherwise, return to step 1, using the residual vector computed in step 2 as the new CG source.

We observe that this algorithm outperforms straight double precision CG by approximately a factor of 2 — as one would expect if the calculation is truly memory bandwidth-limited — provided the local lattice volume on each node is sufficiently large to avoid communications bottlenecks.

In Figure 5 we plot the CG residual as a function of the wall clock running time of the inverter for a series of benchmark inversions of Equation (51) on the 24ID ensemble. These benchmarks show the inverter performance as we sequentially introduce even-odd preconditioning, Cayley-form preconditioning, and finally, mixed precision CG. We also plot the time required to solve the family of linear systems

[TABLE]

using multishift CG for the same set of poles $\{\beta_{k}\}$ used in the rational approximation to $x^{-1/2}$ in the RHMC evolution that generated the 24ID ensemble. This allows a baseline estimate of the cost of evaluating the EOFA Hamiltonian or pseudofermion force against the cost of evaluating the RHMC Hamiltonian or pseudofermion force at the same quark mass. We observe a factor of 3.9 speed-up for fully preconditioned EOFA over the even-odd preconditioned RHMC system. In both cases the underlying operator being inverted is $\mathscr{D}_{\rm DWF}$; the slower RHMC benchmark demonstrates the overhead associated with multishift CG relative to solving a single system with standard CG, both due to the inability to fully utilize mixed precision methods and due to the additional linear algebra required at each iteration.

VI.2 Heatbath

Achieving the full performance improvement suggested by the inversion benchmarks in Section VI.1 is complicated by the form of the EOFA heatbath, which is expected to be more expensive than the RHMC heatbath, even with efficient EOFA code. Applying $\mathcal{M}_{\rm EOFA}^{-1/2}$ (Eqn. (29)) requires two independent CG inversions per pole used in the rational approximation to $x^{-1/2}$, since multishift CG is not applicable: we use two algorithmic techniques to reduce this cost. The first is a forecasting technique initially proposed by Brower et al. Brower et al. (1997) in the context of more general HMC simulations, and later used successfully by TWQCD in the context of the EOFA heatbath Chen and Chiu (2014a). The idea is the following: given a set of solutions $\{\psi_{k}\}_{k=1}^{N}$ to Equation (51) for $N$ different poles $\{\beta_{k}\}_{k=1}^{N}$, one can use the linear combination

[TABLE]

minimizing the functional

[TABLE]

as the initial CG guess for the next inversion with pole $\beta_{N+1}$. The coefficients $c_{k}$ satisfy

[TABLE]

and can be computed explicitly using e.g. Gauss-Jordan elimination. Since Equation (54) is the same functional minimized by the conjugate gradient algorithm itself, accurate initial guesses can be computed for modest $N$ provided the $\{\beta_{k}\}_{k=1}^{N+1}$ are similar in magnitude. In Figure 6 we test this forecasting technique using the 24ID ensemble and a rational approximation with 8 poles, and find that the iteration count required to solve Eqn. (51) to a tolerance of $10^{-10}$ is more than halved for the last few poles.

The second technique we have used to accelerate the heatbath is motivated by observing that the coefficients entering into Equation (29) span several orders of magnitude for a typical rational approximation to $x^{-1/2}$. We find typical values $k\alpha_{l}\gamma_{l}^{2}/\alpha_{0}\sim\mathcal{O}(10^{-3}-10^{-5})$, suggesting that the inversions can be performed with reduced stopping tolerances relative to the desired accuracy of $\mathcal{M}_{\rm EOFA}^{-1/2}\psi$, since the solution vectors are ultimately multiplied by small coefficients when the result is formed. We have explored the following simple optimization scheme to relax the stopping conditions for each pole:

Choose a desired tolerance for the heatbath, $\varepsilon_{\rm tol}$, where $\varepsilon$ is defined by Equation (30). 2. 2.

Choose one of the inversions required to compute $\mathcal{M}_{\rm EOFA}^{-1/2}$ according to Equation (29), and relax the stopping tolerance until the overall error in the heatbath $\varepsilon$ reaches $\varepsilon_{\rm tol}$. 3. 3.

Iterate over each inversion until all stopping conditions have been tuned.

We report results for the 24ID ensemble in Table 9. Using a rational approximation with 6 poles, and $\varepsilon_{\rm tol}=10^{-10}$, we observe that the total heatbath time is more than halved while only slightly increasing the error. We have also checked that the final error and heatbath running time after tuning is insensitive to the exact order in which the stopping tolerances are tuned.

VII Large-Scale EOFA Calculations

In this section we turn to two ongoing ensemble generation calculations currently being performed by the RBC/UKQCD collaboration. The first is a strong-coupling $N_{f}=2+1$ $24^{3}\times 64\times 24$ Iwasaki+DSDR lattice (24ID) intended for exploratory studies and calculations requiring high statistics Mawhinney (2016). The second (32ID) has been used for a first-principles calculation of the ratio of Standard Model $CP$-violation parameters $\epsilon^{\prime}/\epsilon$ from $\Delta I=1/2$ $K\rightarrow\pi\pi$ decays in Ref. Bai et al. (2015a). RBC/UKQCD is currently generating more gauge field configurations to reduce the statistical errors in the $\Delta I=1/2$ decay amplitudes. Both ensembles have physical quark masses and large volumes, allowing for tests of the performance of EOFA in the context of state-of-the-art domain wall fermion calculations.

Tables 10 and 11 summarize the details of the integrator parameters and nesting for these evolutions. The ensembles labeled RHMC correspond to the evolutions of Ref. Mawhinney (2016) (24ID) and Ref. Bai et al. (2015a) (32ID-G). For the ensembles marked EOFA, we have changed the strange quark (light quark) action to EOFA for the 24ID (32ID-G) ensemble and retuned the details of the integrator as described in the remainder of the section. For the 32ID-G ensemble — where, due to the $G$-parity flavor doubling, the EOFA action naturally describes the degenerate light quark pair — we have also switched from an Omelyan integrator to a force gradient integrator, and inserted additional Hasenbusch preconditioning determinants.

VII.1 24ID Ensemble

We use the 24ID ensemble as a straightforward benchmark of RHMC against an equivalent EOFA simulation to describe a physical heavy quark flavor. Here this is the strange quark, but $N_{f}=2+1+1$ simulations with dynamical strange and charm quarks are another obvious target of EOFA. We make no serious attempt to retune the integrator after switching to EOFA beyond tuning the heatbath step with the following procedure:

Compute the largest and smallest eigenvalues of $\mathcal{M}_{\rm EOFA}$ (Eqn. (10)) for a few thermalized configurations of the gauge field, and use these measurements to inform the bounds of the rational approximations to $x^{-1/2}$ constructed via the Remez algorithm. 2. 2.

Add poles to the rational approximation, with all CG stopping tolerances set to $r_{\rm MC}$, until $\varepsilon<r_{\rm MC}$ (Eqn. (30)) is reached. 3. 3.

With the rational approximation now fixed from step 2, tune the CG stopping tolerances corresponding to each pole, following the procedure described in Section VI.2, and keeping $\varepsilon<r_{\rm MC}$.

After tuning the heatbath, we then ran a single trajectory of the RHMC evolution and the EOFA evolution on a 256-node Blue Gene/Q partition. For the EOFA ensemble, we compare two schemes. The first (“dense”) is a straightforward implementation of Möbius DWF as proposed in Ref. Chen and Chiu (2014a): we invert Equation (51) directly, where $H=\gamma_{5}R_{5}\mathscr{D}_{\rm EOFA}$ and the other dense 5D operators appearing in the EOFA action are listed explicitly in Appendix A.2. We also do not apply the final step in our heatbath tuning procedure, leaving all CG stopping tolerances in the heatbath fixed at $r_{\rm MC}=10^{-10}$. In the second EOFA scheme (“preconditioned”) we fully tune the heatbath step and apply the Cayley-form preconditioning detailed in Appendix B.2 to inversions of Equation (51). Timing breakdowns for the strange quark part of the evolution are reported in Table 12.

We observe that the dense EOFA formalism is actually somewhat slower than RHMC: the additional complexity of the EOFA heatbath, together with the more expensive inversions of the dense 5D operator $\mathscr{D}_{\rm EOFA}$, negate the expected performance gains from the simpler forms of the Hamiltonian and force evaluations. We emphasize, however, that we have made no attempt to retune the integrator details to optimize for EOFA. After introducing Cayley-form preconditioning — so that we are inverting the tridiagonal operator $\mathscr{D}_{\rm DWF}$ rather than $\mathscr{D}_{\rm EOFA}$ when we solve Equation (51) — we find that EOFA outperforms RHMC by a significant factor of 3.5.

VII.2 32ID-G Ensemble

One particularly promising feature of EOFA in the context of $G$-parity ensembles is the potential for aggressive Hasenbusch mass preconditioning of the light quark determinant; this makes the 32ID-G ensemble a particularly interesting case study since the EOFA formalism is used to describe a physical mass light quark pair. In Ref. Bai et al. (2015a) the RBC/UKQCD collaboration observed that mass preconditioning is not particularly effective for the RHMC light quark determinant, since each molecular dynamics step requires one multishift inversion of $\mathscr{D}^{\dagger}\mathscr{D}$ evaluated at the numerator quark mass and two multishift inversions of $\mathscr{D}^{\dagger}\mathscr{D}$ evaluated at the denominator quark mass. The latter two solves become prohibitively expensive if many intermediate masses are introduced, negating the expected gain from integrating the preconditioned pseudofermion forces with larger step sizes. The EOFA force, on the other hand, is no more expensive to evaluate than the force associated with the standard quotient action (Eqn. (38)), so it is natural to expect better performance from Hasenbusch preconditioning.

In Table 13 we list details of the tuning runs we have used to explore potential schemes for evolving the 32ID-G ensemble with EOFA light quarks. We started by switching from an Omelyan integrator, for which the leading errors are $\mathcal{O}(\delta\tau^{2})$, to a force gradient integrator, for which the leading errors are $\mathcal{O}(\delta\tau^{4})$, and studied the effects of inserting mass preconditioning determinants one at a time (runs 1-7). We then identified two promising mass preconditioning schemes — one with four intermediate masses (runs 8-10), and the other with five intermediate masses (runs 11-14) — and continued tuning the step size, CG stopping conditions, and heatbath, to optimize the job time per trajectory and Monte Carlo acceptance. The initial RHMC scheme used in Ref. Bai et al. (2015a) corresponds to run 1, and the final EOFA scheme we have adopted for our continuing ensemble generation corresponds to run 12.

We find, in practice, that Hasenbusch mass preconditioning is extremely effective for the EOFA light quark determinant. In addition to reducing the size of the pseudofermion force, we also observe that the largest eigenvalue of the EOFA action, Equation (10), decreases rapidly as $m_{2}\rightarrow m_{1}$. As a consequence, the heatbath is also less expensive with Hasenbusch preconditioning, since, as we increase the number of intermediate masses, we can simultaneously decrease the range and number of poles entering into the rational approximation used for each determinant. Table 14 summarizes the measured spectral range, the heatbath error, and the total heatbath cost for each of the runs 2-7. For this ensemble the first Hasenbusch mass reduces the cost of the heatbath by more than a factor of two, and subsequent Hasenbusch masses essentially leave the cost fixed.

For each of the runs 2-7 we generated ten trajectories, beginning from the same seed configuration, and analyzed the resulting distributions of $F_{\rm RMS}$ and $F_{\rm max}$. In panel (a) of Figure 7 we plot distributions of $F_{\rm max}$ from 850 trajectories of the production RHMC ensemble generation calculation (run 1). Since we are using exactly the same RHMC action for the strange quark on the RHMC and EOFA ensembles, we tune by adjusting the number and magnitude of the intermediate light Hasenbusch masses such that the forces associated with each of the light quark determinants are comparable to the strange quark force. This allows us to simplify the integrator layout to a three-level scheme, with the light and strange quark determinants updated on the same level. We find that four intermediate Hasenbusch masses are sufficient to ensure that the strange quark force is dominant in the sense of $F_{\rm RMS}$, and that five intermediate Hasenbusch masses are sufficient in the sense of $F_{\rm max}$. Panel (b) shows an analogous force distribution for the latter mass preconditioning scheme.

In runs 8-10 we explore further tuning of a scheme with four light Hasenbusch masses, and in runs 11-14 we explore further tuning of a scheme with five light Hasenbusch masses. We note that the Monte Carlo acceptance is relatively poor in runs 8-10 — as we argued in Section V, this is consistent with the view that the acceptance should be controlled by the largest integration errors accrued during the trajectory, which are proportional to $F_{\rm max}$ rather than $F_{\rm RMS}$ — and thus have abandoned this mass preconditioning scheme in favor of the scheme used in runs 11-14. We have then tuned the step size of the outermost integrator level ($\delta\tau$) and the CG stopping tolerance used in the molecular dynamics evolution ($r_{\rm MD}$) to minimize the mean time required to generate an accepted gauge field configuration, resulting in the scheme of run 12. In addition, we have applied the heatbath tuning procedure described in Section VI.2 in all of the runs 8-14, allowing us to relax CG stopping tolerances for the individual solves in the heatbath, while keeping the overall error bounded by $r_{\rm MC}=10^{-10}$. For the final scheme (12) this optimization further reduced the cost of the light quark heatbath from approximately 2300 s, as reported in Table 14, to approximately 850 s after tuning.

Comparing the fully tuned EOFA scheme (12) to the original RHMC scheme (1) in Table 13, we find that we are able to generate EOFA trajectories a factor of 4.2 times faster than RHMC trajectories, while maintaining a slightly higher acceptance rate of 93%. We emphasize, however, that this improved performance is only partially attributable to the simpler form of the EOFA Hamiltonian and force evaluations: we have also switched from an Omelyan integrator to a force gradient integrator, retuned the step sizes and integrator layout, and, in some cases, applied optimizations to the EOFA simulation that are not applicable to RHMC simulations (e.g. mixed precision CG). Figure 8 briefly summarizes the respective techniques used in the RHMC and EOFA evolution schemes. We have now adopted the EOFA scheme tested in run 12 for ensemble generation in our ongoing $\Delta I=1/2$ $K\rightarrow\pi\pi$ calculation Kelly (2016). We expect the resulting performance gain to enable up to four times as many measurements in our current production run as we would have been able to generate using the initial RHMC evolution scheme, enabling a significantly more precise first-principles determination of the Standard Model ratio $\epsilon^{\prime}/\epsilon$.

VIII Conclusion

In this work we have explored the viability of the exact one flavor algorithm (EOFA) as an alternative to the rational Hybrid Monte Carlo (RHMC) algorithm in molecular dynamics simulations of lattice QCD with domain wall fermions and periodic or $G$-parity boundary conditions. We have verified the formal equivalence of EOFA to RHMC through statistical tests of the EOFA action (Section IV), and checked, using a series of inexpensive, small volume ensembles with heavy pions, that physical observables such as the plaquette, quark condensates, topological susceptibility, and low energy spectrum are consistent between ensembles generated using EOFA and ensembles generated using RHMC (Section V). We have then discussed preconditioning and tuning techniques for EOFA simulations (Section VI and Appendix B), and finally, demonstrated that EOFA can substantially outperform RHMC for state-of-the-art lattice QCD simulations with large volumes and physical quark masses (Section VII). In particular, we find that we are able to generate gauge field configurations for the ongoing RBC/UKQCD calculation of the $\Delta I=1/2$ $K\rightarrow\pi\pi$ decay amplitudes a factor of 4.2 times faster with EOFA. The keys to this dramatic speed-up are a novel preconditioning technique which relates inversions of the EOFA Dirac operator ($\mathscr{D}_{\rm EOFA}$) to cheaper inversions of the standard domain wall fermion Dirac operator ($\mathscr{D}_{\rm DWF}$), and the ability to apply mixed precision defect correction solvers and extensive Hasenbusch mass preconditioning in the context of EOFA.

Future work will explore further physics applications of EOFA. We intend to generate variants of the 24ID ensemble with non-degenerate up and down quark masses in the near future. These ensembles will enable exploratory studies of isospin breaking effects in the meson and baryon spectra, as well as in other precision lattice calculations such as the extraction of the CKM matrix element $V_{us}$ from semileptonic kaon decays Boyle et al. (2015a). Other potential applications include domain wall QCD simulations with dynamical charm quarks in the sea, and simulations with light, $SU(3)$-symmetric quarks. The latter simulations could be used, for example, to better constrain the strange quark dependence of our $SU(3)$ chiral perturbation theory studies Mawhinney and Murphy (2016), or to probe the location of the critical point separating the crossover and first-order phase transition regions in three-flavor domain wall QCD at finite temperature.

Acknowledgments

The authors are grateful to Peter Boyle, Ting-Wai Chiu, and Norman Christ for helpful discussions in support of this work. Calculations were performed using the Blue Gene/Q computers of the RIKEN-BNL Research Center and Brookhaven National Lab. The software used includes the CPS QCD code (https://github.com/RBC-UKQCD/CPS) C. Jung (2014) (RBC and UKQCD Collaborations), supported in part by the USDOE SciDAC program, the BAGEL (http://www2.ph.ed.ac.uk/~paboyle/bagel/Bagel.html) assembler kernel generator for high-performance optimized kernels and fermion solvers Boyle (2009), and the Grid data parallel C++ QCD library (https://github.com/paboyle/Grid) Boyle et al. (2015b). R.D.M. and D.J.M. are supported in part by U.S. DOE grant #DE-SC0011941. C.J. is supported in part by U.S. DOE Contract #AC-02-98CH10886 (BNL). C.K. is supported by the Intel Corporation.

Appendix A EOFA Operators for Shamir and Möbius DWF

In this appendix we list the operators which enter into $\mathscr{D}_{\rm EOFA}$ (Eqn. (6)) and the EOFA action (Eqn. (10)). The more general case of DWF with Zolotarev-type domain wall fermions is constructed implicitly in Ref. Chen and Chiu (2014a); we explicitly list these operators for the more restrictive cases of Shamir and Möbius DWF used in our simulations. We use $\Theta_{s}$ to denote the discrete Heaviside theta function

[TABLE]

and assume even $L_{s}$. The operators $\Omega_{\pm}$ and $\Delta_{\pm}$ are related by the identity

[TABLE]

and the Möbius operators reduce to the corresponding Shamir operators in the limit $c=d=1/2$. We note that the dense Möbius expressions listed here are not used inside the inverter; we instead invert the preconditioned system discussed in Appedix B.2.

A.1 Shamir Kernel

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A.2 Möbius Kernel

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix B Four-Dimensional Even-Odd Preconditioning

The inversions required to compute the exact one flavor Hamiltonian can be accelerated using a standard checkerboarding technique: we label lattice sites as “even” if $\left(x+y+z+t\right)\equiv 0\pmod{2}$ or “odd” if $\left(x+y+z+t\right)\equiv 1\pmod{2}$. This naturally induces a block structure in the Dirac operator $\mathscr{D}$, which can be $LDU$ decomposed as

[TABLE]

Left-multiplying the linear system $\mathscr{D}\psi=\phi$ by

[TABLE]

results in the equivalent system

[TABLE]

leading to the following trick: assuming $\mathscr{D}_{ee}^{-1}$ is available in an explicit form, it suffices to invert

[TABLE]

with $\widetilde{\phi}_{o}\equiv\phi_{o}-\mathscr{D}_{oe}\mathscr{D}_{ee}^{-1}\phi_{e}$. This system only involves the odd sublattice, and is thus substantially cheaper to invert than $\mathscr{D}$ using an iterative algorithm like CG. The solution on the even sublattice can then be reconstructed for a trivial additional cost as

[TABLE]

This technique is already well understood in the context of RHMC with Shamir or Möbius DWF; in this appendix we describe how to generalize the method to the exact one flavor algorithm.

In the context of EOFA, the generic linear system one needs to invert takes the form

[TABLE]

We choose to multiply by an overall factor of $\gamma_{5}R_{5}$, rewriting the system as

[TABLE]

for the following reasons: first, we wish to re-use the existing high-performance implementation of the Wilson $\not{D}$ kernel in the BAGEL library without modification, and second, overall factors of $\gamma_{5}R_{5}$ will cancel inside the inverter since we use CG applied to the normal equations and $(\gamma_{5}R_{5})^{\dagger}(\gamma_{5}R_{5})=1$. Since $\left(\mathscr{D}_{\rm EOFA}\right)_{eo}=\left(\mathscr{D}_{\rm DWF}\right)_{eo}$ and $\Delta_{\pm}\propto\delta_{xx^{\prime}}$ in the 4D bulk, only the operators coupling sites of the same parity need to be modified to implement even-odd preconditioned EOFA. We take somewhat different approaches for the Shamir and Möbius cases.

B.1 Shamir Kernel

Recall that for the Shamir kernel $\mathscr{D}_{\rm DWF}=\mathscr{D}_{\rm EOFA}$, so the extension of an inverter for the even-odd preconditioned $\mathscr{D}_{\rm DWF}$ operator to instead solve Eqn. (82) is straightforward. With $\mathscr{D}=\mathscr{D}_{\rm DWF}$, the same parity fermion matrix has the tridiagonal block structure

[TABLE]

One can check by explicit calculation that the shift operators have the form

[TABLE]

so one can consider the operator appearing in Eqn. (82) as a slight generalization of Eqn. (83) to

[TABLE]

with

[TABLE]

and

[TABLE]

where the index $i$ denotes the chirality of the shift operator. $\mathscr{D}_{ee}^{-1}$ can be efficiently applied using the $LDU$ decomposition of $\mathscr{D}_{ee}$, again as a slight generalization of the standard Shamir DWF case.

B.2 Möbius Kernel and Cayley-Form Preconditioning

Using Eqn. (7) we can write $\mathscr{D}_{\rm EOFA}$ in the form

[TABLE]

The action of the operator appearing in Eqn. (82) on lattice sites of the same parity, then, is given by

[TABLE]

with $M_{+}$, $M_{-}$, and $\Delta_{\pm}$ as defined in equations (70)-(73). The matrix elements of $\mathscr{D}_{ee}^{-1}=\mathscr{D}_{oo}^{-1}$ can be found by explicit numerical inversion as part of the setup cost; this is a trivial overhead since it suffices to invert only the $ss^{\prime}$ subblock of $\mathscr{D}_{ee}$. In this form, the exact factorization of the fermion determinant in Eqn. (2) comes at the cost of dense $L_{s}\times L_{s}$ matrix operations. We argue that it is possible to do significantly better by introducing an additional preconditioning step.

We note that the system defined by Eqn. (82) can be more efficiently inverted for the case of Möbius DWF by using the operator $\widetilde{\mathscr{D}}^{-1}$ as a right preconditioner, resulting in an equivalent system in terms of $\mathscr{D}_{\rm DWF}$. For the special case $\beta=0$ this is straightforward: observing that the relationship between $\mathscr{D}_{\rm EOFA}$ and $\mathscr{D}_{\rm DWF}$ (Eqn. (6)) can be used to manipulate

[TABLE]

it suffices to solve $\mathscr{D}_{\rm DWF}\psi^{\prime}=\gamma_{5}R_{5}\phi$, from which $\psi=\widetilde{\mathscr{D}}\psi^{\prime}$ can be recovered at the cost of a single additional matrix multiplication. While we observe that $\mathscr{D}_{\rm DWF}^{\dagger}\mathscr{D}_{\rm DWF}$ has a slightly larger condition number than $\mathscr{D}_{\rm EOFA}^{\dagger}\mathscr{D}_{\rm EOFA}$, leading to a modest increase in the total number of CG iterations required to invert the system, $\mathscr{D}_{\rm DWF}$ also has a tridiagonal stencil in the fifth dimension, and can thus be applied in $\mathcal{O}(L_{s})$ operations — unlike the $\mathcal{O}(L_{s}^{2})$ operations required for the dense $\mathscr{D}_{\rm EOFA}$ — leading to a substantial reduction in wall clock time for the inversion.

The $\beta\neq 0$ case is more involved, but can be treated in a similar manner. Right preconditioning Eqn. (82) with $\widetilde{\mathscr{D}}^{-1}$ leads to

[TABLE]

where we have used $\gamma_{5}P_{\pm}=\pm P_{\pm}$. We define a new, preconditioned, shift operator $\widetilde{\Delta}_{\pm}$ by

[TABLE]

and note that since $(\widetilde{\Delta})_{eo}=(\widetilde{\Delta})_{oe}=0$, Eqn. (91) can be inverted efficiently even with $\beta\neq 0$ provided we can apply the operator $(\mathscr{D}_{\rm DWF})_{ee}\pm\beta\widetilde{\Delta}_{\pm}$ and its inverse in $\mathcal{O}(L_{s})$ operations. This turns out to be possible after observing that $\widetilde{\Delta}_{\pm}$ is rank-one, i.e. it can be written as a vector outer product

[TABLE]

To see this, we start by decomposing $\widetilde{\mathscr{D}}$ into its chiral components — $\widetilde{\mathscr{D}}=\widetilde{\mathscr{D}}_{+}P_{+}+\widetilde{\mathscr{D}}_{-}P_{-}$ — in terms of which we can also decompose

[TABLE]

The 5D structure of these operators can be worked out by direct calculation, leading to Eqn. (93), with

[TABLE]

Matrix-vector products involving the preconditioned shift operator and a pseudofermion field can be computed from this decomposition as

[TABLE]

and

[TABLE]

The inverses can be applied using the Sherman-Morrison formula:

[TABLE]

In terms of

[TABLE]

which can be constructed numerically using the tridiagonal matrix algorithm Datta (2010), the necessary factors can be written as

[TABLE]

and

[TABLE]

which allow Eqn. (98) to be applied to a pseudofermion vector in $\mathcal{O}(L_{s})$ operations.

In Figure 9 we benchmark representative even-odd preconditioned inversions of Eqn. (82) on the 24ID ensemble, with and without additional preconditioning by $\widetilde{\mathscr{D}}^{-1}$, at the physical strange quark mass. In addition to observing a substantial improvement in terms of wall clock time for the inversion, we note that this preconditioning scheme also has the advantage that it requires little new code — assuming an existing high-performance implementation of $\mathscr{D}_{\rm DWF}$ — since $\mathscr{D}_{\rm EOFA}$ is never applied directly in the preconditioned formalism.

Appendix C Additional Plots for Small Volume Reproduction Tests

C.1 Evolution of the Plaquette, Quark Condensates, and Topological Charge

C.2 Effective Mass Plots