The Convergence of Markov Chain Monte Carlo Methods: From the Metropolis Method to Hamiltonian Monte Carlo

From its inception in the 1950s to the modern frontiers of applied statistics, Markov chain Monte Carlo has been one of the most ubiquitous and successful methods in statistical computing. The development of the method in that time has been fueled by not only increasingly difficult problems but also novel techniques adopted from physics. Here, the history of Markov chain Monte Carlo is reviewed from its inception with the Metropolis method to the contemporary state‐of‐the‐art in Hamiltonian Monte Carlo, focusing on the evolving interplay between the statistical and physical perspectives of the method.     

New Monte Carlo method for the self-avoiding walk

We introduce a new Monte Carlo algorithm for the self-avoiding walk (SAW), and show that it is particularly efficient in the critical region (long chains). We also introduce new and more efficient statistical techniques. We employ these methods to extract numerical estimates for the critical parameters of the SAW on the square lattice. We find=2.63820 ± 0.00004 ± 0.00030=1.352 ± 0.006 ± 0.025v=0.7590 ± 0.0062 ± 0.0042 where the first error bar represents systematic error due to corrections to scaling (subjective 95% confidence limits) and the second bar represents statistical error (classical 95% confidence limits). These results are based on SAWs of average length 166, using 340 hours CPU time on a CDC Cyber 170–730. We compare our results to previous work and indicate some directions for future research.     

Hierarchical Monte Carlo methods for fractal random fields

Two hierarchical Monte Carlo methods for the generation of self-similar fractal random fields are compared and contrasted. The first technique, successive random addition (SRA), is currently popular in the physics community. Despite the intuitive appeal of SRA, rigorous mathematical reasoning reveals that SRA cannot be consistent with any stationary power-law Gaussian random field for any Hurst exponent; furthermore, there is an inherent ratio of largest to smallest putative scaling constant necessarily exceeding a factor of 2 for a wide range of Hurst exponentsH, with 0.30<H<0.85. Thus, SRA is inconsistent with a stationary power-law fractal random field and would not be useful for problems that do not utilize additional spatial averaging of the velocity field. The second hierarchical method for fractal random fields has recently been introduced by two of the authors and relies on a suitable explicit multiwavelet expansion (MWE) with high-moment cancellation. This method is described briefly, including a demonstration that, unlike SRA, MWE is consistent with a stationary power-law random field over many decades of scaling and has low variance.     

Potential-Decomposition Strategy in Markov Chain Monte Carlo Sampling Algorithms

We introduce the potential-decomposition strategy （PDS）, which can be used in Markov chain Monte Carlo sampling algorithms. PDS can be designed to make particles move in a modified potential that favors diffusion in phase space, then, by rejecting some trial samples, the target distributions can be sampled in an unbiased manner. Furthermore, if the accepted trial samples are insumcient, they can be recycled as initial states to form more unbiased samples. This strategy can greatly improve efficiency when the original potential has multiple metastable states separated by large barriers. We apply PDS to the 2d Ising model and a double-well potential model with a large barrier, demonstrating in these two representative examples that convergence is accelerated by orders of magnitude.     

Potential-Decomposition Strategy in Markov Chain Monte Carlo Sampling Algorithms

We introduce the potential-decomposition strategy (PDS), which can be used in Markov chain Monte Carlo sampling algorithms. PDS can be designed to make particles move in a modified potential that favors diffusion in phase space, then, by rejecting some trial samples, the target distributions can be sampled in an unbiased manner. Furthermore, if the accepted trial samples are insufficient, they can be recycled as initial states to form more unbiasedsamples. This strategy can greatly improve efficiency when the original potential has multiple metastable states separated by large barriers. We apply PDS to the 2d Ising model and a double-well potential model with a large barrier, demonstrating in these two representative examples that convergence is accelerated by orders of magnitude.     

The pivot algorithm: A highly efficient Monte Carlo method for the self-avoiding walk

The pivot algorithm is a dynamic Monte Carlo algorithm, first invented by Lal, which generates self-avoiding walks (SAWs) in a canonical (fixed-N) ensemble with free endpoints (hereN is the number of steps in the walk). We find that the pivot algorithm is extraordinarily efficient: one effectively independent sample can be produced in a computer time of orderN. This paper is a comprehensive study of the pivot algorithm, including: a heuristic and numerical analysis of the acceptance fraction and autocorrelation time; an exact analysis of the pivot algorithm for ordinary random walk; a discussion of data structures and computational complexity; a rigorous proof of ergodicity; and numerical results on self-avoiding walks in two and three dimensions. Our estimates for critical exponents are=0.7496±0.0007 ind=2 and= 0.592±0.003 ind=3 (95% confidence limits), based on SAWs of lengths 200N10000 and 200N 3000, respectively.     

Nonergodicity of local,length-conserving Monte Carlo algorithms for the self-avoiding walk

It is proved that every dynamic Monte Carlo algorithm for the self-avoiding walk based on a finite repertoire of local,N-conserving elementary moves is nonergodic (hereN is the number of bonds in the walk). Indeed, for largeN, each ergodic class forms an exponentially small fraction of the whole space. This invalidates (at least in principle) the use of the Verdier-Stockmayer algorithm and its generalizations for high-precision Monte Carlo studies of the self-avoiding walk.     

Coupling control variates for Markov chain Monte Carlo

We show that Markov couplings can be used to improve the accuracy of Markov chain Monte Carlo calculations in some situations where the steady-state probability distribution is not explicitly known. The technique generalizes the notion of control variates from classical Monte Carlo integration. We illustrate it using two models of nonequilibrium transport.     

High-precision Monte Carlo test of the conformai-invariance predictions for two-dimensional mutually avoiding walks

Let _l be the critical exponent associated with the probability thatl independentN-step ordinary random walks, starting at nearby points, are mutually avoiding. Using Monte Carlo methods combined with a maximum-likelihood data analysis, we find that in two dimensions ₂=0.6240±0.0005±0.0011 and ₃=1.4575±0.0030±0.0052, where the first error bar represents systematic error due to corrections to scaling (subjective 95% confidence limits) and the second error bar represents statistical error (classical 95% confidence limits). These results are in good agreement with the conformal-invariance predictions ₂=5/8 and ₃=35/24.     

Auxiliary-field Monte Carlo methods

I discuss Monte Carlo algorithms for quantum many-body systems that employ an auxiliary field to linearize a two-body interaction. These reduce the evaluation of the partition function to sampling many one-body evolutions in a fluctuating field. Fermions and bosons are treated on an equal footing. Applications to potential models and to quantum spin systems are discussed. This work was supported in part by the National Science Foundation, grants PHY82-07332 and PHY85-05682. The potential-model studies were done in collaboration with G. Sugiyama, while A. Khan and T. Troudet were responsible for the work on the quantum spin systems.     

A myopic random walk on a finite chain

We solve analytically the problem of a biased random walk on a finite chain of ‘sites’ (1,2,…,N) in discrete time, with ‘myopic boundary conditions’—a walker at 1 (orN) at timen moves to 2 (orN − 1) with probability one at time (n + 1). The Markov chain has period two; there is no unique stationary distribution, and the moments of the displacement of the walker oscillate about certain mean values asn → ∞, with amplitudes proportional to 1/N. In the continuous-time limit, the oscillating behaviour of the probability distribution disappears, but the stationary distribution is depleted at the terminal sites owing to the boundary conditions. In the limit of continuous space as well, the problem becomes identical to that of diffusion on a line segment with the standard reflecting boundary conditions. The first passage time problem is also solved, and the differences between the walks with myopic and reflecting boundaries are brought out.     

求解中子动力学方程的加权蒙特卡罗方法

为了实现基于蒙特卡罗方法的中子动力学计算,在传统的直接蒙特卡罗动力学方法的基础上,提出了一种加权蒙特卡罗动力学方法。该方法通过引入粒子权重的概念,隐式考虑中子俘获反应和裂变反应过程中中子数目的变化,避免了模拟粒子的数目随时间的变化,降低了统计偏差,消除了程序计算过程中粒子的存库操作,提高了计算精度。基于单能点堆模型,开发了中子动力学计算程序NECP-Dandi,进行了大量数值验证与分析,包括无缓发中子、单组缓发中子、六组缓发中子、正阶跃反应性引入、负阶跃反应性引入、正脉冲反应性、负脉冲反应性和正线性反应性引入等情况。数值结果表明,相比于直接蒙特卡罗动力学方法,加权蒙特卡罗动力学方法在计算结果的精度和计算效率上有较为明显的改进,程序结构更为简洁。     

Note on Monte Carlo methods without necessary detailed balance

Simple criteria for convergence of Monte Carlo algorithms not necessarily requiring detailed balance for any specified transition probability are derived and it is shown that it is possible to view the algorithm as a superimposition of a Brownian motion on configurational space coupled to the transition probabilities. As such, the error contributions due to a particular Monte Carlo algorithm and the integration limits in configuration space must be distinguished from those due to the nonuniform sampling of the Brownian motion, and criteria related to the number of steps required to distinguish these errors are provided for the simplest cases involving one dimension and symmetrical probability distributions.     

Bilinear quantum Monte Carlo: Expectations and energy differences

We propose a bilinear sampling algorithm in the Green's function Monte Carlo for expectation values of operators that do not commute with the Hamiltonian and for differences between eigenvalues of different Hamiltonians. The integral representations of the Schrödinger equations are transformed into two equations whose solution has the form _a(x) t(x, y)_b(y), where _a and _b are the wavefunctions for the two related systems andt(x, y) is a kernel chosen to couplex andy. The Monte Carlo process, with random walkers on the enlarged configuration spacex y, solves these equations by generating densities whose asymptotic form is the above bilinear distribution. With such a distribution, exact Monte Carlo estimators can be obtained for the expectation values of quantum operators and for energy differences. We present results of these methods applied to several test problems, including a model integral equation, and the hydrogen atom.     

Critical exponent for the loop erased self-avoiding walk by Monte Carlo methods

A Monte Carlo simulation was performed for loop-erased self-avoiding walks (LESAW) to ascertain the exponentv for the Z² and Z³ lattices. The estimated values were 2v=1.600±0.006 in two dimensions and 2v=1.232±0.008 in three dimensions, leading to the conjecturev=4/5 for the two-dimensional LESAW. These results add to existing evidence that the loop-erased self-avoiding walks are not in the same universality class as self-avoiding walks.     

Random walk with persistence

The exact analytic result is obtained for the Fourier transform of the generating functionF(R,s)= _n=0 s ⁿ P(R,n), whereP(R,n) is the probability density for the end-to-end distanceR inn steps of a random walk with persistence. The moments R ²(n), R ⁴(n), and R ⁶(n) are calculated and approximate results forP(R,n) and R ^–1(n) are given.     

模拟回火马尔可夫链蒙特卡罗全波形分析方法

针对传统的全波形分析方法不能快速自动处理全波形数据的缺点,提出了一种模拟回火马尔可夫链蒙特卡罗全波形分析法,用于求解全波形数据中的波峰数和峰值位置等参量.该方法采用Metropolis更新策略求解波峰数量和噪声两个参量,以达到快速求解的目的;而峰值位置和波峰幅值则采用改进的模拟回火策略求解,通过添加的主动干预回火步骤实现对参量更新过程的有效探测,以满足对速度或运算收敛性的要求.模拟回火马尔可夫链蒙特卡罗全波形分析方法以马尔可夫算法为基础,仍保持马氏链的收敛性,从而保证本方法具有良好的鲁棒性,实现对全波形数据的自动化处理.     

Transition Matrix Monte Carlo Method

We present a formalism of the transition matrix Monte Carlo method. A stochastic matrix in the space of energy can be estimated from Monte Carlo simulation. This matrix is used to compute the density of states, as well as to construct multi-canonical and equal-hit algorithms. We discuss the performance of the methods. The results are compared with single histogram method, multi-canonical method, and other methods. In many aspects, the present method is an improvement over the previous methods.     

Monte Carlo methods in radiative transfer and electron-beam processing

Monte Carlo methods (MCMs) are the most versatile approaches in solving the integro-differential equations. They are statistical in nature and can be easily adapted for simulation of the propagation of ensembles of quantum particles within absorbing, emitting, and scattering media. In this paper, we use MCM for the solution of the Boltzmann transport equation, which is the governing equation for both radiative transfer and electron-beam processing. We briefly outline the methodology for the solution of MCMs, and discuss the similarities and differences between the two different application areas. The focus of this paper is primarily on the treatment of different scattering phase functions.     

输运问题蒙特卡罗模拟方法回顾及展望北大核心CSCD

蒙特卡罗(MC)方法具有复杂几何处理能力强,方法通用灵活,核数据完备,模拟忠实于物理过程等特点,成为中子学数值模拟的首选方法之一。在核能领域,MC方法得益于计算机的快速发展,在辐射屏蔽、反应堆堆芯临界安全分析、乏燃料后处理、放射性废物处置、核设施退役、核事故应急、放射性石油测井、核医学等领域均有广泛应用。对MC方法及软件输运计算做简要回顾,并对未来发展进行展望。