Stochastic resolution-of-the-identity auxiliary-field quantum Monte Carlo: Scaling reduction without overhead

The accurate ab initio simulation of the ground state properties of molecules and solids is essential for the understanding of major swaths of chemistry and physics. All known techniques for the reliable calculation of electronic structure are challenged by large systems containing strongly interacting electrons. Here, mean field approaches fail, and the exponential scaling of a brute force solution of the Schrödinger equation renders the description of all but the smallest systems impossible. A large array of powerful methods have been developed, which have been successfully employed in the study of large correlated systems.^1–4 Among the many available electronic structure methodologies with which to attack such systems, quantum Monte Carlo (QMC) stands out as a unique tool for the calculation of ground state energies due to its attractive combination of scalability and accuracy.² 

While there are many flavors of QMC,^2,5–7 our focus in this work centers on the auxiliary-field QMC (AFQMC) approach.^8–10 AFQMC is a variant of projector QMC, which is formulated in a second-quantized determinant space. The unbiased version of AFQMC, often referred to as free-projection AFQMC, is exact, in principle, but has an exponential scaling with system size due to noise growth caused by the fermionic sign (or phase) problem.¹⁰ Zhang and co-workers developed the phaseless approximation to AFQMC (ph-AFQMC), which has become a practical and accurate tool for the ab initio simulation of molecules and materials. In ph-AFQMC, the phase of walker wavefunctions during imaginary-time propagation is constrained by the a priori chosen trial wavefunction, |Ψ_T⟩. While this approach is biased, it can yield answers that systematically approach the exact solution either by the release of the constraint¹¹ or via the choice of progressively more sophisticated trial functions.¹² 

In standard AFQMC simulations, walker propagation scales as $O(OM2)$ where O is the number of electrons and M is the number of single-particle basis functions. Such a cubic-scaling propagation per sample makes ph-AFQMC attractive for general ab initio problems. While there have been many successful ph-AFQMC applications to ab initio problems,^13–27 there are several remaining challenges that still need to be addressed before ph-AFQMC can become a universal tool for the study of large scale correlated electronic structure problems. The challenge that we address in this work is the steep cost of the local energy evaluation that is necessary for estimating the ph-AFQMC energy. The local energy of a walker with a wavefunction |ψ⟩ is given by

The cubic-scaling of the walker propagation mentioned above is asymptotically irrelevant in standard ph-AFQMC calculations because the local energy evaluation scales quadratically with system size without exploiting the sparsity or low-rank (LR) structure. This quartic complexity arises from the two-electron repulsion integral (ERI) tensor,

where {ϕ_μ} are the underlying single-particle basis functions. In ph-AFQMC, this integral is most commonly factorized via the Cholesky decomposition (CD) into

where L denotes a Cholesky vector. Even with CD, the evaluation of the local energy remains a quartic-scaling task in ph-AFQMC. For systems with translational symmetry, it is possible to work with plane waves and achieve cubic-scaling with¹⁵ or without²¹ fast Fourier transform. In this work, however, we focus on general basis functions that do not exploit such symmetry.

We mention two notable previous approaches that may be used to address this problem. The first is the tensor hypercontraction decomposition (THC)^28–30 proposed by Martínez, Sherill, and co-workers to factorize the ERI tensor into

Malone et al. successfully applied the THC factorization of the ERI tensor to ph-AFQMC simulations and showed that the local energy evaluation can be brought down to $O(cTHC2OM2)$ where the value of c_THC was found to be ∼8 for accurate local energy evaluation.²⁶ While the improved asymptotic scaling is satisfying, the resulting algorithm has a steep overhead such that the actual crossover between the conventional algorithm and the THC variant occurs for very large system sizes in practice.^26,31

An alternative approach has been proposed by Motta et al.³² In this formulation, the nested diagonalization of Cholesky vectors proposed by Peng and Kowalski³³ is used to exploit the underlying low-rank structure of the ERI tensor. We refer to this factorization as the low-rank (LR) factorization. In the LR factorization, one writes

which arises from

This factorization achieves an asymptotically cubic-scaling local energy evaluation due to the fact that the number of terms in the summation over α in Eq. (6) is limited to $O(logN)$⁠, where N is the system size. Again, due to a large overhead, this asymptotic scaling will generally not be achieved until the system reaches ≳1000 electrons.³² 

The THC and LR factorization schemes enable a cubic-scaling local energy evaluation algorithm, but the overhead associated with both approaches is large enough that the actual efficiency crossover will generally not be observed for medium-sized molecules without obvious sparsity or low-rank structure. We note, however, that these algorithms offer quadratic scaling memory requirements, which results in a reduction of the usual cubic to quartic memory requirement of standard ph-AFQMC. In fact, such a reduction in the memory cost may be the biggest benefit gained from using the THC and LR factorizations. We will return to a discussion of the memory reduction afforded by these techniques before concluding.

Given the above facts, it is clear that there is a need for an overhead-free cubic-scaling algorithm, which can accelerate ab initio ph-AFQMC simulations for medium-sized systems. Based on the earlier work of Baer and Neuhauser,³⁴ Neuhauser and co-workers developed the use of stochastic orbitals as a means of reducing scaling in a variety of electronic structure methodologies, including density functional theory (DFT),^35–37 time-dependent DFT (TDDFT),³⁸ second-order Møller–Plesset perturbation theory (MP2),^39–41 second-order perturbative Green’s functions (GF2),^42–44 the random phase approximation,⁴⁵ GW,⁴⁶ and the Bethe–Salpeter equation.⁴⁶ 

Drawing from this body of work, we explore the combined use of the stochastic resolution-of-the-identity (sRI) formulation of Takeshita et al.⁴¹ with ph-AFQMC. The central idea underlying sRI is a stochastic resolution of the Coulomb operator, which was first developed in Ref. 36. We will demonstrate that sRI can be naturally incorporated into ph-AFQMC and used to reduce the formal scaling of the local energy evaluation to cubic-scaling when combined with CD and to quadratic-scaling when combined with THC or LR factorization. Most notably, we will numerically demonstrate that the combined use of Cholesky decomposition and the stochastic resolution of the identity (CD–sRI) with a simple variance reduction (VR) technique allows for a cubic-scaling algorithm without overhead. Given this attractive feature, we expect the CD–sRI approach to become the standard method for the computation of such quantities for medium-sized systems and THC–sRI and LC–sRI formulations to be the preferred approach for very large ones.

This paper is organized as follows: (1) we give a brief review of the ph-AFQMC algorithm, (2) we analyze the formal scaling of the existing local energy evaluation strategies, (3) we examine scaling reduction in these strategies when combined with sRI, and (4) we present numerical examples (hydrogen chains and water clusters) to show speed-up and scaling reduction of the CD approach with sRI.

To set precise notation for our discussion of scaling, in the following, we use O to denote the number of occupied orbitals, M to denote the number of single-particle basis functions, X to denote the number of Cholesky vectors (or auxiliary basis functions), and N to denote the system size in general.

We briefly review the basic algorithmic details of AFQMC here and refer interested readers to Ref. 47 for a more extended review of modern AFQMC development. We start from an ab initio Hamiltonian Ĥ given by

where Ĥ₁ and Ĥ₂ are the one-body and two-body part of Ĥ, respectively. As with other projector QMC methods, AFQMC starts from

where Ĥ is the Hamiltonian, τ denotes imaginary time, |Ψ₀⟩ is the exact ground state, and |Φ₀⟩ can be any wavefunction with a non-zero overlap with the true ground state. A key feature of AFQMC is that one employs the Hubbard–Stratonovich transformation along with the Trotter approximation to simplify the many-body propagator in Eq. (8) and reduce the problem to that of a series of one-body problems in a fluctuating auxiliary field. With the Trotter approximation for a given time step Δτ, the propagator is expressed as

The Hubbard–Stratonovich transformation allows for the rewriting of the two-body propagator as an integration over auxiliary fields x,

where the one-body operator $v^$ is obtained from

Equation (11) is usually achieved via the CD of the ERI tensor. Alternatively, density-fitting can be used. With this decomposition, the total propagator for given fields x and a given time step Δτ is written as

By the virtue of the Thouless theorem,^48,49 the application of this propagator to a single determinant remains a single determinant. This allows for an efficient walker propagation in AFQMC.

One can reduce statistical fluctuations greatly by employing a mean-field subtraction,⁴⁷ which is closely related to the shifted contour technique developed by Rom and co-workers.⁵⁰ The mean-field subtraction technique redefines the one-body and two-body using the expectation value of $v^P$ with a trial wavefunction,

We then define

which maintains Ĥ = $Ĥ1′$ + $Ĥ2′$⁠. With this, the propagator in Eq. (12) uses $Ĥ1′$ and $v^−⟨v^⟩T$ instead of Ĥ₁ and $v^$⁠. This simple subtraction has been shown to be effective in reducing statistical fluctuations,^47,50 and this is what is used throughout this work.

The imaginary-time equation-of-motion is dictated by the propagator in Eq. (12). With importance sampling via a trial wavefunction |Ψ_T⟩, we write the global AFQMC wavefunction at imaginary-time τ as

where $Nw$ is the number of walkers. Within the phaseless approximation,¹⁰ we update the ith walker weight and determinant via

where the force-bias, $x¯i$⁠, is defined as

The phaseless importance function is defined as

with the hybrid importance function given by

the overlap ratio of the ith walker S_i given by

and the phase θ_i(τ) given by

This completes our compact description of the ph-AFQMC algorithm.

In practical AFQMC calculations, one usually targets the ground state energy estimated via

where i indexes the ith walker, $wi$ is the weight of the ith walker, and E_L[ψ_i] is the local energy of the ith walker, as defined in Eq. (1). Without any approximations, the two-body contribution to the local energy is, via Wick’s theorem, given by

where the Green’s function G is a walker-dependent quantity defined as

Here, C_T is the molecular orbital (MO) coefficient matrix for the trial wavefunction and Θ is defined as

with C being the MO coefficient matrix for the walker determinant |ψ⟩.

Without exploiting any structure in the ERI tensor, the local energy evaluation scales quadratically with system size. For instance, a standard way to evaluate Eq. (1) is to use the so-called “half-rotated” (HR) ERI tensor,^14,26

while storing (ir| js) in memory. This algorithm requires $O(O2M2)$ memory and scales as $O(O2M2)$⁠. Each walker additionally requires $O(O2M)$ memory to save intermediates that appear when evaluating Eq. (28). Furthermore, the formation of (ir| js) at the beginning of the QMC run scales as $O(OM4)$⁠, which is more expensive than the local energy evaluation, although it needs to be performed only once. Despite these relatively steep scaling behaviors, HR is perhaps the most widely used local energy algorithm when the memory cost is affordable. We refer this algorithm to as the “half-rotated” (HR) local energy evaluation.

Another standard evaluation procedure is to work with density-fitted integrals or Cholesky vectors directly. This may be achieved by using Eq. (3) within the local energy expression, yielding

which scales as $O(O2MX)$⁠. In general, X is 4–5 times larger than M, and therefore, Eq. (29) is slower than Eq. (28). However, Eq. (29) requires only $O(OMX)$ memory storage, which makes it better suited for large-scale simulations. To save intermediates, each walker additionally requires $O(O2X)$ memory. We refer this algorithm to as the “Cholesky decomposition” (CD) local energy evaluation.

Finally, we note that some recent advances in integral factorization techniques have allowed for a cubic-scaling local energy evaluation with a quadratic memory storage.²⁶ We first discuss the THC strategy^28–30 invented by Sherill, Martínez and, co-workers, which was applied to AFQMC by Malone et al.²⁶ Specifically, we use Eq. (4) within the local energy evaluation, yielding

All of the tensor contractions in Eq. (30) are simple matrix–matrix multiplications, leading to cubic scaling with system size. More precisely, the local energy calculation scales as $O(cTHC2OM2)$⁠, where we define the THC rank to be c_THCM. The memory requirement is set by storing $ηpP^$ and M that have only a quadratic number of non-zero values. Furthermore, each walker additionally requires $O(cTHC2M2)$ memory to save intermediates that naturally appear in Eq. (30). While these asymptotic properties are highly attractive, the THC algorithm has been shown to have a large prefactor compared to Eqs. (28) and (29) because c_THC can be quite large in practice.²⁶ For example, c_THC was found to be ∼8 for the accurate calculation of the local energy of diamond within the double-zeta basis.

Similarly, for the LR factorization of Motta et al.,³² we use Eq. (6) with the AFQMC local energy expression, yielding

which scales as $O(nrOMX)$ where the rank n_r sets the upper limit on the summation over α and β and scales like log(N). Note that this approach also effectively achieves a cubic-scaling local energy evaluation algorithm. The memory requirement is $O(nrOX+nrMX)$⁠, and each walker additionally requires $O(nrOX)$ to save intermediates. However, as with the use of THC, n_r is usually large, and therefore, for practical applications, useful cubic-scaling may not be observed,³² and only memory saving may be practically useful.

We follow the explication of the stochastic resolution-of-the-identity (sRI) technique proposed by Takeshita et al.,⁴¹ which has been shown to lower the scaling of MP2 (RI-MP2)⁴¹ and the second order Green’s function method (RI-GF2).^43,44 sRI is based on the simple mathematical observation that one may represent the Kronecker delta function with stochastic functions as

where N_ξ is the number of stochastic samples and θ^ξ is an sRI basis function whose entry is randomly chosen to be ±1. In practice, we will limit the number of samples to a finite number, which does not scale with system size. This restriction has been shown to be sufficient to achieve a fixed statistical error per particle in MP2 and GF2. The key feature of this approach is that it does not assume any structure (either sparsity or low-rank) in the underlying ERI tensor while still reducing the cost. Due to this fact, the overhead of this approach is almost negligible, especially when looking at size-intensive quantities for which ph-AFQMC applications are also well-suited.

We emphasize that the scaling of a QMC algorithm should not be discussed without the assessment of the underlying statistical error. For instance, one may argue that propagation within ph-AFQMC scales cubically, but the number of samples associated with a fixed statistical error scales linearly with system size. Then, this makes the ph-AFQMC propagation scale quadratically with system size for a fixed statistical error. The usual cubic scaling of the ph-AFQMC propagation quoted in the literature²¹ is the scaling for a fixed statistical error per particle. Similarly, the lower scaling of sRI-MP2 and sRI-GF2 approaches cited above is only justified when considering observables for a fixed statistical error per particle. This makes the combination of these two methods very natural. We also note that it is still possible to obtain the total energies for a fixed statistical error in both the ph-AFQMC and sRI methods, albeit with an increased cost.

We first re-write Eq. (28) with two Kronecker deltas,

Next, one may employ two sets of stochastic orbitals to represent these Kronecker deltas and write the local energy as

where

and

The formation of Eq. (35) is the bottleneck, scaling as $O(NξO2M2)$⁠. The additional memory usage of each walker scales as $O(NξO2M)$⁠. Thus, there is no reason to employ this algorithm as the asymptotic scaling is not improved by sRI and no memory saving is obtained.

We next insert Eq. (32) into Eq. (29) and show that the scaling of the local energy evaluation is reduced from quartic to cubic. The combination of the sRI expression to the CD-based local energy algorithm is most natural in the auxiliary basis function space. Specifically,

We further define an intermediate tensor R,

The formation of R scales as $O(NξOMX)$⁠, which is cubic-scaling for a size-intensive quantity N_ξ. The Coulomb term in Eq. (29) is thus expressed as

which scales as $O(OMX)$⁠. The extra memory cost for storing walker-specific intermediates scales as $O(X)$⁠. Using R,

where

The summation over P is the term now replaced by the summation over ξ, which lowers the scaling from $O(OMX)$ to $O(OMNξ)$⁠. The walker-specific extra memory cost is reduced from $O(X)$ to $O(Nξ)$⁠. The exchange term in Eq. (29) is

and scales as $O(O2MX)$ with $O(O2X)$ walker-specific extra memory cost. This term now becomes

where

which scales as $O(O2MNξ)$⁠. Its walker-specific memory cost scales as $O(NξO2)$⁠. This completes the demonstration of a scaling reduction of Eq. (29) to cubic scaling. In summary, the asymptotic scaling of the CD–sRI algorithm is $O(NξOMX+NξO2M)$⁠, which includes the cost of the formation of R as well. The additional memory cost due to storing R and other intermediates scales as $O(NξOM+NξO2)$⁠. This is also an improvement over the conventional CD algorithm.

We next apply sRI to Eq. (30) and show that the overall scaling can be reduced to quadratic. With the use of two sets of sRI insertions, we write

with

where we define

With a proper set of contractions, one can show that this expression scales as $O(NξcTHCM2+NξOM+NξcTHC2M2)$⁠. The additional walker memory usage due to storing intermediates comes at a cost of $O(NξcTHCM+Nξ2cTHCM)$⁠. THC–sRI improves both scaling and memory usage over the conventional THC algorithm.

The same conclusion can be deduced for the LR–sRI algorithm by employing two sets of sRI expressions. With sRI, one may write Eq. (31) as

with

Here, we have defined

With an appropriate series of contractions, it can be shown that this local energy evaluation scales as $O(NξnrMX+NξOM+NξnrOX)$⁠. The extra memory required for storing walker-specific intermediates is $O(NξnrX)$⁠. As with THC–sRI, we observe an improvement using LR–sRI in both scaling and memory usage over the LR algorithm.

The implementation of the sRI local energy estimator can be achieved with only minor modifications to the existing AFQMC programs. There are two stochastic samplings to complete the evaluation of the sRI global energy: one is the standard AFQMC walker local energy sampling for each walker in Eq. (24) and the other is the sRI stochastic orbital sampling in Eq. (32). A simple algorithm for this is as follows:

1. Each walker samples a set (or two sets depending on the choice of algorithms) of N_ξ sRI orbitals θ.

2. Evaluate the corresponding local energy expression derived in Sec. II C.

3. The global mixed estimator energy is then estimated in the standard manner as in Eq. (24) as before.

Minimizing N_ξ is critical in order to achieve the aforementioned lower scaling with as little overhead as possible. Fortunately, within the QMC setup, it is natural to set N_ξ = 1 per sample because the representation power of sampled sRI orbitals increases as stochastic samples are accumulated throughout the imaginary-time propagation. In other words, at each time step, each walker samples one sRI orbital, and the global estimator is obtained by stochastically averaging over many such samples. In the limit of infinite statistical sampling, this algorithm converges to the exact ph-AFQMC energy (i.e., the one without sRI).

The use of a variance reduction (VR) technique^51,52 for sRI-ph-AFQMC has been found to be very effective. For any of the algorithms outlined in this work, we write

where ψ_T is the trial wavefunction. As suggested by this equation, we compute the “correlation” contribution to the local energy via sRI and the mean-field energy is reused. By expressing the energy in this form, the statistical fluctuations are found to be greatly reduced. This technique is often referred to as the “control variate” approach in Monte Carlo.^51,52 We note that the efficacy of the control variate approach in this work will ultimately depend on the quality of the trial wavefunction as well. As long as the trial wavefunction is of good quality, the variance reduction will be effective.

We have discussed a total of four local energy evaluation strategies and how the sRI can reduce the scaling of three of them. We summarize this in Table I. With the sRI, we have formulated one cubic-scaling (CD–sRI) and two quadratic-scaling (THC–sRI and LR–sRI) algorithms.

For each algorithm, we have also discussed the additional memory requirement of each walker for storing intermediates. This is summarized in Table II. With sRI, we have one algorithm (CD–sRI) that requires additional quadratic-scaling memory and two algorithms that require additional linear-scaling (THC–sRI and LR–sRI) memory. These are all improvements over their conventional counterparts.

Among these algorithms, we choose CD–sRI to demonstrate the behavior of sRI-ph-AFQMC methods due to its simple implementation, reasonable leading speedup of O/N_ξ, and leading walker-specific memory saving of $O(OX/(NξM))$⁠. The implementation of THC–sRI and LR–sRI is straightforward if THC and LR are already implemented. However, in PAUXY,⁵³ these are currently unavailable, which led us to choosing CD–sRI for the purpose of demonstration. Moreover, we apply the sRI only to the exchange contribution (i.e., E_K) because the Coulomb contribution (i.e., E_J) can be evaluated at cubic cost within the standard algorithm. More specifically, we write (with the VR technique)

where E_K,CD[ψ_T] is evaluated only once in the beginning and used throughout the simulation. Since E_K,CD–sRI is evaluated twice (once for each of ψ and ψ_T), this adds an additional prefactor of 2. This prefactor increase is negligible compared to the variance reduction that we gain, as we shall see. The resulting algorithm is overall cubic-scaling per sample.

All calculations were performed with an open-source python-based AFQMC program called PAUXY.⁵³ The pair-branching algorithm was employed for the population control.⁵⁴ The time step of 0.005 a.u. was used in all computations. We set N_ξ = 1 throughout all sRI calculations. Increasing N_ξ provided almost no improvement in statistical efficiency for systems considered in this work. This enables a speedup of the order of O (the number of electrons) compared to the conventional CD algorithm. We use the cc-pVDZ⁵⁵ basis set in all examples presented below. We used restricted Hartree–Fock as a trial wavefunction throughout because systems considered in this work are dominated mainly by dynamical correlation. The truncation threshold for the Cholesky decomposition was set to 10⁻⁵. A total of 160 walkers were used in every calculation presented below. QMCPACK^56,57 was used to cross-check the total energies reported in this work along with the underlying population bias associated with the pair-branching algorithm. The reported total energies have negligible population bias and time step error.

For numerical examples, we picked a series of one-dimensional (1D) hydrogen chains (H-chains) and water clusters, which clearly illustrates the following:

1. The VR technique discussed in Sec. II D is very effective. The VR was found to be so effective with restricted Hartree–Fock trial wavefunctions for systems considered here, which allowed for a very small number of stochastic orbitals to sample (i.e., N_ξ = 1).

2. The overhead of using CD–sRI with VR for computing ph-AFQMC total energies as opposed to using CD is negligible, unlike for the THC and LR approaches.

3. CD–sRI is a cubic-scaling local energy algorithm if implemented correctly (i.e., with a proper tensor contraction ordering).

(1) and (2) are rather difficult to show analytically, and therefore, we provide numerical examples to support our claims. While (3) was formally shown in Sec. II C 2, we provide timing benchmarks to bolster support for this analysis. For 1D H-chains, the inter-hydrogen distance was fixed to be 1.6 bohr, and the water cluster geometries were taken from Ref. 31.

We first focus on the effect of the VR technique discussed in Eq. (52). In Fig. 1, we present the energy per H atom as a function of imaginary time. In both H₁₀ and H₄₀, we observe the same trend. Without VR, CD–sRI fluctuates significantly around the CD results, while CD–sRI with VR behaves statistically very similarly to CD. These results indicate that the overhead of CD–sRI with VR is almost negligible, which becomes more evident when looking at statistically averaged total energies.

In Table III, we compare CD and CD–sRI with VR for the total energies of H-chains. The difference in total energies between two algorithms is statistically insignificant since they are all within error bars. This is not surprising since CD–sRI is not expected to introduce any biases into the estimator and should recover the CD result with enough samples. Furthermore, the magnitude of statistical error is similar when comparing the two algorithms for a similar number of statistical samples.

We have also carried out similar numerical experiments for water clusters, (H₂O)_N, with N = 2, 4, 8. A representative imaginary-time propagation trajectory is given in Fig. 2. The conclusions drawn from the H-chain results hold for these cases as well. With VR, the statistical fluctuations of CD–sRI are close enough to CD. This is true for all values of N = 2, 4, 8, and we expect this to hold for larger clusters. Finally, we present the absolute energies of (H₂O)_N computed via CD and CD–sRI in Table IV. The total energies from two different algorithms lie within their respective error bars. Similarly to Table III, two methods exhibit comparable magnitude of error bars, given a similar number of statistical samples. This strongly suggests that CD–sRI with VR performs as well as CD without overhead.

We finish our discussion with timing benchmarks of walker memory saving of H-chains and water clusters. For a consistent increase in the number of Cholesky vectors, we employ the density-fitting approximation instead. The auxiliary basis set used here is that of Weigend et al.⁵⁸ The timing benchmark here was obtained from an Intel(R) Xeon(R) Platinum 8268 central processing unit (CPU) 2.90 GHz processor, and only a single thread was used. Furthermore, a single local energy evaluation of a single walker is considered for the purpose of demonstration. The timing results are reported in Fig. 3. Panel (a) shows the results for 1D H-chains up to N = 300. The observed scaling for CD is $O(N3.13)$⁠, whereas it is $O(N2.12)$ for CD–sRI (estimated with R² greater than 0.99). For the obtained data points, there is no crossover between CD and CD–sRI. CD–sRI is always faster than CD, and there is no overhead for using CD–sRI compared to CD for a single sample. Similarly, panel (b) shows the timing benchmark for water clusters up to N = 64. The observed scaling for CD is $O(N3.63)$⁠, whereas it is $O(N2.94)$ for CD–sRI (estimated with R² greater than 0.999). The CD algorithm does not exhibit asymptotic scaling for these system sizes, which could be due to the efficient use of Basic Linear Algebra Subroutines (BLAS). CD–sRI appears to be quite close to its theoretically predicted asymptotic scaling. The memory results are obtained based on Table II and are shown in Fig. 4. Even for the smallest systems examined, CD–sRI requires an order of magnitude less memory for each walker. Due to the scaling differences between CD and CD–sRI, the difference in the memory usage becomes only larger as the system size increases. This highlights the utility of CD–sRI in reducing the additional memory usage for each walker.

In summary, based on simple benchmarks on 1D H-chains and water clusters, we have successfully demonstrated the utility of the CD–sRI algorithm. With the variance reduction technique introduced here, the CD–sRI algorithm achieves a cubic-scaling local energy evaluation without overhead.

In this work, we have shown how the stochastic resolution-of-the-identity (sRI) technique developed by Takeshita et al. can be seamlessly integrated with the local energy evaluation of phaseless auxiliary-field quantum Monte Carlo (ph-AFQMC). We considered four different existing local energy evaluation strategies: a half-rotated (HR) electron repulsion integral tensor approach, a Cholesky decomposition-based approach (CD), tensor hypercontraction (THC), and low-rank (LR) factorization approaches. We have carefully analyzed their formal scalings and discussed possible scaling reduction as well as walker-specific memory reduction when combined with sRI.

We found that HR–sRI achieves neither scaling reduction nor walker-specific memory reduction. It scales quadratically with system size like the original HR approach with the same memory usage. On the other hand, the CD approach, which formally scales quadratically with system size, can be reduced to cubic-scaling when combined with sRI (CD–sRI). Furthermore, the walker-specific memory usage is reduced to quadratic from cubic. Similarly, the cubic-scaling approaches, THC and LR, can be reduced to quadratic-scaling with sRI. The additional walker-specific memory requirement is also reduced from quadratic to linear. Without sRI, previously available algorithms using THC and LR achieve cubic scaling, but the overhead associated with them limits the applicability of these algorithms. Therefore, we have examined a new cubic-scaling algorithm, CD–sRI, with a particular focus on characterizing the overhead associated with the algorithm.

We applied the CD–sRI algorithm to one-dimensional hydrogen chains (1D H-chains) and water clusters. With a variance reduction (VR) technique developed for CD–sRI, CD–sRI with VR exhibits no overhead compared to CD up to H₈₀ and up to (H₂O)₈. By no overhead, we mean that CD–sRI and CD exhibit a similarly large error bar for a similar number of statistical samples. Therefore, CD–sRI adds no additional cost for computing the total ph-AFQMC energy for a fixed statistical error. Furthermore, we have performed a timing benchmark for H-chains up to H₃₀₀ and water clusters up to (H₂O)₆₄. We observed reduced scaling from CD–sRI compared to CD. Specifically, the observed scaling of CD was $O(N3.13)$ for H-chains and $O(N3.63)$ for water clusters, while we observed $O(N2.12)$ for H-chains and $O(N2.94)$ for water clusters in the case of CD–sRI. Based on these observations, we argue that the CD–sRI algorithm should be the standard cubic-scaling local energy algorithm with an expected generic speed-up compared to the conventional CD algorithm.

The remaining challenges for large-scale ph-AFQMC simulations are the memory bottleneck and the cubic-scaling walker propagation that occurs every time step. Both THC and LR approaches have successfully addressed the memory bottleneck^26,32 as they require only quadratic-scaling storage for each walker. More theoretical development is needed for accelerating the walker propagation. With THC–sRI or LR–sRI, the local energy evaluation scales quadratically with system size. In these cases, the bottleneck for ph-AFQMC is the walker propagation. It has been shown that a significant speed-up of the propagation is possible by utilizing graphical processing units (GPUs).^59,60 Furthermore, the significant reduction in the requirements for walker-specific extra memory will be particularly useful with GPUs since in this case, the available memory is highly limited. With sRI, the propagation is the bottleneck in AFQMC for all system sizes. Therefore, it will be interesting to explore ideas such as sRI to design new algorithms for walker propagation, which scale quadratically (or less) with system size. We reserve this and the exploration of additional projects centered on THC–sRI or LR–sRI with GPUs for large-scale ph-AFQMC applications for future work.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

We are indebted to Fionn Malone for stimulating discussions and pointing out the potential usefulness of control variate in reducing the variance. We also thank Fionn Malone for help with QMCPACK,^56,57 which was used to verify the lack of population bias of reported numbers in this work. We thank Evgeny Epifanovsky for useful discussions on efficient implementations of the local energy evaluation. We thank Miguel Morales and Shiwei Zhang for helpful comments. D.R.R. acknowledges support from Grant No. NSF-CHE 1954791.