# Preconditioning the non-relativistic many-fermion problem

###### Abstract

Preconditioning is at the core of modern many-fermion Monte Carlo algorithms, such as Hybrid Monte Carlo, where the repeated solution of a linear problem involving an ill-conditioned matrix is needed. We report on a performance comparison of three preconditioning strategies, namely Chebyshev polynomials, strong-coupling approximation and weak-coupling expansion. We use conjugate gradient (CG) on the normal equations as well as stabilized biconjugate gradient (BiCGStab) as solvers and focus on the fermion matrix of the unitary Fermi gas. Our results indicate that BiCGStab is by far the most efficient strategy, both in terms of the number of iterations and matrix-vector operations.

## 1 Introduction

The Hybrid Monte Carlo (HMC) HMC_Duane (); HMC_Gottlieb () algorithm provides an efficient way of calculating the properties of many-body Fermi systems. HMC achieves high efficiency by performing global updates by means of Molecular Dynamics (MD) trajectories that are integrated with a finite time step. At each step, the solution of a linear problem

(1) |

is needed. In this problem, is a large but typically sparse matrix which provides a real-space representation of the fermion operator that defines the dynamics of the system. While is typically neither symmetric nor positive definite, it is always possible to consider solving Eq. (1) via conjugate gradient (CG) iteration, as satisfies these properties CG (). However, as CG is a Krylov-type method based of successive matrix-vector (MV) operations using , the efficiency of CG depends critically on the condition number of that matrix. Unfortunately, in most cases of interest may be extremely ill-conditioned, leading to a rapid growth of the number of CG iterations required to reach a preset convergence criterion. Such a situation manifests itself in many problems at low enough temperatures or strong enough couplings. Apart from the obvious disadvantages of a rapidly growing computational cost, any iterative method may eventually become unstable if the number of iterations grows without bounds. For these reasons, CG and similar Krylov methods are preferentially used together with a suitable preconditioner. It should be pointed out that while direct methods such as LU decomposition avoid the problems related to ill-conditioning, problems of physical interest are typically large enough such that the storage of is not feasible. It is thus desirable to resort to “matrix-free” methods based on MV operations only, which require significantly less memory and scale much more favorably with system size.

Whenever is ill-conditioned, the problem is much more severe for . An attractive option is then to perform the solve the problem in two steps, by considering the linear systems given by

(2) |

However, since is typically neither symmetric nor positive definite, these equations require more sophisticated algorithms such as biconjugate gradient (BiCG) BiCG () or stabilized biconjugate gradient (BiCGStab) vanderVorst (), which is a modern development of the former. One is thus dealing with a much less ill-conditioned problem, at the price of solving two successive problems. Nevertheless, preconditioning remains instrumental in accelerating and stabilizing the solution process. In its simplest form, preconditioning amounts to considering the linear problem

(3) |

where represents an approximation to . A similar procedure applies to Eq. (2) given approximations to the inverses of and . The procedure in Eq. (3) is referred to as “left preconditioning”. So-called “right preconditioning” is also possible by insertion of between and , followed by multiplication with once the solution has been obtained. In this work, we shall be mainly concerned with the former approach.

While the choice of is in principle arbitrary, preconditioning is only useful if the solution of Eq. (3) to a given accuracy requires fewer CG iterations than the solution of Eq. (1). Moreover, as iterative methods tend to accumulate roundoff error, an iterative solution may not even be possible without preconditioning. In such cases preconditioning also serves to stabilize the algorithm by enabling convergence. Two obvious extreme choices are those of setting and setting , however no benefits are achieved in either case, as the first one amounts to no preconditioning at all and the second one requires the solution of the original problem in order to determine . Therefore, what one seeks is compromise between the computational simplicity of the first option and the effectiveness of the second one.

Over the years, many preconditioning strategies were developed and tested on a wide range of problems (see e.g. Chapter 11 in Ref. BroydenVespucci (), or Chapters 9 and 10 in Ref. Saad ()). In general, preconditioners specifically designed for a particular problem tend to be superior to generic approaches that build in no knowledge about the physics of the system in question. With this in mind, we will explore several approaches to preconditioning which combine some knowledge of the features of , and , in order to find a strategy which is both computationally efficient and effective in accelerating the convergence of the linear problem at hand.

In this work, we are mainly concerned with problems that arise in the context of HMC simulations of non-relativistic, (3+1)-dimensional many-fermion systems. The HMC approach forms a central part of state-of-the-art Lattice QCD calculations LQCDReview (), which involve a relativistic system of fermions. The popularity of HMC stems from the great algorithmic efforts of the Lattice QCD community to enable simulations on large lattices. On the other hand, HMC remains little known in other areas such as condensed-matter physics CondMatt () and nuclear structure Leeetal (), where many of the problems are non-relativistic and determinantal Monte Carlo (DMC) DMC () remains the method of choice. The most serious drawback of DMC is the poor scaling of the computation time with the lattice volume , namely . In contrast, the use of global updates and iterative solvers enables the HMC algorithm to scale as . In spite of this obvious advantage, the widespread use of HMC implementations for non-relativistic fermions has been hampered by problems related to the operator , which is severely ill-conditioned at low temperatures (see e.g. AbeSeki ()). It is therefore of great interest to explore various preconditioning strategies, as these hold the key to finding a practical HMC implementation which is competitive with current DMC calculations.

The class of problems at hand is characterized by a zero-range interaction, such that the Hamiltonian is given by

(4) |

where the kinetic energy is

(5) |

and the potential energy

(6) |

such that and are creation and annihilation operators for particles of spin at point , respectively. We next seek to recast the grand canonical partition function

(7) |

in a form amenable to a Monte Carlo calculation. We proceed by discretizing imaginary time into slices of length using a Trotter-Suzuki decomposition Suzuki () followed by a continuous Hubbard-Stratonovich transformation (HST) HS () of the form recently proposed by Lee Lee (). It should be noted that these last two steps are, to a certain extent, a matter of choice. For instance, continuous-time formulations exist Rombouts () although they have as yet not been combined with HMC. The choice of HST is also dictated by computational preferences Hirsch (), although it is not clear at this time whether an HMC implementation is compatible with a discrete HST. The end result is a path integral representation of the partition function

(8) |

(9) |

where the explicit form of the block matrices

(10) |

and , which results from the specific choice of HST. In momentum space, the kinetic energy is given by

(11) |

where we have chosen the dispersion relation for our calculations. The integer vectors assume values on a three-dimensional momentum lattice given by , where . Eq. (9) provides a representation of the fermionic operator referred to earlier, which is real and manifestly not symmetric.

One of the most popular and successful implementations of HMC is known as the -algorithm HMC_Gottlieb (), which combines the stochastic evaluation of the fermion determinant in Eq. (8) with the MD evolution of the auxiliary Hubbard-Stratonovich field (originally the gauge field in the case of Lattice QCD). To this end, the pseudofermion representation

(12) |

is applied, giving

(13) |

for the corresponding action, where . The next step is to introduce an auxiliary field that will play the role of a conjugate momentum to the field in the molecular dynamics evolution. This is accomplished by multiplying the partition function by a constant in the form of a path integral over ,

(14) |

where the MD Hamiltonian is given by

(15) |

In the pseudofermion formulation, field configurations are sampled by evolving and according to the classical MD equations of motion that follow from the MD Hamiltonian. It should be noted that the pseudofermion field may be sampled exactly from a Gaussian heat bath and is kept constant during the MD evolution. In practice, the MD evolution should be performed using a reversible symplectic integrator in order to ensure detailed balance. At each step in the MD evolution the “fermion force”

(16) |

is calculated, which requires frequent computation of the vector . This represents the most time-consuming part of the HMC algorithm. One of the most obvious advantages of the -algorithm is that the direct calculation of a large determinant is avoided, in favor of the solution of a much smaller linear problem a fixed number of times. Specifically, in our case this involves the repeated solution of Eq. (1) in its preconditioned form given by Eq. (3), where the structure of the matrix is constant but the auxiliary field varies at each MD step. This should be contrasted with the cost of computing the full inverse, which is much more expensive than solving a single linear problem. By avoiding the calculation of the determinants and inverses, this algorithm enables global updates of the auxiliary field via MD evolution. Moreover, all HMC approaches can be rendered free of systematic errors associated with a finite MD integration step size by means of a Metropolis accept/reject step at the end of each MD trajectory.

In Sec. 2, we study the effect of preconditioning via Chebyshev polynomials, which provide an approximation to . Chebyshev polynomials provide a generic approach which is easy to apply, but computationally expensive. In Sec. 3, we proceed to study the weak-coupling approximation (WCA) method of preconditioning and its transpose separately. In Sec. 4 we explore the strong-coupling approximation (SCA) to which was originally proposed in Ref. Scalettaretal (). Our findings are summarized in Sec. 5, followed by a concluding discussion in Sec. 6.

## 2 Chebyshev polynomials

Chebyshev polynomials have frequently been applied in the preconditioning of relativistic quantum field theory problems. In fact, they have led to the development of a whole new class of HMC-type algorithms commonly referred to as Polynomial Hybrid Monte Carlo (PHMC) FrezzottiJansen (), in which preconditioners are used to separate modes that evolve at different rates in an MD trajectory. The Chebyshev preconditioning technique is based on approximating with a Chebyshev polynomial of degree ,

(17) |

where the coefficients are given by

(18) |

and the roots by

(19) |

such that the parameter provides a lower bound for the range of validity of the approximation. Indeed, in the interval the Chebyshev polynomials provide a good approximation, in the sense that the relative error

(20) |

where

(21) |

is exponentially suppressed with the degree of the polynomial.

The preconditioner we seek is obtained from Eq. (17) upon replacing , assuming that the eigenvalue spectrum of is normalized in the range . One of the most appealing features of the Chebyshev preconditioning method is that its speed and effectiveness can be tuned by varying the degree of the polynomial . Moreover, since is applied in its factorized form, the speed of the preconditioner is directly linked to the speed with which is applied. The order in which the factors of Eq. (17) are applied turns out to be critical for numerical stability, due to cancellations between the various terms that occur naturally when evaluating a polynomial in finite precision arithmetic. In practice it becomes crucial beyond to address this issue. Various orderings of the factors were explored in Ref. Bunk (), as well as a recursive Clenshaw algorithm which largely eliminates the accumulation of round-off error. Our experience indicates that the performance of the “bit-reversal” ordering of Ref. Bunk () is nearly identical to that of the Clenshaw algorithm, and it is slightly faster and requires less memory, in addition to being much simpler to implement. The bit-reversal ordering will therefore be our method of choice.

A serious disadvantage of the Chebyshev approach is that it is only applicable to the preconditioning of the normal equations, as it requires the matrix to have positive definite eigenvalue spectrum. Also, as becomes extremely ill-conditioned at low temperatures, polynomial orders larger than are needed for effective preconditioning, which makes the Chebyshev approach rather computationally intensive. As will become evident in Sec. 5, this problem can largely be avoided by separately inverting and , which however requires a different type of preconditioner. We shall now turn to the description of such strategies.

## 3 Weak coupling expansion

In this section, we demonstrate how the weak-coupling limit of can be used as the starting point of an efficient preconditioning strategy for and . It is useful to note that the representation of given in Eq. (9) can be decomposed as

(22) |

where corresponds to the non-interacting case, and

(23) |

where the block matrices are given by

(24) |

The inverse can then be formally expanded in powers of , giving

(25) |

where . One can then determine a posteriori whether the truncated series represents a good approximation to , which is a reasonable expectation for . The most computationally intensive part of such an approximation is the application of . However, this matrix is diagonal in frequency-momentum space. It is thus convenient to apply after a four-dimensional Fourier transform of the relevant vector, followed by an anti-Fourier transform to return to the original basis. This approach is referred to as Fourier acceleration. In frequency-momentum space, is given by

(26) |

up to a constant normalization factor, with

(27) |

Care must be taken that the boundary conditions in the temporal direction are antiperiodic, as befit fermions, rather than periodic as commonly assumed in FFT routines. A possible complication arises due to the fact that the radius of convergence of Eq. (25) can be quite limited depending on the spectrum of , which is unknown for all practical purposes. We therefore turn to techniques which may accelerate the convergence of the expansion.

### 3.1 Convergence acceleration methods

Here, we will briefly review the practical aspects of the Euler and van Wijngaarden methods for convergence acceleration of series expansions, without concerning ourselves with the underlying theory. The interested reader is referred to the standard literature, see e.g. Refs. VanWijngaarden (); Lauwerier (); NR (). Given a sequence of partial sums (in this case of sign-alternating terms)

(28) |

where , the Euler method consists of defining a set of new sequences given by

(29) |

where . According to this method, instead of using as an estimate of the infinite sum, one can arrive at a much better estimate as follows. Use the original sequence to define new ones using Eq. (29) until all the terms available have been exhausted. The last sum, namely will consist of only one term, which is the estimate we seek. In practice one usually takes , and it is not uncommon to achieve convergence improvements by several orders of magnitude with a handful of terms as a starting point. It has also been pointed out in the literature that is often an even better approximation than . We shall apply this method to the series that results of applying Eq. (25) to a given vector.

The sequence is commonly referred to as the Euler transform of the original sequence , and it can be shown to be given by

(30) |

where

(31) |

The idea of van Wijngaarden was to first multiply each term in the original series by non-vanishing arbitrary constants and then perform an Euler transformation. It was assumed that a moment generating function exists such that

(32) |

and it was shown that the sum of the original series is given by

(33) |

where

(34) |

In practice, van Wijngaarden’s transformation can turn a slowly convergent or even divergent series into a rapidly convergent series. It can also be shown that the so-called special van Wijngaarden transformation, for which and , where is a free parameter, does not change the Borel sum of the original series and corresponds to the Laplace transform of the Euler transformation Lauwerier ().

## 4 Strong coupling approximation

The idea of constructing a preconditioner based on a strong-coupling approximation was originally developed by Scalettar et al. in Ref. Scalettaretal (). It consists of constructing the inverse of a strong-coupling approximation to , which we shall denote by , where has the same structure as but with the substitution

(35) |

where it should be noted that is diagonal in coordinate space. In this approach, all the blocks will have this property, which makes it extremely inexpensive from a computational point of view. In order to find the inverse, one defines a factorization

(36) |

where

(37) |

and

(38) |

for which it is straightforward to show that the individual blocks are given by

(39) | |||||

(40) |

for and

(41) | |||||

(42) |

where the parameter is tuned such that the number of CG iterations is minimized. The matrix

(43) |

then plays the role of the preconditioner in this approach. Notice that the matrices and are trivial to invert, in the sense that the linear problems and can be solved in a fast and straightforward fashion. For example, is given by the recursive expression

(44) |

While Ref. Scalettaretal () found this approach promising, it presents at least two disadvantages. Firstly, it is not easy to improve in a systematic fashion, in contrast with the Chebyshev approximation. Indeed, one of the key properties of this preconditioner is that it only involves diagonal matrices. Any improvement involving the kinetic energy operator will eliminate this property and make the approach much more computationally expensive. Secondly, just like the Chebyshev method, this preconditioner aims at approximating the inverse of , whereas inverting and separately is preferable due to the dramatic increase in the condition number of relative to . However, it should be pointed out that this preconditioner is computationally very inexpensive and can furthermore be combined with the Chebyshev preconditioner, using the latter as a secondary filter. It is not yet known whether such a combined approach is useful in practice.

## 5 Results

The number of CG or BiCGStab iterations required to solve the linear problems of Eqs. (1) and (2) is given in Table 1 for the various preconditioners, with a more comprehensive study of the scaling of the second-order weak-coupling preconditioner (WCE2) in Fig. 1. The convergence criterion, referred to as the tolerance parameter, was taken to be a reduction in the norm of the residual by a factor of . The achieved absolute accuracy is approximately constant for each space-time volume , but deteriorates slowly as and are increased. However, it should be noted that in typical HMC simulations, a tolerance parameter of is considered adequate. Our test runs have been performed using a imaginary time step of , with couplings of and . The corresponding values of are and . The results shown in Tables 1 and 2 correspond to , which is close to the unitary limit.

We have used auxiliary field configurations with randomly distributed entries in the interval in order to characterize the performance of the preconditioners. While this is the proper interval for the chosen HST, the distribution differs from that of a thermalized HMC simulation, where the values of tend to cluster around and . However, in practice the net effect of using thermalized instead of random configurations is to increase the number of CG iterations by a factor of across all parameter values. We have confirmed this behavior in a number of cases, and we therefore conclude that the relative merits of the preconditioners remain unaffected by the issue of thermalized versus random configurations. We have also found that 10 random field configurations suffices to provide an estimate of the average number of CG iterations with an uncertainty of less than in all cases. The uncertainly can be further suppressed by increasing the number of sample configurations. So far, we have made limited use of Euler–Van Wijngaarden acceleration techniques. For the second-order weak-coupling expansion (WCE2), the simplest Euler method yielded a moderate but definite improvement over the non-accelerated implementation. For the fourth-order expansion (WCE4), the Euler method was found to reduce the number of iterations by up to a factor of 2. Finally, we have varied the right-hand side from a constant vector to one with randomly distributed entries and found no impact on the results.

The number of MV operations consumed by each method is shown in Table 2. Per iteration, CG on the normal equations requires one application of as well as the preconditioner . On the other hand, BiCG and BiCGStab require the application of both and (as well as their respective preconditioners) for each iteration. Thus, in the absence of a preconditioner, the cost per iteration of CG, BiCG and BiGStab is the same. The main advantage of using the normal equations with CG is that one avoids the need to successively solve two linear problems. The MV operations performed per CG iteration when using Chebyshev preconditioning can be estimated as follows: the application of a polynomial of degree to a vector requires MV operations; adding the application of the matrix itself yields a total of operations. Estimating the computational cost of the weak-coupling preconditioner is somewhat more involved. The BiCG or BiCGStab iterations contribute two MV operations from the application of and . For each power of , applying contributes one MV operation to the total cost. However, the non-interacting matrix is applied using a four-dimensional FFT, which scales as , which differs from the scaling of the MV operations. The overall computational cost of the approach using the weak-coupling preconditioner of degree in is given by MV operations, where is the factor in front of the scaling law for the 4D FFT, relative to that of the 3D FFT. While we have not attempted to estimate directly (it is implementation-dependent), we find that in practice the gain in CPU time can be lower by a factor of compared with the gain factors quoted in Table 2, where we set . Notice however that the gain is still substantial, especially as the total number of MV operations increases only mildly as a function of .

Lattice size | CG | BiCGStab | |||||
---|---|---|---|---|---|---|---|

Cheb8 | Cheb16 | Cheb32 | SCA | WCE0 | WCE1 | WCE2 | |

294 | 167 | 104 | 1,790 | 128 | 77 | 55 | |

386 | 226 | 130 | 2,388 | 147 | 84 | 62 | |

489 | 259 | 163 | 2,810 | 152 | 93 | 67 | |

534 | 316 | 187 | 3,470 | 162 | 97 | 71 | |

875 | 460 | 288 | 5,223 | 183 | 114 | 78 | |

1,102 | 610 | 382 | 7,104 | 196 | 118 | 88 | |

864 | 530 | 276 | 5,987 | 208 | 119 | 87 | |

3,301 | 998 | 604 | — | 246 | 148 | 105 | |

— | 1,306 | 735 | — | 279 | 162 | 121 |

While the strong-coupling approximation is extremely inexpensive from a computational point of view, its effect on the solution process is rather mild, at least for the values of considered here. Thus, the number of CG iterations remains rather high, which is an undesirable feature as the iterative solution process is prone to accumulate numerical round-off error after seveal hundred iterations, which is especially problematic for the BiCG and BiCGStab solvers. As a rule of thumb, in our HMC calculations we aim to keep the iterations below the range. Though computationally cheap, a preconditioner which is unable to reduce the iteration count below is of limited practical usefulness, especially at large where the number of iterations may increase dramatically. The effective cost of the strong-coupling preconditioner in terms of MV operations is difficult to estimate, as it involves no FFTs. We have found, empirically, that the number of MV operations required per CG iteration can be conservatively estimated to be . In light of the exceedingly high iteration count shown in Table 1 for this preconditioner, however, this MV estimate is irrelevant.

Lattice size | CG | BiCGStab | ||||||
---|---|---|---|---|---|---|---|---|

Cheb8 | Cheb16 | Cheb32 | SCA | WCE0 | WCE1 | WCE2 | Gain | |

5,292 | 5,678 | 6,864 | 5,370 | 768 | 924 | 990 | 7x | |

6,948 | 7,684 | 8,580 | 7,164 | 882 | 1,008 | 1,116 | 8x | |

8,802 | 8,806 | 10,758 | 8,430 | 912 | 1,116 | 1,206 | 9x | |

9,612 | 10,744 | 12,342 | 10,410 | 972 | 1,164 | 1,278 | 10x | |

15,750 | 15,640 | 19,008 | 15,669 | 1,098 | 1,368 | 1,404 | 14x | |

19,836 | 20,740 | 25,212 | 21,313 | 1,176 | 1,416 | 1,584 | 17x | |

15,552 | 18,020 | 18,216 | 17,961 | 1,248 | 1,428 | 1,566 | 12x | |

59,418 | 33,932 | 39,864 | — | 1,476 | 1,776 | 1,890 | 23x | |

— | 44,404 | 48,510 | — | 1,674 | 1,944 | 2,178 | 27x |

We have also tested the weak-coupling approximation with CG and found that it provides little benefit, in contrast to the dramatic improvement when used with BiCGStab. This situation arises as we have only attempted left preconditioning with CG, and thus the full potential of the preconditioning strategy is not realized. For this to be the case, one would need to precondition with on the left and with on the right. In the case of BiCGStab this complication does not appear, as and are preconditioned separately. Finally, in Table 3 we show the behavior of the weak-coupling preconditioner as a function of the coupling, for and . We note that the case of (), typically failed to converge in less than 15,000 MV operations.

Lattice size | ||||||
---|---|---|---|---|---|---|

WCE2 | WCE4 | WCE2 | WCE4 | WCE2 | WCE4 | |

630 | 600 | 1,206 | 1,320 | 2,358 | 4,110 | |

702 | 690 | 1,584 | 1,530 | 3,816 | 9,330 | |

864 | 780 | 2,178 | 2,190 | 6,444 | 12,720 |

## 6 Summary and Conclusions

In this work, we have evaluated three different strategies to precondition the non-relativistic many-fermion problem: Chebyshev polynomials, a strong-coupling approximation and a weak-coupling expansion. We have argued that such preconditioning forms a central part of HMC calculations, as the frequent solution of a linear problem involving the ill-conditioned fermion matrix is at the heart of the HMC algorithm. We have considered the cases of normal equations which is tractable using the CG algorithm, as well as the two-step approach of solving followed by , using the BiCG or BiCGStab algorithms. Our results indicate that both the Chebyshev polynomials (Sec. 2) and the weak-coupling expansion (Sec. 3) can be effective preconditioners, especially when high orders are used. However, from the point of view of performance, Chebyshev polynomials represent an expensive choice. Additionally, the Chebyshev approach requires a matrix with a positive definite eigenvalue spectrum, which forces us to work with the CG solver on the normal equations. This is unfortunate, as tends to be a rather ill-conditioned matrix (especially at low temperatures), which makes the normal system extremely ill-conditioned. The result is that the number of MV operations required to solve the problem with a preset accuracy grows rapidly with the size of the problem, in particular with , which controls the temperature.

We have shown that solving followed by using the BiCGStab algorithm provides an alternative, extremely promising approach. The use of a weak-coupling preconditioner supplemented with a simple convergence acceleration technique solves the problem elegantly and provides dramatically enhanced performance in terms of both number of iterations and MV operations, which translates into significant savings of CPU time and improved scaling of the HMC algorithm at low temperature. While we have only provided a first, sketchy comparison of the two methods, the difference between using CG with Chebyshev preconditioning and BiCGStab with weak-coupling preconditioning is so striking that the latter technique is the obvious choice. The weak-coupling preconditioner represents a good example in which knowledge of the structure of the matrix, as well as the physical problem at hand, allows for the construction of an effective and efficient strategy to accelerate the solution process.

The strong-coupling approximation of Ref. Scalettaretal () provides an alternative to Chebyshev preconditioning. While computationally inexpensive, this preconditioner turned out to be less efficient in terms of its ability to reduce the number of CG iterations. Apparently, this approach fails to capture the relevant physics at the range of couplings in the vicinity of the unitary limit. In an effort to shed more light on this problem, we explored the behavior of the weak-coupling preconditioner as a function of the coupling , and found that its effectiveness breaks down somewhere between and . This situation might be remedied by using a strong-coupling expansion in powers of , although our attempts to apply such an expansion have so far not been successful. However, we believe this to be a topic worthy of further investigation. Alternatively, one can reduce the value of the imaginary time step , which serves to extend the usefulness of the weak-coupling expansion to larger values of , given that the relevant expansion parameter is . However, this has the drawback of increasing the value of required to achieve a given temperature, which can be taxing on the method.

Various cases of interest that could and should be studied lie beyond the scope of this work. Among these are the case of finite effective range, relevant for nuclear and neutron matter calculations, which would make the potential energy operator more dense in real space (it is diagonal for the zero-range interaction considered here). Other dispersion relations, such as that of the conventional Hubbard model should be studied as well, particularly since this might open the door to large-scale HMC simulations in the fields of solid-state and atomic physics. Finally, we have focused here on a problem in dimensions, although the methods employed in this study carry over directly to applications in lower dimensions.

## Acknowledgements

We would like to thank Dick Furnstahl, Richard Scalettar and Nandini Trivedi for useful discussions and encouragement. We acknowledge support under U.S. DOE Grants No. DE-FG02-00ER41132 and DE-AC02-05CH11231, UNEDF SciDAC Collaboration Grant No. DE-FC02-07ER41457 and NSF Grant No. PHY–0653312. This study was supported in part by the Academy of Finland through its Centers of Excellence Program (2006 - 2011), the Vilho, Yrjö, and Kalle Väisälä Foundation of the Finnish Academy of Science and Letters, and the Waldemar von Frenckell and Magnus Ehrnrooth Foundations of the Finnish Society of Sciences and Letters. Part of this work was performed using an allocation of computing time from the Ohio Supercomputer Center.

## References

- (1) S. Duane, A. D. Kennedy, B. J. Pendleton, D. Roweth, Phys. Lett. B 195, 216 (1987);
- (2) S. A. Gottlieb, W. Liu, D. Toussaint, R. L. Renken, R. L. Sugar, Phys. Rev. D 35, 2531 (1987).
- (3) M. R. Hestenes, E. Stiefel, J. Res. Nat. Bureau of Standards, 49, 409 (1952).
- (4) C. Lanczos, J. Res. Nat. Bureau of Standards, 49, 33 (1952); R. Fletcher,“Conjugate gradient methods for indefinite systems”; G. A. Watson (Ed.), Numerical Analysis, Lecture Notes in Mathematics, Vol. 506, (Springer, Berlin-Heidelberg, 1976).
- (5) H. A. van der Vorst, SIAM (Soc. Ind. Appl. Math.) J. Sci. Stat. Comput. 13(2), 631 (1992).
- (6) C. G. Broyden, M. T. Vespucci, Krylov solvers for linear algebraic systems, (Elsevier, Amsterdam, 2004).
- (7) Y. Saad, Iterative methods for sparse linear systems, (SIAM, Philadelphia, 2003).
- (8) T. DeGrand and C. DeTar, Lattice Methods for Quantum Chromodynamics, (World Scientific: Singapore, 2007); K. Orginos, J. Phys.: Conf. Ser. 46, 132 (2006); T. A. DeGrand, The presence and future of lattice QCD in Computational Methods in Field Theory, edited by H. Gausterer and C.B. Lang, Lecture Notes in Physics 409 (Springer, 1992) 159; U. Wolff, High precision simulations with fast algorithms in Computational Methods in Field Theory, edited by H. Gausterer and C.B. Lang, Lecture Notes in Physics 409 (Springer, 1992) 127.
- (9) Quantum Simulations of Condensed Matter Phenomena, edited by J. D. Doll and J. E. Gubernatis (World Scientific, Singapore, 1990); Quantum Monte Carlo Methods in Condensed Matter Physics, edited by M. Suzuki (World Scientific, Singapore, 1993); W. von der Linden, Phys. Rep. 220, 53 (1992).
- (10) E. Epelbaum, H. Krebs, D. Lee, Ulf-G. Meiner, Eur. Phys. J. A 40, 199 (2009); Phys. Rev. Lett. 104, 142501 (2010).
- (11) R. Blankenbecler, D. J. Scalapino, R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
- (12) T. Abe, R. Seki, Phys. Rev. C 79, 054002 (2009); ibid. 054003 (2009).
- (13) H. F. Trotter, Proc. Am. Math. Soc. 10, 545 (1959); M. Suzuki, Commun. Math. Phys. 51, 183 (1976); Phys. Lett. A 146, 319 (1990).
- (14) R. L. Stratonovich, Sov. Phys. Dokl. 2, 416 (1958); J. Hubbard, Phys. Rev. Lett. 3, 77 (1959).
- (15) D. Lee, Phys. Rev. C 78, 024001 (2008); Prog. Part. Nucl. Phys. 63, 117, (2009).
- (16) S. M. A. Rombouts, K. Heyde, N. Jachowicz, Phys. Rev. Lett. 82, 4155 (1999).
- (17) J. Hirsch, Phys. Rev. B 28, 4059(R) (1983).
- (18) R. T. Scalettar, D. J. Scalapino, R. L. Sugar, D. Toussaint, Phys. Rev. B 36, 8632 (1987).
- (19) B. Bunk, S. Elser, R. Frezzotti, K. Jansen, Comput. Phys. Commun. 118, 95 (1999).
- (20) R. Frezzotti, K. Jansen, Phys. Lett. B 402, 328 (1997); Nucl. Phys. Proc. Suppl. 63, 943 (1998); R. Frezzotti, K. Jansen, Nucl. Phys. B 555, 395 (1999); ibid. 555, 432 (1999).
- (21) A. van Wijngaarden, Indag. Math., 15, 522 (1953).
- (22) H. A. Lauwerier, SIAM (Soc. Ind. Appl. Math.) J. Math. Anal. 6, 96 (1975).
- (23) W. H. Press et al., Numerical Recipes in FORTRAN, (2 Ed., Cambridge University Press, Cambridge, England, 1992).