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

TL;DR

This paper evaluates and compares different error estimation methods for energy estimators in quantum Monte Carlo simulations, proposing a hybrid approach for reliable error analysis across various series lengths.

Contribution

It introduces a hybrid error estimation method and heuristic schemes for identifying the equilibrated phase in QMC energy series analysis.

Findings

Hybrid analysis provides reliable error estimates for all series lengths.

MSER heuristic effectively identifies the start of the equilibrated phase.

Comparison of three error estimation methods demonstrates their strengths and limitations.

Abstract

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.

Equations

$${\sigma }_{\bar{X}}^2$$

$$\tau$$

$$c_{\lambda }\approx \frac{1}{n-\lambda }\sum _{i=1}^{n-\lambda }\left\{\left(X_i-\bar{X}\right)\left(X_{i+\lambda }-\bar{X}\right)\right\},(\lambda <n).$$

$$X_i$$

$$\bar{X}$$

$$\tau =\left(1-\sum _{\lambda =1}^p{c}_{\lambda }{\pi }_{\lambda }\right)/{\left(1-\sum _{\lambda =1}^p{\pi }_{\lambda }\right)}^2,$$

$$W_j(m)=\frac{1}{m}\sum _{l=1}^m{X}_{m\left(j-1\right)+l}.$$

$$\mathrm{MSER}\left(d\right)=\frac{s_X^2\left(d\right)}{n-d},s_X^2\left(d\right)\equiv \frac{1}{n-d}\sum _{i=1}^{n-d}{\left(X_{i+d}-\bar{X}\right)}^2$$

$$\text{WREE}\left(d\right)=\text{REE}\left(d\right)/\sqrt{n-d}$$

Full text

I Introduction

With the increase in availability of large-scale computers, Quantum Monte Carlo (QMC) methods have spread rapidly owing to the embarassing parallelizability of such algorithms. 1; 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.8

Diffusion Monte Carlo (DMC)9 is one such QMC method with a computational scaling of $\mathcal{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. 5 The drawback of DMC is the requirement to use the fixed-node approximation10 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 full-configuration interaction QMC (FCIQMC) method,13; 14; 15 which stochastically solves the equations of full-configuration interaction (FCI), by sampling with discrete particles. Although the scaling of calculation cost is still exponential15 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.16 Secondly is Coupled Cluster Monte Carlo (CCMC), which stochastically solves Coupled Cluster (CC) equations.17; 18 Since 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 (warm-up) phase.19 Underestimation 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; 26 Our 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 $\left\{{{X_{i}}}\right\}_{i=1}^{n}$ is given as follows without any assumptions:20

[TABLE]

Here, $\tau$ is the estimation of time-correlation length (steps).000 The correlation length $\tau$ is differently defined in two papers, 20; 21 but they give the same definition of error. We took the newer definition of $\tau$.21

${{c_{\lambda}}}$ is the autocorrelation function with lag $\lambda$, which is approximately given by the finite number $n$ of samples as:

[TABLE]

Here, the approximation of $c_{\lambda}$ is inaccurate when the number of terms $n-\lambda$ to be summed up is small. Thus, in equation 2 we limit the summation over $\lambda$ to values before $c_{\lambda}$ becomes negative for the first time, since later $c_{\lambda}$ oscillates around zero as shown in Figure 1. The resulting $\tau$ is then used to calculate $\sigma_{\bar{X}}^{2}$.

AutoRegressive (AR) Model

The AR model assumes that the random process of Markov chain sampling $\left\{{{X_{i}}}\right\}_{i=1}^{n}$ can be reasonably described by 21

[TABLE]

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 $\sigma_{a}^{2}$. The coefficients $\left\{{{\pi_{i}}}\right\}_{i=1}^{p}$ and the variance $\sigma_{a}^{2}$ are fitted to the given time-series using Yule–Walker equation.28 The number of coefficients $p$ is decided based on Akaike’s Information Criterion (AIC). 29 Large 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 $\tau$ is calculated by

[TABLE]

where $c_{\lambda}$ are defined by equation 3.

von Neumann blocking

Von Neumann blocking takes into account the non-normality of the distribution of the sampling, in contrast to the two above-mentioned methods. The given time-series $\left\{{{X_{i}}}\right\}_{i=1}^{n}$ is divided into blocks with block size $m$, and a new time-series produced:

[TABLE]

Both the non-normality of the distribution and the correlation length are reduced in the new time-series $\left\{{{W_{j}}}\right\}_{i=1}^{n}$. 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.22 After blocking, the variance of the statistical average is given by $\sigma_{\bar{X}}^{2}=\frac{1}{{k\left({k-1}\right)}}\sum\limits_{j=1}^{k}{{{\left[{{W_{j}}\left(m\right)-\bar{W}}\right]}^{2}}}$, $\left({\bar{W}\equiv\bar{X}}\right)$.

Here, we also note the very commonly used blocking approach by Flyvbjerg and Petersen.30 This 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, MSER24 and min-WREE, which are implemented in HANDE code.25; 26

Mean 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:

[TABLE]

Here, ${s_{X}^{2}\left(d\right)}$ is a constant independent of $d$ corresponding to the case where $\left\{{{X_{i}}}\right\}_{i=d}^{n}$ does not include a warm-up phase. When a warm-up phase is present, ${s_{X}^{2}\left(d\right)}$ increases from the constant value according to how much warm-up phase remains in $\left\{{{X_{i}}}\right\}_{i=d}^{n}$. 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, ${\text{REE}}\left(d\right)$, weighted by $1/\sqrt{n-d}$:

[TABLE]

We evaluated ${\text{REE}}\left(d\right)$ using von Neumann blocking 23 in this work.

IV QMC calculation details

We employed CASINO31 for DMC calculations. We used a Slater-Jastrow trial wave-function.9 The determinant is generated by the Hartree–Fock method using a STO-6G Gaussian basis set.32 The 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.34 Each sample of the energy time-series is given by averaging the local energies9 over all of the walkers for every QMC iteration. The influence of the population fluctuation and the population control31 is not considered in this work.

We performed FCIQMC and CCMC calculations with HANDE.25; 26 The reference Slater determinant is prepared by Hartree–Fock method with cc-pVDZ Gaussian basis set35 using Psi4.36 The target number of walker population is 500 and the time step is 2.0$\times 10^{-5}$ a.u.-1. Each sample of the energy time-series is given as the instantaneous projected energy, which is a ratio between $\left\langle{{{\rm{D}}_{0}}\left|{\hat{H}}\right|{\Psi}}\right\rangle$ and ${N_{0}}\equiv\left\langle{{{\rm{D}}_{0}}\left|{\Psi}\right.}\right\rangle$ 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\sigma$ confidential interval (CI) as shown in Figure 2. The concordance rate for a 1$\sigma$ 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 Figure 2bc). All of the error estimation methods give lower concordance rate than the ideal value for $1\sigma$ confidential interval, 68.27%, when the length of time-series $n$ is small. This is typically observed for error estimation.21 For 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/\sqrt{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 $\tau$ fully depends on the autocorrelation function $c_{\lambda}$ through equations 1 and 2, so we examined how the shape of the autocorrelation function $c_{\lambda}$ changes according to the length $n$ of time-series in Figure 4: The autocorrelation function $c_{\lambda}$ becomes negative more quickly for smaller $n$ by the $c_{\lambda}$ oscillating since it is estimated by insufficient number of terms, $n-\lambda$, 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_{\lambda}$ does not much affect the AR model, although its estimation also depends on autocorrelation functions $c_{\lambda}$ through equation 6. This is because taking a product with ${\pi_{\lambda}}$ drastically reduces the contribution of $c_{\lambda}$ with large lag $\lambda$:  Figure 5 shows the expansion of the parameters ${{\pi_{\lambda}}}$, where they are terminated or converged to zero only within a few terms. Therefore, just a few terms of $c_{\lambda}$ from small lag $\lambda$ is used to calculate $\tau$ 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) is the most reliable for small (middle/large) length time-series, 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 the larger of the errors estimated by both methods: ${\sigma_{\bar{X}}}\left({{\rm{hybrid}}}\right)=\max\left\{{{\sigma_{\bar{X}}}\left({{\rm{AR\,model}}}\right),{\sigma_{\bar{X}}}\left({{\rm{Straatsma}}}\right)}\right\}.$

The concordance rate for the hybrid method is shown as ‘AR model + Straatsma’, in Figure 2 and always comparatively close to 68.27 %.

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 non-plateau 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.

VI Conclusion

We compared the reliability of three kinds of error estimation methods, Straatsma,20 AR model,21 and von Neumann blocking,22; 23 in a statistical manner, when they are applied to the energy time-series given by applying DMC9, FCIQMC13, and CCMC17; 18 to the neon atom. In the case of FCIQMC/CCMC, it is shown that Straatsma (the AR model) is the most reliable for comparatively long (short) length of time-series, respectively. We established that the assumption in the AR model significantly reduces the influence of the oscillation of autocorrelation functions for short time-series lengths, but we concluded that the systematic error from the assumption is pronounced for long lengths. We devised a hybrid scheme, which takes the larger error from the ones estimated by Straatsma and AR model, and established it works for any length of time-series. In contrast, for DMC, we showed that von Neumann blocking is the most reliable for any length of time-series. This method has an advantage of considering the non-normality of the distribution of time-series, and we showed an strong evidence that the advantage is essential to analyse the DMC results. We have also tested two kinds of warm-up steps estimation schemes, MSER24 and min-WREE25, and established that MSER gives constant and reasonable estimations of warm-up steps, independent of the length and for time-series made by any of the QMC methods.

Acknowledgments

The computation in this work has been performed using the facilities of the Research Center for Advanced Computing Infrastructure (RCACI) at JAIST. T.I. is grateful for financial suport from Grant-in-Aid for JSPS Research Fellow (18J12653). K.H. is grateful for financial support from a KAKENHI grant (JP17K17762), a Grant-in-Aid for Scientific Research on Innovative Areas “Mixed Anion” project (JP16H06439) from MEXT, PRESTO (JPMJPR16NA) and the Materials research by Information Integration Initiative (MI2I) project of the Support Program for Starting Up Innovation Hub from Japan Science and Technology Agency (JST). R.M. is grateful for financial supports from MEXT-KAKENHI (17H05478 and 16KK0097), from Toyota Motor Corporation, from I-O DATA Foundation, and from the Air Force Office of Scientific Research (AFOSR-AOARD/FA2386-17-1-4049). R.M. and K.H. are also grateful to financial supports from MEXT-FLAGSHIP2020 (hp170269, hp170220). A.J.W.T. thanks the Royal Society for a University Research Fellowship (UF110161).