Constructing auxiliary dynamics for nonequilibrium stationary states by variance minimization

Large deviation theory (LDT) is a framework to extend the formalism of equilibrium statistical mechanics to nonequilibrium systems.¹ Much of LDT is concerned with summarizing the dynamics of the system, as expressed via the fluctuations of typical trajectories. Ensembles of rare trajectories display a fascinating behavior reminiscent of phase transitions and criticality. Recent work illustrates such a dynamical behavior both in lattice systems, such as simple exclusion processes,^2–9 constrained kinetic models,^10,11 models of self-assembly,^12–14 and dissipative hydrodynamics,^15,16 and in continuum systems in the form of driven or active Brownian particles,^14,17–21 as well as open quantum systems.^22–24 Recently, LDT has also been shown to offer a route to calculating nonlinear transport coefficients in model and realistic systems.^25,26

Accessing properties of interest, in all but the simplest systems, requires numerical tools. In the context of LDT, this takes the form of sampling techniques such as the cloning algorithm or diffusion Monte Carlo (DMC)^27–32 and transition path sampling (TPS) or path integral Monte Carlo,^33,34 or representing the non-equilibrium distribution with an explicit ansatz, for example, matrix product states.^9,11,35 The primary challenge in sampling methods is the problem of exponential variance in the estimator for the large deviation function as a function of the bias. Several approaches have been suggested to ameliorate this variance problem that can be interpreted as different forms of importance sampling. Nemoto et al. demonstrated the use of control forces that can be parameterized and optimized via a feedback loop.^10,32 Klymko et al. formulated a probability-conserving reference model to umbrella sample trajectories for certain growth processes.¹³ A different approach has been formulated by Jacobson and Whitelam to directly compute the rate function via an optimizable variational ansatz.³⁶ In previous work, we showed that guiding distribution functions (GDFs), as introduced for quantum DMC calculations,³⁷ define an auxiliary dynamics that can be used to importance sample computations of large deviation functions under bias²⁰ and highlighted the connection to the generalized Doob’s transform.^38–40 In the current work, we describe a practical numerical technique to generate good guiding distribution functions in both lattice and continuum simulations of large deviation functions using the idea of variance minimization. This again draws from the quantum field and, in particular, the methods of variational Monte Carlo (VMC).⁴¹ There has been other recent work on optimizing auxiliary dynamics by using a variational estimator to determine a GDF for a continuum system;⁴² our work provides a different perspective based on different techniques that are less computationally expensive than performing a DMC based optimization.

In the current work, the quantity of interest in large deviation theory is the cumulant generating function (CGF) ψ(λ) which is analogous to the free energy of equilibrium statistical mechanics. It is computed via an ensemble average over trajectories given by

where λ is a field conjugate to the observable $O=∑t=1tNo(Ct+,Ct−)$⁠, with o being an arbitrary function of configurations (⁠$Ct$⁠) at adjacent times, t+ and t−, and t_N is the final trajectory time. $P[C(tN)]$ is the likelihood of a given trajectory $C(tN)={C0,C1,…,CtN}$⁠. The trajectories are generated by the master equation $∂tpt(C)=Wpt(C)$⁠, where $pt(C)$ is the probability of a configuration of the system, $C$⁠, at time t, and $W$ is a linear operator. It can be shown that ψ(λ) is the largest eigenvalue of a “tilted” operator $Wλ$⁠, i.e., $Ξ|Wλ=Ξ|ψ(λ)$⁠, where ⟨Ξ| is the corresponding dominant left eigenvector.⁴³ The effect of the tilt is to reweight the transition probabilities of $W$⁠. In the discrete case, $Wλ(C,C′)=W(C,C′)e−λo(C,C′)(1−δC,C′)−R(C)δC,C′$⁠, where $R(C)=∑C≠C′W(C,C′)$ is the exit rate.

$Wλ$ is not a Markovian operator (i.e., the sum of transition probabilities is not normalized). Consequently, when calculating ψ(λ) via Monte Carlo techniques, it is necessary to track this additional normalization constant or weight, whose variance grows exponentially with |λ|. It is desirable, thus, to instead consider an auxiliary dynamics generated by a modified operator $W̃λ=Ξ^WλΞ^−1$ (related to the generalized Doob’s transformation^18,38) from which the CGF can be obtained as

where the diagonal matrix $Ξ^=∑CΞ̃(C)|CC|$ is constructed from the dominant left eigenvector of $Wλ$ [i.e., $Ξ̃(C)=⟨Ξ|C⟩$], $1|$ is the uniform left vector, and |p₀⟩ is the initial distribution of configurations. The normalization of $W̃λ$ is completely independent of the configuration, and the variance due to the bias is removed; thus, the transformation by $Ξ^$ carries out a form of importance sampling. Although propagation with $W̃λ$ requires the knowledge of the exact eigenvectors, in the GDF approach, we simply approximate these eigenvectors with guiding functions of our own construction and carry out dynamics with $W̃λ$ using the diffusion Monte Carlo algorithm.²⁰ The quality of the GDF importance sampling then depends on the degree of overlap between the approximate and the exact left eigenvector of $Wλ$⁠. The problem is therefore reduced to finding a high quality GDF in order to compute the CGF and its associated cumulants with good statistical efficiency.

As mentioned earlier, we can determine GDFs using ideas that originate from quantum diffusion Monte Carlo calculations, where the analogous problem is to determine a guiding function that best approximates the ground-state of a quantum Hamiltonian.⁴¹ This is termed a variational Monte Carlo (VMC) calculation. While energy minimization is commonly used for this purpose,^44,45 $Wλ$ is not Hermitian, and thus, its spectrum is not necessarily bounded. However, we can use an associated property of eigenstates, viz., the variance of the quantity,

must vanish for eigenstates. This quantity, which we call the local CGF, is the analog of the local energy for which variance minimization has previously been explored in quantum Monte Carlo.⁴¹ The nonequilibrium variance minimization problem thus corresponds to minimizing

where $⟨Λ({p})⟩≡∑Cw(C)Λ({p},C)$ is the average of the local CGF, {p} are the parameters used to characterize the GDF, $Ξ̃({p},C)$⁠, and $w(C)$ is a normalized sampling distribution. The calculation involves minimizing Eq. (4) with respect to {p} over a fixed set of configurations ${C}$⁠. Note that (4) can be sampled from any $w(C)$⁠, but for $w(C)=Ξ̃({p},C)‖Ξ^({p})‖$⁠, we obtain

which is an estimator for ψ(λ) in the sense that the approximate sign is replaced by an equality for a GDF, that is, the exact eigenstate. Despite the absence of a bounded variational principle, this estimator of ψ(λ) in practical terms is often useful, and we refer to this as the VMC estimator for ψ(λ), in complete correspondence to its quantum counterpart. Alternate forms of $w(C)$ could be used to control the statistical error of the VMC estimator σ²({p}, λ), but we have not explored this dependence.

Whereas a strictly zero variance is guaranteed for eigenstates, we emphasize that a smaller (non-zero) variance of the local CGF does not strictly imply a better GDF. A more rigorous metric is the actual reduction of the standard deviation in the estimate of the CGF [ε(ψ)] from DMC (or indeed TPS, which can easily be adapted to use a GDF). Decreasing ε(ψ), thus, defines a better GDF, and this is reflected in statistical independence of samples (the trajectories which are being generated). For the DMC algorithm, the indicator of statistical independence is the fraction of independent walkers f_I, which is the ratio of the number of unique walkers at a given simulation time (t) to the total number of walkers.^20,32 These walkers are ancestors for all walkers at simulation time t′ > t. Generally, walkers at later times are descended from a smaller set of (unique) walkers at previous times, leading to a growth in correlated walkers as the simulation proceeds. Empirically, f_I is seen to be less susceptible to statistical noise than ε(ψ) but (for a given λ) still bears a monotonic relationship with it, and we will, therefore, primarily use f_I as the measure of quality of a GDF. In the case of perfect importance sampling, using the exact auxiliary dynamics, f_I is equal to 1 at all times. In the other limit, if all walkers are correlated, f_I = 0. Because f_I(t) measures the correlation among walkers as a function of time, it must be smallest at t = 0,^14,20 which is what we will report. It is important to note that although improvements in f_I yield improvements in ε(ψ), the relationship between the two is not linear.

To demonstrate our procedure, we carry out simulations on a continuum system and a lattice model. For the continuum system, we consider the prototypical driven Brownian walker, where the observable of interest is the entropy production. This system consists of N = 10 particles (at location R = {r_i}) moving in a periodic potential [V(r) = $v$₀ cos(2πr)] on a ring of size L = 1.0 under the influence of an external driving force (f) and a fluctuating field represented by Gaussian white noise. These particles also interact via a pairwise repulsive force [$f^(ri,rj)=−f^(rj,ri)$ and $|f^|=α⁡exp[−(|ri−rj|/rc)2]$]. The transformed tilted propagator that includes auxiliary dynamics due to the GDF $Ξ̃(R)$ is given by²⁰ 

where $Fi(R)=fx^−∂iV(ri)x^+∑j≠if^(ri,rj)$⁠. The adjoint operator $Wλ†=∑i∂i2+(Fi(R)⋅x^+2fλ)∂i+fλ(fλ+Fi(R)⋅x^)$ represents the norm-breaking term that is handled via branching. Trajectories for $W̃λ$ are generated from Langevin dynamics $∂tri=Fi(R)⋅x^+2fλ+2∂i⁡lnΞ̃(R)+ηi$⁠, where the random force, η_i, satisfies ⟨η_i(t)⟩ = 0 and ⟨η_i(t)η_i(t′)⟩ = 2δ(t − t′). The entropy production s(t) is reflected in the biasing term exp[λt_Ns(t)] with $s(tN)tN=∑i∫0tNfri̇(τ)dτ$ which is absorbed into the dynamics of the tilted Fokker–Planck operator.^14,20,46

For this system, we parameterize the GDF as $Ξ̃(R)=∏iϕ(ri)∏j<i[J(ri,rj)]$⁠, where the single-particle function, ϕ(r_i), is obtained from the non-interacting eigenstate of $W̃λ$ [$f^(ri,rj)=0$] generated with M = 101 plane wave modes²⁰ and $J(ri,rj)=∑k1,k2J̃(k1,k2)eik1r1+ik2r2$⁠. The parameters $J̃(k1,k2)$ are estimated by minimizing (4).

Shown in Fig. 1(a) is the large deviation function and f_I computed using the no guiding function, the non-interacting guiding function (J = 1), and the variational form above (a similar form has been explored independently⁴²). We parameterized J(k₁, k₂) with 241 parameters (21 plane waves per particle), and we minimized the variance such that σ²({p}, λ) < 4.0 × 10⁻³ for all λs we considered. The minimization procedure is started at small |λ| using an initial guess for parameters that produces a uniform state since we know that for λ = 0, the exact ⟨Ξ| is uniform. The optimized parameters are used as an initial guess for the next (nearby) λ. A trace of the minimization is shown in Fig. 2 at λ = −0.5 obtained using a simple simplex algorithm. For this system, the minimization required 10 000 configurations [sampled from the guessed GDF as the distribution $w(C)$] for −0.3 ≤ λ ≤ 0.3 (λ ≠ 0) in order to avoid getting stuck in local minima (which produced a poor GDF). For |λ| > 0.3, 2000 configurations were sufficient.

From the reduction in f_I [inset of Fig. 1(a)], it is evident that continuum calculations can be made much more efficient with appropriate GDFs. As noted earlier, the improvement in f_I estimated at t = 0 can imply greater efficiency gains when the full trajectory space is considered. In Fig. 1(b) (inset), we show the improvement in the corresponding standard deviation in the subsequent DMC calculation, which can be reduced by a large factor, although statistical noise in this measure means that it is difficult to give a precise estimate of the factor. Additionally, although we used the same observation time (t_N) to compute the CGF and cumulants for all types of sampling, we find that the results converge much faster with t_N for simulations done with auxiliary dynamics.

We now consider an interacting non-equilibrium problem on a lattice, namely, the current fluctuations of a periodic weakly asymmetric simple exclusion process (WASEP).⁴⁷ The WASEP models the transport of N particles on a lattice with L sites. Here, N is chosen as 0.3L (rounded to the nearest integer). The configuration of the particles is defined by a set of occupation numbers with hard-core constraints, n_i = {0, 1}, e.g., $C={0,1,…,1,1}$⁠. The tilted propagator for a current bias λ,

yields particles hopping to the right with rate $p=12eE/L$ and to the left with rate $q=12e−E/L$⁠, where the 1/L factor in hopping gives the weakly asymmetric limit, whose large scale behavior (universality class) is described by the Edwards–Wilkinson equation.^4,48 Here, $b^i†$ (⁠$b^i$⁠) creates (destroys) particles on site i, and n_i (⁠$1−ni$⁠) counts the number of particles (holes). For subsequent calculations, we set E = 10 and use periodic boundary conditions, i = L + 1 → 1 and i = 0 → L.

Unlike the continuum system, where the soft-core interaction means the non-interacting solution is a sensible starting point to construct the GDF, the hard-core interaction requires a different treatment. Here, we consider a GDF that is a product purely of n-particle correlation factors, e.g., for n = 2, $Ξ̃(C)=∏p<qJ2(rp,rq)$⁠, where r_p, r_q denote the positions of the particles in the configuration $C$ and for n = 3, we use J₃(r_p, r_q, r_s), and the variational parameters are the values [i.e., J₂(r_p, r_q)]. This form is sometimes referred to as a correlator product state (CPS) in quantum systems.^50,51 To enforce periodic boundary conditions (PBC), we use J₂(r_pq), J₃(r_pq, r_qr, r_qs) where the inter-particle distances are defined with a minimum image convention, i.e., r_pq = min(|r_p − r_q|, L − |r_p − r_q|), and J₂, J₃ are symmetric under cyclic permutations of their arguments (cyclic permutation symmetry, rather than full symmetry, was used to reflect the handedness of hopping around in the model). To minimize the number of variational parameters in the large calculations below, for r_pq > R_cut, we set J₃ = 1, where R_cut is a cutoff distance. For comparison, we have also considered a GDF of the cluster mean-field (CMF) form described in the work by Ray et al.²⁰ 

Figure 3(a) shows a comparison of f_I as a function of λ for a L = 16 model with different GDFs. In the WASEP, the short-range CMF is unable to capture the long-range correlations present in the system, and therefore, we do not get much improvement in efficiency using this GDF. However, with the correlator product state, we can obtain large improvements [e.g., we see from Fig. 3(b) that using J₃, ε(ψ) is improved by two orders of magnitude].

In order to illustrate the flexibility of this ansatz, we have further performed calculations for different system sizes L = 24–96. Due to the reduction in standard deviation, despite the large system size, we needed only a modest number of walkers in the DMC procedure N_w = 10 000–120 000, which was sufficient to estimate the susceptibility $χ(λ/L)=d2ψ(λ/L)d(λ/L)2$ (computed as a correlation function). For these calculations, in order to reduce the number of parameters, we used a cutoff distance of R_cut = 16 for J₃. The variance minimization was carried out using 4000–8000 fixed configurations sampled from the GDF to optimize 85–155 parameters depending on the system size. The error of the VMC estimator for ψ(λ), as a measure of the GDF quality, is shown in Fig. 4(a), where we see that the relative error grows with |λ| as more particles become correlated. Note that for these calculations, we did not spend a lot of effort to optimize the more extreme λ values as we were interested in the range over which the susceptibility peaks, where the VMC error is <5%. In this model, the trajectories undergo a continuous phase transition from a uniform state to a traveling wave at some λ_c, provided that $E>Ec=π/ρ(1−ρ)$⁠, where ρ = N/L.^3,49 This can be detected from the growing susceptibility of the system shown in Fig. 4(b). Macroscopic Fluctuation Theory (MFT) suggests that $λc=−E+E2−Ec2=−2.7198$ (for ρ = 0.3). However, quantitative accuracy of MFT is not guaranteed, and although a Bethe ansatz formulation exists, the precise value of λ_c for L → ∞ is not explicitly known.⁴ Our DMC calculations suggest a critical point near λ_c ∼ −3.9 for the largest system sizes (L = 64 and 96). We note that our results for ψ(λ/L) are in excellent agreement with the largest system size (L = 64) considered in the work of Espigares et al.,⁴⁹ also plotted in the inset of Fig. 4(a).

In this article, we showed how to compute guiding distribution functions using the technique of variance minimization originating in variational Monte Carlo calculations of quantum systems. This provides a systematic route to statistically efficient Monte Carlo computation of large deviation functions using the diffusion Monte Carlo or cloning algorithm, which we demonstrated in the continuum for the problem of the Brownian walker, as well as on the lattice for the WASEP model. The very general nature of variance minimization means that the possibilities for different guiding distribution functions are limited only by one’s imagination and they can be adapted to very complex systems. The form of the GDF should typically be motivated by properties of the system. For instance, we have demonstrated that for the continuum system, the repulsive nature of the interactions can be accommodated by a Jastrow type of GDF. Similarly, the long-range order present in the WASEP case can be accounted for by the CPS form. Other types of forms may be motivated by strategies employed in quantum ground state calculations, as the Backflow ansatz,⁵² and, additionally, different forms can be combined. Finally, we note that obtaining a good form for the guiding distribution function, much like obtaining a compact wavefunction in a quantum problem, is likely to provide important analytic insights into the behavior of the non-equilibrium system of interest.

The authors would like to thank Rob Jack, Vivien Lecomte, Juan P. Garrahan, and David Limmer for fruitful and engaging discussions. U.R. was supported by the Simons Collaboration on the Many-Electron Problem and the California Institute of Technology. G.K.-L.C. is a Simons Investigator in Theoretical Physics and was supported by the California Institute of Technology and the U.S. Department of Energy, Office of Science, via Grant No. DE-SC0018140. These calculations were performed with CANSS, available at https://github.com/ushnishray/CANSS.

In the main text, we described how the computation of the cumulant generating large deviation function (CGF) can be greatly enhanced with the appropriate form of auxiliary dynamics. The GDF structure that we have introduced in the 1D case is of the form $Ξ̃(R)=∏iϕ(ri)∏j<i[J(ri,rj)]$⁠. To carry out simulations with this GDF, we require the adjoint of the transformed Fokker–Planck operator. For the Brownian walker in the main text, we find

where the non-interacting single-particle operator is given by

where γ_i(r_i) = f − ∂_iV(r_i) is the effective single-particle force acting on a particle. This expression provides the norm-breaking term for a subsequent DMC calculation and can be generalized straightforwardly to higher dimensions.

In the main text, we have mentioned that an improvement of the GDF is indicated by the improvement in the fraction of independent walkers (f_I) estimated from DMC (this would be reflected in a corresponding decrease in the autocorrelation time for TPS), and we reported f_I evaluated for t = 0. Here, we show the full f_I(t) in Fig. 5 for calculations using different types of dynamics at λ = −5.0 for the 1D WASEP system. Note that f_I(1) = 1, but as one goes back in simulation time, the walkers are descended from a smaller and smaller set of ancestors. The GDFs J₂ and $J̃2$ use the same J₂ ansatz in the main text but illustrate the effect of optimizing for more steps in the variance minimization of J₂, leading to an improved f_I(t). This is also reflected in Table I, where the VMC estimator of the CGF appears to converge with an increase in the flexibility of the GDF, while the standard deviation is systematically improved.