Faster and lower scaling orbital-space Variational Monte Carlo
Iliya Sabzevari Sandeep Sharma





Department of Chemistry and Biochemistry, University of Colorado Boulder, Boulder, CO 80302, USA

sanshar@gmail.com
Department of Chemistry and Biochemistry, University of Colorado Boulder, Boulder, CO 80302, USA

In this work, we introduce three algorithmic improvements to reduce the cost and improve the scaling of orbital space variational Monte Carlo (VMC). First, we show that by appropriately screening the one- and two-electron integrals of the Hamiltonian one can improve the efficiency of the algorithm by several orders of magnitude. This improved efficiency comes with the added benefit that the resulting algorithm scales as the second power of the system size O(N^2), down from the fourth power O(N^4). Using numerical results, we demonstrate that the practical scaling obtained is in fact O(N^1.5) for a chain of Hydrogen atoms, and O(N^1.2) for the Hubbard model. Second, we introduce the use of the rejection-free continuous time Monte Carlo (CTMC) to sample the determinants. CTMC is usually prohibitively expensive because of the need to calculate a large number of intermediates. Here, we take advantage of the fact that these intermediates are already calculated during the evaluation of the local energy and consequently, just by storing them one can use the CTCM algorithm with virtually no overhead. Third, we show that by using the adaptive stochastic gradient descent algorithm called AMSGrad one can optimize the wavefunction energies robustly and efficiently. The combination of these three improvements allows us to calculate the ground state energy of a chain of 160 hydrogen atoms using a wavefunction containing ∼ 2× 10^5 variational parameters with an accuracy of 1 mE_h/particle at a cost of just 25 CPU hours, which when split over 2 nodes of 24 processors each amounts to only about half hour of wall time. This low cost coupled with embarrassing parallelizability of the VMC algorithm and great freedom in the forms of usable wavefunctions, represents a highly effective method for calculating the electronic structure of model and ab initio systems.

Improved speed and scaling in orbital space variational monte carlo
    Sandeep Sharma
    December 30, 2023
===================================================================



Quantum Monte Carlo (QMC) is one of the most powerful and versatile tools for solving the electronic structure problem and has been successfully used in a wide range of problems<cit.>.

QMC can be broadly classified into two categories of algorithms, variational Monte Carlo (VMC)<cit.> and projector Monte Carlo (PMC). In VMC one is interested in minimizing the energy of a suitably chosen wavefunction ansatz. The accuracy of VMC is limited by the functional form and the flexibility of the wavefunction employed. PMC on the other hand is potentially exact and several variants exist, such as diffusion Monte Carlo (DMC)<cit.>, Auxiliary field quantum Monte Carlo (AFQMC)<cit.>
, Green's function Monte Carlo (GFMC)<cit.> and full configuration interaction quantum Monte Carlo (FCIQMC)<cit.>. Although exact in principle, in practice all versions of PMC suffer from the fermionic sign problem when used for electronic problems, except in a few special cases. A commonly used technique for overcoming the fermionic sign problem is to employ the fixed-node or fixed-phase approximation that stabilizes the Monte Carlo simulation by eliminating the sign problem at a cost of introducing a bias in the calculated energies. The fixed-node or the fixed-phase is often obtained from a VMC calculation and thus the accuracy of the PMC calculation depends to a large extent on the quality of the VMC wavefunction. Thus VMC wavefunction plays a pivotal role and to a large extent determines the final accuracy obtained from a QMC calculation. (FCIQMC is unique in its ability to control the sign problem self-consistently by making use of the initiator approximation.) Traditionally the most commonly used variant of QMC has been the real space VMC followed by the fixed-node DMC calculation which is able to deliver results at the complete basis set limit. Although extremely powerful, one of the shortcomings of real space QMC methods is that there is less error cancellation than is typically observed in basis set quantum chemistry. This fact, coupled with the recent success of both AFQMC and FCIQMC (both of which work in a finite basis set), has led to renewed interest in development of new QMC algorithms that work in a finite basis set<cit.>. The present work is an attempt in this direction and presents algorithmic improvements for performing orbital space VMC calculations. 

To put the orbital space VMC in the broader context of wavefunction methods, it is useful to classify the various wavefunctions used in electronic structure theory into three classes. The wavefunctions in the first class allow one to calculate the expectation value of the energy, ⟨ψ|H|ψ⟩, deterministically with a polynomial cost in the number of parameters, examples of these include the configuration interaction wavefunction<cit.> and matrix product states<cit.>. The second class of wavefunctions only allow one to calculate the overlap of the wavefunction with a determinant (or a real space configuration), ⟨|ψ⟩, deterministically with a polynomial cost in the number of parameters, examples include the Slater-Jastrow (SJ)<cit.>, correlator product states<cit.>, Jastrow-antisymmetric geminal power (JAGP)<cit.>. Finally, the third class of wavefunction are the most general and do not allow polynomial cost evaluation of either the expectation value of energy or an overlap with a determinant, examples of these include the coupled cluster wavefunction<cit.> and the projected entangled pair states<cit.>. The first class of wavefunction are the most restrictive, but are the easiest to work with and often efficient deterministic algorithms are available to variationally minimize their energy. For the second class of wavefunctions one has to resort to VMC algorithms to both evaluate the energy and optimize the wavefunctions. Finally, the third type of wavefunctions are the most general and are often quite accurate but one needs to resort to approximations to evaluate their energies. In this work we will focus on the second class of wavefunction, in particular the wavefunction comprising of a product of the correlator product states and a Slater determinant (CPS-Slater). 

The rest of the paper is organized as follows. In Section <ref>, we will briefly outline the main steps of the VMC algorithm, followed by the new algorithmic innovation introduced to improve the efficiency and reduce the scaling of the algorithm. Next, in Section <ref>, we will describe how the algorithm is used to calculate the energy and optimize the CPS-Slater wavefunction. We will also discuss in detail the cost and computational scaling of the algorithm
. Finally in Section <ref> we will present the results obtained by this algorithm for model systems including the 1-D hydrogen chain and 2-D Hubbard model. 



§ ALGORITHMIC IMPROVEMENTS

In VMC, the expectation value of the Hamiltonian for a wavefunction Ψ(), where  are the set of wavefunction parameters, is calculated by performing a Monte Carlo integration.

    E =   ⟨Ψ() |H|Ψ()⟩/⟨Ψ() |Ψ()⟩ = ∑_ |⟨Ψ() |⟩|^2 ⟨|H|Ψ()⟩/⟨|Ψ()⟩/∑_ |⟨Ψ()|⟩|^2
    
     =   ∑_ρ_ E_L[] 
    
    =   ⟨ E_L[]⟩_ρ_

where, E_L[] = ⟨|H|Ψ()⟩/⟨|Ψ()⟩ is the local energy of a determinant
, the last expression denotes that the energy is calculated by averaging the local energy over a set of determinants 
that are sampled from the probability distribution ρ_ =  |⟨Ψ() |⟩|^2 /∑_ |⟨Ψ()|⟩|^2. 
With Equation <ref> as the background, we can describe the VMC algorithm by splitting it into three main tasks


  * For a given determinant |⟩ we have to calculate the local energy E_L[].

  * We have to generate a set of determinants |⟩ with a probability distribution ρ_ for a given wavefunction Ψ.

  * And finally, we need an optimization algorithm that can minimize the energy of the wavefunction by varying the parameters . 

We introduce algorithmic improvements in each of these tasks which are described in detail in the next three sections.


 
 

 §.§ Reduced scaling evaluation of the local energy
 The local energy E_L[] of a determinant |⟩ is calculated as follows
 
    E_L[] =   ⟨|H|Ψ()⟩/⟨|Ψ()⟩
    
    =   ∑_⟨|H|⟩⟨|Ψ()⟩/⟨|Ψ()⟩
    
    =   ∑_ H_,⟨|Ψ()⟩/⟨|Ψ()⟩,

 where the sum is over all determinants  that have a nonzero transition matrix element (H_, = ⟨|H|⟩) with . The number of such non-zero matrix elements H_, are on the order of n_e^2 n_o^2, where n_e is the number of electrons and n_o is the number of open orbitals. This number increases as a fourth power of the system size and a naive use of this formula results in an algorithm that scales poorly with the size of the problem. 
 
 To reduce the cost of calculating the local energy we truncate the summation over all  to just a summation over those , that have a Hamiltonian transition matrix element H_, > ϵ, 
 
    E_L[, ϵ] =   ∑_^ϵ H_,⟨|Ψ()⟩/⟨|Ψ()⟩,

 where ϵ is a user defined parameter. Note that in the limit that ϵ→ 0, we recover the exact local energy, E_L[, ϵ]  → E_L[]. It is useful to note that when a local basis set is used the number of elements H_, that have a magnitude larger than a fixed non-zero ϵ scale quadratically with the size of the system. To see this, let's consider double excitations, the argument for the single excitations is similar. For a given determinant (), we can obtain another determinant by a double excitation of electrons from say orbitals i and j to orbitals a and b to give another determinant () with a matrix element 
 |H_,| = |⟨ a b|i j⟩ - ⟨ a b | j i⟩|. Because of the locality of orbitals, all orbitals decay exponentially fast and thus the integral ⟨ a b|i j⟩ is negligible unless a is close to i and b is close to j. In other words this integral is non-zero for only O(1) instances of a and b and because there are an O(N^2) number of possible i and j, the total number of non-negligible integrals is also O(N^2). Thus if we are able to efficiently screen the transition matrix elements for a given ϵ≠ 0, no matter how small the ϵ is, we are guaranteed to obtain a quadratically scaling evaluation of the local energy. 
 
 To see how the parameter ϵ can be effectively used to dramatically reduces the cost of the local energy evaluation, let us examine the one electron and the two electron excitations separately. 
 
 First let's see how the more numerous two electron excitations can be screened. The value of the Hamiltonian transition matrix element between a determinant, , contributing to the summation, that is related to determinant, , by the double excitation of electrons from orbitals i and j to orbitals a and b is equal to 
 |H_,| = |⟨ a b|i j⟩ - ⟨ a b | j i⟩|. The magnitude of the Hamiltonian matrix element only depends on the four orbitals that change their occupation and does not depend on the rest of the orbitals of the determinants  or . To use this fact efficiently, we stores the integrals in the heat bath format whereby, for all pairs of orbitals i and j we store the tuples {a,b,⟨ a b|i j⟩ - ⟨ a b | j i⟩} in descending order by the value of |⟨ a b|i j⟩ - ⟨ a b | j i⟩|. Now for generating all determinants  that are connected to  by double excitation that is larger than ϵ, one simply loops over all pairs of occupied orbitals i and j in determinant . For each of these pairs of orbitals one loops over the tuple {a,b,⟨ a b|i j⟩ - ⟨ a b | j i⟩}, until the value of |⟨ a b|i j⟩ - ⟨ a b | j i⟩| > ϵ and one exits this inner loop as soon as this inequality is violated. The CPU cost of storing the integrals in the heat bath format is O(N^4ln(N^2)), but this is essentially the same cost as reading the integrals O(N^4) and only needs to be done once for the entire calculation. The algorithm described here is essentially exactly the same as the one used to perform efficient screening of two electron integrals in the Heat Bath Configuration interaction algorithm<cit.>, which is a variant of the selected configuration interaction algorithm<cit.> and is much more efficient than other variants.
 
It is important to recall that although it might seem that there are only O(N^2) single excitations, and thus one does not need to screen them, this is in fact not true. The reason being the Hamiltonian transition matrix element between two determinants,  and , that are connected by the single excitation of an electron from orbital i to orbital a is given by H_, = ⟨ i|a⟩ + ∑_j∈ occ⟨ ij|aj⟩ - ⟨ ij|ja⟩. There are O(N^2) such connections and the cost of calculating each of these matrix elements scales as O(N) thus making the entire cost O(N^3). We have observed that without appropriate screening of these one electron integrals their cost starts to dominate the cost of calculating the local energy.  To screen the singly excited determinants we first calculate the quantity S_i,a = |⟨ i|a⟩| + ∑_j |⟨ ij|aj⟩ - ⟨ ij|ja⟩|, for all pairs of orbitals i and a. Notice, that S_i,a > |H_,|. Thus, while generating singly excited determinants from a determinant , by exciting an electron from an occupied orbital i to an empty orbital a, we discard all excitations where S_i,a <ϵ.
 
 Thus for a given ϵ, only determinants, , that make a non-zero contribution to the local energy are ever considered. As we will demonstrate in Table <ref>, this algorithm can be used to discard a very large fraction of the determinants from the summation of Equation <ref> providing orders of magnitude improvement in the efficiency of the algorithm while incurring only a small error in the overall energy. In Section <ref> we will show that the contribution to the local energy (H_,⟨|Ψ()⟩/⟨|Ψ()⟩) for each of these O(N^2) determinants  can be evaluated at a cost O(1), allowing one to evaluate the local energy with a cost that scales quadratically with the size of the system. Recently, two other algorithms have appeared in literature that sample the matrix elements H_, stochastically<cit.> and semistochastically<cit.> instead of screening the summation deterministically as we have proposed here. Both these algorithms should have the same asymptotic scaling as the current algorithm, however, in practice the benefits of reduced scaling only manifest themselves for large system sizes.
 
 


 §.§ Continuous time Monte Carlo for sampling determinants

The usual procedure for generating determinants  according to a probability distribution ρ_ uses the Metropolis-Hastings algorithm in which a Markov chain is realized by performing a random walk according to the following steps 


  * Starting from a determinant , a new candidate determinant |⟩ is generate with a proposal probability P(←).

  * The candidate determinant is then accepted if a randomly generated number between 0 and 1 is less than min(1, ρ() P(←)/ρ()P(←)), otherwise we stay at determinant  for another iteration.		

The algorithm guarantees that with sufficiently large number of steps one generates determinants according to the probability distribution ρ as long as the principle of ergodicity is satisfied. This states that it should be possible to reach any determinant  from a determinant  in a finite number of steps. The draw back of the algorithm is that if the one chooses the proposal probability distribution poorly, such that determinants  that have a small ρ() are proposed with a high probability, then several moves will be rejected leading to long autocorrelation times. Although the algorithm places very few restriction on the proposal probability distribution, and one can in principle devise a probability distribution that leads to very few rejections, in practice it is far from easy to do this. 

In this work we bypass the need for devising complicated proposal probability distributions, by using the continuous time Monte Carlo (CTMC), also known as the kinetic Monte Carlo (KMC), Bortz-Kalos-Lebowitz (BKL) algorithm<cit.> and Gillespie's algorithm<cit.>. This is an alternative to the Metropolis-Hastings algorithm and has the advantage of being a rejection free algorithm and every proposed move is accepted. Just as in the case of the Metropolis-Hastings algorithm, there is considerble flexibility in how the algorithm is implemented but in this work we use the following steps:


  * Starting from a determinant  calculate the quantity 

    r(←) = (ρ()/ρ())^1/2 = |⟨|Ψ()⟩/⟨|Ψ()⟩|

 for all determinants  that are connected to  by a single excitation or a double excitation. 

  * Calculate the residence time t_ = 1/∑_ r(←) and stay on the walker  for the time t_ (the residence time can also be viewed as a weight).

  * Next, a new determinant is selected, without rejection, out of all the determinant  with a probability proportional to r(←).		


To get a better understanding of the algorithm, it is useful to think of each determinant () as a reactive species in a unimolecular reaction network that can be transformed into another determinant () through a unimolecular reaction with the ratio of the forward and reverse rates determined by the equilibrium constant (in our case k(←)/k(←) = K_eq =  | ⟨|Ψ()⟩/⟨|Ψ()⟩ |^2). As long as all determinants are participating in the reaction network, it will reach an equilibrium state in which a determinant () will have a concentration   

    c() ∝|⟨|Ψ()⟩ | ^2.

The CTMC algorithm can be viewed as a way of driving this unimolecular reaction network stochastically to equilibrium. 

Notice, that the same equilibrium state will be reached irrespective of the relative magnitude of the rate constants k(←) versus k(←).

Within the reaction network picture, this is akin to catalyzing a reaction which can lead to faster equilibration but does not change the equilibrium state itself. This allows considerable freedom in how the algorithm is implemented. The rates chosen by us in Equation <ref> are only one example of a choice that one can make. In fact, as long as the algorithm satisfies the condition of detailed balance r(←)/r(←) = ρ()/ρ(), and the principle of ergodicity, it is guaranteed to sample the determinants with the correct probability.


Although the CTMC can be used in virtually all Monte Carlo simulations, to the best of our knowledge it has never been used in VMC before. This is partly because one has to calculate and store the quantities r(←) for a potentially large set of connected configurations and this is usually quite expensive. However, it is interesting to note that in the VMC algorithm, the quantities |⟨|Ψ()⟩/⟨|Ψ()⟩| are already used in the calculation of the local energy (see Equation <ref>) and just by storing those quantities the CTMC algorithm can be used with almost no overhead. We notice that using the CTMC algorithm leads to a shorter autocorrelation time and results in a more efficient VMC calculation.

A potential difficulty that might arise due the use of CTMC is that one could encounter a situation in which although a set of determinants have a large overlap with the current wavefunction, they are never reached because we start our Monte Carlo simulation from a determinant that is not connected to these determinants through a set of significant Hamiltonian transition matrix elements (ergodicity of the simulation is broken). For example let us imagine a situation in which two states |Ψ_1⟩ and |Ψ_2⟩ are nearly degenerate but have different irreducible representations. If we introduce a small perturbation in the Hamiltonian that breaks the symmetry, then the resulting ground state will be a linear combination of the two states |Ψ_1⟩ and |Ψ_2⟩. In such a situation, the CTCM algorithm when performed with a screening parameter ϵ, that is larger than the magnitude of the perturbation, will not be able to sample both states and will most likely be stuck in one or the other state depending on the initial determinant chosen to start the Monte Carlo calculation. This situation is difficult to diagnose because it is likely that the energy obtained will be reasonably accurate, but the properties such as correlation function, etc. will be quite inaccurate. Although it is difficult to be certain that such an eventuality is eliminated, it can be avoided to a large extent by fully localizing the orbital basis, and trying to break all possible spatial symmetries. This will ensure that the approximate symmetries of the problem will not influence the Monte Carlo walk.



 §.§ AMSGrad algorithm for optimizing the energy

The optimized wavefunction (Ψ()) is obtained by minimizing its energy with respect to its parameters . The gradient of the energy of the wavefunction with respect to   can be evaluated as follows

    _i =   ∂ E/∂_i
    
    =   ∂⟨Ψ() |H|Ψ()⟩/⟨Ψ() |Ψ()⟩/∂_i
    
    =    2⟨Ψ_i() |H - E|Ψ⟩/⟨Ψ() |Ψ()⟩
    
    =    2 ∑_ |⟨Ψ() |⟩|^2  ⟨Ψ_i()|⟩/⟨Ψ()|⟩⟨|H - E|Ψ()⟩/⟨|Ψ()⟩/∑_i |⟨Ψ()|⟩|^2
    
    =     2 ⟨⟨Ψ_i()|⟩/⟨Ψ()|⟩ (E_L[]-E)⟩_ρ_

where, |Ψ_i()⟩ = |∂Ψ()/∂_i⟩ and in going from the first line to the second we have assumed that all parameters are real. Note that during the calculation of the energy the determinants  are generated according to the probability distribution ρ_ and for each of these determinants a local energy E_L[] is evaluated. Thus to obtain the gradient of the energy, the only additional quantity needed is a vector of gradient ratios ⟨Ψ_i()|⟩/⟨Ψ()|⟩. We will show that these can be calculated at a cost that scales linearly with the number of parameters in the wavefunction. In the wavefunctions that we have used in this work, the number of parameters themselves scale quadratically with the size of the system. Thus the gradient can be calculated at a cost that scales no worse than the local energy.

The two most common algorithms for minimizing the energy of the wavefunctions in VMC are the linear method (LM)<cit.> and the stochastic reconfiguration (SR)<cit.>. The former is a second order method and can be viewed as an instance of the augmented Hessian method in which one repeatedly solves the generalized eigenvalue equation

    ([ E_0  ^T;       H ]) ([ 1; Δ ]) =   ([ 1 0; 0 S ]) ([ 1; Δ ])

where, S_ij = ⟨Ψ_i|Ψ_j⟩,  H_ij = ⟨Ψ_i|H|Ψ_j⟩, Δ is the update to the wavefunction parameters and |Ψ_i⟩ = 1/√(⟨Ψ|Ψ⟩)(|Ψ_i⟩ - ⟨Ψ|Ψ_i⟩/⟨Ψ|Ψ⟩|Ψ⟩). (Here we won't go into the details of how the matrices S and H are calculated.)

The SR algorithm can be thought to be performing projector Monte Carlo by repeated application of the propagator (I - Δ t H) to the wavefunction |Ψ(t)⟩ and at each step projecting the resulting wavefunction on to the tangent space constructed from the wavefunction gradients |Ψ_i⟩. It can be shown that this results in a linear equation


    S_ijΔ = -(Δ t)

which has to be solved at each step to obtain Δ, the update to the wavefunction parameters.
  Note that in both these methods, one has to solve an equation at each iteration which has a CPU cost of O(N_p^3) and a memory cost of O(N_p^2), where N_p is the number of wavefunction parameters, which we will show scales quadratically with the size of the system. The CPU cost can be reduced to O(n_s N_p) and O(n_s N^4) respectively for SR and LM method by using a direct method which avoid building the matrices, where n_s is the number of Monte Carlo samples used in each optimization iteration<cit.>.

  
  Both these methods are quite effective at optimizing the wavefunction, in particular, the LM requires fewer iterations and is often able to find the optimized wavefunction containing about 1000 parameters in less than 10 iterations. The effectiveness of the method is somewhat adversely affected while using direct method because one often needs to include level shifts to remove the singularity in the overlap matrix S, which can result in slower convergence.
  
  In this work we instead use a flavor of the adaptive stochastic gradient descent (SGD) method called AMSGrad<cit.>. The use of SGD methods have become popular in machine learning, where one is often interested in optimizing a cost function that depends non-linearly on a set of parameters and which can only be evaluated with a stochastic error by using batches of finite samples. This problem is of course very similar to the one we are interested in solving in VMC. The use of these adaptive stochastic gradient descent in the context of VMC was first proposed by Booth and coworkers<cit.>.
  
  The stochastic gradient descent methods have the advantage that, (1) their CPU and memory cost scales linearly with the number of wavefunction parameters and (2) the update in parameters depends strictly linearly on the gradient and thus does not introduce a systematic bias. The disadvantage of these methods is that traditionally they have been thought to be too slow and needing several thousand iterations to reach convergence, thus rendering them relatively ineffective. We show that by choosing appropriate parameters in the AMSGrad algorithm, we obtain convergence in a few tens to hundreds of iterations, which coupled with the fact that each iteration is much cheaper than LM and SR algorithm, makes it a powerful algorithm to solve challenging non-linear optimization problems in VMC.

A general algorithm for adaptive methods is as follows

    𝐦^(i) =    f(^(0),⋯, ^(i))
    𝐧^(i) =    g(^(0),⋯, ^(i))
    Δ_j =    -α^(i)𝐦_j^(i)/√(𝐧^(i)_j)

where f and g are some functions that take in all the past gradients and generate vectors 𝐦^(i) and 𝐧^(i), that are then used to calculate the updates Δ to the parameters. The various flavors on adaptive methods, ADAGrad<cit.>, RMSProp[T. Tieleman and G. Hinton. RmsProp: Divide the gradient by a running average of its recent magnitude. COURSERA: Neural Networks for Machine Learning, 2012.], ADAM<cit.> and AMSGrad<cit.> differ in the function f and g. In AMSGrad f and g respectively calculate the exponentially decaying moving average of the first and second moment 

    𝐦^(i) =    (1-β_1)  𝐦^(i-1) + β_1^(i)
    𝐧^(i) =   max(𝐧^(i-1), (1-β_2) 𝐧^(i-1) +  β_2(^(i)·^(i))

of the gradients () respectively, with the caveat that the second moment at iteration (i) is always greater than at iteration (i-1), i.e. 𝐧^(i)≥𝐧^(i-1). This ensures that the learning rate, α^(i)/√(𝐧^(i)), is a monotonically decreasing function of the number of iterations. This was shown to improve the convergence of AMSGrad for a synthetic problem over ADAM, which is essentially identical to AMSGrad with the only difference being that g calculates the exponentially decaying moving average of the second moment of the gradient. In our experiments with various VMC calculations we have found AMSGrad to always be superior to ADAM and RMSProp. For most of the systems studied in this article, unless otherwise specified, we have used the parameters α^(i) = 0.01, β_1 = 0.1, β_2=0.01, which are significantly more aggressive than the recommended values<cit.> of α^(i) = 0.001, β_1 = 0.1, β_2=0.001.



§ COMPUTATIONAL CONSIDERATIONS

In this section we will first briefly describe the wavefunction consisting of the product of a correlator product state<cit.> and a Slater determinant and then show that the cost of each stochastic sample is dominated by operations that cost O(N^2) for ab-initio Hamiltonians and O(N) for the Hubbard model.



 §.§ Correlator product states and Slater determinant
 
The wavefunction (|Ψ⟩) used in this work consists of a product of the correlator product states (CPS) (Ĉ) and a Slater determinant (|Φ⟩)

    |Ψ⟩ =   Ĉ|Φ⟩
    
    |Φ⟩ =   ∏_x<N(∑_iθ_ix a_i^†) |0⟩
    Ĉ =   ∏_λ( ∑_𝐧_λ c_𝐧_λP̂_𝐧_λ)

where, the subscripts (i, j, ⋯) represent the local orbital basis and subscripts (x, y, ⋯) represent the delocalized molecular orbitals, N is the number of electrons, θ is the matrix of molecular coefficients, λ represents the correlators (see Figure <ref>), 𝐧_λ represents all the Fock states in the correlator λ and finally the P̂_𝐧_λ = |𝐧_λ⟩⟨𝐧_λ| is the projector onto the Fock states of the correlator. CPS is an instance of the tensor network state, and is closely related to several other wavefunctions such as Huse-Elser wavefunction, the Jastrow factors used in the resonating valence bond states and Laughlin wavefunction. In all our calculations we use local correlators containing at most 5 sites. Although the number of parameters increase exponentially with the size of the correlator, the number of local correlators themselves only increase linearly with the size of the system. In this work we always use local correlators and thus the number of parameters in CPS c_𝐧λ only increases linearly with the size of the system, however the number of parameters in the matrix θ increases quadratically. 



 §.§ Local energy and Gradient ratios
 At each optimization step in VMC, one needs to evaluate the local energy and the gradient of the energy with respect to the wavefunction parameters. To evaluate these, one needs to be able to calculate the overlap ratio ⟨|Ψ()⟩/⟨|Ψ()⟩ for all the O(N^2) determinants, , that are connected to the current determinant, , with a non-negligible Hamiltonian transition matrix element, as well as the gradient ratio ⟨Ψ_i()|⟩/⟨Ψ()|⟩ for all the O(N_p) = O(N^2) parameters in the wavefunction.

The ratio of the overlap is equal to the product of the ratios of the overlaps with the CPS and the Slater determinant

    ⟨|Ψ()⟩/⟨|Ψ()⟩ = Ĉ[]/Ĉ[](θ[])/(θ[])

where, θ[] is the square matrix constructed by only retaining those rows and columns from θ that are occupied in  and Φ.

The CPS overlap with a determinant, Ĉ[], is just the product of the coefficients (c_𝐧_λ) of all correlators (λ) present in . Therefore, the ratio of the CPS overlap (Ĉ[]/Ĉ[]) is equal to the ratio of the coefficients for only those correlators that contain at least one of the four orbitals that differ between  and , as all correlators in common cancel. This ratio can be evaluated with an O(1) cost as follows. For each orbital in the system we maintain a list of all the correlators that contain that orbital. Because of the locality of the correlators, there are only a maximum of O(1) such local correlators per orbital and thus the ratio of the correlators can be calculated in a time that scales as O(1).

The cost of calculating the ratio of the Slater determinants ((θ[])/(θ[])) can be reduced from O(N^3) to O(N) by using the Woodbury lemma[M. Brookes, The Matrix Reference Manual, (2011); see online at http://www.ee.imperial.ac.uk/hp/staff/dmb/matrix/intro.html.] 

    (θ[])/(θ[]) = R[]_ai = ∑_x θ_ax(θ[]^-1)_xi

which holds when  is obtained from  by exciting an electron from orbital i to orbital a. Similar expressions exist for the ratio of the determinants that are related by a double excitation, for example, when two α electrons are excited from orbitals i and j to orbitals a and b, then the ratio is given by R[]_a,iR[]_b,j- R[]_b,iR[]_a,j. Thus, by precalculating and storing the matrix R[]_a,i the ratios of determinants that are related by single and double excitations can also be calculated with an O(1) cost. To calculate R one needs to have access to the inverse of the matrix θ[]. Performing direct inversion is expensive and has an O(N^3) cost. We avoid the expensive direct inversion by simply updating the inverse of the matrix using the Sherman-Morrison formula 

    θ[]^-1 = θ[]^-1 - (θ[]^-1 u) (v^Tθ[]^-1)/I + v^Tθ[]^-1u

that has a cost O(N^2), where  and  are related by a single excitation from orbital i to orbital a, u_x1 =δ_xi and v_x1 = θ_ax - θ_ix. Thus one only needs to calculate the inverse once at the beginning of the calculation with an O(N^3) cost and all subsequent updates can be calculated at an O(N^2) cost. It is also worth noting that calculating the matrix R[] itself has an O(N^3) cost, however, this can be reduced to O(N^2) by updating the matrix R[] at each iteration

    R[] = R[] - θ(θ[]^-1 u) (v^Tθ[]^-1)/I + v^Tθ[]^-1u

which has a cost O(N^2).

The elements of the vector of gradient ratios ⟨Ψ_i()|⟩/⟨Ψ()|⟩ for CPS and Slater determinant parameters can be calculated as 

    ⟨Ψ_c_𝐧λ()|⟩/⟨Ψ()|⟩ =   ⟨Ψ()|⟩/c_𝐧λ
    ⟨Ψ_θ_xi()|⟩/⟨Ψ()|⟩ =   ⟨Ψ()|⟩(θ[]^-1)_ix,
 
where each equation can be implemented in O(1) time as long as matrix θ[]^-1 is available.





 §.§ Computational scaling of the algorithm
 
The Table <ref> contains the formal computational scaling of the various steps of the algorithm. The cost of the algorithm is O(n_sN^2) down from the usual algorithm that has a cost of O(n_sN^4). The essential reason for the reduction in the cost is the introduction of the screening parameter ϵ, and the use of a first order gradient descent method. 



To reason about the scaling of n_s with the system size, let us define the scaling of a method as the cost of performing a calculation to obtain a constant error per electron as the size of the system increases. This is in line with the usual definition used to define a linear scaling method in electronic structure theory. With this definition the scaling of the method is closely linked to the size-extensivity of a method, for example, with this definition the scaling of a method that is not size-extensive  is not a very useful concept because the error per particle keeps on increasing and in the limit of a large enough system size the correlation energy obtained goes to zero. Let us consider the wavefunction of a large system composed of independent, unentangled, and identical subunits. The CPS-Slater wavefunction where both the correlators and Slater determinant are optimized has enough flexibility to describe the wavefunction of such a system by factorizing the overall wavefunction into a product of wavefunctions of individual subunits. For such a system the variance (σ^2) of the total wavefunction scales linearly with the number of such subunits or the size of the systems i.e. σ^2 ∝ N. The error estimate (e) of a Monte Carlo calculation that uses n_s^i uncorrelated samples is thus e ∝(N/n_s^i)^1/2. Note that n_s^i is not equal to the number of stochastic samples n_s, because in our calculations each stochastic update usually moves a single or at most two electrons. This indicates that there is likely a serial correlation of length N in the Monte Carlo samples and thus n_s ≈ n_s^i N. Thus e ∝(N^2/n_s )^1/2 and if we want a constant error per particle (e/N), then we need to perform a constant number of Monte Carlo iterations. In other words, n_s needed is independent of the size of the system if a constant error per particle is needed.

The analysis of scaling becomes somewhat more complicated because often the optimization algorithms are quite sensitive to the noise in the calculated energy and gradient. Thus, one might envision a scenario in which given a wavefunction one can calculate an accurate relative energy with an n_s that is independent of the system size, however, to optimize the wavefunction to minimize the energy a system size independent n_s is not sufficient. In Section <ref> we give numerical evidence using the hydrogen chain of increasing size to show that the adaptive SGD is stable with respect to the stochastic noise and displays only a weak system size dependent convergence rate with a constant n_s.




§ RESULTS

In this section we give more details about the saving in the CPU cost and error incurred due to the use of the screening parameter ϵ; we will also discuss the effectiveness of the AMSGrad in optimizing the wavefunction and will end with results on benchmark systems.



 §.§ Accuracy and scaling

Table <ref> shows the effect on the accuracy and CPU cost as the parameters ϵ is varied. Notice that there is almost a factor of 300 reduction in the CPU cost accompanied by an error of 41 milliHartree as we go from ϵ=0 to ϵ = 10^-3. However, interestingly, the optimized wavefunction itself is quite accurate as is indicated by the energies in the final column entitled “Final Energy". This energy is calculated by evaluating the energy using a tighter ϵ = 10^-8 with a wavefunction that was obtained by optimizing using a looser ϵ. Thus a useful strategy for calculating the energies is to perform the optimization with a loose ϵ and then calculate the final energy using a single shot calculation with a tighter ϵ. It should be noted that the calculations in Table <ref> were performed on a H_50 at a bond length of 2.2 a_0, and the improvements in CPU cost will be larger when one goes to larger systems or longer bond lengths. 


Next let us examine the scaling of the algorithm with the size of the system for a 1-D hydrogen chain and the 2-D Hubbard model. The results of the calculations are plotted in Figure <ref>. The figure shows the impressive reduction in the CPU time even for a chain of H_20 which only contains 20 electrons. Further, this speed up increases substantially with the size of the system because of the N^4 versus the N^1.5 scaling of the standard algorithm and the new algorithm respectively. The N^1.5 scaling of the Hydrogen chain is surprising given the fact that the lowest scaling operation in the Table <ref> is O(N^2). The reason for the lower scaling is as follows, the largest CPU time is consumed in calculating the local energy. In calculating the local energy one has to loop over determinants  in Equation <ref>, and we expect N^2 scaling because of the double excitation from a doubly occupied orbital i to an empty orbital a with the Hamiltonian transition matrix element equal to ⟨ a a| i i⟩. Because of the long range nature of the Coulomb interaction this integral decays slowly and we expect it to be important for even very large systems. However, this excitation often does not apply because most determinants  and  contain several singly occupied orbitals; the double occupation is suppressed due to the locality of the orbitals. For the hydrogen chain in the minimal basis with a bond length of 2.2 a_0, the integrals for the type ⟨ ab|ij⟩, where a,b and i,j are spatially close together are almost always very small.


The almost linear scaling of the Hubbard model is due to the fact that only a linear number of local single excitations are present. The scaling of the local energy calculations is thus almost exactly linear, but with the relatively small quadratic scaling expense of updating the determinant, the overall scaling gradually increases with the size of the system and asymptotically will approach quadratic scaling.



 §.§ Performance of the optimizer

Figure <ref> shows the energy per electron as the AMSGrad optimizer is used to minimize the energy of the wavefunction consisting of a product of 5-site CPS and a UHF determinant for H_20, H_40 , H_80 and H_160 molecules. The number of stochastic samples n_s per optimization step was equal to 240,000 and was the same for all hydrogen chain lengths. It is interesting to note that the rate of convergence of the different hydrogen chains is almost independent of the size of the system. The slightly slower convergence in case of H_160 and H_80 was due to the fact that in the initial few iterations we had to use a smaller step size of α = 0.001. After performing 5 iterations the step size was increased to α=0.01 and was subsequently kept constant. This was necessary to build up reasonable estimates of 𝐧^(i) and 𝐦^(i), the exponentially decaying moving average of the first and second moment of the gradient (see Equations <ref> and Equation <ref>) without which the energy tends to fluctuate wildly for these larger systems. Finally, Table <ref> shows that the same energy/electron down to three decimal places is obtained by carrying out the optimization using virtually the same setting in all four cases. This demonstrates that the AMSGrad optimizer is fairly tolerant to the absolute magnitude of the noise and is likely to deliver the same relative accuracy as long as the relative noise is kept constant with the size of the system. 

Perhaps the most impressive aspect of the calculation is that the CPU cost of performing the optimization for the H_160 molecule was merely 25 CPU hours, which when split between 2 nodes containing 24 processors each amounted to a wall time of just over half hour.







 §.§ 2-D Hubbard Model Benchmark

Here we present benchmark results on the 2-D Hubbard model  with periodic boundary conditions that is tilted by 45^o as shown in Figure <ref>. The tilted lattice is used because one can obtain a restricted Hartree Fock  solution, which simplifies the subsequent optimization of the CPS-UHF wavefunction. All results are calculated using 5-site overlapping correlators, that are centered on each site of the lattice as shown in Figure <ref>. The largest lattice considered here consisted of 162 orbitals and contained ∼ 2× 10^5 variational parameters. All optimizations took less than 250 iterations with each iteration using between 7× 10^6-8× 10^6 Monte Carlo samples, embarrassingly parallelized over 100 cpu cores. The overall cost of the largest calculation was about 700 CPU hours which equaled 7 hours of wallclock time. This is quite satisfactory for a quantum Monte Carlo calculation of this size. It is worth pointing out that the converged results that we have obtained for the 98 site Hubbard model with U/t=8 is 0.004 /t lower than the one obtained using the RMSprop algorithm<cit.>. This highlights the importance of using aggressive settings in the adapted SGD methods, but more importantly points to the fact that it is often difficult to decide when the optimization has converged. This is because one often observes a relatively long tail at the end of optimization where the energy very gradually decreases with the number of iterations. We have noticed such an effect in our calculations as well and thus it is difficult to estimate the accuracy of the calculations presented here. This provides one motivation for using the second order methods such as LM, which can be used to greatly speed up the convergence when the wavefunction parameters are already close to their optimal values. 





§ CONCLUSIONS

In this work we have introduced several innovations including the efficient screening of integrals, use of continuous time Monte Carlo to sample determinants and finally the use of AMSGrad which displays robust and efficient convergence of the wavefunction in the model systems studied here.

Out of the three innovations introduced above, the first two can be straightforwardly used for any system, however, more work is needed to study the effectiveness of AMSGrad and other such adaptive SGD methods in optimizing the wavefunction of more realistic systems of interest in quantum chemistry. Such work is already underway and we are also exploring a suitable combination of SGD methods in conjunction with direct second order methods; using the former for the bulk of the optimization and switching to the latter at the end of the optimization when the wavefunction is already nearly converged. 

Finally, we are also exploring two different approaches to go beyond the VMC framework and correct the shortcomings of the wavefunction. First, is the use of stochastic perturbation theory<cit.> that was recently generalized to correct any variational wavefunction. Second, is the use of fixed node orbital-space projector Monte Carlo<cit.> that is guaranteed to provide variational energies which are at least as good as the VMC energies.




§ ACKNOWLEDGEMENTS

We would like to thank Cyrus Umrigar, Eric Neuscamman and George Booth for several helpful discussions. The funding for this project was provided by the national science foundation through the grant CHE-1800584.




@ifundefinedendmcitethebibliography
  
58








f
subitem(mcitesubitemcount)

  
  

[Foulkes et al.(2001)Foulkes, Mitas, Needs, and
  Rajagopal]RevModPhys.73.33
Foulkes, W. M. C., Mitas, L., Needs, R. J., and Rajagopal, G. (2001) Quantum
  Monte Carlo simulations of solids. Rev. Mod. Phys. 73,
  33–83
[Nightingale and Umrigar(1999)Nightingale, and Umrigar]NigUmr-BOOK-99
Nightingale, M. P., and Umrigar, C. J., Eds. Quantum Monte Carlo Methods
  in Physics and Chemistry; NATO ASI Ser. C 525; Kluwer: Dordrecht, 1999
[Toulouse et al.(2015)Toulouse, Assaraf, and
  Umrigar]TouAssUmr-AQC-15
Toulouse, J., Assaraf, R., and Umrigar, C. J. Introduction to the variational
  and diffusion Monte Carlo methods. Advances in Quantum Chemistry. 2015; p
  285
[Kolorenc and Mitas(2011)Kolorenc, and Mitas]KolMit-RPP-11
Kolorenc, J., and Mitas, L. (2011) Applications of quantum Monte Carlo
  methods in condensed systems. Rep. Prog. Phys. 74
[Becca and Sorella(2017)Becca, and Sorella]sorellabook
Becca, F., and Sorella, S. Quantum Monte Carlo Approaches for Correlated
  Systems; Cambridge University Press, 2017
[McMillan(1965)]Mcmillan
McMillan, W. L. (1965) Ground State of Liquid He^4. Phys.
  Rev. 138, A442–A451
[Ceperley et al.(1977)Ceperley, Chester, and
  Kalos]CepCheKal-PRB-77
Ceperley, D., Chester, G. V., and Kalos, M. H. (1977) Phys. Rev. B
  16, 3081
[GRI()]GRIMM1971134
Monte-Carlo solution of Schr
[Anderson(1975)]Anderson75
Anderson, J. B. (1975) A random walk simulation of the Schrodinger equation:
  H+3. The Journal of Chemical Physics 63, 1499–1503
[D. M. Ceperley and B. J. Alder(1980)]CeperleyAlder80
D. M. Ceperley and B. J. Alder, (1980) Ground state of the electron gas by a
  stochastic method. 45, 566
[Zhang and Krakauer(2003)Zhang, and Krakauer]ZhaKra-PRL-03
Zhang, S., and Krakauer, H. (2003) Quantum Monte Carlo method using phase-free
  random walks with slater determinants. Phys. Rev. Lett. 90,
  136401
[Mario and Shiwei()Mario, and Shiwei]Mario-WIREs
Mario, M., and Shiwei, Z. Ab initio computations of molecular systems by the
  auxiliary field quantum Monte Carlo method. Wiley Interdisciplinary
  Reviews: Computational Molecular Science 0, e1364
[Runge(1992)]Runge1992
Runge, K. J. (1992) Quantum Monte Carlo calculation of the long-range order in
  the Heisenberg antiferromagnet. Physical Review B 45
[Trivedi and Ceperley(1990)Trivedi, and Ceperley]Trivedi1990
Trivedi, N., and Ceperley, D. M. (1990) Ground-state correlations of quantum
  antiferromagnets: A Green-function Monte Carlo study. Physical Review
  B 
[Sorella(1998)]Sorella98
Sorella, S. (1998) Green function Monte Carlo with stochastic
  reconfiguration. Physical Review Letters 80, 4558–4561
[Booth et al.(2009)Booth, Thom, and Alavi]BooThoAla-JCP-09
Booth, G. H., Thom, A. J. W., and Alavi, A. (2009) Fermion Monte Carlo
  without fixed nodes: A game of life, death, and annihilation in Slater
  determinant space. J. Chem. Phys. 131
[Cleland et al.(2010)Cleland, Booth, and
  Alavi]CleBooAla-JCP-10
Cleland, D., Booth, G. H., and Alavi, A. (2010) Communications: Survival of
  the fittest: Accelerating convergence in full configuration-interaction
  quantum Monte Carlo. J. Chem. Phys. 132
[Petruzielo et al.(2012)Petruzielo, Holmes, Changlani,
  Nightingale, and Umrigar]PetHolChaNigUmr-PRL-12
Petruzielo, F. R., Holmes, A. A., Changlani, H. J., Nightingale, M. P., and
  Umrigar, C. J. (2012) Semistochastic Projector Monte Carlo Method.
  Phys. Rev. Lett. 109, 230201
[Holmes et al.(2016)Holmes, Changlani, and Umrigar]Holmes2016
Holmes, A. A., Changlani, H. J., and Umrigar, C. J. (2016) Efficient Heat-Bath
  Sampling in Fock Space. Journal of Chemical Theory and Computation
  12, 1561–1571
[Zhao et al.(2017)Zhao, Ido, Morita, and
  Imada]PhysRevB.96.085103
Zhao, H.-H., Ido, K., Morita, S., and Imada, M. (2017) Variational Monte Carlo
  method for fermionic models combined with tensor networks and applications to
  the hole-doped two-dimensional Hubbard model. Phys. Rev. B 96,
  085103
[Tahara and Imada(2008)Tahara, and Imada]Tahara2008
Tahara, D., and Imada, M. (2008) Variational Monte Carlo method combined with
  quantum-number projection and multi-variable optimization. Journal of
  the Physical Society of Japan 77, 1–16
[Neuscamman(2013)]Eric2013
Neuscamman, E. (2013) The Jastrow antisymmetric geminal power in Hilbert
  space: Theory, benchmarking, and application to a novel transition state.
  The Journal of Chemical Physics 139, 194105
[Ferrari et al.(2018)Ferrari, Parola, Sorella, and
  Becca]PhysRevB.97.235103
Ferrari, F., Parola, A., Sorella, S., and Becca, F. (2018) Dynamical structure
  factor of the J_1J_2 Heisenberg model in one
  dimension: The variational Monte Carlo approach. Phys. Rev. B
  97, 235103
[Schwarz et al.(2017)Schwarz, Alavi, and
  Booth]PhysRevLett.118.176403
Schwarz, L. R., Alavi, A., and Booth, G. H. (2017) Projector Quantum Monte
  Carlo Method for Nonlinear Wave Functions. Phys. Rev. Lett.
  118, 176403
[Sorella et al.(2007)Sorella, Casula, and Rocca]Sorella07
Sorella, S., Casula, M., and Rocca, D. (2007) Weak binding between two
  aromatic rings: Feeling the van der Waals attraction by quantum Monte Carlo
  methods. The Journal of Chemical Physics 127, 14105
[Neuscamman(2016)]Neuscamman2016
Neuscamman, E. (2016) Improved optimization for the cluster Jastrow
  antisymmetric geminal power and tests on triple-bond dissociations.
  Journal of Chemical Theory and Computation 12, 3149–3159
[Kurita et al.(2015)Kurita, Yamaji, Morita, and
  Imada]PhysRevB.92.035122
Kurita, M., Yamaji, Y., Morita, S., and Imada, M. (2015) Variational Monte
  Carlo method in the presence of spin-orbit interaction and its application to
  Kitaev and Kitaev-Heisenberg models. Phys. Rev. B 92,
  035122
[Knowles and Handy(1984)Knowles, and Handy]knowles1984new
Knowles, P., and Handy, N. (1984) A new determinant-based full configuration
  interaction method. Chem. Phys. Lett. 111, 315–321
[White(1992)]White1992
White, S. R. (1992) Density matrix formulation for quantum renormalization
  groups. Phys. Rev. Lett. 69, 2863
[White(1993)]White1993
White, S. R. (1993) Density-matrix algorithms for quantum renormalization
  groups. Phys. Rev. B 48, 10345
[Östlund and Rommer(1995)Östlund, and Rommer]Ostlund1995
Östlund, S., and Rommer, S. (1995) Thermodynamic limit of density matrix
  renormalization. Phys. Rev. Lett. 75, 3537
[Gutzwiller(1963)]PhysRevLett.10.159
Gutzwiller, M. C. (1963) Effect of Correlation on the Ferromagnetism of
  Transition Metals. Phys. Rev. Lett. 10, 159–162
[Changlani et al.(2009)Changlani, Kinder, Umrigar, and
  Chan]Changlani2009
Changlani, H. J., Kinder, J. M., Umrigar, C. J., and Chan, G. K. L. (2009)
  Approximating strongly correlated wave functions with correlator product
  states. Physical Review B - Condensed Matter and Materials Physics
  80, 1–8
[Neuscamman et al.(2011)Neuscamman, Changlani, Kinder, and
  Chan]Neuscamman2011
Neuscamman, E., Changlani, H., Kinder, J., and Chan, G. K. L. (2011)
  Nonstochastic algorithms for Jastrow-Slater and correlator product state
  wave functions. Physical Review B - Condensed Matter and Materials
  Physics 84, 1–9
[Bartlett and Musiał(2007)Bartlett, and
  Musiał]RevModPhys.79.291
Bartlett, R. J., and Musiał, M. (2007) Coupled-cluster theory in quantum
  chemistry. Rev. Mod. Phys. 79, 291–352
[Verstraete et al.(2008)Verstraete, Murg, and Cirac]peps
Verstraete, F., Murg, V., and Cirac, J. (2008) Matrix product states, projected
  entangled pair states, and variational renormalization group methods for
  quantum spin systems. Advances in Physics 57, 143–224
[Holmes et al.(2016)Holmes, Tubman, and Umrigar]Holmes2016b
Holmes, A. A., Tubman, N. M., and Umrigar, C. J. (2016) Heat-Bath
  Configuration Interaction: An Efficient Selected Configuration Interaction
  Algorithm Inspired by Heat-Bath Sampling. 
[Huron et al.(1973)Huron, Malrieu, and Rancurel]Huron1973
Huron, B., Malrieu, J.-P., and Rancurel, P. (1973) Iterative perturbation
  calculations of ground and excited state energies from multiconfigurational
  zeroth-order wavefunctions. The Journal of Chemical Physics
  58, 5745
[Wei and Neuscamman()Wei, and Neuscamman]ericScaling
Wei, H., and Neuscamman, E. https://arxiv.org/pdf/1806.08778.pdf
[Bortz et al.(1975)Bortz, Kalos, and Lebowitz]BORTZ197510
Bortz, A., Kalos, M., and Lebowitz, J. (1975) A new algorithm for Monte Carlo
  simulation of Ising spin systems. Journal of Computational Physics
  17, 10 – 18
[Gillespie(1976)]GILLESPIE1976403
Gillespie, D. T. (1976) A general method for numerically simulating the
  stochastic time evolution of coupled chemical reactions. Journal of
  Computational Physics 22, 403 – 434
[Umrigar et al.(2007)Umrigar, Toulouse, Filippi, Sorella, and
  Hennig]UmrTouFilSorHen-PRL-07
Umrigar, C. J., Toulouse, J., Filippi, C., Sorella, S., and Hennig, R. G.
  (2007) Alleviation of the fermion-sign problem by optimization of many-body
  wave functions. Phys. Rev. Lett. 98, 110201
[Toulouse and Umrigar(2007)Toulouse, and Umrigar]Toulouse2007
Toulouse, J., and Umrigar, C. J. (2007) Optimization of quantum Monte Carlo
  wave functions by energy minimization. Journal of Chemical Physics
  126
[Toulouse and Umrigar(2008)Toulouse, and Umrigar]TouUmr-JCP-08
Toulouse, J., and Umrigar, C. J. (2008) Full optimization of Jastrow-Slater
  wave functions with application to the first-row atoms and homonuclear
  diatomic molecules. J. Chem. Phys. 128, 174101
[Sorella(2001)]PhysRevB.64.024512
Sorella, S. (2001) Generalized Lanczos algorithm for variational quantum Monte
  Carlo. Phys. Rev. B 64, 024512
[Neuscamman et al.(2012)Neuscamman, Umrigar, and
  Chan]Neuscamman2012
Neuscamman, E., Umrigar, C. J., and Chan, G. K. L. (2012) Optimizing large
  parameter sets in variational quantum Monte Carlo. Physical Review B -
  Condensed Matter and Materials Physics 
[Reddi et al.(2017)Reddi, Kale, and Kumar]Reddi2017
Reddi, S. J., Kale, S., and Kumar, S. (2017) on the Convergence of Adam and
  Beyond. 1–23
[C. Duchi et al.(2011)C. Duchi, Hazan, and Singer]adagrad
C. Duchi, J., Hazan, E., and Singer, Y. (2011) Adaptive Subgradient Methods for
  Online Learning and Stochastic Optimization. 12, 2121–2159
[Note1()]Note1
T. Tieleman and G. Hinton. RmsProp: Divide the gradient by a running average of
  its recent magnitude. COURSERA: Neural Networks for Machine Learning,
  2012.
[Kingma and Ba(2014)Kingma, and Ba]adam
Kingma, D., and Ba, J. (2014) Adam: A Method for Stochastic Optimization.
  
[Note2()]Note2
M. Brookes, The Matrix Reference Manual, (2011); see online at
  http://www.ee.imperial.ac.uk/hp/staff/dmb/matrix/intro.html.
[LeBlanc et al.(2015)LeBlanc, Antipov, Becca, Bulik, Chan,
  Chung, Deng, Ferrero, Henderson, Jiménez-Hoyos, Kozik, Liu, Millis,
  Prokof'ev, Qin, Scuseria, Shi, Svistunov, Tocchio, Tupitsyn, White, Zhang,
  Zheng, Zhu, and Gull]LeBlanc2015
LeBlanc, J. P. F. et al.  (2015) Solutions of the Two Dimensional
  Hubbard Model: Benchmarks and Results from a Wide Range of Numerical
  Algorithms. 041041, 1–28
[Qin et al.(2016)Qin, Shi, and Zhang]Qin2016
Qin, M., Shi, H., and Zhang, S. (2016) Benchmark study of the two-dimensional
  Hubbard model with auxiliary-field quantum Monte Carlo method.
  Physical Review B 
[Sharma et al.(2017)Sharma, Holmes, Jeanmairet, Alavi, and
  Umrigar]ShaHolUmr-JCTC-17
Sharma, S., Holmes, A. A., Jeanmairet, G., Alavi, A., and Umrigar, C. J. (2017)
  Semistochastic Heat-bath Configuration Interaction method: selected
  configuration interaction with semistochastic perturbation theory. J.
  Chem. Theory Comput. 13, 1595–1604
[Guo et al.(2018)Guo, Li, and Chan]doi:10.1063/1.5031140
Guo, S., Li, Z., and Chan, G. K.-L. (2018) Communication: An efficient
  stochastic algorithm for the perturbative density matrix renormalization
  group in large active spaces. The Journal of Chemical Physics
  148, 221104
[Sharma()]sharma18
Sharma, S. https://arxiv.org/abs/1803.04341
[van Bemmel et al.(1994)van Bemmel, ten Haaf, van Saarloos, van
  Leeuwen, and An]PhysRevLett.72.2442
van Bemmel, H. J. M., ten Haaf, D. F. B., van Saarloos, W., van Leeuwen, J.
  M. J., and An, G. (1994) Fixed-Node Quantum Monte Carlo Method for Lattice
  Fermions. Phys. Rev. Lett. 72, 2442–2445
[ten Haaf et al.(1995)ten Haaf, van Bemmel, van Leeuwen, van
  Saarloos, and Ceperley]Haaf1995
ten Haaf, D. F. B., van Bemmel, H. J. M., van Leeuwen, J. M. J., van
  Saarloos, W., and Ceperley, D. M. (1995) Proof for an upper bound in
  fixed-node Monte Carlo for lattice fermions. Phys. Rev. B 51,
  13039


