Rare events in stochastic dynamical systems tend to cluster around their most likely realization. As a result they have predictable features that can be calculated via some optimization problem. This profound observation has been made in numerous fields, and used for example, to explain phase transitions in statistical mechanics [18], derive Arrhenius' law in chemical kinetics [1], or use semiclassical trajectories in quantum field theory [49]. Large deviation theory (LDT) [52] gives a mathematical justification to these results and provide us with an action, or rate function, to minimize in order to calculate paths of maximum likelihood, also known as instantons. The theory also gives exponential asymptotic estimates of rare event probabilities. While this information is already useful in many cases, more refined estimates are often desirable. These 'prefactor' calculations attempt to quantify the effect of Gaussian fluctuations around the instanton, a notion that has also been separately rediscovered in the literature through various means [5]. For example, in the context of chemical reaction rates, next order refinements of the exponential reaction rate are known as the Eyring-Kramers law [19,38]. Similarly in quantum field theory, perturbing around the semiclassical trajectory, the second order variations leads to a Gaussian path-integral, which ultimately results in an additional contribution in the form of a ratio of functional determinants [49].
Over the last 2 decades, several computational methods have been developed to calculate instantons. Among others, we refer to the string method in the context of gradient flows [15,17], the minimum action method [16,22], the adaptive minimum action method (aMAM) [53] and the geometric minimum action method (gMAM) [34,35,50]. These methods are now efficient enough to be used in the context high-dimensional systems, including stochastically driven partial differential equations arising in fluid dynamics [28,30].
In contrast, surprisingly little work has been done on the numerical side of prefactor calculations (see however [45]). The main objective of this paper is to show how to extend methods such as MAM or gMAM to efficiently estimate prefactors in the context of the calculations of expectations, probabilities, and mean first passage times.
Consider a family of stochastic differential equations (SDEs) for ∈ℝ , with drift vector field ∶ℝ →ℝ and diffusion matrix = ⊤ , where ∈ℝ × ,
Here, is an -dimensional Wiener process, and we have introduced a small parameter >0 to characterize the strength of the noise. For simplicity we will assume that the diffusion matrix ∈ℝ × is positive-definite (hence invertible) and constant (but not necessarily diagonal), that is the case of additive Gaussian noise-the generalization of the methods presented below to a covariance matrix that depends on , that is multiplicative Gaussian noise, is straightforward. Large deviations theory [24,52] indicates that, in the limit as →0 , the solutions to (1.1) that contribute most to the probability of an event or the value of an expectation are likely to be close to the minimizer of the Freidlin-Wentzell rate function subject to appropriate boundary conditions. This action functional is given by
with the Lagrangian
Here ⟨ , ⟩ stands for the Euclidean scalar product between the vectors and and we introduced the norm induced by ,
. If the diffusion matrix is the identity, this norm becomes simply the Euclidean length.
The minimizer of the action (1.2) is referred to as path of maximum likelihood or instanton,and it can be found by solving the corresponding Euler-Lagrange equations ̇ = (1.4) with boundary conditions appropriate to the event under consideration. Thus, in the context of LDT, the leading order estimation of probabilities or expectations can be reduced to the solution of the deterministic system (1.4). Alternatively, there is a Hamiltonian formulation to the problem. Taking the Legendre transform of the Lagrangian and introducing the momentum = ∕ ̇ , we obtain the Hamiltonian
and the Euler-Lagrange equations for the instanton become ̇ = ( )+ , ̇ = -(∇ ( )) ⊤ .
(1.6)
The instanton gives the leading contribution to the exponential decay of the probability for observing a rare event. In order to obtain sharp estimates, one needs to furthermore consider prefactor contributions.
Intuitively these prefactors can be calculated by accounting for the effects of the fluctuations around the instanton , which can be done by linearizing the solution of the SDE (1.1) around and considering = ∇ ( ( )) + .
(1.7)
The solution to this equation defines a Gaussian process, and the prefactor contribution to expectations, probabilities, or mean first passage times can be calculated as specific expectations over this process. In turns, these expectations are ratios between determinants of specific positivedefinite matrices or operators that can be expressed in terms of the solutions of deterministic Riccati equations, as can be intuited by analogy with results from optimal control theory. Our objectives here are to: (i) formulate these Riccati equations, including their boundary conditions, in the specific cases of the calculation of expectations, probabilities, and mean first passage time; and (ii) develop efficient numerical methods for their solution.
From a theoretical point of view, our approach builds on a large corpus of works dealing with expansion beyond the exponential estimate of LDT and the evaluation of quadratic path integrals or Wiener integrals, as initiated by Schilder [44]. Results from probability theory and stochastic analysis in this direction include for example the pioneering works by Kifer [37] and Azencott [2]. On the analytical side, formal asymptotic expansions were used for example in [40,41], the WKB expansion in [21], optimal control theory in [3], and more recently potential theory in [5,12]. This last approach allows one to establish rigorously the Eyring-Kramers law for reversible system [11], and has also been extended to other problems, such as lattice gas models [13]. Similarly, prefactor calculations recently included non-reversible systems [9]. From a mathematical perspective it is much harder to treat the infinite dimensional case, even though recent breakthroughs have been made at least for systems in detailed balance [4,6,7,20]. From a computational viewpoint, none of the references above deal with explicit calculation of the instanton or the prefactors, and the equations they derive for these objects, while often very general, do not lend themselves automatically to numerical implementation. Explicit (especially numerical) computation of prefactors is usually confined to quantum field theory, where the evaluation of Gaussian path integrals is a classical result [23,26,39,43], see also [14] for a recent review. However, the Lagrangian in quantum mechanics is usually assumed to take a special form with no first order time derivative in the Euler-Lagrange equation (corresponding to a stochastic process in detailed balance in the stochastic interpretation). The processes we are interested in are considerably more general from that perspective, and require the more general methods we develop here. These methods built on formulas derived by asymptotic analysis of backward Kolmogorov equations (BKEs), and are complementary to those recently proposed in [45]using a path integral approach.
Our main results can be summarized as follows: (1) We provide formal but short and simple proofs of propositions establishing sharp estimates for expectations, probability densities, probabilities, exit probabilities, and mean first passage times. While most of these results can be found in some form in the literature, they are scattered in many different papers, and we believe it is useful to collect and summarize them in one place. The methods used in these proofs can also be extended to more general situations not covered here. (2) We phrase each of our theoretical statements in a way geared towards numerical implementation, unlike what is usually found in the literature on the topic. (3) We use the geometric approach from gMAM to provide statements that are valid at infinite time (i.e. on the invariant measure of the process), by explicitly computing this infinite time limit via mapping of ∈ (-∞, 0] onto the normalized arc-length ∈( 0 ,1 ]along the instantons involved. (4) We formally generalize our results to the infinite-dimensional setup, with applications to stochastic partial equations of reaction-advection-diffusion type. (5) We also generalize our result to examples driven by non-Gaussian noise, specifically Markov jump processes (MJPs) in specific limits. (6) We illustrate the applicability of our method through tests cases in finite and infinite dimension, in which we discuss how to perform the numerical calculation involved.
In terms of limitations, our work focuses on situations where the stochastic system at hand has a single attracting point in the limit of vanishing noise. This setup is of interest in several situations, but it excludes the important problem of analyzing rare transitions between metastable states. Some of the tools we develop here, in particular the Riccati equations whose solutions enter the expressions for the prefactors, may be useful to analyze these transitions, but we will not consider them here.
As stated before, we are interested in obtaining sharp asymptotic estimates for expectations, probabilities, exit probabilities, and mean first passage times. We will do so under the generic assumption that the LDT optimization problem associated with each of these questions, that is the minimization of the action in (1.2) subject to appropriate constraints and/or boundary conditions, is strictly convex. This simplifying assumption guarantees that the solution of the equations presented below exists and is unique, and therefore allows us to avoid dealing with local minimizers, flat minima, etc. that may require to generalize/amend some of the statements below. While this is not necessary difficult to do, at least formally, it leads to a zoology of subcases that we want to avoid listing.
To make (1.1) have unique solutions, we will also make:
Assumption 1.1. The vector field is 2 (ℝ ) and such that:
and the matrix is such that:
This assumption guarantees [42] that the solution to the SDE (1.1) exists for all times and is ergodic with respect to a unique invariant measure with a probability density function ∶ℝ → (0, ∞). This density is the unique solution to
where here and below the colon denotes the trace, that is ∶ ∇∇ = tr( [∇ ⊗ ∇]). We will also make a stronger assumption:
Assumption 1.2. The ODE ̇ = ( )has a single fixed point located at * (i.e. * is the only solution to ( ) = 0), which is linearly stable locally (i.e. the real part of all eigenvalues of the matrix ∇ ( * ) are strictly negative), and globally attracting (i.e. any solution to ̇ = ( )approaches * asymptotically).
This assumption implies that becomes atomic on = * as →0 , a property we will need when looking at expectations or probability on the invariant measure.
The remainder of this paper is organized as follows: In Section 2 we will first consider the problem of calculating sharp estimates of expectations, both at finite time (Sections. 2.1 and 2.2) and on the invariant measure of (1.1) (Section 2.3). In Section 3,w ew i l ls h o w how to calculate sharp estimates of probabilities densities at finite (Section 3.1) and infinite times (Section 3.2). In Section 4 we then build on these results to calculate probabilities at finite times (Section 4.1) and on the invariant measure of the process (Section 3.2). Finally, in Section 5 we consider the problem of mean exit time calculation. These results are illustrated on finite dimensional examples throughout the paper, and on infinite dimensional examples in Section 6, where we consider linear and nonlinear reaction-advection-diffusion equations with random forcing. The results in this paper are mostly for diffusions, but they can be generalized to other set-ups, in particular MJPs: this is discussed in Section 7. We end the paper with some conclusions in Section 8 and defer some technical results to Appendices.
Given the observable ∶ℝ →ℝwith ∈ 2 (ℝ ), consider
where the expectation is taken over samples of the SDE (1.1) conditioned on 0 = and evaluated at the final time = <∞ . Expectations of the form (2.1) are an important first step to compute all other probabilistic variables, as they correspond to the moment generating function of the probability distribution of the random variable ( ). Thus, they represent a problem dual to the one of finding the associated probability densities. This will become clear in Section 3, where we use the results presented here to derive similar sharp limits for probability densities.
We have the following proposition:
Proposition 2.1. Let ( ( ), ( )) solve the instanton equations
and ( ) be the solution to the Riccati equation
integrated backwards in time from = to =0along the instanton ( ). Then the expectation (2.1) satisfies
where
)
)
.6) Remark 2.1. The function ( , ) is typically referred to as the prefactor. LDT gives a rougher estimate lim →0 log ( , ) log ̄ ( , ) =1, (2.7) which would be unaffected if we were to neglect the prefactor ( , ) in ̄ ( , ). Of course, this prefactor is key to get the more refined estimate in (2.4). Remark 2.2. If the SDE (1.1) is modified into = + ̃ + √ . (2.8) with ̃ ∶ℝ →ℝ is 1 (ℝ ) with bounded derivatives, Proposition 2.1 can be generalized by replacing (2.6)with
while leaving unchanged all the other equations in the proposition. The statements in the propositions below can similarly be straightforwardly amended to apply to (2.8), but for the sake of brevity we will stick to (1.1).
It is well-known that satisfies the BKE
where is the generator of the process (1.1):
Look for a solution of (2.12) of the form
where ( , ) satisfies the Hamilton-Jacobi equation
If ( ( ), ( )) solves the instanton equations in (2.2), we have ( )=∇ ( , ( )) and a direct calculation shows that
implying
As a result, to show that (2.4) holds, it remains to establish that the factor ( , ) has a limit as →0with lim →0 (0, ) = ( , ). To this end, notice that ( , ) satisfies
Formally taking the limit as →0on this equation, we deduce that lim →0 ( , ) = ( , ), where ( , ) solves
(2.20)
Setting ( ) = ( , ( )) on the instanton path, so that (0) = (0, (0)) = (0, ), and using the instanton equation ̇ = ( )+ ,wefindasevolutionequationfor
where the matrix ( ) is defined as the Hessian of ( , ) evaluated along the characteristics (i.e. the instanton path ( )):
( )=∇∇ ( , ( )). In order to solve the equation (2.21)f o r ( ), we need an equation for ( ). Differentiating the Hamilton-Jacobi equation (2.15)twice with respect to and evaluating the result at = ( ), it is easy to show that solves the Riccati equation in (2.3). Therefore, by integrating (2.21), we deduce
)
which terminates the proof. □
An alternative justification of the Riccati equation (2.3) can be given by introducing a stochastic process that samples the Gaussian fluctuations around the instanton. This approach opens up the possibility to alternatively compute the prefactor contribution as an expectations via sampling techniques. This result is well-known (see e.g. [24]) and can be phrased as:
Proposition 2.2. The prefactor (2.6) satisfies
where ( ( ), ( )) are defined as in Proposition 2.1, solves
and the expectation 0 is taken over realizations of (2.24) conditioned on 0 =0.
where ( ) is the instanton solution to (2.2). By invoking Girsanov's theorem, we can write ( ) as
where is the Radon-Nikodym density
.
(2.27) Since = ( ) + √ , after expanding both and in , it is easy to see that the leading order contribution is precisely given by exp
)
The next order vanishes due to the criticality of the minimizer. The first correction term in the exponential is therefore ( 0 ), that is it gives the prefactor ( , ),andreads
where = . Noticing that the term exp
is a itself Radon-Nikodym density for the change of measure from the random process defined in (2.24)t o = , we can therefore alternatively write the right hand-side of (2.28)a s in (2.23). □
Note that the formula (2.23) immediately tells us that the prefactor, defined as the limit as →0of the ratio between ( , ) and ̄ ( , ), is unity when both the drift and the observable are linear. Note also that computing the expectation by using the change of measure from the original process to one representing fluctuations around the instanton can be seen as an approximated way (up to terms of order ( )) of performing importance sampling via Monte-Carlo method: how to do this importance sampling exactly is harder in general, as discussed for examplein [51].
Finally, note that another, more direct, proof of Proposition 2.2 goesasfollows.Itiseasytosee that
where ( , ) solves
with
The solution of (2.31) can be expressed as
where ( ) and ( ) solve (2.3)a nd(2.21), respectively. Therefore (0, 0) = (0), consistent with the statement in Proposition 2.1. In Appendix B we list a few more expressions that relate the solution to Riccati equations like (2.3) to expectations over the solution of a linear SDE like (2.24). These expressions are useful as they give a possible route to solve the Riccati equation (2.3) via sampling, which is less accurate in general but simpler than solving the Riccati equation itself.
Given the observable ∶ℝ →ℝwith ∈ 2 (ℝ ), consider the expectation
where ( ) denotes the invariant density solution to (1.8). (2.34) can also be written as
) Assumption 1.1 guarantees that this limit exists for appropriate (e.g. if | | is bounded), and is independent on . The next proposition shows how to calculate a sharp estimate for (2.34) using the procedure of gMAM that allows us to compactify the physical time interval [0, ∞) onto [0,1]. Proposition 2.3. Let ( ̂ ( ), ̂ ( )) solve the geometric instanton equations
with =| ( ̂ )| ∕| ̂ | ,and ̂ ( ) be the solution to the Riccati equation
= ∇∇ ( ̂ (1)), (2.37) integrated backwards in time from =1to =0along the instanton ̂ ( ). Then the expectation (2.34) satisfies lim →0 ̄ =1, (2.38)
where
with
)
.
(2.40)
For more details about the geometric instanton equations we refer the reader to [29,31,34]as well as Appendix A. Note that the parametrization of ̂ can be chosen arbitrarily: In the calculations, it is convenient to use normalized arc-length, that is impose that | ̂ ( )| = , where the constant is the length of the instanton.
Proof. In order to speak meaningfully about the limit →∞ , we introduce a reparametrization of time, ( ) ∶ [0, 1] → ℝ + , to compactify the infinite "physical" time interval. In particular, we choose ( ) such that
and denote ̂ ( ) = ( • )( ). We then have
where we introduced ( ) = ∕ . As a consequence, ̂ ( ) fulfills the geometric instanton equations (2.36). These equations can be used to show that, in the limit as →∞ , the initial condition becomes irrelevant, and we can consider
, so that in the limit →∞ , we must have ̇ = ( ) f o ra ni n f i n i t e amount of time if the action is to remain finite. As a consequence, for →∞ , the global minimizer will decay towards the unique fixed point * , as all initial points are attracted towards it.
In order to find the global minimum, we can consider the two separate problems of first approaching the fixed point, and subsequently leaving again. For this purpose, consider the trajectory ( ), and corresponding reparametrized trajectory ̂ ( )
where
(2.44)
It follows that ̂ ( ) corresponds to the trajectory that deterministically decays into the fixed point * starting from , and then corresponds to the minimizer ̂ , solution to equations (2.36)fr om then on.
Since * is the unique fixed point and all ∈ℝ are attracted to it, such an ̂ 1 exists and is unique for all . On the other hand, since ( ̂ 1 , ̂ 1 ) fulfill instanton equations on ∈ [-1, 0], and ̂ 1 = (̂ 1 ), we have that ̂ 1 =0 and the corresponding action vanishes on the interval ∈ [-1, 0]. Therefore, the action associated with ̂ ( ) is equal to the action associated with ̂ ( ). The problem of finding a global minimizer starting at reduces to the problem of solving (2.36).
Another consequence of taking the limit →∞is that the Hamiltonian is zero along the instanton, that is ( ̂ ( ), ̂ ( )
(2.45)
Since ̇ | ( ) = ( ) ̂ ( ), this allows us to deduce the following expression for :
This is the expression stated in the proposition. It remains to be shown that the reparametrization of the Riccati equation, (2.37), is well posed. To this end, notice that as →0 ,wehave ̂ ( ) → * and ̂ ( ) → 0. Therefore,
which leads to the conclusion that, since ∇ ( * ) is negative definite by definition, we have ̂ ( ) → 0 as →0 . □
An easy example that can be computed explicitly is the case of diffusion in a potential landscape,
where ∈ℝ ,a n d ∶ℝ →ℝ is a potential with a unique minimum * , such that ∫ ℝ --1 ( ) < ∞ for all >0 -we can set ( * )=0without loss of generality. Given an observable ∶ℝ →ℝsuch that ( ) grows sufficiently slower than ( ) as | | →∞,wewant to obtain a sharp estimate of the expectation
where ( ) is the density of the invariant measure of (2.48), given by
Combining (2.49)and(2.50) we see that is given by
and both integrals can be estimated by Laplace's method. The result is that lim →0 ∕ ̄ =1with ̄ given by
det ( )-∇∇ ( )
where ( ) = ∇∇ ( ) is the Hessian of the potential and the point ∈ℝ is the solution of the optimization problem =argmin ∈ℝ ( ( ) -( )).
(2.53)
We will now show that Proposition 2.3 yields the same result. First, we know explicitly that the instanton fulfills
where primes again denote derivatives with respect to arc-length. Further, the endpoints of the instanton are * = ̂ (0) and = ̂ (1): the first is by definition, and the second since the final condition for the equation for ̂ ( ) gives
where the second equality follows from the second equation in (2.54). Since (2.55) is also the Euler-Lagrange equation for the minimization problem in (2.53), we deduce = ̂ (1). Therefore, using
we deduce that the exponential term in the expectation is given by exp
))
=exp -1 ( )-( ) .
(2.57)
Second, we have
This means that the Riccati equation (2.37) is here given by
= ∇∇ ( ).
(2.59)
Let us look for a solution of the form
for a matrix ∈ℝ × ,w i t h (1) = Id,a n d (1) (1)=∇∇ ( ) to fulfill the boundary conditions for ̂ .E q u a t i o n( 2.60)c a nb ew r i t t e na s = -1 ̂ which, from Liouville's formula,
(2.61)
From Equation (2.60) it also follows that
which we can use in conjunction with the Riccati equation to get
or equivalently
(2.64)
Using again Liouville's formula (this time for Ψ= -) yields the relation
into which we can insert the boundary conditions (1) (1) = ∇∇ ( ), (1) = Id and (0) = 0, as well as det (0) from Equation (2.61)toobtain
This gives eventually the prefactor contribution,
det ( )-∇∇ ( )
.
(2.67)
Therefore we recover ̄ given in (2.52), as needed.
10970312, 2024, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/cpa.22177 by College Of Staten Island, Wiley Online Library on [26/02/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Note that the above results can be generalized to a diffusion in a potential landscape with mobility matrix, that is systems of the form
where ∈ℝ × is the symmetric, positive definite mobility matrix, and 1∕2 ∈ℝ × is a symmetric matrix with 1∕2 1∕2 = . The procedure above can be repeated by replacing ̂ with -1∕2 ̂ -1∕2 and with -1∕2 -1∕2 . The result (2.52) is unchanged.
Here we estimate the probability density function of in the limit of small . Denote this density by ( , ), so that, given any suitable ∶ℝ →ℝ,
The next proposition shows how to get a sharp estimate of ( , ) at any by purely local considerations:
Proposition 3.1. Let ( , ( ), , ( )) solve the instanton equations
and let , ( ) be the solution to the forward Riccati equation
In addition let ( , ) be defined as
Then the probability density function ( , ) of satisfies
where ̄ ( , ) is given by
)
3)for , is structurally equivalent to the equation (2.3)for , = -1 , except for the boundary conditions. The boundary conditions necessitate solving , forward in time. Notably, this direction of integration is also the one in which the equation for , is numerically stable, as can be intuited for example by considering the Ornstein-Uhlenbeck process ( ) = -,for >0(see Section 3.1.1 below). This feature will become even more apparent in the context of stochastic partial differential equations, see the discussion after Proposition 6.1 in Section 6.
Remark 3.2. In Appendix B we discuss how to solve (3.3) via sampling of the solution of a nonlinear (in the sense of McKean) SDE.
Intuitively, the local estimation of of ( , ) implied by Proposition 3.1 is possible because in the limit →0 , this probability density is dominated by the instanton, and the prefactor contributions can again be estimated from the Gaussian fluctuations around this instanton. More concretely, we will see in the proof of Proposition 3.1 below (see (3.23)) that , ( ) is the inverse of the Hessian of the action ( , ), that is
This implies that we can also write (3.7)as
) .
(3.9)
This form shows that ̄ ( , ) is normalized in the limit as →0 , a result we state as Lemma 3.2. The density ̄ ( , ) defined in (3.7) is asymptotically normalized, that is
Proof. Starting from expression (3.9)f o r ̄ ( , ), which as we show in the proof of Proposition 3.1 is equivalent to (3.7), let us evaluate the integral ∫ ℝ ̄ ( , ) by Laplace's method. To this end, we must first identify the point where ( , ) is minimal. It is easy to see that this point is = ( ), where ( ) is the solution to the ODE ̇ = ( ), (0) = . Indeed, ( ) is also the instanton obtained when , ( ) ( )=0, and from (3.5)i tg i v e s ( , ( )) = 0, which implies that the minimum at is necessarily a global minimum. Similarly, , ( ) ( )=0implies that , ( ) ( ) = 0. Therefore
10970312, 2024, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/cpa.22177 by College Of Staten Island, Wiley Online Library on [26/02/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
where to get the second equality we used =( -( ))∕ √ as new integration variable and only keep the zeroth order term in that contributes to the limit. □
In the proof of this lemma, only the behavior of ̄ ( , ) near its maximum at ( ) matters, but let us stress that the expression for ̄ ( , ) given in (3.7) offers a much finer approximation of this density valid arbitrary far away from ( ).T h i sw i l l allow us to use ̄ ( , ) to estimate expectations dominated by tail events, as shown in Section 3.3.
Proof of Proposition 3.1. Consider the expectation
We can alternatively express this expectation as
We will show that ̄ ( , ) allows us to estimate this expectation for all ∈ℝ . Assuming that ( , ) is of the form
we can evaluate the integral in (3.13) by Laplace's method to obtain
where
From Proposition 2.1 we also know that
)
where ( ( ), ( )) solve the instanton equations (2.2)and ( ) solves the Riccati equation (2.3); note that, since ( ) = ⟨ , ⟩ here, the boundary conditions at = reduce to ( ) = and ( ) = 0. Comparison between (3.15)a n d( 3.17) implies that ( ) = , that is the instanton equations for ( ( ), ( )) ≡ ( , ( ), , ( )) reduce to the form in (3.2). In addition, ( , ) is givenby(3.5)and
) .
(3.18)
Using Proposition C.1 with ( ) = , ( ), we can identify ( )= ( ),and ( )= , ( ) to deduce that
since ( ) = 0 and , (0) = 0 and so det(Id -( ) , ( )) = det(Id -(0) , (0)) = 1. Therefore, to finish the proof it remains to evaluate det[∇ ∇ ( , )].T ot h i se n d ,n o t i c efi r s t that, since ( , ) isgivenby(3.5), it satisfies
Indeed, the solution to this equation can be expressed as in (3.5), and the boundary condition follows from the fact that, for small , to leading order the solution to the instanton equations (3.2) reads
This implies that
consistent with the boundary condition in (3.20). Therefore, if we introduce
so that ∇ ∇ ( , ) = -1 , ( ), then an equation for , ( ) can be derived by (i) differentiating (3.20) twice with respect to and evaluating the result at = , ( ) to obtain an equation for ∇ ∇ ( , , ( )), and (ii) using this equation to derive one for , ( ). The result is the Riccati equation (3.3), in which the boundary condition follows from lim →0 ∇ ∇ ( , ) =
1 2 This one-dimensional and linear case is the easiest possible non-trivial scenario, and in particular we know explicitly that
) .
(3.25)
We want to compare this analytical result to the approximation ̄ ( , ) giveninEquation(3.7).
In fact it turns out that ̄ ( , ) of proposition 3.1 is exact in this case, since there is no higher order contribution in in the prefactor. To show this, we have the instanton equations
which are solved by
(3.27)
The exponential estimate therefore yields exp
For the prefactor, we have the Riccati equation
which is solved by
and since ∇∇ ( ) = 0 in this example, we obtain
which is precisely the analytical result.
We can generalize Proposition 3.1 to calculate the probability density of the invariant measure, using again the procedure of gMAM to compactify physical time:
Proposition 3.3. Let ( ̂ ( ), ̂ ( )) solve the geometric instanton equations and let ̂ ( ) be the solution to the forward Riccati equation
wherewedenote
Finally, let ̂ ( ) be given by
Then the probability density function ( ) of satisfies
in which ̄ ( ) is given by either one of the two equivalent expressions
The function ̂ ( ) defined in (3.37) is important: it is the quasipotential of with respect to * . Interestingly ̂ ( )=[∇∇ ̂ ( ̂ ( ))] -1 and so (3.7) can also be written as
where we defined
The second expression in (3.40) is consistent with the one derived in [9]. Note that we can find yet another illuminating form to write (3.40) by use of the transverse decomposition that decomposes the drift ( ) into a gradient of the quasipotential, and a transverse portion ( ), defined as
Such a decomposition is guaranteed to exist under the assumption that the quasipotential is continuously differentiable [24, chap. 3.4]. In terms of ( ) (3.41) reduces to ̂ ( , ) = ∇ ⋅ ( ̂ ( )). Note that for a reversible diffusion, that is a gradient flow, we have ( ) = 2 ̂ ( ), and so ( ) = 0 and ̂ ( , ) vanishes. More generally, in flows that have ∇⋅ ( )=0, the invariant measure remains the Gibbs-measure,
leading the authors of [9] to call the quantity exp(-
the non-Gibbsianness of the flow. It is important to remark, though, that this simplification is only available in the infinite time case on the invariant measure, as for finite times the corresponding quantity ∇ ∇ ( , ) becomes singular at =0 . Further, since ∇⋅ is generally not available explicitly, the first form of (3.40) is the one most easily used for numerical computation.
We can use (3.40) to show that lim →0 ∫ ℝ ̄ ( ) = 1, that is the equivalent of Lemma 3.2 holds in the infinite time limit using the geometric expressions above. Basically, this is because the dominating point in this integral is = * , for which ̂ * ( ) = 0 and ̂ * ( ) = * =∇∇ ̂ -1 ( * ). This also shows that, ( √ ) away from * , ( ) can be approximated by the Gaussian density with mean * and covariance * given by
and this density is normalized. However, for locations that are (1) away from * , this Gaussian approximation is no longer valid, and the full expression (3.39)mustbeusedasestimateofthe probability density on the invariant measure, as shown in Section 3.
that is
) )
This shows that the two expressions for ̄ ( ) in (3.39) are identical. □ Remark 3.3. Note that this is consistent with the intuitive interpretation that ̂ quantifies the covariance of fluctuations around the instanton. For →∞the fluctuations will "thermalize" around the fixed point, which corresponds to considering the linearized (Ornstein-Uhlenbeck) dynamics around * , the covariance of which solves the Lyapunov equation (3.34).
Remark 3.4. We seemingly encounter a practical problem with the forward Riccati equation (3.35) at =0 , as also discussed in [10]: Since the instanton remains at the fixed-point * foraninfinite amount of time, we have (0)=0,a sw e l la s ̂ (0) = 0, and thus the forward Riccati equation (3.35) reads 0=0at =0 -a similar issue arises with the equation for ̂ in (3.33)a n di s discussed in Appendix A. This problem is only apparent however and the limit of ̂ ( ) as →0 can be straightforwardly obtained by writing (3.35)as
and sending →0at both sides using l'Hôpital rule to compute the limit of the right hand side.
Accounting for the fact that ̂ (0) = * and ̂ (0) = 0 since ̂ (0) = 0, the result is the following equation for ̂ (0): Since (0) > 0 and -∇ ( * ) is positive-definite as * is a stable fixed point of ̇ = ( ) by assumption, the matrix ℭ is positive-definite and invertible, meaning that (3.50)h a sau n i q u e solution ̂ (0).Once ̂ (0) has been calculated, one can initialize the integration of ̂ ( ) with arc-length stepsize ∆ > 0 using for example
For the examples presented in Sections 3.2.2 or 4.2.1, it turns out that ∇∇ ( * )=0and thus (0) = 0. Therefore, approximating ̂ (∆ ) by * is correct to (∆ 2 ), which is precise enough for the numerical scheme we employ, and is therefore what we used in practice. However, we note that it is straighforward to go to higher order if necessary: for example, the derivation to obtain ̂ (0) when ̂ (0) = 0 is given in appendix A.2.
Similar to the expectations for gradient diffusion, discussed in Section 2.3.1, we can also derive a formula for the small approximation of the invariant density of the gradient diffusion process (2.48). The result is known to be
where the only approximation made is on the prefactor that can be evaluated by Laplace's method (using ( * )=0): =(2 ) ∕2 ((det ( * )) -1∕2 + ( )). Let us show that the small approximation of Proposition 3.3 recovers this result, including the normalization factor.
First, the instanton contribution yields
which is solved by
so that the exponential large deviation contribution is given by
which recovers the exponential part of ̄ ( ). If we additionally notice that for a gradient system we know explicitly the transverse decomposition to be trivial, ( ) = 0 and thus ̂ ( , )=0, the second form of (3.39)inproposition3.3 immediately yields the limiting density of (3.53).
It is instructive to recover the same result by explicitly solving the Riccati equation. Notice that ̂ ( )=∇∇⟨ ( ̂ ( ), ̂ ( )⟩ =-( ) ( ̂ ( )), so that forward Riccati equation (3.35) reduces to Combining these results we indeed recover exactly the limiting density given in (3.53).
Example: Invariant measure of a nonlinear, irreversible process in ℝ 2
For all above examples, analytical results were available from case-specific calculations, since the systems are very simple. The easiest example of a system where no analytical results can be easily derived is a system which is nonlinear, and further the drift is not given by a gradient of a potential. Since the latter is always the case in 1D, we need to consider a state space of at least dimension two. In this case, in order to compute expectations, probability densities, or probabilities with our approach, we need to solve the corresponding equations, both instanton equation and Riccati equation, by numerical means. As a concrete example, consider the non-gradient drift given by
For positive , ∈ ℝ, this drift corresponds to a nonlinear attractive force towards the fixedpoint * = (0, 0).T h ep a r a m e t e r ≠ 0 adds a swirl on the drift, so that the total dynamics are 1 2 -1 ̂ ( ) ,at = (1, 1) as a function of the nonlinearity parameter , obtained by numerically integrating the Riccati equation along the instanton using proposition 3.3. Clearly, the strength of the nonlinearity influences the prefactor. The other parameters are again =0 . 5 , =1 .
no longer gradient. The behavior of this drift can be seen by the streamlines in Figure 1 (left). Note in particular that the system is not rotationally symmetric in the presence of cubic terms (i.e. ≠ 0).
In order to estimate numerically the density ̄ ( ) defined in (3.39), we need to:
(1) compute the instanton ̂ ( ) connecting (0,0) to , as well as the parametrization ( ),using gMAM; and (2) solve the forward Riccati equation (3.35).
The results of this procedure are shown in Figure 1. For a given arbitrary endpoint (1,1), we can obtain the invariant density by sampling the process for long times, and approximating the density by a normalized histogram. The result is depicted as shaded contours on the left. Note that for this procedure, not only do we need many samples hitting close to the point (1,1), but furthermore statistics everywhere else in state space, because these determine the normalization constant. The instanton connecting (0,0) to (1,1) is depicted here as a solid line, while the drift is given by the flowlines.
We now want to establish how the prefactor contribution, ̄ ( ) 1 2 -1 ̂ ( ) , changes when changing the parameter , which determines the strength of the nonlinear term. In particular, for the linear case =0 , we know that ̄ ( ) 1 2 -1 ̂ ( ) =1∀ ∈ℝ 2 , since the invariant measure is Gaussian, and its covariance can be computed by solving a Lyapunov equation. For the nonlinear case, >0 , this is no longer possible, and we need to resort to numerical results instead. In the right panel of Figure 1 we show the computation of the prefactor through proposition 3.3, computing the instanton and solving the Riccati equation. For =0 , this reproduces the known linear result, but we obtain values for finite values of as well. These results are in agreement with results from sampling the whole process, and looking at the normalized histogram at the point (1,1) for different ∈[ 0 ,2 ] . The numerical parameters here are = 0.5, = 1, ∈ [0, 2], = 0.25, and the histogram binning is done with square bins of side-lengths ∆ = 0.02 and with 10 7 samples per value of .
We can use Proposition 3.1 to give an expression alternative to that in Proposition 2.1 for finite-time expectations:
Proposition 3.4. Let ( ( ), ( )) solve the instanton equations in (2.2), and ( ) be the solution to the forward Riccati equation
wherewedenote ( ) = ∇∇⟨ ( ( )), ( )⟩. Then the factor ( , ) defined in (2.6)ofProposition 2.1 can also be expressed as
)
.64) Proof. We can arrive at expression (3.64)f o r ( , ) in two ways. We can use Proposition C.1 with ( ) solution to (2.3)a n d ( ) solution to (3.63) to deduce from (C.2) that
The left hand side is (3.64) and the right hand side is (2.6), showing that these two equations contain identical expressions. Alternatively, we can evaluate
by Laplace's method, using the expression in (3.7)for ̄ ( , ). This calculation show that this integral is asymptotically equivalent to ̄ ( , ) in (2.5)with ( , ) given by (3.64) instead of (3.64), thereby establishing again that these two equations give the same ( , ) since this factor is independent on . □
Similarly, we can use Proposition 3.3 to give an expression alternative to that in Proposition 2.3 for expectations on the invariant measure: Proposition 3. . Let ( ̂ ( ), ̂ ( )) solve the instanton equations in (2.36), and ̂ ( ) be the solution to the forward Riccati equation
where * solves (3.34) and we denote ̂ ( ) = ∇∇⟨ ( ̂ ( )), ̂ ( )⟩. Then the factor ̂ defined in (2.40) of Proposition 2.3, can also be expressed as
) .
(3.68)
The proof of this proposition is similar to that of Proposition 3.4.
Having access to pointwise estimates of the probability density allows us to estimate probabilities by integration. For example, let ∶ℝ →ℝbe some observable, and assume we want to compute the probability that ( ) exceeds some value ∈ℝ . This probability can be expressed as
and calculating for various values of gives the complementary cumulative distribution function (aka tail distribution) of the random variable ( ). For concreteness, we will focus on the calculation of this tail probability for one value of which, without loss of generality, we can set to zero. In this case, (4.1) can also be interpreted as the probability that the stochastic process (1.1) hits the set
at time = . This is a problem interesting in its own right, and we will denote this probability as
To make the problem interesting, we will work under: This assumption implies that the event is noise-driven, and as a result ( , ) → 0 as →0 , which is the nontrivial case. We also need is unique, and so is the instanton leading to it. In addition is not a critical point of .
Intuitively we demand that the set does not admit multiple distinct points that can be reached by competing instantons of identical action. Then, we have: and let , ( ) solve the forward Riccati equation in (3.3). Let also ( , ) and ( , ) be defined as
where , ( ) = ∇∇⟨ ( , ( )), , ( )⟩. in which
) . (4.8)
In addition , the unit normal ̂ can be expressed as ̂ = , ( )∕| , ( )|.
The situation of Proposition 4.3 is depicted in Figure 2.
Remark 4.1. As we will see in the proof, the factor involving ( , ) describes the effect of the geometry (specifically, the curvature) of the set around the maximum likelihood hitting point : If the set is planar, then ∇∇ ( ) = 0 and the factor evaluates to unity. Equivalently, this factor disappears for linear observables . As will also be clear from the proof, in the definition of ( , ) in (4.6) we can replace |∇ ( )| -1 ∇∇ ( ) by ∇ ̂ ( ), that is in (4.8) we can replace ( , ) with This is because ∇ ̂ =|∇ | -1 ∇∇ + |∇ )| -1 ̂ (∇∇ ̂ ) and the extra factor |∇ | -1 ̂ (∇∇ ̂ ) does not contribute to the ratio ⟨ ̂ , ( , ) ̂ ⟩| det ( , )| -1 in (4.8)
. This shows that this expression is intrinsic to the set , that is it does not depend on the way we parametrize its boundary by the zero level-set of , as it should be. For computational purposes, using the parametrized versionin(4.6) is more convenient, however.
Proof. Let us evaluate the integral in (4.3) by Laplace's method using the ansatz (3.14)for ( , ).
To this end notice that the point specified in (4.5) is also given by =a r g m i n ∈ ( , ) ∈ . (4.10)
Notice also that , ( )=∇ ( , ), and that ̂ = , ( )∕| , ( )| is the inward pointing unit normal to at . Therefore, to perform the integral
we split the integration variable into components parallel to ̂ and perpendicular to ̂ ,as
where ∥ ∈ℝand ⟂ ∈Ω ⟂ ̂ with
For any ∈ we have ( ) ≥ 0. Further using ( ) = 0,wecanexpand around to obtain
Since ⟨ ⟂ ,∇ ( )⟩ =0by definition, we have for points in around that
Effectively, to the relevant order in , can be approximated by a paraboloid
where ∇∇ ( )∕|∇ ( )| is the (normalized) curvature of the set at .
We can now evaluate the integral in (4.3) by Laplace's method,
Here
)
and
where det ⟂ is defined as follows: Given an invertible, positive definite, symmetric ∈ℝ × and a unit vector ̂ we define det ⟂ via
In other words, det ⟂ is the determinant evaluated only in the space Ω ⟂ ̂ perpendicular to a vector ̂ . As shown in Appendix E we have
and thus, inserting in (4.17), this equation gives (4.7). Finally, note that by definition of det ⟂ we can replace |∇ ( )| -1 ∇∇ ( ) by ∇ ̂ ( ) in (4. 19), which justifies the alternative expression in (4.9) for ( , ). □
The equivalent of (4.3) at infinite time is
where ( ) is the invariant density defined by (1.8) and the set is defined as before. Then we have: and let ̂ ( ) solve the forward Riccati equation in (3.35). Let also ̂ and ̂ be given by
) where In addition , the unit normal ̂ can be expressed as
The proof is similar to that of Proposition 4.3.
To show the applicability of these propositions to an actual system, we re-use the nonlinear, irreversible process in ℝ 2 , as defined in Equation (3.62), as a simple example for which the solution is not readily accessible by analytical considerations. Instead of computing the invariant density, as in Section 3.2.2 we use proposition 4.4 to compute the probability of hitting the set on the invariant measure, where is the half-space defined by
with ̂ ≈ (0.6304, 0.7762) chosen specifically such that the point =( 1 ,1 )is the maximum likelihood point on . Because of this, the solutions to the instanton equations and the Riccati equation do not need to be re-computed, and we can insert their results into equation (4.27) to obtain the instanton and prefactor estimate for the hitting probability on the invariant measure.
The results of this experiment are shown in Figure 3: While the dynamics and instanton in Figure 3 (left) look identical to the original problem in Section 3.2.2, we are now trying to estimate the probability of hitting the set denoted by the light green shading. The result shown in the right panel of Figure 3 confirms that the asymptotic prediction of proposition 4.4 agrees with the Monte-Carlo simulations for different values of . The parameters are = 0.5, = 1, = 0.5,and we are taking samples =10 6 samples each value of .
⊂ℝ be the complement of the set defined in (4.2) that is
Under Assumption 4.1, we know from the argument given after Proposition 3.3 that
We wish to estimate the exit probability
in the limit as →0 , since this surface integral enters the limiting expression for the mean first passage time (MFPT) of the process to set .Wehave Proposition .1. Let ( ̂ ( ), ̂ ( )) solve to the instanton equations in (3.33) with specified as
and let ̂ and ̂ be given by (4.25). Denote by ̂ =∇ ( )∕|∇ ( )| the inward pointing unit normal vector on at . Then the exit probability satisfies
where
Proof. The proof has the same ingredients as the proof of proposition 4.4, except instead of integrating over the paraboloid approximation ̃ , we integrate over its boundary ̃ . No normal integral is performed and therefore the scaling in and | ̂ (1)| differs. □
With this knowledge, we can now compute exit times. Consider the entrance into the domain for a particle starting at ∈ . The mean exit time ∶ →ℝ + is given by the solution of the problem
where is the generator of our process defined in (2.13). In the following, we quickly review the standard approach to use a boundary layer expansion in order to approximate the exit time .
Multiplying the first equation by the invariant measure and integrating we find-after applying Green's theorem and making use of the fact that * =0 [25]:
Here, ̂ ( ) = ∇ ( )∕|∇ ( )| denotes the inward pointing normal unit vector of the surface at ∈ . An approximation of ⟨ ̂ ( ), ∇ ( )⟩ for ∈ can be obtained by boundary layer analysis. For this purpose, we first expand ( ) = ∕ ( ) (5.9) suchthatwefindfrom(5.7) the corresponding equation for = -∕ ≈0.
(5.10)
Since we assumed that ⟨ ̂ ( ), ( )⟩ <0for all ∈ , the appropriate scaling of the boundary layer is = -̂ . In the scaled variables, the equation for becomes
with the solution = ̃ 1-⟨ ̂ ( ), ( )⟩ (5.12)
This means that we obtain for the exit time the expression = ̃ ∕ 1-⟨ ̂ ( ), ( )⟩ (5.13) and therefore
which we can use in the solvability condition (5.8). From there, for away from the boundary layer, we obtain
The following proposition shows how to estimate ( ) in the limit as →0 :
Proposition .2. Let ( ̂ ( ), ̂ ( )) solve the instanton equations in (3.33) with as specified in (5.4), and let ̂ and ̂ be given by (4.25).Denote by ̂ =∇ ( )∕|∇ ( )| the inward pointing unit normal vector on at .Thenforany ∈ , the mean first passage time ( ) satisfies lim →0 ( ) ̄ =1 (5.16)
where
) .
(5.17)
Proof. From (5.15) it follows with the use of proposition 5.1 that
where we additionally used that lim →0 ∫ ( ) = 1. □ Remark 5.1. Another interesting case is when we demand ⟨ ̂ ( ), ( )⟩ =0everywhere on , such that corresponds exactly to the basin of attraction of the process. In this case, the situation is more complicated. Generally, one expects a different scaling of the form
The underlying assumptions need to make sure that the quasipotential is twice differentiable at the exit point, which is true for gradient systems, but not generally true for an arbitrary drift. We refer to [9,40] for details.
For the 1D Ornstein-Uhlenbeck process,
∈ℝ, the formula for the MFPT is known exactly [48]. Concretely, the expected time for the process (5.20) to leave the set = [-∞, ] is given by
) .
(5.21)
In the limit →0 , truncating the prefactor at ( 3∕2 ), this yields the limiting result
Since necessarily = , and there is no perpendicular direction in 1D, Proposition 5.2 tells us that ̄ = exp( -1 2 ) with given by
Using
-1∕2
we obtain
in agreement with the analytical result in (5.22).
Consider instead the situation of ∈ℝ 2 , but still
We are interested in the expected time to enter the complement of the translated circle [33]o f radius around ( -, 0), that is exit
We want to compare this result against the translated half-plane,
where 1 is the first component of =( 1 , 2 ). The most likely exit point, located at =( ,0 ) , is identical in both cases, as is the instanton, but we expect to find different prefactors due to the difference in curvature at between the two sets. In the left panel of Figure 4 we illustrate the problem setup. We will focus on the case > when the instanton lies on the -axis: at = the instanton becomes degenerate and we obtain caustics, that is due to rotational symmetry every point on the circle is equally likely to be the exit point.
When > , we have, basically identically to Section 5.0.1,
. for different , in comparison to the theoretical prediction (5.28). For →∞ , this ratio should converge to 1, but for small radii the curvature has a measurable effect. For example, at a radius = 1 4 , the measured exit time is roughly 25% smaller than for the planar exit.
The only difference between the two cases (5.25)a n d( 5.26) are the curvature contributions ̂ and ̂ | , respectively. For the half-space | , the curvature at is 0, and thus ̂ | =Id. For the circle, instead, we can choose
is the zero level-set of . Then,
and thus
Defining the scalar contribution of the curvature to the prefactor = (det ⟂ ̂ ) -1∕2 as | for the planar case and for the case of a circle with radius , respectively, we obtain
In order to measure this curvature prefactor experimentally, we performed the following numerical experiment: For different radii from = 0.25 to = 100, we performed =2⋅1 0 5 Monte-Carlo simulations each, and measured the mean time DNS to exit the set .W ep e rformed further 2⋅10 5 Monte-Carlo simulations in the planar case | to obtain DNS | . The ratio of the measured passage times of the planar case to the circular case yields a numerical estimate of the curvature component only:
which can be compared against the analytical prediction
This comparison is shown in Figure 4.For →∞ , indeed converges to 1, but for small radii the curvature has a measurable effect, and the measured discrepancy to the planar prediction agrees with the correction term given in (5.28). The other parameters are =1 , =0.1,and =5⋅1 0 -3 .
In the previous sections, we derived prefactor estimates and sharp limits in finite dimension. All of these estimates have a counterpart in infinite dimension, that is when applied to stochastic partial differential equations.
For concreteness we will focus on reaction-advection-diffusion equations of the type (for more general equations see Appendix F)
Here ∈[ 0 ,∞ ) , ∈Ω⊂ℝ ,w i t hΩ compact, and ∶[0,∞)×Ω→ℝ; the vector field ∶Ω→ℝ is in 2 (Ω); the diffusion tensor ∶Ω→ℝ ×ℝ is in 2 (Ω), symmetric, ( ) = ( ), and positive-definite for all ∈Ω ; ∶Ω→ℝis in 2 (Ω); the reaction term ∶ℝ→ℝ ,is in 2 (Ω), nonlinear in general; and the noise is white-in-time Gaussian with covariance
where ∶Ω×Ω→ℝgiven and in 2 (Ω × Ω).IfΩ is a rectangular domain in ℝ , we can impose periodic boundary condition on , assuming that all the other functions in (6.1) are also periodic; otherwise, denoting by ̂ ( ) the unit normal to Ω,weimposethat ( ) = 0 and ̂ ( ) ⋅ ( ) = 0 on Ω, ( , ) = 0 for all ∈ Ωand ∈Ωor ∈ Ωand ∈Ω , and also that ̂ ( ) ⋅ ( )∇ ( , ) = 0 for all ∈ Ω and ≥ 0. (
In this infinite-dimensional setting, the determinants must be replaced by functional determinants, the instanton equations become PDEs, and the Riccati equation must be replaced by a functional variant as well. While these changes require a whole new set of techniques to rigorously prove validity, our results and methods remain formally applicable in infinite dimension.
We will omit the proof of this proposition since it essentially amounts to translating the equations in Proposition 4.3 to the infinite dimensional setting, and using the fact that the term involving the tensor ( , ) disappears since the equivalent of ∇∇ is zero for a linear observable as in (6.4). Similar reformulations of our other propositions are straightforward as well, and for the sake of brevity we will therefore omit writing them down.
Note the instanton equation for ( , ) in (6.5) as well as the Riccati equation in (6.7)f o r ( , , )are well-posed as formulated, that is forward in time; similarly the instanton equation for ( , ) in (6.5) is well-posed as formulated, that is backward in time.
In Sections 6.1 and 6.2, we confirm numerically that Proposition 6.1, or its counterparts for other quantities, produces results that agree with those obtained via direct sampling of the SPDE in (6.1). As example problems, we consider two special cases of (6.1): First, in Section 6.1 we study a linear advection-reaction-diffusion equation with spatially non-homogeneous forcing, for which some analytical results can be derived. The probability that the concentration exceeds a threshold at a given location can then be estimated by a mixed analytical-numerical approach. Second, in Section 6.2, we study a reaction-advection-diffusion equation with cubic nonlinearity, where we need to resort to fully solving numerically all involved instanton and Riccati equations.
Here we consider the stochastic one-dimensional advection-diffusion-reaction equation,
where ∈[ 0 ,1 ]periodic. This a special case of (6.1)with = , ( ) = -,and ( ) = -. For the covariance of the white-in-time forcing ( , ) we take
where 1 ( ) is a mollifier of length <1concentrated around 1 ∈ [0, 1]. For the observable, we take
where 2 ≠ 1 and the expectation is taken on the invariant measure of the solution to (6.12).
Intuitively, the scenario we are investigating is therefore that of a pollutant, the density of which is described by ( , ), along a one-dimensional periodic channel. The pollutant is randomly emitted into the environment at a spatial location 1 , and gets advected and diffused conservatively, but decays over time with rate . We are interested in measuring extreme concentrations of the pollutant around the location 2 somewhere else in the channel.
6.1.1 Finite-dimensional analogous case Equation (6.12) is a linear SPDE, an infinite-dimensional generalization of the (non-normal,dimensional) Ornstein-Uhlenbeck process,
for Γ∈ℝ × , where Γ ≠ Γ ⊤ (non-symmetric), ΓΓ ⊤ ≠ Γ ⊤ Γ (non-normal), is an -dimensional Wiener process, not necessarily invertible, and = ⊤ . In analogy to the above scenario, we can ask for probabilities on the invariant measure of the form =ℙ( ∈ ), where ={ ∈ℝ |⟨ , ⟩ ≥ }. (6.16)
Intuitively, we want to estimate the probability that the process reaches the far-side of a plane with normal ,distance away from the origin.
If we define the symmetric, positive semi-definite matrix as the solution of the Lyapunov equation
the quasi-potential of (6.15)isgivenby
Since is not necessarily invertible, we interpret = -1 to be the solution of = if it exists, and otherwise the quasi-potential is set to infinity. The final point of the geometric instanton ̂ ( ) for observable value must be given by =a r g m i n ∈ ( ( )) = ⟨ , ⟩
Here, we must assume that is not in the kernel of , which is the same as saying that is in the support of the invariant measure of (6.15). The action is
which follows from evaluating ( ) at = given by (6.19).
To estimate the prefactor, we need to solve the -equation with appropriate boundary condition. Because the Equation (6.15)islinear , ̂ ( ) =
1 2 . As a result lim →0 ̄ =1 with ̄ = √ ⟨ , ⟩ 4 2 exp ( -2 ⟨ , ⟩ -1 ) (6.21)
where we used
Due to the linearity of the system, only the end location of the instanton at =1plays a role.
Coming back to the infinite-dimensional case, we can proceed similarly, noticing that (6.12)can be written as
where
is a linear differential operator, acting on functions in 2 . Notably, isnotnormal,andweneed to solve the Lyapunov equation
where ( , ) is symmetric and positive semi-definite in 2 and ⊤ =( ) ⊤ , that is it is the differential operator acting on the second variable of ( , ). Explicitly (6.26) reads
By extending the argument for the finite dimensional case to the functional setting, we also deduce that the endpoint at =1of the geometric instanton ̂ ( , ) for hitting the set is given by
in analogy to Equation (6.19) by replacing the inner product with the 2 inner product. The corresponding action is (compare (6.20))
Concerning the prefactor, we have (compare (6.23))
6.30) 10970312, 2024, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/cpa.22177 by College Of Staten Island, Wiley Online Library on [26/02/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License and thus (compare (6.22)) = 1 2 ( ∫ 1 0 ∫ 1 0 2 ( ) ( , ) 2 ( ) ) 1∕2 , (6.31) so that in total, lim →0 ℙ ( ∫ 1 0 ( ) 2 ( ) ≥ ) ( ) =1 (6.32)
with
)
as in Equation (6.21). Clearly, this probability is Gaussian, as expected when we consider a linear system with Gaussian input.
In order to obtain , we numerically solve the operator Lyapunov equation (6.26) for a finite difference representation of . As an example, we take the advection-diffusion-reaction equation (6.12), with a velocity field
that is with negative divergence at the center of the domain, =1 ∕ 2 . We expect (if we were not to condition on any outcome, and would force homogeneously in space) that in a typical configuration the concentration has higher variance at =1 ∕ 2and lower around =0 , where it is depleted by the velocity. We also choose 1 =1∕32and 2 =1∕4, that is the forcing is localized on the very left of the domain, where concentration is released randomly, and we are sensing concentration further down the channel, where it is transported by advection and diffusion. We also check the results above numerically by comparing the instanton prediction derived here with the result of a direct simulation of the advection reaction diffusion equation (6.12) with stochastic forcing localized according to (6.13), and counting the number of times the observable exceeds a given threshold.
In the left panel of Figure 5 we show the corresponding instanton: It has a non-trivial dependence on both the location of the forcing as well as the location of the observation. Additionally shown are the localized functions defining the forcing and the observation, as well as the velocity field. In the right panel of Figure 5 we show a comparison against direct simulation results of the SPDE. In particular we note that the mean of all observed events within resembles the instanton. Further, the highest variance in the direct simulation is clearly observed where noise is injected, as well as close to =1 ∕ 2 , at the sink of the velocity field.
Depending on and , the probability can be computed via (6.32). In Figure 6 we show a comparison of these probabilities against numerics, obtained by integrating the SPDE (6.12) for a long time, and observing how often each threshold was exceeded. We note that the prediction indeed captures not only the scaling, but also the correct prefactor.
Endpoint of the instanton configuration (rescaled, red), that is solution of the minimization problem (6.28) for the advection-reaction-diffusion equation (6.12) with velocity field ( ) (green), conditioning on observing a high concentration in the region specified by 2 (blue), where concentration is released randomly around location 1 (yellow). Here, =1 0 -foot_1 , = 1, ( ) = 0.4 sin(2 ), =5⋅1 0 -2 , = 0.0025, 1 =1∕32, FIGURE 6 Probabilities for the advection-diffusion-reaction problem with localized forcing obtained by direct numerical simulations (blue), from formula (6.32) (black). Also shown is the instanton prediction (black dotted), that is only the exponential contribution from the large deviation scaling. It is rescaled by a factor 10 -2 to obtain the best fit in the tail. The difference between the instanton prediction and ( ) is noticeable.
For the linear reaction-advection-diffusion equation of Section 6.1,wewereabletoharnesslinearity to avoid explicitly using Proposition 6.1 and numerically solving the instanton and Riccati equations. This is no longer possible if we choose the system to be nonlinear, for example by taking a cubic reaction term. In this case, we no longer have access analytically to the invariant measure, the instanton trajectory, or the solution to the prefactor terms. Concretely, we consider the nonlinear stochastic partial differential equation
where the spatial variable is periodic on the domain [-∕2, ∕2] and ∈[ -,0 ] . This equation is a special case of (6.1)with ( ) = , ( ) = 0,and ( ) = --3 : the parameters and control the linear and nonlinear part of the reaction term, respectively. As velocity field we pick
which always advects to the right, but with spatially varying speed. As a consequence, the whole equation is no longer translation invariant. We will assume that the noise is white in space and time, with covariance ( ( , ) ( , ))= ( -) ( -). (6.37) This is a rougher noise than the one in the SPDE (6.5), which is allowed because (6.35) is wellposed in 1D with this forcing [20]. We are interested in the probability that a sample on the invariant measure of (6.35) exceeds the threshold at the location =0 , that is
Note that this problem can either be viewed as a nonlinear version of the reaction-advectiondiffusion equation of Section 6.1, or as an infinite dimensional version of the ℝ 2 process given in Section 4.2.1.
To apply our method we need to first solve the geometric instanton equations s i n g the appropriate boundary conditions. Since we are focusing in this example on the limit → ∞, an efficient way of numerically solving these equations using an iterative scheme and arclength parametrization has been discussed in previous work [29]. Since for Equation (6.35) the Hamiltonian is
we immediately obtain the (geometric) instanton equations in this case as
In this case, the boundary conditions ̂ (0) = 0 and ̂ (1) = ( ) allow direct application of the previously developed geometric techniques [29] to solve the above system iteratively in order to find the instanton solution.
Once the instanton is found, we need to solve the corresponding (geometric variant of the) forward Riccati equation for ( , , )given by (6.7).
For the specific SPDE (6.35), the equation for ( , , ) = ( , , )is
( , , )̂ ( , ) ̂ ( , ) ( , , ) +2 ( -) (6.41) to be solved with the initial condition, (0) = * , solution to the Lyapunov equation -1 2 2 * + 2 * + * + 1 2 ( ) * + ( ) * = ( -).
One can numerically integrate the instanton equations (6.40) to obtain the most likely configuration that achieves a given amplitude at =0at final time =0 . Fourier transforms were used to calculate the spatial derivatives. When solving the associated Riccati equation for , however, we chose an equally distant point grid to discretize the time interval [-, 0] for stability reasons. Here, the time can be found from the geometric parametrization [34] and the instanton solution needs to be interpolated onto this new point grid. For the solution of the Riccati equation on the equidistant grid, exponential time differencing [36] was employed and the diffusive term ( 2 + 2 ) can be treated for numerical efficiency using the two-dimensional fast Fourier transform.
In the left panel of Figure 7 we compare this numerically computed instanton with the filtered instanton from direct simulations of the SPDE (6.35) using the method described in [28]. As can be seen, the mean realization with =0.7at =0in the SPDE is very similar to the instanton. In particular, both instanton and the typical sample from the SPDE show an asymmetry around =0coming from the spatially inhomogeneous velocity field (6.36). Also shown is one standard deviation of the fluctuations around the instanton as shaded red region. In the right panel of Figure 7 we show the whole evolution of the instanton in arc-length parameter that compactifies the infinite time interval ∈[ 0 ,∞ )into ∈[ 0 ,1 ] . Clearly visible is the (inhomogeneous) movement of the peak as it is advected with the positive velocity ( ) given in Equation (6.36).
In the left panel of Figure 8 shows the conjugate momentum ̂ ( , ) along the arc-length parameter , and in the right panel we show the quantity ( ̂ ( , )) ̂ ( , ) ̂ ( , , ): the exponential of its integral along and enters the expression in (6.11) for the prefactor component ( , 0 ).
Comparison of the predicted probabilities ( ) of exceeding the threshold . The dark blue dots shows the probability estimated via direct numerical simulations for the nonlinear case ( =2 ), while the light blue line is the theoretical prediction from instanton and Riccati equation. Clearly, the instantons, together with the corresponding prefactors, approximate these probabilities well. For comparison, we show the analytical prediction that can be obtained for the linear case ( =0 ), highlighting the fact that the nonlinearity indeed plays a role for the tail probabilities in particular. Similarly, we show the situation without advection ( ( ) = 0), demonstrating that the advection term modifies the probabilities.
The combination of both results inserted into Proposition 6.1 allows us to determine the prefactor. In Figure 9 we compare the result from direct simulations of the SPDE (6.35) to the predictions of the instanton equations together with the prefactor, showing clear agreement especially for large values of , as expected. A comparison is further made to the analytical result available for the linear case, =0 , which clearly shows that the nonlinear term affects the probability. Similarly we compare against the case without advection ( ( ) = 0), which is gradient, demonstrating that the advective term also has a considerable (opposite) effect on the probabilities. In this numerical example, we chose = 256 grid points for the discretization in space for a domain of size =16,and = 5000 grid points in time for the direct simulations for a temporal domain of size =2 5 . The instanton is computed with = 2000 discretization points in the arc-length parameter. Additional parameters are =1 , =0 . 6 ,and =2 . The noise level is set to =0 . 1 ,andwe collected samples =10 5 samples in the direct simulations.
It is possible to generalize the above results to continuous time MJPs that cannot be represented by an SDE with additive Gaussian noise like in (1.1). The intuition is similar: considering a WKB approximation to the BKE allows us to obtain a Hamilton-Jacobi equation that defines the behavior of the instanton. If we furthermore keep track of the prefactor in the WKB, we will obtain an additional equation for it on the next order, which will yield a generalization of the Riccati equations to obtain sharp prefactor estimates.
Concretely, consider a continuous-time MJP on the state space ℝ with the generator
which encodes a reaction network for ∈ℕspecies with ∈ℕreactions, each reaction leading to a change in species defined by the vector ∈ℝ , and happening with rate ( ) ∈ ℝ + possibly depending on the state ∈ℝ . Depending on the microscopic model at hand, we can expand ( ) in terms of orders of , ( ) = (0) ( )+ (1) ( ) + ( 2 ). tr( ) + (1) , ( ) = 1 , (
with (1) =2 ∑ (1) ( )( ⋅∇ -1). Note that, similar to the derivations for Gaussian SDEs, ( ) plays the role of prefactor along the instanton (as ( ) = ( , ( )) and is used to compute the eventual prefactor correction (which is merely (0)).
As before, we obtain an evolution equation for = ∇∇ by differencing the HJB equation twice,
to be solved with boundary conditions that depend on the scenario, and where subscripts of the Hamiltonian denote differentiation. We can also derive an equation for , .12) This allows one to compute sharp estimates for expectations, probability densities, hitting probabilities and exit times in a similar way as before, replacing the instanton equations (1.6) and its variants with (7.6), and the forward and backward Riccati equations with (7.11) and (7.12).
Consider the following continuous-time MJP, inspired by [8], in which a particle hops on a grid with spacing , ∈ ℤ(see Figure 10), that is the spatial coordinate becomes continuous for and the generator is given by = -1 ( + ( )( ( + ) -( )) + -( )( ( -) -( ))). (7.14) By construction, this process is in detailed balance with respect to the Gibbs distribution ∞ ( )= -1 -2 -1 ( ) , where = ∑ ∈ ℤ -2 -1 ( )
that is ( ) plays the role of the free energy and that of the temperature. In the continuum limit, →0 , the rates (7.13) can be expanded as for the MJP with generator (7.1), with the particle starting at 0 =0at =0 . Solving this problem numerically by our approach amounts to performing the following steps: ̄ ( ) = ( 2 ( ) 2 ( ) ) 1∕2 exp ( 1 2 ∫ 0 ( ( ), ( )) ( ) + (1) ( ( )) ) ×e xp ( --1 ( ∫ -( , )
)) (7.22) This procedure can be carried out for various ∈ℝto, for example, investigate the probability ( ) forrareevents,thatislarge . Displayed as black solid line in Figure 11 is the estimate ̄ ( ) at = 5.12. Here, we choose ( ) = 1 4 4 . For the numerical computation of the instanton, we employed a simple forward Euler scheme integrating the instanton equations (7.6)f o r w a r d and backward multiple times until convergence, with = 512 time discretization points, and ∆ = 10 -2 temporal resolution. The blue markers compare the instanton prediction against a numerical computation of the actual stochastic process for finite (but small) ≪1by using the Gillespie algorithm [27]. Numerical parameters are =5⋅1 0 -4 , and we took samples =10 4 samples to estimate ( ).
In this paper, we have proposed explicit formulas to calculate the prefactor contribution of expectations, probabilities, probability densities, exit probabilities, and mean first passage times for stochastic processes, both at finite time and over the invariant measure. The approach gives sharp estimates in the small noise limit, corresponding to the next order (pre-exponential) term to the Freidlin-Wentzell large deviation limit. This allows us to compute these probabilistic quantities in absolute terms, that is including normalization constants, for finite noise values, instead of merely producing the exponential scaling. This feature makes the approach valuable whenever full quantitative estimates of probabilities are required, as is the case in almost all applications, for example in physics, chemistry, biology, and engineering. From a physical point of view the prefactor formulas given above represent an explicit evaluation of the fluctuation determinant. In principle, these ideas have been formulated multiple times in various contexts (compare Section 1.3), from quantum field theory over linear-quadratic control to calculus of variations. As noted by Schulman [46]: "Methods for handling the quadratic Lagrangian are legion and have been well developed since the earliest work on path integrals. Oddly enough, papers on the subject continue to appear and may give some historian of science material for a case history on the nondiffusion of knowledge." Importantly, though, numerical methods for the calculation of prefactors in the general setup we consider here remain, to the best of our knowledge, mostly nonexistent. Concretely, for the prefactor terms we (i) formulate algorithms suitable for their explicit numerical computation, (ii) phrase them for the more general class of Lagrangians encountered in LDT, and (iii) treat the case of the invariant measure and infinite time horizon. Computations on stochastic partial differential equations highlight the fact that our results produce correct results even in the irreversible infinite dimensional case, and are efficient enough to be used for quantitative estimates of probabilistic quantities in regimes where direct sampling is inaccessible.
where we used = + ̂ to change integration variable with
10970312, 2024, 4, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/cpa.22177 by College Of Staten Island, Wiley Online Library on [26/02/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
=1∕4,and =1 0 -2 ,aswellas = 128. Right: Comparison between instanton and the result from direct numerical simulations, conditioning on hitting the set . Notably, the mean realization recovers the instanton, but variances are very high around = 1∕32, where the concentration fluctuations are inserted, and also high around =1 ∕ 2 , where the velocity field compresses the fluctuations of the concentration field.