Making the most of data: Quantum Monte Carlo Post-Analysis Revisited

In quantum Monte Carlo (QMC) methods, energy estimators are calculated as the statistical average of the Markov chain sampling of energy estimator along with an associated statistical error. This error estimation is not straightforward and there are several choices of the error estimation methods. We evaluate the performance of three methods, Straatsma, an autoregressive model, and a blocking analysis based on von Neumann's ratio test for randomness, for the energy time-series given by Diffusion Monte Carlo, Full Configuration Interaction Quantum Monte Carlo and Coupled Cluster Monte Carlo methods. From these analyses we describe a hybrid analysis method which provides reliable error estimates for series of all lengths. Equally important is the estimation of the appropriate start point of the equilibrated phase, and two heuristic schemes are tested, establishing that MSER (mean squared error rule) gives reasonable and constant estimations independent of the length of time-series.


I. INTRODUCTION
2][3] QMC is one of the most accurate ab initio methods, and it is often used for systems which cannot be sufficiently accurately described by Density Functional Theory (DFT) [4][5][6][7] or which are too large to apply post Hartree-Fock methods. 8ffusion Monte Carlo (DMC) 9 is one such QMC method with a computational scaling of O(N 3 ) for a system of N electrons, and, as such, it can be applied even to large-sized systems including more than 1000 electrons. 5he drawback of DMC is the requirement to use the fixed-node approximation 10 to avoid the sign problem, which introduces a systematic error dependent upon the quality of the nodes of a trial wavefunction, and, although there are ways to suppress this error, 11,12 they make calculations considerably more expensive.
Two newer QMC methods in quantum chemistry have attracted interest of late, as they are not constrained by the fixed-node approximation: Firstly the fullconfiguration interaction QMC (FCIQMC) method, [13][14][15] which stochastically solves the equations of fullconfiguration interaction (FCI), by sampling with discrete particles.Although the scaling of calculation cost is still exponential 15 in the number of electrons like FCI, the prefactor of scaling curve is significantly reduced.Thus, this method can be applied to medium-sized systems. 16Secondly is Coupled Cluster Monte Carlo (CCMC), which stochastically solves Coupled Cluster (CC) equations. 17,18Since the parameter space of a truncated CC calculation is smaller than that of FCI, CCMC will in general have a smaller memory cost than FCIQMC.
QMC methods commonly provide an energy estimator as the statistical average of sampling a Markov chain, also producing an estimate of the statistical error.It is difficult to estimate the error reliably due to the following reasons: 19 (i) The samples are not independent of each other but correlated along the simulation time evolution.(ii) When the distribution of sampling is non-normal, the probability distribution of the mean value is also non-normal unless the number of sampling is large enough to satisfy the central limit theorem.In this work, we examine the performance of three characteristic automatic error estimation methods, Straatsma, 20 , the AutoRegressive (AR) model 21 , and blocking analysis based on von Neumann's ratio test for randomness 22,23 (von Neumann blocking) for the energy time-series obtained by applying DMC, FCIQMC, and CCMC to the neon atom.From these data we establish recommendations for the most reliable error estimation method for different lengths of time series, and devise a new hybrid scheme applicable to any length of time series.
Another important issue on the post-analysis of QMC is to determine the length of the pre-equilibration (warmup) phase. 19Underestimation of the length gives a systematic error in the energy but its overestimation also increases the statistical error.In this work, we tested two heuristic methods to estimate the warm-up steps.One is MSER (mean squared error rule) 24 and the other is min-WREE (minimization of weighted relative error of the error) inspired by the post-analysis implemented in HANDE QMC code for stochastic quantum chemistry. 25,26Our analysis establishes that the estimation of warm-up steps by min-WREE changes depending on the length of time-series.On the other hand, MSER makes reasonable and constant estimation of warm-up steps, independent of the length of time-series.

II. ERROR ESTIMATION METHOD
In this section we elucidate the error estimation methods used: Straatsma, 20 AR model, 21 and von Neumann blocking. 23

Straatsma
Straatsma et al. show that the variance of the statistical average of stationary time-series including n samples {X i } n i=1 is given as follows without any assumptions: 20 Here, τ is the estimation of time-correlation length (steps). 27c λ is the autocorrelation function with lag λ, which is approximately given by the finite number n of samples as: Here, the approximation of c λ is inaccurate when the number of terms n − λ to be summed up is small.Thus, in equation 2 we limit the summation over λ to values before c λ becomes negative for the first time, since later c λ oscillates around zero as shown in Figure 1.The resulting τ is then used to calculate σ 2 X .

AutoRegressive (AR) Model
The AR model assumes that the random process of Markov chain sampling {X i } n i=1 can be reasonably described by 21 The i-th sample is given as a linear combination of the previous steps and a random Gaussian noise a i with average 0 and variance σ 2 a .The coefficients {π i } p i=1 and the variance σ 2 a are fitted to the given time-series using Yule-Walker equation. 28The number of coefficients p is decided based on Akaike's Information Criterion (AIC). 29arge number of coefficients are needed to accurately describe the stochastic process of the given time-series, but, if it is too large, it becomes over-fitting.AIC aims to provide an appropriate compromise p value.
The estimation of the correlation length τ is calculated by where c λ are defined by equation 3.
von Neumann blocking Von Neumann blocking takes into account the nonnormality of the distribution of the sampling, in contrast to the two above-mentioned methods.The given timeseries {X i } n i=1 is divided into blocks with block size m, and a new time-series produced: Both the non-normality of the distribution and the correlation length are reduced in the new time-series The block size m is decided such that each W j (m) can be regarded to be sampled from independent and identical normal distributions, based on von Neumann's rate test for randomness. 22After blocking, the variance of the statistical average is given by σ Here, we also note the very commonly used blocking approach by Flyvbjerg and Petersen. 30This method performs integration of autocorrelation functions almost equivalent to equation 1 yet utilizing blocking, aiming to reduce the computational cost of analysis: It is mathematically similar to Straatsma 20 but it suffers a systematic bias stemming from blocking as well as von Neumann blocking.

III. ESTIMATION SCHEMES OF WARM-UP STEPS
We introduce two schemes for warm-up steps estimation in this section, MSER 24 and min-WREE, which are implemented in HANDE code. 25,26an squared error rule (MSER) MSER aims to give an adequate estimate of warm-up steps d as minimizing a sum of the systematic error from the warm-up phase and the statistical error.The number of warm-up steps, d, is determined by minimizing the following quantity: Here, s 2 X (d) is a constant independent of d corresponding to the case where {X i } n i=d does not include a warmup phase.When a warm-up phase is present, s 2 X (d) increases from the constant value according to how much warm-up phase remains in {X i } n i=d .Meanwhile 1/(n−d) monotonically increases according to d, and the value of d minimizing their product gives an appropriate estimate of the number of warm-up steps.

Minimization of weighted relative error of the error (min-WREE)
This scheme estimates the warm-up steps as minimizing the relative error of error of the statistical average, REE (d), weighted by 1/ √ n − d: We evaluated REE (d) using von Neumann blocking 23 in this work.

IV. QMC CALCULATION DETAILS
We employed CASINO 31 for DMC calculations.We used a Slater-Jastrow trial wave-function. 9The determinant is generated by the Hartree-Fock method using a STO-6G Gaussian basis set. 32The Jastrow factor consists of one-and two-body terms and includes 42 parameters in total.For the DMC calculations, the target population of walkers is set to be 1024 and the time step is 0.005 a.u.−1 .We also used the cusp correction scheme, 33 which replaces the shape of orbitals nearby ionic cores with Slater functions to satisfy the Kato cusp conditions. 34Each sample of the energy time-series is given by averaging the local energies 9 over all of the walkers for every QMC iteration.The influence of the population fluctuation and the population control 31 is not considered in this work.
We performed FCIQMC and CCMC calculations with HANDE. 25,26The reference Slater determinant is prepared by Hartree-Fock method with cc-pVDZ Gaussian basis set 35 using Psi4. 36The target number of walker population is 500 and the time step is 2.0×10 −5 a.u.−1 .Each sample of the energy time-series is given as the instantaneous projected energy, which is a ratio between D 0 Ĥ Ψ and N 0 ≡ D 0 |Ψ for every QMC iteration.The influence of the population fluctuation and the population control 31 is not considered in this work.

V. RESULTS AND DISCUSSION
We prepared one thousand different energy time-series for the neon atom, with the same calculation settings but with the different random seeds, using DMC, FCIQMC, and CCMC methods, for different lengths of time-series, respectively.We applied the error estimation methods to them and surveyed the concordance rate between the energy means and the reference mean value within the estimated errors with 1σ confidential interval (CI) as shown in Figure 2. The concordance rate for a 1σ CI is ideally 68.27%.Thus, when the observed rate is closer to this value, the error estimation is regarded to be more reliable.Here, the reference mean value is given by taking an average of very long length of time-series, and the error is just less than 3% of those of the energy means of 1000 time-series.First, we discuss the case of FCIQMC/CCMC (see Fig- ure 2bc).All of the error estimation methods give lower concordance rate than the ideal value for 1σ confidential interval, 68.27%, when the length of time-series n is small.This is typically observed for error estimation. 21or comparatively short lengths of time-series, the AR model shows the highest concordance rate among them.The comparison of the estimated errors shown in Figure 3bc further distinguishes the AR model from the others: Only the AR model reproduces that the estimated error normally decreases in proportion to 1/ √ n.It clearly proves the advantage of taking AR model of equation 4. On the other hand, the lowest concordance rate is measured for von Neumann blocking.The von Neumann's criteria to check randomness and normality tends to be not effective for small numbers of data, 22 so it underestimates the correlation length for small length of time-series.
Straatsma gives the intermediate concordance rates for small lengths n.The calculated τ fully depends on the autocorrelation function c λ through equations 1 and 2, so we examined how the shape of the autocorrelation function c λ changes according to the length n of timeseries in Figure 4: The autocorrelation function c λ becomes negative more quickly for smaller n by the c λ oscillating since it is estimated by insufficient number of terms, n − λ, through equation 3. The truncation of the sum in equation 2 is so drastic that Straatsma underestimates the correlation time.In contrast, the oscillation of c λ does not much affect the AR model, although its estimation also depends on autocorrelation functions c λ through equation 6.This is because taking a product with π λ drastically reduces the contribution of c λ with large lag λ: Figure 5 shows the expansion of the parameters π λ , where they are terminated or converged to zero only within a few terms.Therefore, just a few terms of c λ from small lag λ is used to calculate τ in AR model.
When the length of time-series n is comparatively large, the concordance rate of AR model converges to 68.27 % the most slowly.This is because the assumption of equation 4 in AR model cannot fully describe the target random process, and therefore its reliability is reduced.To summarize, the AR model (Straatsma)    Autocorrelation functions for different length n of time-series, given by CCMC method.For comparatively small n, the autocorrelation function c λ apparently includes a noise and it cannot be seen that c λ gradually converges to zero along with the lag λ increasing.
is the most reliable for small (middle/large) length timeseries, respectively.We have therefore devised a hybrid scheme of AR model and Straatsma, which works reasonably for any length of time-series.It simply adopts Parameters π λ in AR model fitted to different lengths of CCMC time-series.The number of the parameters is determined by AIC. 29 The parameters are terminated or converged around zero within a few terms, regardless of the length of time-series.
the larger of the errors estimated by both methods: σ X (hybrid) = max {σ X (AR model) , σ X (Straatsma)} .The concordance rate for the hybrid method is shown as 'AR model + Straatsma', in Figure 2  We performed the same test of the error estimation methods for DMC time-series.Von Neumann blocking gives the closest concordance rate to the ideal value, 68.27%, for any length of time-series, and only the concordance rate of this method reaches 68.27 %.The difference from the case of FCIQMC/CCMC comes from that the distribution of DMC energy time-series tends to be non-normal: Figure 6 clearly shows that the distributions of the time-series given by DMC have larger skewness and kurtosis than those of FCIQMC/CCMC.As mentioned in section II, only von Neumann blocking can take into account non-normality for error estimation, which would be the reason why von Neumann blocking the most works in the case of DMC.This difference in non-normality also explains why the AR model gives the lowest concordance rate: the AR model assumes that the randomness between the neighboring steps is expressed by normally distributed noise, so it would not be possible to make a description when the distribution of time-series is non-normal.
Finally, we discuss the performance of the estimation schemes of warm-up steps, MSER and min-WREE.We apply these schemes to the time-series including nonplateau part and removed the estimated warm-up steps d.Then, we applied 'von Neumann blocking' to obtain the concordance rate and the statistical error.Figure 2 shows that MSER basically gives higher concordance rate and smaller statistical error especially in the case of DMC: MSER is superior to min-WREE on the whole.The advantage is further distinguished comparing the average of warm-up steps.They are shown in Figure 7 with the standard errors.It clearly shows that the estimations of min-WREE strongly depends on the length of time-series and largely scattered.On the other hand, MSER estimates constant warm-up steps with small variances, independent of the length of time-series.

FIG. 1 .
FIG.1.A typical autocorrelation function of energy timeseries generated by CCMC calculation.It rapidly decreases and oscillates around zero as the lag, λ, increases, showing that the correlation between Xi and X i+λ decreases.

FIG. 2 .
FIG.2.The concordance rates between the energy means and the reference mean value within the errors with 1σ confidential interval estimated by the error estimation methods.The concordance rates are surveyed for the time-series generated by (a)DMC, (b)FCIQMC, and (c)CCMC, for different lengths of time-series.The black horizontal line shows 68.27 %, which is the ideal value for 1σ confidence interval.When the measured concordance rate is closer to this value, the error estimation is regarded to be more reliable.

FIG. 3 .
FIG. 3. The statistical errors estimated by the error estimation methods for time-series given by (a)DMC, (b)FCIQMC, and (c)CCMC for different lengths of time-series.
FIG. 4.Autocorrelation functions for different length n of time-series, given by CCMC method.For comparatively small n, the autocorrelation function c λ apparently includes a noise and it cannot be seen that c λ gradually converges to zero along with the lag λ increasing.
FIG. 5.Parameters π λ in AR model fitted to different lengths of CCMC time-series.The number of the parameters is determined by AIC.29 The parameters are terminated or converged around zero within a few terms, regardless of the length of time-series.

FIG. 7 .
FIG. 7.The warm-up steps estimated by MSER and min-WREE for the time-series generated by DMC, FCIQMC, and CCMC.The error bars correspond to the statistical errors.
and always com-paratively close to 68.27 %.