Accelerated Markov Chain Monte Carlo Simulation via Neural Network–Driven Importance Sampling

@article{kim2026acceleratedmcmc,
  title={Accelerated Markov Chain Monte Carlo Simulation via Neural Network-Driven Importance Sampling},
  author={Kim, Michael and Cai, Wei},
  journal={arXiv preprint arXiv:2602.12294},
  year={2026}
}

1. Introduction

Atomistic simulations play a crucial role in understanding the fundamental processes governing the behavior of materials at the microscopic scale. However, they are fundamentally limited by the accessible time scale of brute-force methods. Simulations frequently remain trapped within metastable states for extended periods, hindering the observation of transitions between these states. Consequently, capturing these transitions that drive long-term evolution can be prohibitively expensive. This intrinsic limitation is widely recognized as the time-scale or rare event problem[1]–[3].

Various methodologies have been established to overcome the time-scale limitations of standard simulations. Methods such as hyperdynamics[4], steered molecular dynamics[5], and metadynamics[6]–[8] rely on external forces or bias potentials to accelerate the escape from metastable basins. Alternatively, trajectory-based approaches, including the weighted ensemble method[9],[10] and forward flux sampling[11], enhance sampling efficiency by selectively branching trajectories along a reaction coordinate. In parallel, transition path theory[12]–[15] offers a mathematical basis for analyzing the statistics of reactive pathways. Importance sampling methods[16]–[21] represent another established direction, designed to capture dominant transition pathways and estimate dynamical properties of rare events. In this work, we focus on applying importance sampling to extend the time scale of Markov chain Monte Carlo simulations.

The importance sampling method pursues two primary objectives. First, it aims to enhance the sampling of successful transition paths between metastable states without artificially favoring one reaction channel over another: the relative probabilities of competing pathways must be preserved so that the transition mechanism remains unbiased. Second, the applied enhancement must be quantifiable, allowing the true transition rates of the original system to be recovered. In practice, importance sampling modifies the system kinetics by updating the transition probabilities with an importance function[18],[21]. The requirement that these updated probabilities remain normalized serves as a constraint that uniquely determines the optimal importance function. Consequently, importance sampling reformulates the rare event problem into an optimization problem.

Several significant challenges remain. An increase in the size or dimensionality of the system drastically increases the computational burden of finding the optimal importance function. In addition, numerical stability becomes critical at low temperatures, where the magnitude of the optimal importance function near the energy minima of the starting metastable state can become exceedingly small, leading to underflow errors. Most importantly, the accuracy of the estimated transition rates relies heavily on the quality of the importance function: even a small deviation from the optimal solution causes the variance of the estimator to increase drastically, resulting in large errors in the estimated rates.

In this work we propose a robust method to overcome these limitations in Markov chain Monte Carlo (MCMC) simulations of spatially discrete systems. First, to address the challenge of searching for the optimal importance function in high dimensions, we employ a neural network combined with adaptive sampling. Second, to resolve numerical instability at low temperatures, we perform the optimization in logarithmic space by introducing an optimal bias potential, which effectively prevents underflow near energy minima. Finally, we present a general formulation for estimating transition rates, complemented by a branching random walk (BRW) technique that significantly reduces the variance of the estimator.

2. Problem Statement

Consider a domain $\Omega$, a set of allowed microscopic states corresponding to points on a regular grid in a two-dimensional region, as shown in Figure 1. On this domain an energy landscape $E(i)$ is defined, specifying the energy of any microscopic state $i$. The kinetics of the system is defined through the function $r(i\to j)$, the jump rate between any two microscopic states $i$ and $j$ on the grid,

where $\nu_0$ is a constant frequency prefactor and $\mathbb{1}_{\rm nn}(i,j)$ equals $1$ if $i$ and $j$ are nearest-neighbor grid points in $\Omega$ and $0$ otherwise. The jump rate defined in Eq. (1) satisfies the detailed balance condition,

Assume the energy landscape contains two deep basins corresponding to two metastable states. Let A and B be the local energy minima at the bottom of these basins, and assume there are two saddle points, $S_1$ and $S_2$, located on two distinct transition pathways between A and B. At low temperatures, a system starting in the neighborhood of A is likely to spend a long time in the basin containing A, and the transition to the basin containing B is a rare event. When such a rare transition does occur, the transition path is likely to pass through the neighborhood of either $S_1$ or $S_2$.

Importance sampling aims to accelerate the time scale of MCMC simulations while preserving the relative probabilities of distinct pathways. The objectives are twofold: (1) obtain accurate transition rates between metastable states, and (2) reveal the underlying transition mechanisms by sampling transition paths with their correct probabilities. Although the example above is a simple 2-dimensional system, the algorithm is designed to be robust enough to handle high-dimensional systems, where rare events are computationally prohibitive.

3. Methods

3.1 Importance Sampling of Rare Transition Paths

If the system is in microscopic state $i$, the probability of jumping to microscopic state $j$ upon leaving $i$ is

where $r_{\rm tot}(i)$ is the total jump rate out of state $i$ and $\mathcal{N}(i)$ is the set of all microscopic states accessible from (i.e. adjacent to) $i$. We define a path as a Monte Carlo trajectory that starts at A and ends at either A or B. Paths that end at A are called failure paths, and paths that end at B are success paths. The probability of sampling a given success path $\Gamma_{\rm AB} = \{{\rm A}, i_1, i_2, \cdots, i_n, {\rm B}\}$ is

The rare event problem stems from the fact that, under the dynamics defined by $K(i\to j)$, the probability of sampling a success path is extremely small. Direct Monte Carlo simulations are therefore computationally inefficient, as the vast majority of resources are expended sampling failure paths. With importance sampling we construct a new kinetics in which the jump probability $K'$ is modified from the original probability $K$ through an importance function $I(i)$,

where $\mathbb{1}_{\rm A}(j) = 1$ if $j = {\rm A}$ and $0$ otherwise. However, $K'(i\to j)$, being a probability matrix, has to be normalized,

We define the importance function that guarantees the normalization condition, Eq. (6), as the optimal importance function $I_{\rm opt}(i)$. In general, when $I(i)$ does not equal $I_{\rm opt}(i)$ (because the latter is unknown), we must introduce a normalization factor $n'(i)$ to ensure that $K'(i\to j)$ is properly normalized,

The new probability of sampling the same path $\Gamma_{\rm AB}$ under the modified Monte Carlo simulation is

In the event that all the normalization factors equal $1$ (when $I(i) = I_{\rm opt}(i)$), we obtain a remarkable outcome,

This expression implies that the two requirements identified in Section 2 are exactly satisfied: (1) the sampling of all successful paths from A to B is enhanced by the same amount (the enhancement factor is path independent), and (2) the enhancement factor is known, so we can obtain the transition rate of the original system from that of the modified system. In practice $I_{\rm opt}(i)$ is not known, and finding it exactly is at least as difficult as solving the original rare event problem; the well-established approach of importance sampling (widely successful in quantum chemistry[22]) is to develop “good enough” approximations to $I_{\rm opt}(i)$. What has been achieved is the reformulation of the rare event problem into an optimization problem. The optimality condition Eq. (6) can be written equivalently as

In other words, the optimal importance function is the (right) eigenvector of the transition probability matrix $K(i\to j)\left[1 - \mathbb{1}_{\rm A}(j)\right]$ with eigenvalue $1$ (which corresponds to the largest eigenvalue). We set $I({\rm B}) = 1$ without loss of generality. The optimal importance function corresponds to the discrete committor function (Section 3.3). Recalling Eq. (1), modifying the jump probability matrix using $I(i)$ as in Eq. (5) is equivalent to introducing a bias potential,

which is added to the original potential energy $E(i)$ to modify the Monte Carlo trajectory. Searching for the optimal importance function is thus equivalent to searching for the optimal bias potential.

3.2 Generalized Formulation of Importance Sampling of Rare Events

Defining failure and success through single discrete microscopic states (A or B) becomes an artifact when the grid spacing $\Delta x$ is systematically refined: the optimal importance function then exhibits an increasingly sharp gradient around the points A and B. This is analogous to applying boundary conditions at a single point of a continuum satisfying a partial differential equation, which is generally ill-advised. A more robust formulation defines failure and success of a trajectory using a finite region whose size stays constant as the grid is refined.

Instead of selecting two microscopic states A and B, we augment the original system with two additional states, F and S, enlarging the domain to $\tilde{\Omega} = \Omega\cup\{{\rm F},{\rm S}\}$. Multiple microscopic states (around A or B) are given a nonzero jump rate into F or S. If a trajectory enters F it is called a “failure”; if it enters S it is called a “success.” The functions $f(i)$ and $s(i)$ control whether microstate $i$ is allowed to jump into F or S; they are non-zero only in a small area around A and B, respectively, and can be selected to have the form of a truncated Gaussian,

where $\mathbf{x}(i)$ is the position vector of microscopic state $i$ in a $d$-dimensional system, $\|\cdot\|$ is the Euclidean norm, and $\mathbb{1}(\cdot)$ is an indicator function. The prefactor $(\Delta x)^d$ ensures that the flux between state $i$ and F remains nearly constant regardless of the choice of $\Delta x$. After also defining the energies $E({\rm F})$ and $E({\rm S})$ of the two auxiliary states, the jump rate between any two microscopic states $i$ and $j$ ($i \neq j$) in $\tilde{\Omega}$ is

where $\mathbb{1}_{\rm F}$ and $\mathbb{1}_{\rm S}$ are the indicators for states F and S, and $\tilde{r}(i\to j) = 0$ if $i = j$. This jump rate satisfies detailed balance in $\tilde{\Omega}$, exactly as $r(i\to j)$ does in $\Omega$. The goal of this design is to minimize the perturbation of the kinetics of the original system: all we have added are two extra states, whose effect becomes minimal as the total number of states in $\Omega$ becomes very large. Given the jump rates, the probability of jumping to microscopic state $j$ upon leaving $i$ is

where $\mathcal{A}(i) = \mathcal{N}(i)\cup\{{\rm F},{\rm S}\}$ is the set of microscopic states the system is allowed to jump to from $i$. Because $f(i) > 0$ only in a small region around A and $s(i) > 0$ only in a small region around B, the system repeatedly visits F many times before it visits S for the first time. A trajectory that starts and ends at F is a failure path, and one that starts at F and ends at S is a success path; the rare event problem is equivalent to the extremely low probability of sampling success paths using $\tilde{K}(i\to j)$. To overcome this, the importance function $I(i)$, defined on every $i\in\Omega$, modifies the transition probability as

where $\tilde{n}'(i)$ is the normalization factor. The optimality condition is $\tilde{n}'(i) = 1$ for all $i\in\Omega$, which can be equivalently written as

Note that $\mathcal{A}(i)\setminus\{{\rm F}\} = \mathcal{N}(i)\cup\{{\rm S}\}$, i.e. the sum over $j$ does not include state F. Note also that in this formulation importance sampling is not applied to jumps out of state F, i.e. $\tilde{K}'({\rm F}\to j) = \tilde{K}({\rm F}\to j)$, so it is not necessary to define $I({\rm F})$. In the following we set $I({\rm S}) = 1$ without loss of generality.

3.3 Transition Rate Estimation

With the auxiliary states F and S, the transition rate between the two metastable states is very close to the transition rate between the microscopic states F and S, $r_{\rm FS}$, which equals the inverse of the mean first passage time $T_{\rm FS}$ from F to S,

where $\langle\cdot\rangle$ indicates an average over all statistical realizations. Consider a trajectory that starts from F and eventually ends at S, as shown in Figure 2. In a typical situation the system visits F many times before visiting S for the first time, so the mean first passage time reads

where $\langle N_{\rm F}\rangle$ is the average number of times the system returns to F before visiting S, $\langle t_{\rm FF}\rangle$ is the average time of a failure path, and $\langle t_{\rm FS}\rangle$ is the average time of a success path. In evaluating $\langle t_{\rm FF}\rangle$ we adopt the convention of excluding the time accrued while the system stays in state F, so that $r_{\rm FS}$ follows the transition rate of the original system as closely as possible. In the rare event limit, $\langle N_{\rm F}\rangle \gg 1$, and

where $p_{\rm S}({\rm F})$ is the probability of the system starting from F to visit S before returning to F. In general, $p_{\rm S}(i)$ — the probability of visiting S before F starting from $i$ — is the discrete committor function and equals the optimal importance function $I_{\rm opt}(i)$. The transition rate can therefore be written as

where $I_{\rm opt}({\rm F}) = \sum_{j\in\Omega}\tilde{K}({\rm F}\to j)\,I_{\rm opt}(j)$. It is straightforward to compute $\langle t_{\rm FF}\rangle$ using the standard kinetic Monte Carlo algorithm[23]–[25] with $\tilde{r}(i,j)$ and without importance sampling, since failure paths are abundant. Our primary task is thus to determine $p_{\rm S}({\rm F})$ using an approximate importance function $I(i)$.

The transition rate $r_{\rm FS}$ is proportional to the frequency prefactor $\nu_0$ appearing in Eq. (1). Since the grid spacing $\Delta x$ is a numerical artifact and it is of interest to consider the continuum limit $\Delta x\to 0$, we let the frequency prefactor depend on temperature and grid spacing as

where $\mu$ is the mobility constant (we choose $\mu = 1.0\ \mathrm{m^2\,s^{-1}\,eV^{-1}}$), $k_{\rm B}$ is the Boltzmann constant, $T$ is the temperature, and $d$ is the dimension. This allows the transition rate $r_{\rm FS}$ to remain unchanged in the limit $\Delta x\to 0$ (Figure 6a).

3.4 Success Probability Calculation

We now compute the success probability $p_{\rm S}({\rm F})$ given an importance function $I(i)$ that is not necessarily optimal. Define a path $\Gamma$ as a trajectory that starts at F and ends at either F or S, without visiting F or S in between, and let $P(\Gamma)$ be the probability of sampling $\Gamma$ following $\tilde{K}(i,j)$. Let $\{\Gamma_{\rm FF}\}$ denote the set of failure paths and $\{\Gamma_{\rm FS}\}$ the set of success paths. The success probability is

Defining a “reward” function $R(\Gamma) = 1$ if $\Gamma\in\{\Gamma_{\rm FS}\}$ and $0$ otherwise, the success probability can be rewritten as a sum over all paths, and hence as a statistical average over all paths sampled by the Monte Carlo simulation,

Even though Eq. (24) is a simple recipe, for systems exhibiting rare event behavior it is extremely inefficient, because the vast majority of sampled paths are failure paths for which the reward is zero. Using the importance function to modify the transition probability to $\tilde{K}'(i\to j)$ as defined in Eq. (16), the new probability of sampling a specific success path $\Gamma_{\rm FS} = \{{\rm F}, i_1, i_2, \cdots, i_n, {\rm S}\}$ is related to the original probability by

Defining the weight of a path as the product of the normalization factors of all states visited along the path, excluding the starting and final states F and S,

the original probability becomes $P(\Gamma_{\rm FS}) = W(\Gamma_{\rm FS})\,\frac{I(i_1)}{I({\rm S})}\,P'(\Gamma_{\rm FS})$, which substituted into Eq. (24) gives

where $\langle\cdot\rangle'_\Gamma$ means averaging over all paths sampled using the transition probability $\tilde{K}'$, i.e. using importance sampling. The advantage is that if $I(i)$ is close to $I_{\rm opt}(i)$, it is very likely to sample success paths for which the reward is non-zero. If $I(i)$ equals $I_{\rm opt}(i)$, the success probability expression reduces to

meaning the success probability could be computed without sampling any transition paths if $I_{\rm opt}(i)$ were known. However, since in practice the accessible $I(i)$ is only an approximation of $I_{\rm opt}(i)$, the proper way to compute the success probability is through Eq. (27), which requires importance sampling of transition paths. Recall that the first step of every path — the jump out of state F — is not subjected to importance sampling (Section 3.7).

3.5 Efficient Sampling via Branching Random Walk

While Eq. (27) provides an unbiased estimate of $p_{\rm S}({\rm F})$, the number of paths needed for statistical convergence depends strongly on the quality of $I(i)$. In practice $I(i)$ may be so different from $I_{\rm opt}(i)$ that an excessive number of paths is needed, undermining the efficiency of the method. This is addressed by replacing the importance-sampled Monte Carlo simulation with a branching random walk (BRW) process.

The BRW process starts from a single random walker located at state F, with its weight set to unity. At each step the walker updates its weight by multiplying it with the normalization factor of its current state. If the weight $W$ is less than a threshold $W_{\min}$ (e.g. $0.5$), an attempt is made to annihilate the walker with probability $1 - W$; if annihilated, the path is terminated and marked a failure (reward zero), and if it survives, its weight is restored to $1$. If the weight $W$ is greater than a threshold $W_{\max}$ (e.g. $2$), the walker is split into $\lfloor W\rfloor + 1$ walkers with probability $W - \lfloor W\rfloor$, or $\lfloor W\rfloor$ walkers otherwise, with all weights reset to $1$. This can be formulated rigorously with a stochastic rounding function, a technique widely applied in Monte Carlo rendering and particle transport simulations[26]–[28],

At every step, if a walker's weight falls outside the threshold bounds $[W_{\min}, W_{\max}]$, the walker transforms into $\mathcal{R}(W)$ new walkers: $\mathcal{R}(W) = 0$ terminates the walker, while $\mathcal{R}(W)\geq 1$ means the process proceeds with $\mathcal{R}(W)$ walkers. If any walker visits state S it is terminated and its path is marked a success. Once the process is complete, each successful walker receives a base reward of unity, weighted by the product of its current path weight $W(\Gamma)$ and $I(i_1)$, the value of the importance function at the first intermediate state. This weighted summation is a statistical estimate of $\langle R(\Gamma)W(\Gamma)I(i_1)\rangle'_\Gamma$, ensuring the result is unbiased, while the branching mechanism keeps the path weights controlled within $[W_{\min}, W_{\max}]$.

3.6 Neural Network Representation of the Bias Potential

As the grid is refined, the total number of microscopic states increases — drastically so in high dimensions — making it computationally intractable to optimize the importance function $I(i)$, or the bias potential $E_{\rm b}(i)$, for every microscopic state. To overcome this we integrate machine learning techniques[29]–[31] into the framework: the importance function is represented by a neural network $I(i;\theta)$ with trainable parameters $\theta$, optimized towards satisfying the normalization condition $\tilde{n}'(i) = 1$ everywhere through the loss function

In systems exhibiting rare event problems, $I_{\rm opt}(i)$ often yields extremely small values at microscopic states near F, so training a network to approximate the importance function directly can encounter numerical instabilities due to underflow. We address this by optimizing the bias potential instead: by Eq. (11), even when the importance function takes extremely small values the corresponding bias potential still takes values of manageable magnitude. Parameterizing $E_{\rm b}(i;\theta)$ with a neural network keeps it within a numerically stable range, and the loss function is rewritten as

This approach is analogous to a recently proposed method[31]. In practice it is impossible — and unnecessary — to perform the sum in Eq. (30) over all microscopic states: the normalization condition only needs to be well satisfied on states that lie along the transition paths from F to S. We therefore replace the sum with one over the microscopic states sampled by the importance sampling itself, using an adaptive sampling scheme (Algorithm 1). A random walk is initiated from state F, governed by $\tilde{K}'(i\to j)$, with the network parameters held fixed; once enough states are accumulated, the network is updated by minimizing the loss evaluated over the sampled states, and the algorithm returns to the sampling stage with the updated bias potential. Furthermore, we apply a simulated annealing protocol: training starts at an elevated temperature and decreases stepwise to the target temperature, which helps avoid numerical instabilities and benefits convergence.

A useful feature of this framework is that a neural network trained on a coarse grid can be applied directly to accelerate transitions on a finer grid without retraining. This is particularly valuable because the length of the transition paths increases as the grid is refined, raising the cost of bias-potential training (which requires adaptive path sampling). The branching random walk of Section 3.5 brings the variance under control even when the transferred bias potential is not the optimal one.

3.7 Sampling Initial Microscopic States

The first step in path sampling is to sample the microscopic state $i$ jumping out of state F according to $\tilde{K}({\rm F}\to i)$. In low-dimensional systems this can be done by enumerating all states accessible from F, but in high dimensions enumeration is impossible. We achieve a tractable scheme by modifying the control function $f(i)$ of Eq. (12) to

where $\xi$ is a dimensionless constant and $E_{\rm C}(i)$ is a newly introduced confinement potential — an arbitrary convex function having its global minimum at the microscopic state A. For simplicity it is chosen as

with a constant parameter $\eta$; the confinement potential smoothly confines the sampled states to lie near A. According to Eq. (14), the jump rate from F to $i$ then becomes

This significantly simplifies the sampling of initial microscopic states in high dimensions, because Eq. (34) corresponds to the equilibrium distribution governed by the effective potential $E(i) + E_{\rm C}(i) - E({\rm F})$; such states can be readily sampled with a standard Monte Carlo simulation. It is not necessary to compute the total jump rate out of F, because the time spent in F is excluded when calculating the transition rate (Section 3.3). The reverse jump rate follows from detailed balance,

which grows as the microscopic state $i$ gets closer to A, provided the confinement potential is convex with its minimum at A.

4. Results

The importance sampling framework integrated with a neural network is applied to 2-dimensional and 14-dimensional systems to demonstrate its effectiveness and scalability.

4.1 Two-Dimensional System

The 2-dimensional potential energy landscape (in eV) shown in Figure 1 can be written analytically as

The landscape is discretized over the domain $x_1, x_2 \in [-1.5, 1.5]$ using equal grid spacing $\Delta x_1 = \Delta x_2 = 0.1$. The energy minima are located at microscopic state ${\rm A}\ (-1.1, 0.0)$ and microscopic state ${\rm B}\ (1.1, 0.0)$, which correspond to the two metastable states, and the two saddle points at $S_1\ (0.0, 1.0)$ and $S_2\ (0.0, -1.0)$, with energy barriers of 1.02 eV at $S_1$ and 0.98 eV at $S_2$. Transitions from A to B are rare events, particularly at low temperature, and the primary transition pathways proceed through regions proximal to $S_1$ and $S_2$. The truncated-Gaussian control functions of Eqs. (12)–(13) use $\epsilon = 0.1$, $\sigma = 5\times 10^{-3}$, $r_{\rm c} = 0.3$, $E({\rm F}) = -0.5$ eV and $E({\rm S}) = -0.5$ eV.

The neural network architecture comprises two hidden layers, each containing 30 perceptrons with $\tanh$ activation functions. We employ a simulated annealing protocol in which training initially starts at $5800$ K and is subsequently cooled to the target of $500$ K over a sequence of $8$ discrete stages. At each temperature the training is performed for $30$ or $100$ epochs using the ADAM[32] algorithm with a learning rate of $10^{-3}$. For the 2-dimensional system we can enumerate all the states and store the entire transition probability matrix $K(i\to j)$, so the optimal importance function $I_{\rm opt}(i)$ can be solved numerically[18] and the exact optimal bias potential obtained from Eq. (11). The neural network bias potential trained for $100$ epochs at each temperature shows substantially improved agreement with the exact optimal bias potential compared with the one trained for $30$ epochs. Notably, the largest deviations occur at sparsely sampled microscopic states that lie far from the dominant transition pathways (Figure 3).

4.2 Success Probability Estimate

We demonstrate in the 2-dimensional system that the success probability estimator given by Eq. (27) is unbiased even when the importance function is not exactly the optimal solution. The optimal importance function can be used to compute the reference success probability,

which is $2.1899\times 10^{-13}$ at $T = 500$ K ($k_{\rm B}T = 0.043087$ eV). The neural network bias potential $E_{\rm b}(i;\theta)$ is then employed to accelerate the transition and sample 100 successful paths, whose scores are used to estimate the success probability $p_{\rm S}^{\theta}({\rm F})$ via Eq. (27); this accelerated sampling was repeated 100 times to improve statistical accuracy. The estimated success probabilities are $(2.1921\pm 0.0082)\times 10^{-13}$ and $(2.1937\pm 0.0065)\times 10^{-13}$ for models trained for 30 and 100 epochs at each temperature, respectively, both close to the reference value (Table 1).

The deviation of the mean path weights from unity indicates that the neural network importance functions are not optimal; nevertheless, their mean values provide unbiased estimates of the true success probability. The network trained for a longer duration (100 epochs at each temperature) exhibits sharper weight and success-probability distributions, indicating an improved approximation of the optimal importance function (Figure 4).

Table 1. Success probability estimates for the 2-dimensional system at $500$ K from the unbiased estimator Eq. (27), using neural network bias potentials trained for 30 and 100 epochs at each annealing temperature, compared against the reference value computed from the exact optimal importance function. Both estimates are unbiased; the longer-trained model has the smaller variance.

4.3 Branching Random Walk

To study the robustness of the method against approximation errors, we intentionally introduce a discretization mismatch between the training and sampling phases: the neural network importance function trained on a coarse grid ($\Delta x_1 = \Delta x_2 = 0.1$) is used to bias the dynamics of the system on a refined grid ($\Delta x_1 = \Delta x_2 = 0.025$). Sampling without the BRW process over 50 batches of 100 paths each shows that the majority of sampled successful paths have near-zero weights, i.e. they are unimportant paths contributing minimally to the estimate. The estimated success probability is $(1.5198\pm 0.0627)\times 10^{-14}$, while the reference (exact) value is $1.4120\times 10^{-14}$.

Using the BRW process of Section 3.5, the weights of successful paths are maintained within a specified range, and paths with sufficiently small weights are terminated early, enhancing efficiency. Sampling was done with the BRW process over 100 batches of 100 paths each (weights controlled within $[0.5, 1.2]$) to match the variance obtained without BRW, giving an estimated success probability of $(1.5210\pm 0.0690)\times 10^{-14}$. While both approaches yield unbiased estimates, the BRW process is markedly more efficient: the total number of Monte Carlo steps without BRW is $7{,}967{,}789$, whereas with BRW it is $958{,}928$ — an approximately 8-fold speedup (Table 2, Figure 5).

Table 2. Success probability estimates and Monte Carlo cost with and without the branching random walk (BRW) process for the 2-dimensional system, where a bias potential trained on the coarse grid ($\Delta x = 0.1$) is applied to a refined grid ($\Delta x = 0.025$). Both estimators are unbiased, but the BRW process reaches a comparable variance with roughly $8\times$ fewer Monte Carlo steps.

4.4 Transition Rate and Mechanisms

We demonstrate that the transition rates estimated by the method of Section 3.3 are consistent across grid sizes. For systems with different numbers of grid points per dimension (corresponding to grid sizes $\Delta x = 3.0/(N-1)$) at $T = 500$ K, and using the neural network bias potential trained on the $\Delta x = 0.1$ system throughout, the estimated transition rates converge with increasing $N$ to a value in good agreement with the exact transition rate of $6.9191\times 10^{-12}$ for the 2-dimensional system with grid size $\Delta x = 0.025$ (Figure 6a). This confirms that the frequency prefactor of Eq. (22) keeps $r_{\rm FS}$ invariant under grid refinement.

To validate accuracy, we compare against Kramers rate theory[33],[34] for overdamped dynamics. The channel-resolved Kramers rate expressions read

where $\nu_1$ and $\nu_2$ are frequency prefactors associated with the saddle points $S_1$ and $S_2$ (Appendix A). The Arrhenius plot of the estimated transition rates demonstrates good agreement with Kramers rate theory (Figure 6b).

The accelerated sampling also accurately resolves the transition rates associated with the distinct reaction channels. At $T = 500$ K the accelerated Monte Carlo simulation estimates the fraction of the transition rate proceeding through $S_1$ to be $0.2910$, and at $T = 1000$ K the fraction increases to $0.3967$; these agree well with the Kramers predictions of $0.2938$ and $0.3951$, respectively (Figure 6c). Notably, this fraction cannot be determined from the energy difference between the two saddle points alone (i.e. by assuming $\nu_1 = \nu_2$ in Eqs. 38–39); such a simplified analytic theory deviates significantly from the Monte Carlo results. The observed temperature dependence arises from the distinct rate prefactors $\nu_1$ and $\nu_2$ set by the different local curvatures of the potential energy landscape at $S_1$ and $S_2$.

4.5 Fourteen-Dimensional System

A 14-dimensional potential energy landscape (in eV) is constructed by extending the 2-dimensional system of Eq. (36) through the addition of harmonic terms in the remaining coordinates,

with $\zeta = 0.5\ \mathrm{eV^{-1/2}}$, considered at $T = 500$ K. The system is discretized over the domain $x_i \in [-1.5, 1.5]$, a hypercube in 14-dimensional space, using equal lattice spacing $\Delta x_i = 0.1$ for $i = 1, 2, \ldots, 14$. As in the 2-dimensional case, the energy minima are at ${\rm A}: (-1.1, 0.0, 0.0, \ldots, 0.0)$ and ${\rm B}: (1.1, 0.0, 0.0, \ldots, 0.0)$, while the first-order saddle points are at $S_1: (0.0, 1.0, 0.0, \ldots, 0.0)$ and $S_2: (0.0, -1.0, 0.0, \ldots, 0.0)$. The transition pathway of interest depends solely on the first two coordinates of the position vectors.

For the 14-dimensional system we adopt a hybrid parameterization of the neural network bias potential,

the sum of a Gaussian bias potential and a fully connected feed-forward neural network, with $\{A, a_i, c_i, \theta\}$ forming the trainable parameter set. The explicit Gaussian term facilitates rapid exploration of the microscopic states during the early stages of adaptive training, while the feed-forward network serves as a flexible correction capturing the complex, non-Gaussian features of the optimal bias potential. The MLP comprises two hidden layers, each containing 100 neurons with ReLU activation functions, trained directly at the target temperature of $T = 500$ K for $10{,}000$ epochs using ADAM[32] with a learning rate of $10^{-3}$. Projecting the trained bias potential onto the primary 2-dimensional subspace $(x_1, x_2)$ by setting the harmonic coordinates $x_{j\geq 3} = 0$ shows close agreement with the exact optimal bias potential of the 2-dimensional system, with residual deviations remaining minimal within the regions traversed by the dominant transition pathways (Figure 7).

During sampling in the 14-dimensional system, the BRW algorithm restricts the path weights to the interval $[1.0, 1.2]$, and each success probability estimate is derived from a batch of 100 paths (Figure 8). Using the unbiased estimator of Eq. (27), we obtain a mean success probability $\langle p_{\rm S}({\rm F})\rangle = (3.3598 \pm 0.0322)\times 10^{-11}$. Concurrently, unbiased Monte Carlo simulations as described in Section 3.3 yield a mean failure-path duration $\langle t_{\rm FF}\rangle = 4.9805 \pm 0.0564$. Combining these, the transition rate is estimated as $r_{\rm FS} = (6.7459 \pm 0.1429)\times 10^{-12}$, in excellent agreement with the exact transition rate of $6.9191\times 10^{-12}$ calculated for the underlying 2-dimensional system. To highlight the criticality of the reweighting procedure, the estimate obtained by neglecting the path scores (Eq. 28) yields a transition rate of $2.0156\times 10^{-12}$, an underestimation by a factor of three. This discrepancy demonstrates that importance sampling of the transition paths, after the approximate bias potential (or committor function) is trained, is crucial for an unbiased rate estimate.

The accelerated framework also captures the relative transition probabilities of the distinct reaction channels in high dimensions: the accelerated Monte Carlo simulations estimate the fraction of the total transition flux passing through saddle point $S_1$ to be $0.2850$, in close agreement with the prediction of $0.2938$ from Kramers' rate theory for the 2-dimensional system (Table 3). This reaffirms that the method correctly captures the difference in the entropic prefactors ($\nu_1$ and $\nu_2$) arising from the local curvature of the potential energy surface at $S_1$ and $S_2$, in addition to the energy barrier difference between the two saddle points.

Table 3. Transition rate $r_{\rm FS}$ and fraction of transitions through saddle point $S_1$ for the 14-dimensional system at $500$ K, from importance sampling with path reweighting (Eq. 27), the exact transition rate of the underlying 2-dimensional system / Kramers rate theory, and the estimate that neglects path scores (Eq. 28, assuming $I \approx I_{\rm opt}$). Neglecting the reweighting underestimates the rate by a factor of three.

5. Conclusion

We developed a machine learning method for importance sampling of rare events in discrete-domain Monte Carlo simulations. We trained a neural network to learn the optimal bias potential, which is proportional to the logarithm of the optimal importance function. Working with the bias potential allows the neural network to deal with the extremely small values of the importance function inherent to rare event problems. We further demonstrated how this learned bias potential can be used to accurately compute the transition rates between metastable states. The method is rigorously validated on 2D and 14D model potentials, showing excellent agreement with available theoretical predictions. This approach provides a unique combination of importance sampling and trainable bias potentials applied to Markov chain Monte Carlo simulations in the discrete domain, and demonstrates a novel way to accelerate rare event simulations in high-dimensional and low-temperature systems where conventional methods are intractable. The next logical step is to apply this framework to atomistic systems, enabling the simulation of long-term processes such as defect evolution in crystals and protein dynamics in solution.

Appendix A. Kramers' Rate Theory

Kramers' rate theory predicts the transition rate from the basin containing A, through the neighborhood of saddle point $S_1$, to the basin containing B as

where $Z(\cdot)$ denotes the local partition functions,

where $\Omega_{\rm A}$ is the basin containing microscopic state A, $\Omega^{\perp}_{S_1}$ denotes the $(d-1)$-dimensional hyperplane orthogonal to the reaction path at the saddle point $S_1$, and $\Gamma_1$ represents the 1-dimensional minimum energy path passing through $S_1$. Comparing Eq. (42) with Eq. (38) identifies the prefactor

where $\mu$ is related to $\nu_0$ through Eq. (22). If one invokes the harmonic approximation, the prefactor reduces to the more familiar form

where $H(\cdot)$ is the local Hessian matrix around a given microscopic state and $\lambda^{*}(S_1)$ is the negative eigenvalue of $H(S_1)$. In Section 4.4 the more rigorous expression, Eq. (44), rather than the harmonic form Eq. (45), is used to evaluate the Kramers rate prediction.

References

1. A. F. Voter, F. Montalenti, and T. C. Germann. “Extending the Time Scale in Atomistic Simulation of Materials”. Annual Review of Materials Research 32 (2002), 321–346.
2. A. F. Voter. “Introduction to the time scale problem”. In: Proceedings of the International Conference on Computational Nanoscience. 2002.
3. A. F. Voter. “Introduction to the Kinetic Monte Carlo Method”. In: Radiation Effects in Solids. Ed. by K. E. Sickafus, E. A. Kotomin, and B. P. Uberuaga. Dordrecht: Springer Netherlands, 2007, 1–23.
4. A. F. Voter. “Hyperdynamics: Accelerated Molecular Dynamics of Infrequent Events”. Physical Review Letters 78.20 (1997), 3908–3911.
5. S. Izrailev et al. “Steered Molecular Dynamics”. In: Computational Molecular Dynamics: Challenges, Methods, Ideas. Ed. by P. Deuflhard et al. Berlin, Heidelberg: Springer, 1999, 39–65.
6. A. Laio and M. Parrinello. “Escaping free-energy minima”. Proceedings of the National Academy of Sciences 99.20 (2002), 12562–12566.
7. A. Laio and F. L. Gervasio. “Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science”. Reports on Progress in Physics 71.12 (2008), 126601.
8. Y. Wang et al. “Frequency adaptive metadynamics for the calculation of rare-event kinetics”. The Journal of Chemical Physics 149.7 (2018), 072309.
9. G. A. Huber and S. Kim. “Weighted-ensemble Brownian dynamics simulations for protein association reactions”. Biophysical Journal 70.1 (1996), 97–110.
10. D. M. Zuckerman and L. T. Chong. “Weighted ensemble simulation: review of methodology, applications, and software”. Annual Review of Biophysics 46 (2017), 43–57.
11. R. J. Allen, C. Valeriani, and P. R. Ten Wolde. “Forward flux sampling for rare event simulations”. Journal of Physics: Condensed Matter 21.46 (2009), 463102.
12. C. Dellago, P. G. Bolhuis, and P. L. Geissler. “Transition Path Sampling”. In: Advances in Chemical Physics. John Wiley & Sons, 2002, Chap. 1, 1–78.
13. E. Vanden-Eijnden. “Transition Path Theory”. In: Computer Simulations in Condensed Matter Systems: From Materials to Chemical Biology Volume 1. Ed. by M. Ferrario, G. Ciccotti, and K. Binder. Berlin, Heidelberg: Springer, 2006, 453–493.
14. Weinan E. and E. Vanden-Eijnden. “Towards a Theory of Transition Paths”. Journal of Statistical Physics 123.3 (2006), 503–523.
15. P. Metzner, C. Schütte, and E. Vanden-Eijnden. “Transition path theory for Markov jump processes”. Multiscale Modeling & Simulation 7.3 (2009), 1192–1219.
16. H. Kahn and A. W. Marshall. “Methods of Reducing Sample Size in Monte Carlo Computations”. Journal of the Operations Research Society of America 1.5 (1953), 263–278.
17. P. W. Glynn and D. L. Iglehart. “Importance Sampling for Stochastic Simulations”. Management Science 35.11 (1989), 1367–1392.
18. W. Cai et al. “Importance sampling of rare transition events in Markov processes”. Physical Review E 66.4 (2002), 046703.
19. P. L’Ecuyer, M. Mandjes, and B. Tuffin. “Importance Sampling in Rare Event Simulation”. In: Rare Event Simulation using Monte Carlo Methods. John Wiley & Sons, 2009, Chap. 2, 17–38.
20. P. Dupuis and H. Wang. “Importance Sampling, Large Deviations, and Differential Games”. Stochastics and Stochastic Reports 76.6 (2004), 481–508.
21. M. de Koning et al. “Adaptive importance sampling Monte Carlo simulation of rare transition events”. The Journal of Chemical Physics 122.7 (2005).
22. M. H. Kalos and P. A. Whitlock. Monte Carlo Methods. John Wiley & Sons, 2009.
23. A. F. Voter. “Introduction to the kinetic Monte Carlo method”. In: Radiation Effects in Solids. Springer, 2007, 1–23.
24. A. P. J. Jansen. An Introduction to Kinetic Monte Carlo Simulations of Surface Reactions. Vol. 856. Springer, 2012.
25. C. C. Battaile. “The kinetic Monte Carlo method: Foundation, implementation, and application”. Computer Methods in Applied Mechanics and Engineering 197.41–42 (2008), 3386–3398.
26. J. Vorba and J. Křivánek. “Adjoint-driven Russian roulette and splitting in light transport simulation”. ACM Transactions on Graphics (TOG) 35.4 (2016), 1–11.
27. J. Arvo and D. Kirk. “Particle transport and image synthesis”. In: Proceedings of the 17th Annual Conference on Computer Graphics and Interactive Techniques. 1990, 63–66.
28. A. Rath et al. “EARS: efficiency-aware Russian roulette and splitting”. ACM Transactions on Graphics (TOG) 41.4 (2022), 1–14.
29. H. Li et al. “A semigroup method for high dimensional committor functions based on neural network”. In: Mathematical and Scientific Machine Learning. PMLR, 2022, 598–618.
30. J. Yuan et al. “Optimal control for sampling the transition path process and estimating rates”. Communications in Nonlinear Science and Numerical Simulation 129 (2024), 107701.
31. A. R. Mitchell and G. M. Rotskoff. “Committor Guided Estimates of Molecular Transition Rates”. Journal of Chemical Theory and Computation 20.21 (2024), 9378–9393.
32. D. P. Kingma. “Adam: A method for stochastic optimization”. arXiv preprint arXiv:1412.6980 (2014).
33. H. A. Kramers. “Brownian motion in a field of force and the diffusion model of chemical reactions”. Physica 7.4 (1940), 284–304.
34. H. Eyring. “The Activated Complex in Chemical Reactions”. The Journal of Chemical Physics 3.2 (1935), 107–115.