Unbiased Monte Carlo Cluster Updates withAutoregressive Neural Networks

Dian Wudian.wu@epfl.ch

Riccardo Rossi
riccardo.rossi@epfl.ch

Giuseppe Carleo∗
giuseppe.carleo@epfl.ch

Abstract
Efﬁcient sampling of complex high-dimensional probability distributions is a central task in computational science. Machine learning methods like autoregressive
neural networks, used with Markov chain Monte Carlo sampling, provide good
approximations to such distributions, but suffer from either intrinsic bias or highvariance. In this work, we propose a way to make this approximation unbiased
and with low variance. Our method uses physical symmetries and variable-sizecluster updates which utilize the structure of autoregressive factorization. We test
our method for ﬁrst- and second-order phase transitions of classical spin systems,
showing its viability for critical systems and in the presence of metastable states.

1

Introduction

Markov chain Monte Carlo [1] (MCMC) is an unbiased numerical method that allows samplingfrom unnormalized probability distributions, a central task in many areas of computational science.
MCMC is commonly used, for example, in molecular dynamics [2], as well as statistical and quantum physics [3–6]. In addition to fundamental applications, MCMC serves as a physics-inspired
approach to solve a variety of computational problems, including combinatorial optimization [7, 8]
and computer graphics [9]. While MCMC is a generically applicable technique, its implementation
can be plagued by long mixing or autocorrelation time [10]. Various techniques have been proposed
to increase the efﬁciency of MCMC [11], for example, cluster updates [12, 13], parallel tempering [14], the worm algorithm [15], and event-chain Monte Carlo [16]. However, these faster MCMC
algorithms rely on details of the physical system considered, and they cannot be applied generically.
Machine learning (ML) methods, given their intrinsic ﬂexibility in addressing problems in computational physics [17], are being intensively investigated as a way to improve MCMC. Applications inthis direction include, for example, self-learning Monte Carlo methods [18–23], enhanced sampling
driven by neural networks [24, 25], and neural importance sampling [26]. Strongly rooted in theprinciples of statistical physics, variational sampling techniques are among the most promising MLdriven approaches. Generative neural samplers (GNS) [27–29] are a chief example of ML-driven
variational methods. These approaches build on the idea of constructing approximate representations of the original probability distribution at hand. The resulting variational approximations can
efﬁciently perform sampling by construction, thus completely bypassing MCMC. A particularly interesting aspect of this approach is its systematic improvability when using the free energy boundminimization as the guiding principle to gauge the approximation accuracy. The main drawback of
the variational approach, however, is that the estimators of expectation values are intrinsically biasedby the representation error of the approximated distribution. As unbiased estimators are of central
importance in many fundamental applications in physics, recent research has started addressing the
key problem of removing the bias induced by ML variational representations, for example, throughimportance sampling, and incorporating again MCMC strategies [26, 30, 31].
∗

Institute of Physics, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

Fourth Workshop on Machine Learning and the Physical Sciences (NeurIPS 2021).

1.1

Bias in variationally trained generative neural samplers

Consider a system of V classical Ising spins s := (s1 , . . . , sV ), si ∈ {−1, 1}, at inverse temperature β. We use a generative neural sampler (GNS) qθ with parameters θ that variationally approximates the Boltzmann probability distribution p(s) ∝ p̃(s) := e−βE(s) by minimizing a freeenergy upper bound [27], which is equivalent to minimizing the Kullback–Leibler (KL) divergence∑
θ (s)DKL (qθ ∥ p) := s qθ (s) ln qp(s). To construct an expressive qθ , we use an autoregressive neural∏V
network (ARNN) to factorize it into a product of conditional probabilities qθ (s) =: i=1 qθ;i (si |s<i ), where s<i := (s1 , . . . , si−1 ), which also allows us to efﬁciently sample from qθ by samplingfrom {qθ;i } sequentially.
However, the fact that the two distributions are only approximately equal, qθ (s) ≈ p̃(s), implies thatthe samples {s(1) , . . . , s(N ) } drawn from the network carry an intrinsic bias. When these samples
are used to compute the expectation value of a physical observable, the resulting estimator Ō =∑N
1(i)i=1 O(s ) is biased from the true value Ep [O] and it is impossible in general to reliablyNestimate the direction and the magnitude of such bias.
1.2

Neural importance sampling and global updates

Refs. [30, 31] have proposed two closely related solutions to the bias problem. The ﬁrst method,which we denote neural importance sampling (NIS) in the following, consists of using the modi(i)∑Nﬁed unbiased estimator Ō = i=1 w(s(i) )O(s(i) ), where w(s(i) ) := ∑Nw̃(sw̃(s)(j) ) and w̃(s(i) ) :=j=1p̃(s(i) )are the normalized and the unnormalized weights respectively. The second proposed soluqθ (s(i) )tion, which we denote neural global updates (NGU) hereafter, consists of using the GNS as a Markov
chain Monte Carlo (MCMC) proposer: if s is the current Markov chain state, a proposed( state s′ )is′)drawn from the GNS and accepted with the Metropolis probability Pacc (s → s′ ) := min 1, w̃(sw̃(s) .
However, as shown in our numerical experiments, those methods can be plagued by the ergodicity
issue, which produces impractically long autocorrelation time and high variance of the estimators,
due to the generic presence of “exponentially suppressed conﬁgurations”.

2

Methods

2.1

Exponentially suppressed conﬁgurations

We point out an elementary property of the KL divergence: the cost of allowing a single bad approximation qθ (s) scales only logarithmically with the ratio qθ (s)/p(s). Therefore, even when thefree energy is well approximated after the variational training, qθ (s) is still exponentially smaller
than p(s) for a small portion pESC of conﬁgurations, which we call exponentially suppressed conﬁgurations (ESC). A well-trained network has pESC ≪ 1, and those ESC have a limited effect on the
training loss, but rather a strong effect on the autocorrelation time.
Let us consider a Markov chain evolution using NGU, and suppose that the current state s is an
ESC. The ratio w̃(s′ )/w̃(s) in Pacc will be exponentially small for almost any other conﬁgurations′ ; therefore, the Markov chain will be essentially stuck in s for a long time before accepting any
new proposal, and the autocorrelation time of the whole chain will be impractically large. A similarargument applies when considering the variance of the NIS method.
2.2

Neural cluster updates with symmetries

To solve the generic ergodicity problem of neural global update methods, we propose an enhanced
MCMC method with the use of the physical symmetries and the autoregressive structure, which wedenote as neural cluster updates with symmetries (NCUS), as described in Alg. 1 for a 2D latticewith translational and reﬂectional symmetries, and sketched in Fig. 1.
Both the random symmetry operations and the cluster updates help escape from ESC. Assume thatthe current conﬁguration s is an ESC, and we use a random symmetry operation to change s toanother conﬁguration s∗ in the equivalence class C. The probability that all conﬁgurations in C are
2

ESC is on the order of p#C
ESC , so it is extremely unlikely to get stuck within the whole equivalenceclass. Note that the occurrence of ESC does not depend on the physical symmetries, which can bechecked in numerical experiments, but rather on the structure of the network.
′qθ;i (si |s<i ))p̃(s′ ) ∏V
In a cluster update, the weight ratio in Pacc becomes w̃(si=V −k+1 qθ;i (s′i |s′<i ) , whichw̃(s) = p̃(s)is not too far from 1 when k is small. Therefore, the new conﬁguration is closer to the old one
and is easier to be accepted. Although the cluster size k can be sampled from an arbitrary distribution Pcluster (k), numerical experiments have shown that the uniform distribution Pcluster (k) ≡ 1/V
already works better than many other cases we have explored.
Algorithm 1 A step of NCUS on a 2D lattice with Tx × Ty × D4 × Z2 symmetry.
1: Input the current conﬁguration s
2: Sample an integer k ∈ {1, . . . , V } from Pcluster
3: Sample the last k spins and propose the conﬁguration s′
4: Accept s ← s′ with probability Pacc (s → s′ )5: Translate s by a random displacement
6: Reﬂect s along the x axis, the y axis and the diagonal, each with 50% probability
7: Reﬂect s along the z axis (ﬂip all spins) with 50% probability
8: Output s as a sample in the Markov chain

Figure 1: Example of three steps of NCUS applied to a 4 × 4 spin model. The columns correspondto different lines in Alg. 1. The last k spins that can be ﬂipped are highlighted in blue. If a proposal
is accepted, the spins actually ﬂipped are shown in red. For translations, the original borders of thelattice are shown in green. For reﬂections, the plane of reﬂection is shown in yellow, and yellow
borders around the lattice indicate a reﬂection along the z axis (across the xy plane).

3

Numerical experiments

3.1

Ising model

We start to demonstrate the effectiveness of NCUS on the conventional 2D Ising model E(s) :=∑Li,j=1 si,j (si+1,j + si,j+1 ), with periodic boundary conditions sL+1,j = s1,j , si,L+1 = si,1 .
We use a lightweight convolutional ARNN with only 3 convolutional layers and approximately
4 × 103 parameters, and use the same network to compare the different sampling methods. Thanks
to the MCMC bias removal, we do not need the network to approximate the true distribution toan extremely high precision, which will be increasingly difﬁcult for larger lattices. Details of the
network structure, training, and sampling are described in Appendix A.
From Fig. 2 (a), we see that both NGU and NIS have pathologically high autocorrelation times
in the critical region. An inspection of their autocorrelation function in Fig. 2 (b) shows that the3

1.0

2 × 105
NCUS
NCU
NGUS
NGU
NIS
Wolff

1 × 105

105

(b)
1.0

103

0.5

101

r(t)

(a)

0.40

0.45

0.50

0.6

0.8

1.0

0.5

0.0

NCUS
NCU
NGUS
NGU

100

102

104

0.0

0
0.2

0.4

0

2 × 105 4 × 105 6 × 105 8 × 105 1 × 106
t

Figure 2: (a) Integrated autocorrelation time τ as a function of temperature on the 16 × 16 Isingmodel. For NIS, we use the increased variance from the reweighting procedure as the effective autocorrelation time. The inset focuses on their behaviors near the critical point and uses the logarithmic
scale on the y axis. (b) Autocorrelation functions r(t) on the 16 × 16 Ising model at β = 0.44. The
inset uses the logarithmic scale on the x axis to focus on their behaviors at small t.
Markov chain of NGU is essentially non-ergodic in the available simulation time. By contrast,
our proposed method NCUS has no issue in the critical region. A closer inspection of the inset
of Fig. 2 (a) shows that the autocorrelation time of NCUS still increases in the critical region, andthe sampling efﬁciency is improved typically by two orders of magnitude compared with the global
update methods. The performance of NCUS is also comparable to the celebrated Wolff cluster
update method [13], which is speciﬁcally tailored for the Ising model. The ablation studies of neuralcluster updates without symmetries (NCU) and neural global updates with symmetries (NGUS) arealso shown, and the uncertainties of τ are analyzed in Appendix B.
3.2

Frustrated plaquette model

We now study another model that presents a richer physics than the Ising model, and for which,to our knowledge, no traditional cluster update method is applicable. We consider a classical spin-1/2 system with nearest-neighbor J1 , next-next-nearest-neighbor J3 , and plaquette K∑L∑Linteractions E(s) :=J1 i,j=1 si,j (si+1,j + si,j+1 ) + J3 i,j=1 si,j (si+2,j + si,j+2 ) +∑L
K i,j=1 si,j si+1,j si,j+1 si+1,j+1 , with periodic boundary conditions, which we denote as thefrustrated plaquette model (FPM). In this work, we set J1 = J3 = −1, K = 2, where a ﬁrst-ordertransition between the ferrimagnetic (fM) and the paramagnetic (PM) phases is conjectured from ananalogy to the q = 8 Potts model [32, 33]. The comparison of autocorrelation times from different
sampling methods is presented in Fig. 3, which provides numerical evidence that NCUS greatly
alleviates the metastability issue expected near ﬁrst-order phase transitions [34].
4 × 105 (a)

(b)
1.0

1.0
105
NCUS
NCU
NGUS
NGU
NIS

2 × 105
1 × 105

0.5

104
r(t)

3 × 105

NCUS
NCU
NGUS
NGU

103
0.20

0.21

0.22

0.5

0.0

0.23

100

102

104

0.0

0
0.15

0.20

0.25

0.30

0.35

0

2 × 105 4 × 105 6 × 105 8 × 105 1 × 106
t

Figure 3: (a) Integrated autocorrelation time τ as a function of temperature on the 32 × 32 FPM. (b)Autocorrelation functions r(t) on the 32 × 32 FPM at β = 0.2.
4

Theoretically, a ﬁrst-order phase transition occurs when the distribution of energy p(E) ∝N (E) e−βE has two peaks with the same size, as shown in Fig. 4, where N (E) is the numberof conﬁgurations with energy E. A GNS-based sampling method has equal probabilities to generatea sample from the two peaks, and the probability to accept that proposal will be close to 1, if the
network is ideally trained and there is no problem of ESC. Meanwhile, for traditional local-update
MCMC methods, they can only move small horizontal steps in Fig. 4, so it takes more steps (∼ L2 )2and exponentially lower probability (∼ e−βδEL ) for them to walk from the low-energy peak to
the high-energy one, where δE is the typical energy difference in a local update, which does notscale with L. In other words, the exponentially large number of conﬁgurations in the high-energy
peak will not make it easier for local-update MCMC methods to sample from that peak because
it is exponentially hard for the walker to walk between those conﬁgurations in locally connected
paths. NCUS reaches a balance between the two extremes, which solves the problem of ESC andkeeps the autocorrelation time practically low, even if the network is lightweight and cannot ideallyapproximate the true distribution.
L = 16
L = 24
L = 32

0.8

p

0.7
0.6
0.5
3.5

3.0

2.5

2.0

Figure 4: Probability distribution of the energy per site p(ε) for the FPM with different lattice sizes
L at their respective phase transition temperatures, obtained by NCUS.

4

Conclusions

In this work, we have shown a method to systematically remove the bias and keep a low variancein GNS-based sampling, using symmetries of the physical system and cluster updates of the autoregressive factorization, which is viable for critical systems and in the presence of metastable states.
While we have been mainly concerned with the metric of autocorrelation time, we recognize thatthe wall-clock time is another important metric for practical computations. In this respect, when
computing the energy of the system has a negligible computational cost, current neural networkbased methods are not yet competitive with traditional MCMC methods. It can then be argued thatthe ideal application scenario for ML-based methods are those cases where evaluating the integrandis expensive, for example, in determinant quantum Monte Carlo [35] and lattice ﬁeld theory [36].
In future work, computational efﬁciency can be addressed on multiple fronts, for example, by introducing techniques such as hierarchy and sparsity of the neural network models, to reduce thecomputation time and scale up the lattice size by orders of magnitude. After that, we expect that theslow asymptotic growth of the autocorrelation time of GNS will eventually make them outperform
traditional MCMC methods in terms of wall-clock time.

Acknowledgments and Disclosure of Funding
We acknowledge insightful comments and suggestions from Lei Wang and Pan Zhang. Support
from the Swiss National Science Foundation is acknowledged under Grant No. 200021_200336.
The computing power is supported with Cloud TPUs from Google’s TensorFlow Research Cloud(TFRC).

5

References[1] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys.
21, 1087 (1953).
[2] K. Binder, Monte Carlo and Molecular Dynamics Simulations in Polymer Science (Oxford
University Press, 1995).
[3] K. Binder and D. W. Heermann, Monte Carlo Simulation in Statistical Physics (Springer,2010).
[4] W. Krauth, Statistical Mechanics: Algorithms and Computations (Oxford University Press,
2006).
[5] J. Gubernatis, N. Kawashima, and P. Werner, Quantum Monte Carlo Methods (Cambridge
University Press, 2016).
[6] F. Becca and S. Sorella, Quantum Monte Carlo Approaches for Correlated Systems (Cambridge
University Press, 2017).
[7] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Science 220, 671 (1983).
[8] R. Y. Rubinstein and D. P. Kroese, The Cross-Entropy Method: A Uniﬁed Approach to Combinatorial Optimization, Monte-Carlo Simulation and Machine Learning (Springer, 2004).
[9] R. L. Cook, ACM Trans. Graph. 5, 51 (1986).
[10] H. Müller-Krumbhaar and K. Binder, J. Stat. Phys. 8, 1 (1973).
[11] J. S. Liu, Monte Carlo Strategies in Scientiﬁc Computing (Springer, 2004).
[12] J.-S. Wang and R. H. Swendsen, Physica A 167, 565 (1990).
[13] U. Wolff, Phys. Rev. Lett. 62, 361 (1989).
[14] R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 57, 2607 (1986).
[15] N. V. Prokof’ev, B. V. Svistunov, and I. S. Tupitsyn, Phys. Lett. A 238, 253 (1998).
[16] E. P. Bernard, W. Krauth, and D. B. Wilson, Phys. Rev. E 80, 056704 (2009).
[17] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, andL. Zdeborová, Rev. Mod. Phys. 91, 045002 (2019).
[18] D. Levy, M. D. Hoffman, and J. Sohl-Dickstein, in International Conference on Learning
Representations (2018).
[19] J. Song, S. Zhao, and S. Ermon, in Advances in Neural Information Processing Systems, Long
Beach, CA (2017) pp. 5140–5150.
[20] M. Medvidovic, J. Carrasquilla, L. E. Hayward, and B. Kulchytskyy,
arXiv:2012.01442 .

(2020),[21] J. Liu, Y. Qi, Z. Y. Meng, and L. Fu, Phys. Rev. B 95, 041101 (2017).
[22] L. Huang and L. Wang, Phys. Rev. B 95, 035105 (2017).
[23] H. Shen, J. Liu, and L. Fu, Phys. Rev. B 97, 205140 (2018).
[24] L. Bonati, Y.-Y. Zhang, and M. Parrinello, Proc. Natl. Acad. Sci. 116, 17641 (2019).
[25] F. Noé, S. Olsson, J. Köhler, and H. Wu, Science 365, 10.1126/science.aaw1147 (2019).
[26] T. Müller, B. McWilliams, F. Rousselle, M. Gross, and J. Novák, ACM Trans. Graph. 38, 1(2019).
[27] D. Wu, L. Wang, and P. Zhang, Phys. Rev. Lett. 122, 080602 (2019).
6

[28] M. S. Albergo, G. Kanwar, and P. E. Shanahan, Phys. Rev. D 100, 034515 (2019).
[29] S.-H. Li and L. Wang, Phys. Rev. Lett. 121, 260601 (2018).
[30] K. A. Nicoli, S. Nakajima, N. Strodthoff, W. Samek, K.-R. Müller, and P. Kessel, Phys. Rev. E
101, 023304 (2020).
[31] B. McNaughton, M. V. Milošević, A. Perali, and S. Pilati, Phys. Rev. E 101, 053312 (2020).
[32] H. Arisue and K. Tabata, in Non-Perturbative Methods and Lattice QCD (World Scientiﬁc,2001) pp. 233–241.
[33] V. Gorbenko, S. Rychkov, and B. Zan, J. High Energy Phys. 2018 (10), 1.
[34] R. Gheissari and E. Lubetzky, (2016), arXiv:1607.02182 .
[35] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
[36] D. C. Hackett, C.-C. Hsieh, M. S. Albergo, D. Boyda, J.-W. Chen, K.-F. Chen, K. Cranmer,
G. Kanwar, and P. E. Shanahan, (2021), arxiv:2107.00734 .
[37] A. van den Oord, N. Kalchbrenner, and K. Kavukcuoglu, in International Conference on Machine Learning (2016).
[38] S. Elfwing, E. Uchibe, and K. Doya, Neural Netw. 107, 3 (2018).
[39] F. Yu and V. Koltun, in International Conference on Learning Representations (2016).
[40] D. P. Kingma and J. Ba, in International Conference on Machine Learning (2015).
[41] A. Gelman and D. B. Rubin, Stat. Sci. 7, 457 (1992).

7

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reﬂect the paper’s
contributions and scope? [Yes](b) Did you describe the limitations of your work? [Yes](c) Did you discuss any potential negative societal impacts of your work? [No] This israther fundamental physics.
(d) Have you read the ethics review guidelines and ensured that your paper conforms tothem? [Yes]2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [N/A](b) Did you include complete proofs of all theoretical results? [N/A]
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes](b) Did you specify all the training details (e.g., data splits, hyperparameters, how they
were chosen)? [Yes](c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [Yes](d) Did you include the total amount of compute and the type of resources used (e.g., typeof GPUs, internal cluster, or cloud provider)? [Yes]4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [N/A](b) Did you mention the license of the assets? [N/A](c) Did you include any new assets either in the supplemental material or as a URL? [N/A](d) Did you discuss whether and how consent was obtained from people whose data
you’re using/curating? [N/A](e) Did you discuss whether the data you are using/curating contains personally identiﬁable information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if
applicable? [N/A](b) Did you describe any potential participant risks, with links to Institutional Review
Board (IRB) approvals, if applicable? [N/A](c) Did you include the estimated hourly wage paid to participants and the total amountspent on participant compensation? [N/A]8

A

Details of numerical experiments

Our network has 3 convolutional layers, each with kernel size 5. The convolutions are masked to
implement the autoregressive property, as introduced in PixelCNN [37]. The numbers of input,
hidden, and output channels are 1 → 16 → 16 → 1. SiLU activations [38] are applied after the ﬁrst
and the second convolutional layers, which are reported to produce lower loss than ReLU. Sigmoid
activation is applied after the third convolutional layer to restrain the output into (0, 1).
We keep the network to be lightweight, so we can generate a large number of samples and show theeffect of the MCMC bias removal. To capture the long-range correlations in the physical system, thereceptive ﬁeld of the network should be able to approximately cover the whole lattice. We use dilated
convolutions [39] to expand the receptive ﬁeld, and increase the dilation rate in each convolutional
layer by a step size. The receptive ﬁeld radius can be calculated by1s−1
D ((D − 1)d + 2)
,
(1)
2
2
where D = 3 is the number of convolutional layers, s = 5 is the convolution kernel size, and d is thedilation step size. For lattice sizes L = 8, 16, 24, 32, we use d = 1, 2, 3, 4 respectively. The network
has 3, 761 non-masked parameters in total, regardless of the lattice size. It is possible to use larger
networks with the similar structure to obtain better variational free energies before the MCMC bias
removal, as discussed in Refs. [27, 30].
Receptive ﬁeld radius =During training, we use Adam optimizer [40] with a conventional learning rate 10−3 , batch size
64, and take 2 × 104 training steps. To avoid being trapped in local minima, especially at low
temperatures, in the ﬁrst 104 steps we linearly anneal β from 0 to the desired value, which is reported
to produce a lower loss than exponential annealing. We do not use weight regularization or gradientclipping, because the network is shallow and there is no signiﬁcant instability in training.
For NCUS, NCU, and NGUS sampling, we generate 103 Markov chains in parallel, each containing
105 samples. The chains are initialized by samples from the network. The ﬁrst 104 samples ineach chain are discarded to ensure that only the samples after thermalization are taken into account.
For each experiment of NCUS up to L = 32, the Gelman–Rubin diagnostic [41] is less than 1.1,
which conﬁrms that the chains are thermalized. The integrated autocorrelation time (IAT) is lessthan 4 × 103 , which is shorter than the remaining chain length by orders of magnitude. For NGU
and NIS sampling, we generate 102 chains with length 106 , and discard the ﬁrst 105 samples for
NGU. The effective autocorrelation time of NIS is computed individually for each “chain”. On an
NVIDIA V100 GPU, the training for the 32 × 32 FPM takes two hours, and the sampling of NCUS
takes a day.
Our code is available at https://github.com/wdphy16/neural-cluster-update

9

B

Uncertainty analysis(a)

105 (b)

NCUS
NCU
NGUS
NGU
NIS

106
105

NCUS
NCU
NGUS
NGU
NIS

104

104

103

103

102

102

101

101
0.40

0.42

0.44

0.46

0.48

0.50

0.40

0.42

0.44

0.46

0.48

0.50

Figure 5: The uncertainties of the IAT in the inset of Fig. 2 (a). (a) The mean values are computed
using a typical random seed, and the error bars are the standard deviations over all chains that we
generate in parallel. (b) The mean values are computed over 5 random seeds, and the error bars are
the standard deviations of the mean values from each random seed. The IAT of NCUS from different
chains and random seeds are consistent within the same order of magnitude, while those of NGU
and NIS are generally impractical to converge.
106

(a)

(b)

NCUS
NCU
NGUS
NGU
NIS

105

NCUS
NCU
NGUS
NGU
NIS

105

104
104
103
103
0.200

0.205

0.210

0.215

0.220

0.225

0.230

0.200

0.205

0.210

0.215

Figure 6: The uncertainties of the IAT in the inset of Fig. 3 (a).

10

0.220

0.225

0.230