Efficient local energy evaluation for multi-Slater wave functions in orbital space quantum Monte Carlo

Particle-hole excitations from a mean-field reference state have long been used to encode correlation in wave functions. In single-reference wave function methods, one usually starts from a single Hartree–Fock (HF) configuration, and the electron correlation is built on top of this zeroth-order state by adding particle-hole excitations, often variationally or perturbatively.¹ When electron interactions are strong, high rank particle-hole excitations from an HF reference are required for a satisfactory description of the electronic structure. It is computationally infeasible to include arbitrarily high-rank excitations for even moderately sized systems due to the rapid increase in their number with system size. Sometimes, the strong interactions are confined to a small set of orbitals and electrons termed the active space.² For small enough active spaces, it is possible to include all possible excitations or, equivalently, all configurations in the Hilbert space in a variational treatment, and to solve the problem exactly in the active space. This is usually only feasible for active spaces smaller than about 20 electrons in 20 orbitals.³ Various approximate methods have been devised to stretch this restrictive bound on the active space size.^4–8 Recently, selected configuration interaction (CI) methods, first used about 50 years ago,^9,10 have attracted renewed attention because of their ability to target the most important configurations contributing to the wave function at a relatively small cost.^6,11–14 Fast implementations of such techniques are now available that allow calculations with millions of configurations on even small work-stations.

Due to the remarkable accuracy and easy computability of selected CI wave functions, they have been used in quantum Monte Carlo (QMC) approaches as trial wave functions. In real space QMC, they can be used as part of the Jastrow multi-Slater wave functions, where the Jastrow factors serve to capture the dynamic correlation efficiently, thereby greatly truncating the length of the selected CI expansion required for a given accuracy.^15–17 Multi-Slater wave functions have also been used in the second-quantized setting of auxiliary field QMC (AFQMC),^7,18,19 where the error incurred by the phaseless approximation²⁰ can be controlled by using more accurate trial wave functions.^21,22 For calculating properties of wave functions in QMC methods, one needs to evaluate the so-called local quantities for a walker in the MC run. For example, for calculating the energy, one averages the local energy given by

where $Ĥ$ is the Hamiltonian, |ϕ⟩ = ∑_Ic_I|ϕ_I⟩ is the multi-Slater wave function with |ϕ_I⟩ being configurations in the CI expansion obtained by particle-hole excitations from a reference configuration, n_c is the number of configurations, and |n⟩ is the walker given by electronic positions in real space QMC and by orbital occupations in orbital space QMC. At first glance, this expression seems to suggest that the cost of calculating the local energy for the multi-Slater wave function is n_c times the cost of calculating it for a single configuration. However, this cost can be significantly reduced since CI configurations are not independent but are obtained by a small set of excitations from the reference. Various algorithms have been proposed to achieve this speed-up in real space QMC,^13,23–25 most efficient being the one recently proposed by Filippi and co-workers.^24,25 Their algorithm has a cost scaling of O(n³ + n_c), where n is the system size, which has allowed millions of configurations to be included in trial wave functions in real space. While there has been some prior work on using multi-Slater trial wave functions in AFQMC,^26,27 no such drastic speed-up has been reported for calculating energy in orbital space algorithms, thus restricting the number of configurations that can be used. In this article, we present an algorithm for calculating the local energy of a multi-Slater wave function that achieves speed-ups analogous to the real space algorithm for orbital space QMC. We will use the variational Monte Carlo (VMC) treatment of Jastrow multi-Slater wave functions to formulate our technique.

In the following, we start by briefly describing the wave function and defining the required notation. Then, we show how local energy calculations can be sped up by storing some intermediates. Finally, the method is used for polyacetylene ground-state calculations to demonstrate its efficiency.

The Jastrow multi-Slater wave function ansatz is given by

where J_ij are real numbers, $n^i$ is the number operator for spin–orbital i, and |ϕ_I⟩ are electronic configurations forming the multi-configurational state |ϕ⟩. We will suppress spin indices for brevity throughout. Let n and n_c be the number of electrons and number of configurations in the CI expansion, respectively. For scaling considerations, we will assume the number of electrons to be proportional to the number of orbitals, serving as a proxy for the system size.

In general, the single-particle orbital sets used for defining the Jastrow and the CI expansion are different. The second-quantized operators for the spatially local orbitals used in the Jastrow factor²⁸ will be denoted as ${L^μ,L^μ†}$⁠, while those for the canonical orbitals used in the CI expansion will be denoted as ${Ĉν,Ĉν†}$⁠. Indices i, j, a, b will be reserved for local orbitals, and p, q, t, u will be reserved for the canonical ones. We will assume both sets of orbitals to be orthogonal, with the unitary transformation relating them given by

where $Mνμ$ are real numbers forming the coefficient matrix. We define the following slices of the coefficient matrix:

where superscripts show the row indices and subscripts show the column indices used for slicing. Henceforth, all indices will be assumed to be relative to make the notation compact. For example, i_σ will label the i_σth occupied local orbital, t_ν will label the t_νth empty canonical orbital, and so on. All configurations in the CI expansion are related to a single configuration |ϕ₀⟩, termed the reference configuration, via a set of canonical orbital excitations,

The degree of excitation for CI configurations is assumed to be of O(k). In practice, k is usually less than 10. A summary of the notation is shown in Fig. 1.

In this article, we are interested in calculating the energy of this wave function, although similar considerations apply to some other physical properties of interest as well. We calculate the energy using variational Monte Carlo sampling,

where the walkers |n⟩ are configurations in the Hilbert space, and the random walk is performed according to the probability distribution ⟨ψ|ψ⟩. The quantity being averaged, the local energy, is given by

where the configurations |m⟩ are generated by the action of the Hamiltonian on the walker |n⟩. The local energy calculation is the rate limiting step in a VMC calculation, and thus, efficient algorithms for calculating it are essential for feasible VMC simulations. We use the local orbitals used in the Jastrow to define the walkers. This is a sensible choice because the Jastrow is diagonal in this basis and the Hamiltonian, because of its local interactions, has only a quadratically scaling number of two-electron integrals, as opposed to quartic in a general basis. Thus, the application of the Jastrow and the Hamiltonian onto a configuration can be performed efficiently in a local basis.²⁹ This leaves us with the task of calculating the overlap of a walker in a local basis with the multi-Slater wave function defined in a canonical basis. We outline the algorithm for doing this efficiently for local energy calculations in Secs. II B –II D.

For calculating overlap and local energy, we need to evaluate the matrix elements of strings of excitation operators between configurations expressed in local and canonical bases. This section summarizes the identities that will be used to calculate these transition matrix elements. We will use the generalized Wick’s theorem to obtain them.

First, note the following pair contractions:

These can be derived in many different ways.^7,30–33 A brief derivation of the first identity can be found in  Appendix A. Note that given one of the four identities, the rest follow by simply transforming between the $L^$ and $Ĉ$ operators. Using these pair contractions and the generalized Wick’s theorem, we get the following transition matrix element:

where the RHS is the determinant of a block matrix with blocks given by the indicated slices. The RA⁻¹ elements arise from contractions between $L^$ operators, the A⁻¹C elements arise from contractions between the $Ĉ$ operators, and the A⁻¹ and RA⁻¹C − B elements arise from cross contractions between $L^$ and $Ĉ$ operators. The determinant structure results from the fermionic parity factors. This equation will be used often in Secs. II C and II D. For convenience, we define the matrix Δ as

We will not include the Jastrow factor in the following since it can be handled the same way as for the Jastrow Slater wave function. The overlap of a walker with the multi-Slater wave function is given by

The first term in this sum, the overlap of the walker with the reference configuration, is

where |0⟩ is the vacuum state. This follows from the usual (vacuum) Wick’s theorem. Similarly, the overlap of the walker with any other configuration in the CI expansion is also given by a determinant of size n × n. The cost of calculating these determinants scales as O(n³). Thus, the cost of evaluating the total overlap by calculating individual determinant overlaps from scratch has a prohibitive cost scaling of O(n³n_c). We can improve upon this by storing a few matrices at the start of the calculation.

From Eq. (9), we have

where |ϕ_I⟩ is obtained from |ϕ₀⟩ by k excitations [Eq. (5)]. Thus, by calculating the matrix A⁻¹C and det(A) at cost O(n³) once, the overlap with any configuration in the CI expansion can be calculated at cost O(k³). We note that the matrix A⁻¹C can be updated using the Sherman–Morrison formula at cost O(n²) when the walker makes a move and does not need to be calculated from scratch every time. Since k is usually much smaller than n, the total cost scaling of this approach, given by O(n² + n_ck³), proves to be superior in almost all cases.

The local energy of the multi-Slater wave function is given by

where |m⟩ are generated from the action of the Hamiltonian on the walker |n⟩. As mentioned before, an ab initio Hamiltonian generates O(n⁴) excitations, which is reduced to O(n²) when a local basis is used. Because of this large number of excitations, the calculation of local energy constitutes the most demanding part of a VMC calculation. We present two algorithms to evaluate the local energy.

Following the discussion about calculating overlaps, ⟨m|ϕ⟩ is given by a sum of n_c determinants. We can again avoid evaluating these determinants individually from scratch by storing a few matrices. Suppose |m⟩ is obtained from |n⟩ by l excitations as

The overlap of |m⟩ with the excited configuration |ϕ_I⟩ defined in Eq. (5) is given precisely by Eq. (9) as

Thus, by calculating RA⁻¹, A⁻¹C, and Δ once at cost O(n³), each overlap of the form ⟨m|ϕ_I⟩ can be evaluated at cost O(k³) since an ab initio Hamiltonian generates only up to double excitations, i.e., l ≤ 2. Note that all of these matrices can be efficiently updated at cost O(n²) in a Monte Carlo sampling run.

Having calculated these three matrices, we can directly perform the sum in Eq. (14) to get the local energy,

where the overlap ratios can be obtained using Eq. (16). This algorithm has a cost scaling of O(n⁴n_ck³), if screening of the Hamiltonian elements due to locality of orbitals is not considered, and O(n²n_ck³), if it is. This simple algorithm is easy to implement and allows for incorporating Hamiltonian screening. However, the cost scaling is rather steep, and one is restricted to fairly short CI expansions. Note that this algorithm is similar to the real space algorithm proposed by Clark et al.²³ It improves upon their algorithm by obviating the need to calculate extra intermediate matrices.

In the overlap calculation, it was possible to separate the n_c factor from the system size n in the cost scaling. Section II D 2 achieves this for local energy evaluation by storing some intermediates. First, we partition the local energy based on the rank of the Hamiltonian excitation,

where $|nia$ is a configuration obtained from |n⟩ by a single excitation {i → a} and, similarly, $|nijab$ is a configuration obtained from |n⟩ by a double excitation i → a, j → b. The Hamiltonian matrix elements are defined as $Hai[n]=⟨n|H|nia⟩$ and $Habij=⟨n|H|nijab⟩$⁠. In the following, we will suppress the explicit dependence of certain quantities on the walker whenever it is clear from the context.

Consider the calculation of E^S, arising from single excitations. We can further split it into contributions from individual configurations,

$E0S$ can be calculated directly as in Sec. II D 1,

For I ≠ 0, suppose |ϕ_I⟩ is obtained from |ϕ₀⟩ by excitations {p₁ → t₁, …, p_k → t_k} [Eq. (5)]. Using Eq. (16), $EIS$ can be expressed as

Note that the determinants appearing in this sum all have the same k × k block in the right-bottom corner and they only differ in the first row and column. To benefit from this pattern, we Laplace expand the determinants along the first column as

where the ordered set of indices p₁, …, p_k$\$p_μ denotes the set p₁, …, p_k excluding p_μ. In the second line, we have used the linearity of determinants to move the Hamiltonian elements and the sum over Hamiltonian excitations inside the determinants. Evidently, the first term is proportional to $E0S$⁠. To efficiently evaluate the second term, we define the following intermediate:

where the sums have been arranged to show the lowest cost scaling order of contractions given by O(n³). By building this intermediate during the evaluation of $E0S$⁠, the local energy contribution of configuration I can be calculated at cost O(k⁴) as

The total cost scaling of calculating the part of the local energy due to single excitations using Sec. II D 2 is O(n³ + n_ck⁴). Thus, we have successfully separated the n_c factor from the system size in the cost scaling. We note that this part of the algorithm, dealing with single excitations, is similar to the real space algorithm of Filippi et al.,²⁴ albeit formulated differently.

Now, we turn to the calculation of the double excitation part of the local energy, E^D. Again, we split it into contributions due to individual configurations,

Considering |ϕ_I⟩ defined in Eq. (5), its contribution is given by

Note that the determinants appearing in this sum all have the same k × k block in the right-bottom corner and they only differ in the first two rows and columns. To exploit this fact, we Laplace expand the determinants on the RHS along the first two columns,

The first term can again be directly related to the local energy contribution of the reference given as

To efficiently evaluate the second and third terms, we build the intermediates

Both of these can be built at cost O(n⁴) during the evaluation of $E0D$⁠. For example, one can see how D₁ can be calculated at this cost since it involves contractions of the type

obtained by expanding the determinant. Using D₁ and D₂, the second and third terms can be calculated as before due to the linearity of determinants at cost O(k⁴). To calculate the final term, we use the intermediate

This intermediate can be built at cost O(n⁵) since it involves sums like

obtained by expanding the two determinants. Using D₃, the final term can be calculated at cost O(k⁶). Therefore, the total cost scaling of the local energy calculation using Sec. II D 2 is O(n⁵ + n_ck⁶).

Unlike Sec. II D 1, there is no obvious way to use the screening of Hamiltonian elements to reduce the cost scaling of building the intermediates in Sec. II D 2. However, we note that the equation for building the intermediate D₃ resembles a Hamiltonian integral transformation. Techniques like density fitting³⁴ and Cholesky decomposition³⁵ are employed to reduce the cost of such transformations. Tensor hypercontraction³⁶ can be used to reduce the cost scaling to O(n⁴). Another possible way of improving the cost is by updating the intermediates as the walker moves during a Monte Carlo run, instead of building them from scratch every time. Incorporating such techniques into our algorithm will be a task for future research.

With some modifications, these algorithms can be used for calculating local energies of multi-Slater trial wave functions in AFQMC as well. This is outlined briefly in  Appendix B. Finally, we note that the Jastrow factor overlap ratios, which we have ignored in the discussion above, can be absorbed into the intermediates exactly like the Hamiltonian matrix elements. Other details of the VMC procedure can be found in  Appendix C.

In this section, we apply the above formalism to the truncated trans-polyacetylene chain C₂₈H₃₀. A model geometry, with uniform bond lengths given by l(C=C) = 1.34 Å, l(C—C) = 1.45 Å, and l(C—H) = 1.08 Å and 120° bond angles, was used. We used PySCF³⁷ to generate molecular integrals for the 6-31g basis set and the heat bath CI (HCI)^6,38 program Dice to perform an heat bath configuration interaction (HCI) self consistent field³⁹ calculation with the π active space consisting of 28 electrons in 28 orbitals. Intrinsic bond orbitals⁴⁰ obtained by localizing the π orbitals were used in the Jastrow factor. They roughly resemble the p_z atomic orbital on each carbon atom.

We first look at the computational performance of the two algorithms. These calculations were performed on an Intel Xeon Gold 6150 2.7 GHz central processing unit (CPU). Figure 2 shows the computational cost scaling of the two algorithms as a function of the number of configurations in the wave function, showing the cost of calculating 100 local energy samples during a continuous time Monte Carlo run. The configurations were generated using an stochastic HCI (SHCI) calculation with ϵ₁ = 5 × 10⁻⁵, and only those up to quadruply excited from the reference were retained. These configurations constitute the majority of all leading configurations in the SHCI expansion: Of the 5 × 10⁵ leading configurations, only 954 were higher than quadruply excited. Section II D 1 shows a linear scaling with the number of configurations as expected. We store the Hamiltonian in a heat bath format, which allows efficient screening of the two-electron integrals, where integrals below a screening threshold ϵ are ignored. In Fig. 2, we show the cost scaling of Sec. II D 1 for two ϵ values to demonstrate the effect of screening. Because there is essentially only one orbital on each carbon, screening is very efficient in this case and leads to a significant reduction in the cost of local energy evaluation. However, even with aggressive screening, the calculations become very expensive as the number of configurations is increased due to the O(n²n_ck³) scaling.

On the other hand, Sec. II D 2 has a much less severe cost scaling. Up to 1 × 10⁴ configurations, the cost is almost independent of the number of configurations. It is dominated by the O(n⁵) scaling calculation of the intermediate in Eq. (32). Because we have implemented this as a dense tensor contraction, screening does not affect the cost. Beyond 1 × 10⁴ configurations, the linear scaling O(n_ck⁶) starts to become dominant. As more quadruply excited configurations are added to the expansion, the k⁶ factor leads to a change in the slope of the scaling curve. Despite this, the largest calculation with about 4.9 × 10⁶ configurations required 205 s for calculating 100 samples with Sec. II D 2, whereas a linear extrapolation of the Sec. II D 1 times suggests that it would take more than 38.5 h for the same calculation even with aggressive screening. This analysis demonstrates the favorable scaling of Sec. II D 2, making local energy calculations with long HCI expansions feasible.

In Fig. 3, we show the effect of the Jastrow factor on the convergence of the multi-Slater ground state energy. The reference energy was obtained by performing a density matrix renormalization group calculation, which is very accurate for linear systems. SHCI calculations were performed with progressively smaller ϵ₁ values leading to progressively longer expansions. Even after including the perturbation theory correction, SHCI has significant difficulty in converging to the exact energy because of the large number of strongly correlated degrees of freedom in this active space. The ground state is dense in the canonical representation and has significant contributions from highly excited configurations, leading to the slow convergence of SHCI.

For the Jastrow multi-Slater wave function, we used leading configurations from an SHCI calculation performed with ϵ₁ = 5 × 10⁻⁵. Note that these configurations are different from those included in the SHCI wave functions used to examine the convergence of SHCI energies. We made this particular choice to make sure that wave functions with more configurations are variationally superior than those with fewer. This choice also allowed us to use the parameters optimized for the shorter expansions to be used as initial guesses for longer ones, which was crucial for the optimization. We used stochastic gradient descent with momentum⁴¹ to optimize the wave function energies. Jastrow parameters and CI coefficients were optimized together. We used 84 processes with 100 continuous time Monte Carlo (CTMC) samples and 10 burn-in samples each during the optimization and did a final energy calculation with 3000 samples per process to get the statistical error below 0.5 mH. The Jastrow multi-Slater energies can be seen to converge much more rapidly than bare SHCI. The energy error for the wave function consisting of about 1.0 × 10⁴ configurations is 3.5(4) mH, whereas the error in the HCI variational energy with 1.1 × 10⁷ configurations is 48.8 mH. Similar observations have been made for the Hubbard model using transcorrelated Hamiltonians.^42,43

We have presented an algorithm for calculating local energies of multi-Slater wave functions in orbital space Monte Carlo. We used a Jastrow multi-Slater wave function in VMC to demonstrate its efficiency. The algorithm itself is general and should allow the use of longer selected CI expansions as trial wave functions in other orbital space QMC methods like AFQMC and Green’s function Monte Carlo^44,45 as well. It also presents a different formulation of the real space algorithm proposed by Filippi et al. for local energies. We also showed the efficacy of Jastrow factors in vastly improving the convergence of the energy of CI expansions in polyacetylene.

We plan to explore various possibilities for improving the algorithm and using it in other Monte Carlo methods. As mentioned before, methods like density fitting and tensor hypercontraction can be used to reduce the cost of building the intermediates, which will likely dominate the total cost for large systems. We would also like to explore other ways of exploiting the structure in the Hamiltonian. Our pilot implementation uses dense tensor contractions even for the sparse Hamiltonian, which can be improved considerably. As for the Jastrow multi-Slater wave functions, we would like to study the criteria for choosing configurations in the presence of a Jastrow, as well as the possibility of using symmetry projection to improve the convergence of energy with the number of configurations.^46,47 This will be crucial for feasibly optimizing these wave functions because stochastic optimization with a large number of parameters is a difficult task. Our numerical experiments with polyacetylene suggest that an optimization strategy where CI coefficients are progressively optimized in blocks may be more effective than optimizing all of them together. Based on experiences in the community with other nonlinearly parameterized wave functions, it may be worth exploring such optimization techniques.

We thank Joonho Lee for useful discussions about AFQMC. The funding for this project was provided by the National Science Foundation through Grant No. CHE-1800584. S.S. was also partly supported through the Sloan research fellowship.

The data that support the findings of this study are available from the corresponding author upon reasonable request. The code used to generate the results can be found at https://github.com/sanshar/VMC/tree/master.

We start by deriving the first pair contraction identity in Eq. (8). Using the usual (vacuum) Wick’s theorem, we have

where A′ is obtained from A by replacing its ith row with the ath row of R, i.e., $A′i=Ra$⁠. Laplace expanding det(A′) along the ith row, we get

where C_A′ and C_A are the cofactor matrices of A′ and A, respectively. Since A′ and A only differ in their ith row, their cofactor matrices have the same ith rows, i.e., $CA′i=CAi$⁠. The cofactor matrix is related to the inverse by

Thus, we have

leading to the desired result,

As noted in the main text, the rest of the identities in Eq. (8) follow from this one by expressing the $Ĉ$ operators in terms of $L^$ operators.

In this section, we illustrate how Sec. II D 2, explained in the main text within the framework of VMC, can be used in AFQMC with some minor changes. We follow the importance sampled open ended random walk AFQMC procedure described by Motta and Zhang.⁷ While walkers in VMC are represented by different electron configurations in the same one particle orbital set, in AFQMC, walkers are nonorthogonal configurations with different orbitals.

The mixed estimator for ground state energy used in AFQMC is calculated as a weighted average of the local energies of the walkers,

where |ψ_T⟩ is the trial wave function, |ψ₀⟩ is the propagated state, |n⟩ is a walker configuration, $wn$ is the weight, and E_L[n] is the local energy of this walker with respect to the trial wave function,

Suppose the trial wave function used for calculating the mixed estimator is a multi-Slater wave function represented using a canonical orbital set. Green’s function is defined as

where $Ĉ$ denote the canonical basis operators. G is obtained by a simple extension to all orbital indices μ and ν of the second pair contraction identity in Eq. (8). The ab initio Hamiltonian can be written in the canonical basis as

We only focus on the two-body part of the Hamiltonian and analyze the quantity,

which is used to calculate the full local energy as in Eq. (18). Using the definition of |ϕ_I⟩ in Eq. (5) and the generalized Wick’s theorem, we get

Comparing this equation with Eq. (27), we again see the pattern that the determinants in this sum have the same k × k bottom-right corner. This fact allows the use of algorithm II, outlined in Sec. II D 2, by storing the intermediates described in Eqs. (27)–(33).

In order to optimize the wave function energy, we use gradient based methods. The ith component of the energy gradient is sampled according to the equation

where E is the energy of the wave function |ψ⟩ and |ψ_i⟩ is the derivative of the wave function with respect to the ith parameter. Since E_L[n] and E are both available from the energy sampling, we only need to calculate the wave function derivative overlaps to obtain an estimate of the energy gradient.

The Jastrow multi-Slater wave function has three types of parameters that can be optimized: Jastrow elements, CI coefficients, and orbital coefficients. We note that our parameterization is fairly redundant, and it is possible that removing these redundancies can lead to improvements in optimization. We have chosen not to do so because of the simplicity of the redundant parametrization. For Jastrow parameters, the derivative overlap ratio is given by

where n_i and n_j are the occupation numbers of orbitals i and j in the walker. All Jastrow derivative overlaps can be calculated at cost O(n²). For the CI coefficients, the derivate overlap ratio is given by

All CI coefficient derivative overlaps can be calculated at cost O(n_ck³), given the matrices stored during energy sampling. A naive calculation of the orbital coefficient derivatives has a prohibitive cost of O(n²n_ck³), but this can be improved to O(n³ + n_ck³) using the method outlined by Assaraf et al.²⁵ We have not used orbital optimization in this article.

We use continuous time Monte Carlo (CTMC)^48,49 to perform the energy and gradient sampling. We refer the reader to our previous work²⁹ for the details of this sampling technique. To perform efficient CTMC sampling of a wave function |ψ⟩, one needs to calculate the overlap ratios $⟨m|ψ⟩⟨n|ψ⟩$ for all the excitations |m⟩ generated from the walker |n⟩ through the Hamiltonian. In Sec. II D 1, all these overlap ratios are calculated directly and can be stored during the local energy calculation, allowing the CTMC sampling of the full wave function at essentially no extra cost. However, in Sec. II D 2, these overlap ratios are not calculated explicitly, instead they get absorbed in the intermediates. Thus, it is not possible to perform efficient CTMC sampling of the full wave function in Sec. II D 2. Instead, we sample the wave function given by

Since this sampling wave function only contains the reference configuration, its overlap ratios are available in Sec. II D 2. We estimate the energy of the full wave function by taking the ratio of the following quantities obtained by CTMC sampling |ψ₀⟩:

This introduces a bias in the estimate of the energy, which could be severe if |ψ₀⟩ has vanishing contributions from parts of the Hilbert space where |ψ⟩ has significant amplitudes.⁵⁰ We guard against instabilities due to rare events when the walker overlap with the sampling wave function is small by throwing away outlier samples using a threshold for the magnitude of $⟨ψ|n⟩⟨ψ0|n⟩$⁠. In our experiments, this has not lead to significant issues because |ψ₀⟩ is often a good approximation to |ψ⟩. We note that the occurrences of such rare events can be reduced by including more configurations in the sampling wave function, but we have not done so in this study. This will also improve the sampling efficiency at a slightly enhanced cost introduced by the direct calculation of overlap ratios. We are currently investigating possible sampling issues and will address them more thoroughly in a future article. The gradient can also be sampled similarly.