Learning on-the-fly quasi-stationarity diagnostics from low-dimensional surrogates

1 Introduction

• Metastable trapping creates a timescale separation that makes long trajectories computationally challenging, especially when integration steps are much shorter than transition times.The molecular-dynamics version of this challenge is known as the timescale problem.Source: §1 · p. 3

1.1 Local description of metastability

• A quasi-stationary distribution is the long-time law conditioned on remaining in a metastable state, and its exit-time and exit-location properties support accelerated dynamics.Under the QSD, the exit time is exponentially distributed and independent of the exit point.Source: §1.1 · pp. 3–4, §1.2 · p. 4

1.2 Accelerated simulation of metastable dynamics

• Parallel replica dynamics simulates independent QSD-started replicas until the first exit, with expected speedup N under perfect parallel efficiency.Its only approximation is truncation in the Yaglom limit.Source: §1.2 · pp. 4–5

1.3 Challenges in AMD methods

• On-the-fly diagnostics may avoid fixed decorrelation times, but the Gelman–Rubin approach is heuristic, requires a threshold, and can trade efficiency against exit-event bias.Overestimating decorrelation wastes acceleration; a systematically lax diagnostic can bias exit events and long-time kinetics.Source: §1.3 · pp. 5–7

1.4 Contributions of this work and outline

• The method uses simulation-based inference over low-dimensional surrogate processes and an RNN to define an on-the-fly decorrelation stopping time.The RNN replaces the MCMC-inspired diagnostic used with Fleming–Viot particle approximations.Source: §1.4 · p. 7
• The paper proves non-asymptotic estimates for coarse-grained conditioned flows of finite-state killed jump processes, without assuming reversibility.The estimates also yield approximate convergence to the Yaglom limit from convergence of the coarse-grained flow.Source: §1.4 · p. 7
• In an adiabatic slow-fast regime, the coarse-graining estimates hold with explicit rates.The regime is defined by a slow-fast decomposition of the process generator.Source: §1.4 · p. 7
• The RNN diagnostic consistently outperforms the Gelman–Rubin diagnostic when its assumptions hold, with experiments also examining failure modes and a mitigation strategy.The evaluation includes experiments of increasing complexity, including a molecular system.Source: §1.4 · p. 8

• A collective variable maps the full process state to a lower-dimensional space, where a Markov surrogate approximates the conditioned marginals of the generally non-Markovian mapped process.Conditioning on survival makes these probability flows nonlinear, unlike the unconditioned linear-flow setting.Source: §2 · p. 8
• The framework focuses first on finite-state jump processes, while extending it to diffusions requires non-trivial technical adaptations.The true process need not be reversible; detailed balance is not required.Source: §2 · p. 8
• Under a formal slow-variable assumption, non-asymptotic total-variation bounds control coarse-graining bias uniformly in time through a parameter measuring slow-fast timescale separation.The analysis bounds the surrogate-to-mapped-process error and then the mapped-process-to-true-process error.Source: §2 · p. 9

• The processes considered are finite-state jump processes with absorbing killing boundaries, and conditional irreducibility requires connectivity among non-absorbing states and to the boundary.Under this condition, the killing time is almost surely finite.Source: §2.1 · p. 9
• Killing yields a sub-Markovian semigroup, while normalizing its evolution gives the nonlinear flow of the law conditioned on survival.The section introduces the killed semigroup through survival-weighted expectations.Source: §2.1 · pp. 9–10
• For each process, the section identifies a unique principal eigentriplet and names its associated probability measure the quasi-stationary distribution.The stated spectral properties include exponential semigroup decay.Source: §2.1 · pp. 10–11

• Fiber measures define averaging operators that project observables onto functions constant on collective-variable fibers, and compatible initial conditions align the surrogate and true mapped laws.The choice of fiber measures remains a parameter in the proof framework.Source: §2.2 · pp. 11, 12, 14
• Coarse-graining bias measures the total-variation difference between the surrogate’s conditioned flow and the mapped conditioned flow of the true process.The mapped flow is generally non-Markovian, whereas the surrogate flow is Markovian.Source: §2.2 · p. 11
• A backward-error identity expresses the bias through time-dependent probability flows that combine coarse-grained evolution with exact evolution and remain normalized.The identity provides the basis for subsequent bounds on coarse-graining bias.Source: §2.2 · p. 12
• The commutator defect measures intrinsic non-Markovianity of the mapped process, while the surrogate defect measures how well the surrogate intertwines with the true dynamics.The commutator contribution is independent of the surrogate definition; zero surrogate defect corresponds to an intertwining relation.Source: §2.2 · pp. 13–14

2.3 A canonical choice of Z

• The fiber measures are defined by disintegrating the true process’s QSD with respect to the collective variable, and the surrogate is a killed jump process built from that coarse-graining.Its generator gives fiber-to-fiber jump rates by averaging the true process’s rates under the QSD.Source: §2.3 · p. 14
• The construction relates to effective-dynamics definitions based on conditional averages, but uses the local equilibrium measure νX rather than a global invariant measure.The text describes this as an analogous construction with respect to νX.Source: §2.3 · pp. 14–15
• The QSD is invariant under fiber-averaging, and the closure defect is QSD-centered.These are among the stated consequences of the canonical choice.Source: §2.3 · p. 15
• The surrogate’s QSD is explicitly given, and its intertwining defect is controlled by the closure defect.The proof identifies the pushforward of the true QSD as the surrogate QSD, with the same decay eigenvalue.Source: §2.3 · pp. 15–16

2.4 Statement of main results

• Under timescale separation and compatibility, Theorem 5 provides a bias control, quasimode estimates, and a two-sided estimate for all T > 0.The constants are explicit, depend on the initial distribution, and are independent of T.Source: §2.4 · pp. 16–17
• When ∥AX∥/ζ is small relative to the prefactors, Corollary 6 ensures convergence of the surrogate flow toward the QSD as T grows.The diagnostic bound is relevant to quasi-stationarity detection.Source: §2.4 · p. 17
• In the adiabatic regime, ∥AX∥/ζ = O(ε) while the other constants remain uniformly bounded, yielding uniform-in-time O(ε) total-variation convergence and a similarly quantitative diagnostic bound.The result is stated as ε → 0.Source: §2.4 · p. 18
• The main bound favors transient decorrelation accuracy over relative accuracy at long times, while still tending to zero and retaining absolute accuracy as T → +∞.The alternative bound is sharper when T ≫ ζ/∥AX∥ but does not capture timescale separation or the adiabatic limit.Source: §2.4 · p. 18
• The assumptions are difficult to verify in systems with many microstates because they require quantitative spectral information; for complex molecular systems, small ∥AX∥/ζ may remain a working hypothesis.The stated separation condition is described as mild adiabaticity, with the stricter factor of two retained to simplify estimates.Source: §2.4 · p. 18

2.5.2 Slow-fast splitting and proof of (52)

• The bias proof uses timescale separation to split the semigroup into slow and fast variables, combining the resulting bounds through a Duhamel–Dyson identity.The slow component is eliminated from part of the numerator using a linear functional before estimating the fast-variable terms.Source: §2.5 · p. 19, §2.5.2 · pp. 21–24

2.5.1 Preliminary estimates and self-correcting property

• The preliminary estimates show that the error measure’s principal-eigenspace contribution decays faster than its generic spectral bound would suggest.The authors call this faster decay a self-correcting property.Source: §2.5.1 · pp. 20–21

2.5.2 Slow-fast splitting and proof of (52)

• The proof obtains the total-variation bound from the preceding estimate by applying the duality formula and taking a supremum over bounded test functions.The text also notes that the ℓ∞ quasimode estimates do not further control the constants; that would require hY-norm estimates.Source: §2.5.2 · p. 27

2.5.3 Proof of the quasimode estimates (53) and (54)

• The quasimode estimates are derived using the Duhamel–Dyson formula and a Poisson equation on the relevant projected space.The Poisson equation is stated to be well-posed because the principal eigenvalue is simple.Source: §2.5.3 · pp. 27–28

2.5.4 Proof of the two-sided estimate (55)

• The two-sided estimate is proved by using duality, splitting the semigroup, and combining upper and lower bounds.The proof invokes commutation of the generator and principal eigenprojection, along with the Feynman–Kac and Duhamel–Dyson formulas.Source: §2.5.4 · pp. 28–29

2.6 Uniform coarse-graining estimates in the adiabatic limit

• The section studies models whose timescale separation is tuned by a small adiabatic parameter ε.The state space is partitioned into fibers of a collective variable, with dynamics including macroscopic jumps, killing, and fiber-reducible microscopic jumps.Source: §2.6 · p. 29
• Under the stated compatibility condition, Theorem 7 gives ∥AX∥/ζ = O(ε) while the other constants remain uniformly bounded as ε → 0.The theorem verifies the assumptions of Theorem 5 in this regime.Source: §2.6 · p. 30
• The diagnostic estimate requires a positive burn-in time to control a boundary-layer term.The proof bounds this term using a time threshold proportional to ε|log ε|.Source: §2.6 · p. 31

2.7 Particle approximations of non-linear semigroups

• For large systems without a convenient sparse generator, directly propagating the conditioned flow may be intractable, motivating interacting-particle approximations such as Fleming–Viot.The Fleming–Viot process approximates the nonlinear flow using a branching particle system.Source: §2.7 · p. 31
• When a replica reaches the boundary, its state is copied from a uniformly sampled surviving replica, while its past trajectory is retained.Replicas otherwise evolve independently; this branching construction yields the particle approximation of the nonlinear flow.Source: §2.7 · pp. 31–32
• The construction can also be attempted for cadlag one-particle dynamics when the state space is not finite or the generator does not generate a jump process.Source: §2.7 · p. 32
• Propagation-of-chaos results for general one-particle dynamics establish approximation at fixed times, while uniform-in-time results are cited for jump processes.Source: §2.7 · p. 32
• Uniform-in-time approximation for killed diffusion processes remains open in general, although partial results are available.Source: §2.7 · p. 33

3.2 Simulation-based inference of decorrelation times

• The method estimates posterior probabilities of decorrelation from Fleming–Viot observations, training on synthetic surrogate flows and deploying on collective-variable traces without simulating the high-dimensional process during training.The surrogate dynamics are marginalized as nuisance parameters; deployment uses featurized observations of high-dimensional particle approximations.Source: §3 · p. 33, §3.2 · p. 35

3.1 Gelman–Rubin diagnostic

• The Gelman–Rubin stopping rule can be overly conservative because its variance ratio converges algebraically to 1 while the conditioned flow may converge exponentially fast.The proposed method avoids assumption (A2), replacing it with prior knowledge of plausible surrogate processes.Source: §3.1 · p. 34

3.2.1 Honest decorrelation times

• An honest decorrelation time is the first time the conditioned distribution lies within a fixed total-variation tolerance of the quasi-stationary distribution; the surrogate criterion is motivated by coarse-graining results.The honest time depends on the dynamics and initial distribution, not on sample paths, and is generally intractable because it requires the conditioned distribution and QSD.Source: §3.2.1 · pp. 35–36

3.2.2 Synthetic decorrelation ensembles

• The synthetic ensemble places a prior on plausible surrogate generators, while its likelihood is intractable and its finite-particle observations are noisy and biased.The prior should be calibrated with domain knowledge, and the synthetic generators must allow tractable offline computation of ground-truth labels.Source: §3.2.2 · pp. 37, 39

3.2.3 Recurrent model

• Recurrent neural networks support causal, variable-length streaming inference with O(dhist) memory and O(ℓd^2 hist) total inference cost.Their per-step inference cost is constant, making them suitable for long traces.Source: §3.2.3 · pp. 39–40
• CNN features from spatially structured histograms feed LSTM layers, which receive Nτ as their sole design feature; a joint-observation alternative performed poorly.Using Nτ is intended to help filter finite-sample noise while avoiding direct conditioning on elapsed simulation time; that design choice was not tested systematically.Source: §3.2.3 · p. 41

3.2.4 Training objective

• The randomized-horizon objective avoids requiring a fixed maximum diagnostic duration, and its ideal predictor conservatively underestimates successful-decorrelation probabilities.When the horizon exceeds the physical decorrelation time, the training rule yields a fairly balanced set of positive and negative observations, scales simulation cost with that time, and avoids long stationary simulations.Source: §3.2.4 · pp. 42–43
• The ideal predictor keeps prior-averaged false-decorrelation risk below α, conditional on a finite detection time.The randomized horizon must satisfy τL > t0(LZ, Z0) almost surely for the stated conservative-probability guarantee.Source: §3.2.4 · p. 43

3.2.5 On-the-fly diagnostic

• During simulation, the diagnostic accumulates collective-variable bin counts over each lag interval, updates its latent state, and uses a user-defined probability threshold to define decorrelation time.Applying it to multiple collective variables is only heuristic unless the boundary compatibility condition holds simultaneously for all of them.Source: §3.2.5 · p. 45

3.3 Discussion of underlying assumptions

• B1 requires a slow collective variable, while B2 requires Fleming–Viot histogram trajectories to converge uniformly in time as N grows.Source: §3.3 · p. 45
• B3 assumes histogram trajectories predict decorrelation time on average over the synthetic prior, while B4 requires prior mass near the deployment model.Source: §3.3 · p. 45
• For finite-state jump processes, B1's approximation is established, but its extension to boundary-killed diffusions remains open and general-diffusion support for B2 is insufficient.Source: §3.3 · p. 46
• RNN training fits a latent model of finite-N error across N and τ; performance improves with N, but this does not remove the need for B2.Source: §3.3 · p. 46
• Potential remedies are coarse-graining at the Fleming–Viot level or assuming B2 and choosing N large enough for the mean-field description to be accurate.Source: §3.3 · p. 46
• Regardless of B1–B4, the stopping rules are well-defined and can be assessed empirically for exit bias, decorrelation-time overestimation, and parameter sensitivity.Source: §3.3 · p. 46
• When B1–B4 hold, α bounds prior-averaged false-convergence risk and average exit-event bias, enabling a priori target-error selection unlike the variance-ratio threshold.Source: §3.3 · p. 47
• Treating a decorrelated Fleming–Viot population as independent QSD samples remains a caveat, though propagation of chaos implies weak correlations.Source: §3.3 · p. 47

4.1 Validation of theoretical predictions

• Prefactor terms converge to finite values as timescale separation ζ increases, and the observed dependence of the constants on ζ is mild.The supremum-in-time total-variation bias agrees well with the coarse-grained description, with lower bias for fiber-equilibrated initial conditions; long-time decay matches the predicted rate, though the bounds are not sharp.Source: §4.1 · p. 48

4.2 Implementation of the RNN diagnostic

• For molecular applications, the implementation assumes a bounded one-dimensional collective variable and uses killed diffusion surrogates when that variable is slow.The surrogate ansatz motivates construction of the synthetic prior; extending the coarse-graining results to killed diffusions remains to be justified rigorously.Source: §4.2 · pp. 48, 51

4.2.1 Construction of the prior ensemble

• The synthetic prior spans decorrelation times across several orders of magnitude, including metastable decorrelation at low temperatures.The prior uses killed nearest-neighbor jump-process discretizations of one-dimensional surrogate diffusions on a 100-point grid.Source: §4.2.1 · pp. 51–52

4.2.2 Synthetic data generation

• Synthetic Fleming–Viot training trajectories use discretized diffusion dynamics, introducing a discretization bias whose rigorous control is outside this work’s scope.The authors expect this bias to vanish as grid spacing and timestep decrease and to be partly filtered by histogram featurization.Source: §4.2.2 · p. 54

4.2.3 Architecture and training details

• The RNN was trained on 375,000 feature sequences from 37,500 physical landscapes after architecture tuning on a 1,000-sequence dataset.Training used 150 epochs on five CPU cores for approximately 11 hours.Source: §4.2.3 · p. 55

4.3.1 Binary classification performance

• The RNN’s ROC curve dominates Gelman–Rubin’s, with AUCs of 0.96 and 0.82, respectively; the total-variation diagnostic scores 0.72.The ordering holds after aggregating over design parameters; access to the reference QSD and particle histograms alone is insufficient to predict decorrelation accurately under finite-particle noise.Source: §4.3.1 · p. 57
• Across the tested parameter grid, the worst RNN AUC exceeds the best Gelman–Rubin AUC, and both diagnostics are more sensitive to replica count than lag time.RNN and Gelman–Rubin performance improves with more replicas, consistent with a non-negligible finite-particle approximation error.Source: §4.3.1 · p. 57

4.3.2 On-the-fly quasi-stationarity detection

• The on-the-fly comparison tests both diagnostics at α = 0.1 and α = 0.05, counting trajectories with no declaration by 10t0 as failures.The simulations are drawn from the same training-distribution set used in the preceding classification experiment.Source: §4.3 · p. 55, §4.3.2 · p. 57
• Within the synthetic training prior, RNN decorrelation times clustered near the reference, while Gelman–Rubin times were broader, later, and more threshold-sensitive.Gelman–Rubin failed on 9.7% of trajectories at α = 0.1 and 50% at α = 0.05; RNN failure rates were 0.3% and 1%, respectively.Source: §4.3.2 · p. 58, Figure 8 · p. 60

4.3.3 Posterior risk calibration and prior misspecifications

• Across synthetic prior shifts, RNN AUC exceeded Gelman–Rubin AUC, and RNN empirical false-decorrelation risk remained controlled by α.The shifted-prior calibration was generally conservative, with mild miscalibration at low temperatures; excluding metastable substates largely corrected that overconfidence.Source: §4.3.3 · p. 61, Figure 10 · p. 62

4.4.1 Dimer system and collective variables

• The solvated-dimer experiment monitored the intermediate conformation of a 512-atom system using dimer length ξ and solvent coordination ncoord.The system’s potential induces compact, intermediate, and stretched conformations; the free-energy comparison with training-prior samples was described as heuristic evidence of calibration.Source: §4.4.1 · pp. 61, 63, 64

4.3.3 Posterior risk calibration and prior misspecifications

• On the training prior, RNN probabilities were mildly under-confident but generally well-calibrated, and empirical false-decorrelation risk stayed below α on average.The study recorded 1.24 × 10^6 predictions; calibration was especially good at high predicted probabilities.Source: §4.3.3 · p. 60, Figure 9 · p. 62

4.4.2 Molecular dynamics

• Overdamped Langevin dynamics varied dimer mobility through Dd, with Dd = 0.05 representing a slow ξ and Dd = 1 a fast ξ.The authors associate the low-mobility case with a heavy molecule in solution and the high-mobility case with a light solute.Source: §4.4.2 · p. 65

4.4.3 Reference decorrelation times and scale separation

• Reference decorrelation times depended on initial configuration; ξ predicted decorrelation in all tested Dd = 0.05 cases, but converged before ncoord for Dd = 1.Scale separation was at most mild even at Dd = 0.05, so the system was not deep in the adiabatic regime.Source: Figure 13 · p. 67, §4.4.3 · p. 67

4.4.4 On-the-fly decorrelation detection

• For Dd = 1, ξ’s non-slow behavior may produce non-Markovian marginal inputs outside RNN training, making its predictions unreliable.The authors note that Gelman–Rubin’s better performance in this case may reflect system-dependent links between its threshold α and false-decorrelation risk.Source: §4.4.4 · pp. 68–69
• The RNN outscored Gelman–Rubin in ROC AUC for both dimer mobilities, but was miscalibrated at Dd = 1, where ξ is not slow.At Dd = 0.05, the RNN was calibrated for empirical false-decorrelation risk and never failed within the tested horizon; at Dd = 1, its empirical false-decorrelation risks were 57% (α = 0.1) and 41% (α = 0.05).Source: §4.4.4 · p. 68, Figure 14 · p. 69

4.4.5 Bias on the exit event

• Exit-time distributions were indistinguishable from the reference for all diagnostics; at Dd = 0.05, all four diagnostics showed no discernible exit-distribution bias.For Dd = 1, RNN exit probabilities from initial configurations 1 and 4 at α = 0.1 slightly exceeded the target error range, consistent with false decorrelations; Gelman–Rubin showed no discernible bias.Source: §4.4.5 · p. 70, Figure 15 · p. 71
• For the reference, RNN, and Gelman–Rubin populations, exit times showed no dependence on exit point for either Dd value in the tested lax setting α = 0.1.The study expected the more conservative α = 0.05 setting to reduce any bias.Source: §4.4.5 · p. 70, Figure 16 · p. 72

4.4.6 Using multiple collective variables

• With ξ and coordination number, the RNN had zero empirical false-decorrelation risk and a 0.1% failure rate at α = 0.05 for Dd = 1, with no detectable exit-distribution bias.The Gelman–Rubin diagnostic remained more conservative, failed more often, and had slightly lower ROC-AUC; the two-variable mitigation is heuristic because the coordination-number marginals are not well represented by the training prior.Source: §4.4.6 · p. 73, Figure 17 · p. 74

4.4.7 Computational cost and summary statistics

• One RNN inference on a histogram feature costs about one dimer MD step on a laptop, and its cost is independent of atom count.This is anecdotal rather than a systematic benchmark; because relative inference cost decreases with lag time, the authors expect negligible deployment cost for much larger systems.Source: §4.4.7 · pp. 73, 75

Deployability in accelerated molecular dynamics algorithms

• The RNN requires only periodic observations of a Fleming–Viot process in the collective variable, and its inference cost is expected to be negligible relative to interatomic-force evaluation.The cost comparison is made even for the small molecular system considered.Source: §5 · p. 76

Physical information as a modeling input

• The synthetic prior places the method between MCMC-inspired diagnostics and deterministic decorrelation-time definitions, with experiments showing discriminative power and robustness across tested parameter choices.This middle ground may suit applications where a broad range of effective models is known.Source: §5 · p. 76

Interpretability

• The threshold α has a probabilistic interpretation under a well-calibrated prior, but that guarantee is heuristic in practice and is lost under prior misspecification.Under the stated conditions, α controls average false-decorrelation risk over the prior, which is linked to bias in the sampled exit distribution.Source: §5 · p. 76

Failure under misspecification

• Without an out-of-distribution detector, the RNN may confidently declare decorrelation when its assumptions fail; a proposed heuristic takes the most conservative prediction across multiple collective variables.The stated failure example uses a collective variable that is not slow for the system.Source: §5 · p. 76

Choice of collective variables

• The experiments use one-dimensional collective variables, while grid-based training is expected to remain feasible for dimensions 2–4 at most under matrix-free computation and coarser grids.The collective variable need only describe decorrelation within one metastable state, so it may require fewer effective degrees of freedom than the global landscape.Source: §5 · p. 77

Coarse-graining for killed diffusions

• The coarse-graining analysis currently covers finite-state jump processes; extending it to killed elliptic diffusions, then hypoelliptic diffusions, remains future work.The proposed elliptic-diffusion extension could also address a small-parameter adiabatic regime.Source: §5 · p. 77

Clustering of σ(T(ε))

• The proof starts from a reducible operator whose zero eigenvalue has multiplicity equal to the number of fibers, then separates the perturbed spectrum into slow and fast clusters.The reduced operator describes macroscopic jumps between fibers with rates averaged under microscopic measures.Source: A Proof of Theorem 7 · pp. 80–81
• Kato perturbation theory tracks the clustered eigenvalues and projectors, while positivity and conditional irreducibility support the spectral characterization of the reduced dynamics.The proof uses the fiber decomposition to represent the first-order reduced perturbation.Source: A Proof of Theorem 7 · pp. 80–81

Uniform version of Proposition 1

• Semigroup estimates applied to the fast spectral component and the reduced operator establish uniform control of the relevant decay bounds and constants for sufficiently small ε.The proof verifies a uniform spectral gap and prefactor, and notes that the gap estimate could be sharpened.Source: A Proof of Theorem 7 · p. 83

B.1 Synthetic prior construction

• A synthetic dynamics draw samples physical parameters U and D and an inverse temperature β from a uniform range.The energy profile is built from a Gaussian well, internal Gaussian barriers, and a linear correction before rescaling its range to [0, 1].Source: §B.1 · p. 86
• The diffusion coefficient begins with a randomly sampled constant value and is either kept constant or perturbed by Gaussian profiles.The choice between constant and perturbed diffusion is made with probability 1/2.Source: §B.1 · p. 86
• The initial state is sampled uniformly from the grid, while the training and misspecified priors differ in the ranges of β, s, and δ_max.The other parameter distributions are fixed across these priors.Source: §B.1 · p. 86, Table 2 · p. 87

B.2 Network architecture

• The deployed RNN architecture is detailed in Table 3.The passage does not enumerate the architecture components outside the table reference.Source: §B.2 · p. 87, Table 3 · p. 87
• The convolutional layers use a kernel width of 5 and zero-padding size of 2.Source: Table 3 · p. 87