Scalable methods for computing sharp extreme event probabilities in infinite-dimensional stochastic systems

We introduce and compare computational techniques for sharp extreme event probability estimates in stochastic differential equations with small additive Gaussian noise. In particular, we focus on strategies that are scalable, i.e. their efficiency does not degrade upon temporal and possibly spatial refinement. For that purpose, we extend algorithms based on the Laplace method for estimating the probability of an extreme event to infinite dimensional path space. The method estimates the limiting exponential scaling using a single realization of the random variable, the large deviation minimizer. Finding this minimizer amounts to solving an optimization problem governed by a differential equation. The probability estimate becomes sharp when it additionally includes prefactor information, which necessitates computing the determinant of a second derivative operator to evaluate a Gaussian integral around the minimizer. We present an approach in infinite dimensions based on Fredholm determinants, and develop numerical algorithms to compute these determinants efficiently for the high-dimensional systems that arise upon discretization. We also give an interpretation of this approach using Gaussian process covariances and transition tubes. An example model problem, for which we provide an open-source python implementation, is used throughout the paper to illustrate all methods discussed. To study the performance of the methods, we consider examples of stochastic differential and stochastic partial differential equations, including the randomly forced incompressible three-dimensional Navier–Stokes equations.

The estimation of extreme event probabilities in complex stochastic systems is an important problem in applied sciences and engineering, and is difficult as soon as these events are too rare to be easily observable, but at the same time too impactful to be ignored.Examples of such events studied in the recent literature include rogue waves (Dematteis et al. 2018) and wave impacts on an offshore platform (Mohamad and Sapsis 2018), heat waves and cold spells (Gálfi et al. 2019, Ragone et al. 2018), intermittent fluctuations in turbulent flows (Fuchs et al. 2022) and derivative pricing fluctuations in mathematical finance (Friz et al. 2015).A broad perspective on extreme event prediction can be found in Farazmand and Sapsis (2019).Methods to estimate extreme events typically rely on Monte Carlo simulations, including importance sampling (Bucklew 2013), subset simulation (Au and Beck 2001) or multilevel splitting methods (Budhiraja and Dupuis 2019).
A possible theoretical framework to assess extreme event probabilities, which we will follow in this work, is given by large deviation theory (LDT) (Dembo andZeitouni 1998, Varadhan 1984).This approach allows to estimate the dominant, exponential scaling of the probabilities in question through the solution of a deterministic optimization problem, namely finding the most relevant realization of the stochastic process for a given outcome.This realization is sometimes called instanton, inspired by theoretical physics.For stochastic processes described by stochastic differential equations (SDEs), the relevant theory has been formulated by Freidlin and Wentzell (2012), and can be extended to many stochastic partial differential equations (SPDEs).The computational potential of this formulation has been reviewed by Grafke and Vanden-Eijnden (2019).
In addition to the exponential scaling provided by LDT, it is often desirable to obtain asymptotically sharp, i.e. asymptotically exact probability estimates.This requires the evaluation of a pre-exponential factor in addition to the usual leading-order large deviation result, when interpreting LDT as a Laplace approximation.On the theoretical side, there exist multiple results for such precise Laplace asymptotics for general SDEs (Azencott 1982, Ben Arous 1988, Deuschel et al. 2014, Ellis and Rosen 1982, Piterbarg and Fatalov 1995) and certain SPDEs requiring renormalization (Berglund et al. 2017, Friz andKlose 2022), which, however, typically do not include an actual evaluation of the abstract objects in terms of which they are formulated.We concentrate on the case of SDEs or well-posed SPDEs with additive noise here, where computing the leading-order prefactor amounts to evaluating a Fredholm determinant of an integral operator.
Approach.In this paper, we present a sharp and computable probability estimate for tail probabilities P[ f (X T ) ≥ z], i.e. a real-valued function f of a diffusion process (X t ) t∈[0,T ] with state space R n and exceeding a given threshold z at final time T (see Figure 1 for an example of this setup).We demonstrate that in a way to be made precise later on, with real-valued functions C, called the (leading-order) prefactor, and I, called the rate function.The latter is determined through the solution of a constrained optimization problem: where formally η = dB t /dt is the time derivative of the Brownian motion (B t ) t∈[0,T ] , and X T depends on η through (1).The prefactor C is then expressed as a Fredholm determinant of a linear operator which contains the solution of the minimization problem (3), the instanton η z , as a background field and acts on paths δ η∶[0,T ] → R n .We show how to evaluate this operator determinant numerically for general SDEs and SPDEs, and demonstrate through multiple examples that it is possible to do so even for very highdimensional systems with n ≫ 1 arising, for instance, after spatial discretization of an SPDE.Our approach is based on computing the dominant eigenvalues of the trace-class integral operator entering the Fredholm determinant.Related literature.In the physics literature, the leadingorder prefactor computation corresponds to the evaluation of initial position x sample paths of (Xt) t∈[0,T ] final positions XT event set FIG. 1. Visualization of extreme event set (red), a sample path that, from a given initial condition, ends in the extreme event set at final time T (orange) and two typical sample paths that do not end in the event set (blue and green).The gray lines are field lines of the drift vector field b.This is the reason why paths that end in the event set are rare, since the noise must act against the flow of the deterministic vector field b to push the system into the extreme event set.Details of this example problem, which is used throughout the paper as illustration, will be given in Section 1.2 and more rare event paths and the instanton are shown in Figure 2. The implementation of this example is available from a public GitHub repository (Schorlepp, Tong, Grafke and Stadler 2023).
Gaussian path integrals, which is a classical topic in quantum and statistical field theory (Zinn-Justin 2021).There are multiple references dealing with the evaluation of such integrals for the class of differential operators that is necessary for SDEs, such as Corazza and Fadel (2020), Nickelsen and Engel (2011), Papadopoulos (1975).In accordance with these approaches, in the last years, numerical leading order prefactor computation methods for general SDEs and SPDEs via the solution of Riccati matrix differential equations have been established (Bouchet and Reygner 2022, Ferré and Grafke 2021, Grafke et al. 2021, Schorlepp et al. 2021, Schorlepp, Grafke and Grauer 2023).An early example using a similar method is given by Maier and Stein (1996).All of these papers have in common that the leading order prefactor can be evaluated in a closed form by solving a single matrix valued initial or final value problem, thereby bypassing the need to compute large operator determinants directly.We briefly introduce this method in this paper, relate it to the -in some sense complementary -Fredholm determinant prefactor evaluation based on dominant eigenvalues, and discuss possible advantages and disadvantages.We note that for SDEs with low-dimensional state space, it can also be feasible to compute the differential operator determinants, that are otherwise evaluated through the Riccati matrices, directly by discretizing the operator into a large matrix and numerically calculating its determinant, which has been carried out e.g. by Psaros and Kougioumtzoglou (2020), Zhao et al. (2022).
Another perspective on the precise Laplace approximation used in this paper is provided by the so-called secondorder reliability method (SORM), which is used in the engineering literature to estimate failure probabilities, as reviewed e.g. by Breitung (2006), Rackwitz (2001).For example, the asymptotic form of the extreme event probabil-ities in this paper corresponds to the standard form stated by Breitung (1984).In this sense, the method proposed in this paper can be regarded as a path space SORM, carried over to infinite dimensions for the case of additive noise SDEs.The connection of precise LDT estimates to SORM for finite-dimensional parameter spaces has also been pointed out by Tong et al. (2021).
In studies of rare and extreme event estimation, Monte Carlo simulations are commonly used, and various sampling schemes have been designed, some of which have been modified and adapted to systems involving SDEs.These include various importance sampling estimators which can be associated e.g. with the solution to deterministic optimal control problems along random trajectories (Vanden-Eijnden and Weare 2012), with the instanton in LDT (Ebener et al. 2019), or build on stochastic Koopman operator eigenfunctions (Zhang et al. 2022).The method we propose takes a different perspective from these sampling methods-it does not involve sampling, and is only asymptotically exact.
Contributions and limitations.The main contributions of this paper are as follows: (i) Generalizing SORM to infinite dimensions, we introduce a sampling-free method to approximate extreme event probabilities for SDEs (and SPDEs) with additive noise.The method is based on the Laplace approximation in path space and uses second-order information to compute the probability prefactor.(ii) While such precise Laplace asymptotics for SDEs are known on a theoretical level, we show how to evaluate them numerically in a manner that is straightforward to implement and is scalable, i.e. it does not degrade with increasing discretization dimension.We illustrate the method on a high-dimensional nonlinear example, namely estimating the probability of high strain rate events in a three-dimensional stochastic Navier-Stokes flow.(iii) On the theoretical level, we explore the relationship between the proposed eigenvalue-based approach for calculating the prefactor and Riccati methods from stochastic analysis and stochastic field theory.We examine the advantages of each method and provide an interpretation of the involved Gaussian process using transition tubes towards the extreme event, i.e. the expected magnitude and direction of fluctuations on the way to an extreme outcome.
The approach taken in this paper also has some limitations: (i) While we find the probability estimates including the leading-order prefactor to be quite accurate when compared to direct Monte Carlo simulations when these are feasible, these estimates are approximations and only asymptotically exact in the limit as z → ∞.To obtain unbiased estimates, one can e.g.use importance sampling.The instanton and the second variation eigenvalues and eigenvectors can be used as input for such extreme event importance sampling algorithms (Ebener et al. 2019, Tong and Stadler 2022, Tong et al. 2021).(ii) We limit ourselves to SDEs with additive Gaussian noise.For SDEs with multiplicative noise (or singular SPDEs), the leading-order prefactor is more complicated, as the direct analogy to the finite-dimensional case gets lost (Ben Arous 1988).Nevertheless, extensions of the eigenvalue-based prefactor computation proposed here can likely be made, but are beyond the scope of this paper.(iii) The proposed approach assumes that the differential equation-optimization (3) has a unique solution that can be computed.For non-convex constraints, uniqueness may be difficult to prove or may not hold.However, in the examples we consider, we seem to be able to identify the global minimizer reliably by using several different initializations in the minimization algorithm and, if we find different minimizers, by choosing the one corresponding to the smallest objective value.The proposed approach can also be generalized to multiple isolated and continuous families of minimizers (Ellis andRosen 1981, Schorlepp, Grafke andGrauer 2023).
Notation.We use the following notations throughout the paper: The state space dimension is always written as n, a possible time discretization dimension of the interval [0,T ] as n t , and N is exclusively used in section 1.1 for the motivation of our results via random variables in R N .We denote the Euclidean norm and inner product in R N by ∥⋅∥ N and ⟨⋅,⋅⟩ N , respectively, and the L 2 norm and scalar product for R n -valued functions defined on [0,T ] by Convolutions are written as * .The subscript or argument z ∈ R always represents the dependency on the observable value e.g. of the minimizer η z , Lagrange multiplier λ z and projected second variation operator A z , as well as the observable rate function I F (z) and prefactor C F (z).The identity map is in general denoted by Id, and the identity matrix and zero matrix in R N are written as 1 N×N and 0 N×N .The superscript ⊥ always denotes the orthogonal complement, with v ⊥ ∶= (span({v})) ⊥ .Functional derivatives with respect to η ∈ L 2 ([0,T ],R n ) are denoted by δ /δ η.Determinants in R N , as well as Fredholm determinants, are written as det, whereas regularized differential operator determinants are written as Det with the boundary conditions of the operator as a subscript.For two real functions g and h, we write if the functions g and h are asymptotically equivalent as ε ↓ 0. By an abuse of terminology, we use the term "instanton" in this paper to refer to the large deviation minimizer η z for finite-dimensional parameter spaces, and also to both the instanton noise trajectory (η z (t)) t∈[0,T ] and the instanton state variable trajectory (φ z (t)) t∈[0,T ] in the infinitedimensional setup.
We start with a more precise explanation of the concepts described in this introduction in sections 1.1 and 1.2, before summarizing the structure of the rest of the paper at the end of section 1.2.

Laplace method for normal random variables in R N
We start with the finite dimensional setting, following Dematteis et al. (2019), Tong et al. (2021): We consider a collection of N random parameters η ∈ R N that are standard normally distributed, and are interested in a physical observable, described by a function F∶R N → R, that describes the outcome of an experiment under these random parameters.Note that restricting ourselves to independent standard normal variables is not a major limitation as F may include a map that transforms a standard normal to another distribution.To give an example that fits into this setting, η could be all parameters entering a weather prediction model, and F then constitutes the mapping of the parameters to some final prediction, such as the temperature at a given location in the future.Note that the map F may be complicated and expensive to evaluate, e.g.requiring the solution of a PDE.
We are interested in the probability that the outcome of the experiment exceeds some threshold z, i.e.P(z) = P[F(η) ≥ z].Since here z is assumed large compared to typically expected values of F(η), we call P(z) the extreme event probability.To be able to control the rareness of the event, we introduce a formal scaling parameter ε > 0 and consider ε ≪ 1 to make the event extreme by defining . This allows us to treat terms of different orders in ε perturbatively in the rareness of the event and is more amenable to analysis than rareness due to z → ∞.
In the following, we will thus consider z as a fixed constant, while discussing the limit ε → 0. Since η is normally distributed, the extreme event probability is available as an integral, by integrating all possible η that lead to an exceedance of the observable threshold (as identified by the indicator function 1), weighed by their respective probabilities given by the Gaussian densities.Directly evaluating the integral in ( 5) is typically infeasible for complicated sets {η ∈ R N | F(η) ≥ z} and large N.
The central notion of this paper is the fact that in the limit ε ↓ 0, the integral in (5) can be approximated via the Laplace method, which replaces the integrand with its extremal value, times higher order multiplicative corrections.The corrections at leading order in ε amount to a Gaussian integral that can be solved exactly.In effect, the integral (5) is approximated by the probability of the most likely event that exceeds the threshold, multiplied by a factor that takes into account the event's neighborhood.
To make things concrete, we make the following assumptions on F ∈ C 2 (R N ,R) for given z > F(0): )).Necessarily, η z ∈ F −1 ({z}) lies on the boundary, F(η z ) = z, and there exists a Lagrange multiplier λ z ≥ 0 with η z = λ z ∇F(η z ) as a first-order necessary condition.We define the large deviation rate function of the family of realvalued random variables (F( √ εη)) ε>0 at z via ) is positive definite on the (N − 1)dimensional subspace η ⊥ z ⊂ R N orthogonal to the instanton, i.e. we assume a second-order sufficient condition for η z holds.
Then, there is a sharp estimate, in the sense of (4), for the extreme event probability (5) via where the rate function I F determines the exponential scaling, and C F (z) is the z-dependent leading order prefactor contribution that accounts for the local properties around the instanton.Note that the prefactor is essential to get a sharp estimate, which cannot be obtained from mere logasymptotics using only the rate function.The prefactor C F (z) can explicitly be computed via where pr A brief derivation of this result, analogous to the computations of Tong et al. (2021), is included in appendix A1 for completeness.It is also directly equivalent to the standard form of the second order reliability method, as derived e.g. by Breitung (1984).Geometrically, it corresponds to replacing the extreme event set {η ∈ R N | F(η) ≥ z} by a set bounded by the paraboloid with vertex at the instanton η z , the axis of symmetry in the direction of ∇F(η z ) and curvatures adjusted to be the eigenvalues of the −∥∇F∥ −1 -weighted Hessian of F at η z .For the weather prediction example, equations ( 7) and (8) mean the following: We could estimate (5) by performing a large number of simulations of the weather model with a random choice of parameters to obtain statistics on an extremely high temperature event.Instead, we solve an optimization problem over parameters to compute only the single most likely route to that large temperature.When the desired event is very extreme, such a situation can only be realized when all simulated physical processes conspire in exactly the right way to make the extreme temperature event possible.Consequently, only a narrow choice of model parameters and corresponding sequence of events remains that can contribute to the extreme event probability: precisely the instanton singled out by the optimization procedure.The probability of the extreme event is then well approximated by perturbations around that single most likely extreme outcome.
Next, we generalize the statement (7) to the infinitedimensional setting encountered in continuous time stochastic systems.Intuitively, for temporally evolving systems with stochastic noise, there is randomness at every single instance in time, which implies an infinite number of random parameters to optimize over.We generalize the above strategy to the important case of SDEs in R n driven by small additive Gaussian noise, and assemble and compare computational methods to compute I F and C F numerically, even for very large spatial dimensions n stemming from semidiscretizations of multi-dimensional SPDEs.

Generalization to infinite dimensions for SDEs with additive noise
As a stochastic model problem, we consider the SDE on the time interval [0,T ] with a deterministic initial condition and n ∈ N, ε > 0. The drift vector field b∶R n → R n , assumed to be smooth, may be nonlinear and non-gradient.The constant matrix σ ∈ R n×n is not required to be diagonal or invertible.The SDE is driven by a standard n-dimensional  10) that satisfy f (X(T ),Y (T )) ≥ z with z = 3 (red set) and ε = 0.5.Using Euler-Maruyama steps with an integrating factor with step size ∆t = 5 ⋅ 10 −4 , we repeatedly simulated (10) until 100 such rare trajectories were found.The dashed blue line is the state variable instanton trajectory φz, solution of ( 19) with the optimal ηz as forcing.As in Figure 1, the gray lines are field lines of the drift vector field b.
Brownian motion B = (B t ) t∈[0,T ] .We limit ourselves to the estimation of extreme event probabilities (due to small noise ε) of the random variable f (X ε T ), where f ∶R n → R is a smooth, possibly nonlinear observable of the process X ε at final time t = T .
A concrete example of this type of system, already alluded to in the first section, is shown in figure 2. It is given by the SDE The streamlines in the figure show the motion taken by deterministic trajectories of the model at ε = 0. Small magnitude stochasticity in the form of Brownian noise is added, and we ask the question: What is the probability P ε F (z), as defined below in (13), that the system ends up, at time T = 1, in the red shaded area in the top right corner, given by f (x,y) = x + 2y ≥ z = 3?After approximately 1.2 ⋅ 10 7 simulations, 100 such trajectories are found, with some of them shown in light orange in figure 2. These can be considered typical realizations for this extreme outcome, and allow us to estimate P ε F (z) ∈ [6.71 ⋅ 10 −6 ,9.97 ⋅ 10 −6 ] as a 95% confidence interval.While in principle the same approach could be applied to much more complicated stochastic models, such as SPDEs arising in atmosphere or ocean dynamics, it quickly becomes infeasible due to the cost of performing such a large number of simulations.
Instead, we generalize the strategy outlined in the previous section.For the derivation, we make the following, compared to the finite-dimensional case stronger assumptions for technical reasons.To formulate them, we introduce the solution map Then, we assume for all z ∈ R: 1.There is a unique instanton on the z-levelset of F, as a first-order necessary condition.
We define the large deviation rate function for the observable f as 2. The map from observable value to minimizer z ↦ η z is C 1 .In particular 4. The rate function I F is twice continuously differentiable and strictly convex, i.e.I ′′ F > 0.
Under these assumptions and using existing theoretical results on precise Laplace asymptotics for small-noise SDEs, in appendix A2 we sketch a derivation of the following result: For the extreme event probability with z > F(0), the asymptotically sharp estimate (7) holds in the same way as before.The leading order prefactor is now given by where det is now a Fredholm determinant, the second variation δ 2 F/δ η 2 of the solution map F at η = η z is a linear trace-class operator on L 2 ([0,T ],R n ), and pr denotes orthogonal projection in L 2 ([0,T ],R n ).
Applied to the model SDE (10), we must first compute the optimal noise realization η z = (η z (t)) t∈[0,T ] , which has a corresponding optimal system trajectory φ z = (φ z (t)) t∈[0,T ] .This optimal trajectory, shown blue dashed in figure 2, describes the most likely evolution of the SDE (10) from the initial condition (0,0) into the shaded region in the upper right corner, thus leading to an event f (X(T ),Y (T )) ≥ z.Second, through equation ( 14), we can compute the prefactor correction for this optimal noise realization.Inserted into equation ( 7), we obtain P ε=0.5 F (z = 3) ≈ 8.94 ⋅ 10 −6 as an asymptotic, sampling-free estimate, which falls into the estimated interval obtained with direct sampling.The source code to reproduce all results for this example is available in a public GitHub repository (Schorlepp, Tong, Grafke and Stadler 2023).
We add some remarks on the setting: 1. We focus on SDEs with additive noise (9) for simplicity.For the more general case of ordinary Itô SDEs with multiplicative noise σ = σ (X ε t ), the leading order prefactor can still be computed explicitly, but involves a regularized Carleman-Fredholm determinant det 2 (see Simon (1977) for a definition) instead of a Fredholm determinant det, because the second variation of F is no longer guaranteed to be trace-class (Ben Arous 1988).The direct analogy to the finite-dimensional case is only possible for additive noise.
2. We state the theoretical result and computational strategy for ordinary stochastic differential equations, but will also apply them numerically to SPDEs with additive, spatially smooth Gaussian forcing.In this case, we expect a direct generalization of the results for SDEs to hold.
3. Without any additional work, we also obtain a sharp estimate, in the sense of ( 4), for the probability density function From a practical point of view, the remaining question is how to evaluate ( 12) and ( 14), given a general and possibly high-dimensional SDE (9).
Main questions and paper outline.In the remainder of this paper, we will specifically answer the following questions: • How to find the minimizer η z to the differential equation constrained optimization problem ( 12) numerically?This question has been treated in detail in the literature for the setup at hand, and we give a brief summary of relevant references in section 2.1.
• How to evaluate the Fredholm determinant in ( 14) numerically?We show in section 2.2 how to use second-order adjoints to compute the application of the projected second variation operator to functions (or, upon discretization, to vectors), which is the basis for iterative eigenvalue solvers.In section 2.4, we discuss how this allows us to treat very large system dimensions n as long as the rank of σ remains small.
• How does this prefactor computation based on the dominant eigenvalues of the projected second variation operator theoretically relate to the alternative approach using symmetric matrix Riccati differential equations mentioned in the introduction?What are the advantages and disadvantages of the different approaches?We comment on these points in sections 2.3 and 2.4.
• What is the probabilistic interpretation of the quantities encountered when evaluating ( 12) and ( 14)?In how far can they be observed in direct Monte Carlo simulations of the SDE (9)?This is the content of section 3.
After these theoretical sections, illustrated throughout via the model SDE (10), we present two challenging examples in section 4: The probability of high waves in the stochastic Korteweg-De Vries equation in section 4.1, and the probability of high strain events in the stochastic three-dimensional incompressible Navier-Stokes equations in section 4.2.All technical derivations can be found in Appendix A.

NUMERICAL RATE FUNCTION AND PREFACTOR EVALUATION
In this section, we show how the instanton and prefactor for the evaluation of the asymptotic tail probability estimate ( 7) can be computed in practice for a general, possibly high-dimensional SDE (9), and illustrate the procedure for the model SDE (10).Both finding the instanton (section 2.1) and the prefactor (section 2.2) require the solutions of differential equations of a complexity comparable to the original SDE.They therefore become realistic to evaluate numerically even for fairly large problems, provided tailored methods are used, as summarized in section 2.4.Additionally, we compare the adjoint-based Fredholm determinant computation to the approach based on matrix Riccati differential equations in sections 2.3 and 2.4.

First variations and finding the instanton
Here, we discuss the differential equation-constrained optimization problem that determines the instanton noise η z , and briefly review how it can be solved numerically.We reformulate the firstorder optimality condition by evaluating the first variation using an adjoint variable as reviewed by Hinze et al. (2009), Plessix (2006).For any , where the adjoint variable θ (also called conjugate momentum) is found via solving With a = σ σ ⊺ , we recover from (18) the well-known instanton equations, formulated only in term of the state variable φ z and its adjoint variable θ z with optimal noise η z = σ ⊺ θ z : The rate function is given by When formulating the optimization problem in the state variable φ instead of the noise η, the instanton equations (20) are directly obtained as the first-order necessary condition for a minimizer of the Freidlin-Wentzell (Freidlin and Wentzell 2012) action functional S with The numerical minimization of this functional for both ordinary and partial stochastic differential equations is discussed e.g. by E et al. ( 2004), Grafke, Grauer and Schäfer (2015), Grafke, Grauer and Schindel (2015), Grafke and Vanden-Eijnden (2019), Schorlepp et al. (2022).Conceptually, the minimization problem ( 17) is a deterministic distributed optimal control problem on a finite time horizon with a final time constraint on the state variable (Herzog andKunisch 2010, Lewis et al. 2012).The final-time constraint can be eliminated e.g. using penalty methods.Alternatively, for a convex rate function, a primal-dual strategy (Boyd and Vandenberghe 2004) with minimization of 1 2 ∥⋅∥ 2 − λ F at fixed λ can be used.If estimates for a range of z are desired, one can solve the dual problem for various λ , which effectively computes the Legendre-Fenchel transform I * F (λ ), and invert afterwards.If the rate function is not convex, the observable f can be reparameterized to make this possible (Alqahtani and Grafke 2021).To solve the unconstrained problems of the general form min 1 2 ∥⋅∥ 2 − λ (F − z) + µ 2 (F − z) 2 , gradient-based methods with an adjoint evaluation ( 19) can be used, e.g.Schorlepp et al. (2022) use an L-BFGS solver.Simonnet (2022) used a deep learning approach instead.For high-dimensional problems such as multi-dimensional fluids, it may be necessary to use checkpointing for the gradient evaluation, and to use rankσ ≪ n if applicable to reduce memory costs (Grafke, Grauer and Schindel 2015).We comment on this point in more detail in section 2.4.Using second order adjoints as in the next section would also make it possible to implement a Newton solver, cf.Cioaca et al. (2012), Hinze and Kunisch (2001), Hinze et al. (2006), Sternberg and Hinze (2010).
For the model SDE (10), the instanton equations (20) read We implemented a simple gradient descent (preconditioned with a −1 ) using adjoint evaluations of the gradient and an Armijo line search (available in the GitHub repository (Schorlepp, Tong, Grafke and Stadler 2023)) to find the instanton for the model SDE (10).The state equation is discretized using explicit Euler steps with an integrating factor, and the gradient is computed exactly on a discrete level, i.e. "discretize, then optimize".To find the instanton for a given z, we use the augmented Lagrangian method.For each subproblem at fixed Lagrange multiplier λ and penalty parameter µ, gradient descent is performed until the gradient norm has been reduced by a given factor compared to its initial value.All of these aspects are summarized in more detail by Schorlepp et al. (2022).The resulting optimal state vari-able trajectory φ z for z = 3 for the model SDE ( 10) is shown in figure 2.

Second variations and prefactor computation via dominant eigenvalues
Similarly to the previous section, the second variation is also readily evaluated in the adjoint formalism.With this prerequisite, we are able to use iterative eigenvalue solvers to approximate the Fredholm determinant det(Id−A z ).For a comprehensive introduction to the numerical computation of Fredholm determinants, as well as theoretical results on approximate evaluations using integral quadratures, see Bornemann (2010).However, in contrast to Bornemann (2010), we deal with possibly spatially high-dimensional problems, such as the example in section 4.2.Hence, we use iterative algorithms to compute the dominant eigenvalues to keep the number of operator evaluations manageable.
Another application of the adjoint state method shows that applying the second functional derivative of the solution map Here, we use the short-hand notation The trajectories φ and θ in (23) are determined via (19) from η.Note that the second order equations (23) are simply the linearization of ( 19).Together with the projection operator pr η ⊥ z acting as for t ∈ [0,T ], we are now in a position to evaluate the application of the operator A z , as defined in ( 16), to any function δ η∶[0,T ] → R n .Denoting the eigenvalues of the trace-class operator A z by µ ), the Fredholm determinant in the prefactor ( 14) is given by det → 0 in such a way that the product converges.An iterative eigenvalue solver relying solely on matrix-vector multiplication, thus avoiding the explicit storage of the possibly large discretized operator A z as an (n t ⋅ n) × (n t ⋅ n) matrix, can now be used numerically to find a finite number of dominant eigenvalues of A z with absolute value larger than some thresholds, and approximate det(Id−A z ) using these.
For the model example SDE (10), linearizing the state and first order adjoint equations ( 19), the second order adjoint with largest absolute value of Az for the example SDE (10) with z = 3. Discretization of (25) was done with step size ∆t = 5 ⋅ 10 −4 , hence the dimension of the discretized path space variables is 4000 here.Main figure: absolute value of the eigenvalues µ z ) for different m as an approximation for the Fredholm determinant det(Id −Az).We see that the eigenvalues rapidly decay to zero in this example, that similarly, the cumulative product in the inset quickly converges, and that the final estimate det( equations for (10) become We implemented a simple Euler solver for these equations for a given discretized input vector δ η ∈ R 2(n t +1) in the python code (Schorlepp, Tong, Grafke and Stadler 2023) as a subclass of scipy.sparse.linalg.LinearOperator.
To set up this operator, we supply the instanton data (φ z ,θ z ,λ z ) ∈ R 2(n t +1) ×R 2(n t +1) ×R as found using the methods of the previous section 2.1.The LinearOperator class, for which we only need to supply a matrix vector multiplication method instead of having to store the full matrix ∈ R 2(n t +1)×2(n t +1) , can then be used with any iterative eigenvalue solver.Here, we use the implicitly restarted Arnoldi method of ARPACK (Lehoucq et al. 1998), wrapped as scipy.sparse.linalg.eigs in python.Note that in this example, storing the full matrix would be feasible, and the Riccati method discussed in the next section is faster to compute the prefactor.However, we are interested in a scalable approach for large n, where, as discussed in section 2.4 and shown in section 4, the Riccati approach becomes infeasible.We show the results of computing 200 eigenvalues with largest absolute value of the projected second variation operator A z for z = 3 in figure 3.

Alternative: prefactor computation via matrix Riccati differential equations
In Appendix A3, we motivate via formal manipulations that the prefactor ( 14) can also be expressed via the following ratio of zeta-regularized functional determinants (Ray and Singer 1971) of second order differential operators, instead of a Fredholm determinant of an integral operator.This prefactor expression is more natural from the statistical physics point of view, where path integrals in the field variable φ instead of the noise η are typically considered, cf.Zinn-Justin (2021).We obtain in accordance with Schorlepp, Grafke and Grauer (2023), where it was derived directly through path integral computations.Here, Ω is the Jacobi operator of the Freidlin-Wentzell action functional as defined in the appendix A3, and the subscript of the zeta-regularized determinants Det denotes the boundary conditions under which the determinants of the differential operators are computed.Naively evaluating the determinant ratio in ( 26) by numerically finding the eigenvalues of the appearing differential operators is typically not feasible.This is due to the fact that both operators posses unbounded spectra with the same asymptotic behavior of the eigenvalues, which requires computing the smallest eigenvalues of both operators.A threshold for this computation is difficult to set, and while the eigenvalues of both operators should converge to each other as they increase, numerical inaccuracies tend to increase for the larger eigenvalues.Fortunately, there exists theoretical results regarding the computation of such determinant ratios exactly and in a closed form by solving initial value problems (Forman 1987, Gel'fand and Yaglom 1960, Kirsten and McKane 2003, Levit and Smilansky 1977).Using the results of Forman (1987), the prefactor (26) can be computed by solving the symmetric matrix Riccati differential equation for Q z ∶[0,T ] → R n×n and then evaluating This result in terms of a Riccati matrix differential equation is also natural from a stochastic analysis perspective (WKB analysis of the Kolmogorov backward equation (Grafke et al. 2021)), or a time-discretization of the path integral perspective (recursive evaluation method (Schorlepp et al. 2021)).To give intuition for the Riccati differential equation ( 27), note that by letting Q z = γζ −1 with γ(0) = 0 n×n and ζ (0) = 1 n×n , the approach amounts to solving as an initial value problem, whereas the eigenvalue problem δ η 2 δ η = µδ η corresponds to the boundary value problem This means that to evaluate the functional determinant prefactor via the Riccati approach, we consider functions in the kernel of the operator Id−λ δ 2 F/δ η 2 , i.e. eigenfunctions belonging to the eigenvalue 0, but under modified boundary conditions of the operator.In practice, instead of finding the dominant eigenvalues of the integral operator A z of section 2.2 that acts on functions δ η∶[0,T ] → R n , we can integrate a single matrix-valued initial value problem for Q z ∶[0,T ] → R n×n as presented in this section.Even though the Riccati equation ( 27), in contrast to the linear system (30), is a nonlinear differential equation, it is nevertheless advisable to solve (27) instead of (30) numerically because the equation for ζ in (30) has to be integrated in the unstable time direction for the right-hand side term −∇b(φ ) ⊺ ζ .Note also that, depending on the system and observable at hand, the solution of the Riccati equation (27) may pass through removable singularities in (0,T ) whenever ζ (t) in (30) becomes non-invertible, hence direct numerical integration of (27) may require some care (see Schiff and Shnider (1999) and references therein).
For the two-dimensional model SDE (10), the forward Riccati equation for the symmetric matrix where [...] stands for a repetition of the preceding term.We solve the Riccati equation with Euler steps with integrating factor in (Schorlepp, Tong, Grafke and Stadler 2023), and use it to evaluate (28).We do not encounter any numerical problems or singularities in this example.The result for C F (z = 3) agrees with the Fredholm determinant computation using dominant eigenvalues in the previous section 2.2.

Computational efficiency considerations
In this section, we compare the two prefactor computation methods of sections 2.2 and 2.3 using either dominant eigenvalues of the trace-class operator A z evaluated via (23), or the Riccati matrix differential equation ( 27), in terms of their practical applicability as well as computational and memory cost for large system dimensions n ≫ 1.
For the eigenvalue-based approach, we know that ∏ theory, but it is difficult to give bounds on the required number of eigenvalues for an approximation of the Fredholm determinant to a given accuracy.In all examples considered in this paper, at most a few 100 eigenvalues turned out to be necessary for accurate results, even for the three-dimensional Navier-Stokes equations in section 4.2 as a high-dimensional (n = 3 ⋅ 128 3 ≈ 6.3 ⋅ 10 6 ) and strongly nonlinear example.The number of dominant eigenvalues of A z to be computed to achieve a desired accuracy is robust with respect to the temporal resolution and only depends on the (effective, see below) dimension of the system and the level of nonlinearity in the system.In any case, to obtain m eigenvalues of A z with largest absolute value, iterative eigenvalue solvers, either using Krylov subspace methods or randomized algorithms, typically require a number of evaluations of the operator that is equal to a constant times m (Halko et al. 2011).Each evaluation of A z consists of solving two ODEs or PDEs ( 23) with comparable computational complexity to the original SDE.We comment on memory requirements below.
Compared to this, the Riccati approach requires the numerical solution of a single n × n symmetric matrix differential equation as an initial value problem.If n is small, then this is clearly more efficient than computing m > n/2 eigenvalues.However, there may also be problems with the Riccati approach: On the one hand, this approach requires a strictly convex rate function with I ′′ F (z) > 0 at z, as can be seen from ( 26).If this is not satisfied, then a suitable convexification via reparameterization needs to be carried out on a case-by-case basis (Alqahtani and Grafke 2021).While we assumed that the rate function is convex to derive the prefactor (14) in terms of the Fredholm determinant, this assumption is actually not necessary and the eigenvaluebased approach remains feasible regardless of the convexity of the observable rate function I F .Finally, the eigenvalue approach is easier to interpret, while it is not always immediately clear why the Riccati solution may diverge (removable singularities that can be remedied by a suitable choice of integration scheme versus true singularities due to unstable or flat directions of the second variation at the instanton).
We turn to the memory requirements of the prefactor computations strategies, and in particular to their scaling with the system dimension n.Informally, one can think of the Riccati matrix as defined in the (squared) state space of the SDE, in contrast to the eigenvectors of A z that are defined in the noise space that is potentially lower-dimensional.The Riccati equation then integrates a dense n×n array in time by performing n t consecutive times steps of ( 27) and evaluating (28) along the way.This is difficult to achieve directly as soon as (semi-discretizations of) multi-dimensional SPDEs are considered, which are relevant e.g. for realistic fluid or climate models.Usually, large Riccati matrix differential equations, which also arise e.g. in linear-quadratic regulator problems, are solved within some problem-specific lowrank format, see e.g.Stillfjord (2018).In contrast to this, the vectors on which iterative eigenvalue solvers for the Fredholm-determinant based approach need to operate are in general vectors of size n t × n.
As an important class of examples, we now consider sys-tems with large spatial dimension n ≫ 1, for which, however, only a few degrees of freedom are forced, such that the diffusion matrix a = σ σ ⊺ is singular and rankσ ≪ n.Examples for this include fluid and turbulence models with energy injection only on a compactly supported set of either high or low spatial Fourier modes, or climate models with a limited number of random parameters in the model (Margazoglou et al. 2021).In this case, it is straightforward to exploit the small rank of σ within the eigenvalue-based approach to decrease the memory requirements and apply the method even to very high-dimensional models, which we demonstrate for the randomly forced three-dimensional Navier-Stokes equations in section 4.2 in this paper.The idea is that for the eigenvectors δ η of A z , clearly only rankσ many entries are relevant due to the composition with σ and σ ⊺ .Eigenvalue solvers hence act on n t × rankσ vectors, which should fit into memory.This is similar to the computation of the instanton itself, where only the instanton noise η z as a n t × rankσ vector is computed and stored explicitly, as discussed by Grafke, Grauer and Schindel (2015), Schorlepp et al. (2022).The remaining challenge is then to evaluate A z δ η for given δ η ∈ R n t ×rank σ by solving the second order adjoint equations ( 23), without storing the full, prohibitively large n t × n arrays needed for φ z , γ and θ z .Similar to the gradient itself, evaluated via the first order adjoint approach ( 19), this is possible through (static) checkpointing (Griewank and Walther 2000), as illustrated in figure 4. At the cost of having to integrate the first order adjoint equations repeatedly for each noise vector δ η to which A z is applied, and to recursively solve the forward equations for φ z and γ again and again, the memory requirements for the spatially dense fields are only O (logn t ⋅ n) this way.The same problem is encountered and solved similarly in implementations of Newton solvers for high-dimensional PDEconstrained optimal control problems (Cioaca et al. 2012, Hinze and Kunisch 2001, Hinze et al. 2006, Sternberg and Hinze 2010).All in all, in contrast to the Riccati formalism, this permits an easy and controlled strategy that enables to treat very large spatial dimensions within the Fredholmbased prefactor approach, as long as the diffusion matrix possesses a comparably small rank.Note, however, that it is still necessary that the number of eigenvalues needed to approximate det(Id−A z ) remains small for this approach to be applicable in practice.We show numerically in section 4.2 that this is indeed the case for the three-dimensional Navier-Stokes equations as an example.The discussion of this paragraph, with all relevant scalings of computational and memory costs for the different approaches, is briefly summarized in table I.
In conclusion, we recommend using the Riccati equation only in sufficiently "nice" situations for small to moderate system dimensions n.For such systems and diffusion matrices without low-rank properties, and as long as no additional complications such as non-convex rate functions or removable singularities of the Riccati solution are encountered, it is faster than the eigenvalue-based approach, and better suited to analytical computations or approximations since it only involves the solution of initial value problems, in contrast to the boundary value problems that need to be solved to find eigenfunctions of the projected second variation operator A z .On the other hand, the Fredholm determinant computation through dominant eigenvalues is easier to use and implement, requiring only solvers for the origi-TABLE I. Overview of computational and memory costs for finding the prefactor C F (z) either through solving the Riccati equation ( 27), or through determining m dominant eigenvalues of Az.The system's spatial dimension is denoted by n ≫ 1 and the noise correlation has rank σ ≪ n.In the table, c denotes the computational costs of integrating the original S(P)DE once from t = 0 to t = T .Any multiplicative constants were omitted for the costs listed in the table, the eigenvalue-based approach is assumed to use checkpointing as sketched in figure 4, the computational costs of evaluating the quadratic term in (27) were ignored, and the complete instanton data is assumed to be known.The table shows that the eigenvalue-based approach indeed remains feasible for large n.

Riccati Eigenvalues
Memory costs nal SDE, its adjoint, as well as their linearizations.At the cost of introducing numerical errors and a step size parameter h > 0 that needs to be adjusted, one can also approximate the second variation evaluations via or other finite difference approximations, which does not require implementing any second order variations.In this sense, both the numerical instanton and leading-order prefactor computation can quickly be achieved in a black-box like, non-intrusive way when solvers for the state equation and its adjoint are available.Alternatively, the adjoint solver, as well as solvers for the second order tangent and adjoint equation can be obtained through automatic differentiation (Naumann 2011).We also note that in the context of the second order reliability method, there exist further approximation methods that could be used here for the Fredholm determinant prefactor, e.g. by extracting information from the gradient based optimization method that has been used to find the instanton or design point (Der Kiureghian and De Stefano 1991), or through constructing a non-infinitesimal parabolic approximation to the extreme event set (Der Kiureghian et al. 1987).
In any case, for the scenario of possibly multidimensional SPDEs with low-rank forcing, we argue that the eigenvalue approach is to be preferred as it leads to natural approximations and a simpler implementation.However, we remark that in the case of SDEs with multiplicative noise, or SPDEs with spatially white noise that need to be renormalized such as the Kardar-Parisi-Zhang (KPZ) equation, the Riccati approach remains structurally unchanged (Schorlepp, Grafke and Grauer 2023), whereas the Fredholm determinant expression turns into a Carleman-Fredholm determinant and an additional operator trace (Ben Arous 1988), which could potentially be more costly to evaluate.

PROBABILISTIC INTERPRETATION VIA FLUCTUATION COVARIANCES AND TRANSITION TUBES
In this section, we give an intuitive interpretation for some of the quantities encountered in the previous sections.The Step 1: Solve forward for given η z , δη; store fields φ z , γ at all • Step 2: Solve backward and recompute φ z , γ when not available, store at , △, ▽ recursively FIG. 4. Sketch of the checkpointing procedure used to evaluate the second variation operator δ 2 (λ F)/δ η 2 at ηz, applied to δ η, for large system dimensions n ≫ 1 in a memory-efficient way.The instanton noise ηz, the input noise fluctuation δ η, and the return vector σ ⊺ ζ are all stored as dense (nt + 1, rank σ )-arrays for rank σ ≪ n.Given the instanton noise ηz and a noise fluctuation δ η, the first step consists of solving the state equation and linearized state equation for φz and γ simultaneously forward in time from t = 0 to t = T , and storing the fields φz(t i ) ∈ R n and γ(t i ) ∈ R n only at the logarithmically spaced instances t i = •.Afterwards, the first and second order adjoint equations for θz and ζ are simultaneously solved backwards in time from t = T to t = 0 and σ ⊺ ζ (t i ) is stored for each t i .Whenever φz(t j ) and γ(t j ) are needed for the time integration, but not available in storage already, the two forward equations are solved again from the nearest preceding point in time when they are available, and recursively stored at intermediate steps ◻, △, ▽, . . .All fields φz(t i ) ∈ R n and γ(t i ) ∈ R n that are no longer needed during the backwards integration are deleted from memory.
second variation quantifies the linearized dynamics of the SDE (9) around the most likely realization.This implies that dominating eigenfunctions of the second variation correspond to fluctuation modes that are most easily observable.Below, we confirm this with a simple numerical experiment that relates the eigenfunction information with the transition tube along a rare trajectory.The basic object that we consider in this section is the process (X ε t ) t∈[0,T ] as ε ↓ 0, conditioned on the rare outcome f (X ε T ) = z at final time.In other words, we consider only transition paths between the fixed initial state x ∈ R n and any final state in the target set f −1 ({z}) ⊂ R n .The path on which the transition path ensemble concentrates as ε ↓ 0 is given by the state variable instanton trajectory φ z , i.e. the most likely way for the system to achieve f (X ε T ) = z, since deviations from it are suppressed exponentially (Freidlin and Wentzell 2012).One thus has for the mean of the conditioned process.In this sense, by taking conditional averages of direct Monte Carlo simulations of (9) as ε tends to 0, the instanton trajectory φ z is directly observable, and the mean realization agrees with the most likely one for ε ↓ 0. This procedure is sometimes called filtering, and has been carried out e.g. for the one-dimensional Burgers equation (Grafke et al. 2013), the three-dimensional Navier-Stokes equations (Grafke, Grauer andSchäfer 2015, Schorlepp et al. 2022) and the onedimensional KPZ equation (Hartmann et al. 2021).Using the results of the previous sections, we can, however, make this statement more precise and state a central limit-type theorem for the conditioned fluctuations at order √ ε around the instanton: As ε ↓ 0, the process We show in Appendix A4 that C z is fully determined through the orthonormal eigenfunctions δ η z of the projected second variation operator A z with corresponding eigenvalues µ (i) z and associated state variable fluctuations γ (i) z , the solution of the linearized state equation In particular, computing the eigenvalues and eigenfunctions of A z yields a complete characterization of the conditioned Gaussian fluctuations around the instanton.As detailed in the example below, at small but finite ε, C z can be used to approximate the distribution of transition paths at any time t ∈ [0,T ] as multivariate normal N (φ z (t),εC z (t,t)).Effectively, in addition to the mean transition path at small noise, the instanton φ z , we can also estimate the width and shape of the transition tube around it at any t ∈ [0,T ] without sampling within a Gaussian process approximation of the conditioned SDE; see Vanden-Eijnden (2006) for a general introduction to transition path theory, and Archambeau et al. (2007), Lu et al. (2017) for Gaussian process approximations of SDEs based on minimizing the path space Kullback-Leibler divergence, which, in the smallnoise limit and for transition paths, reduce to the Gaussian process considered here.Furthermore, one can show that the forward Riccati approach of section 2.3 recovers the finaltime state variable fluctuation covariance via This directly follows by adapting the forward Feynman-Kac computation used in remark 4 of Schorlepp et al. (2021) to the present calculation of the covariance function ( 35) at final time t = t ′ = T .Note that both, directly from ( 38), as well as from (37) after a short calculation, carried out in Appendix A5, one can see that these results are consistent with the additional final time boundary condition for the state variable fluctuations almost surely, when conditioning on f (X ε T ) = z.In words, the conditioned Gaussian fluctuations at final time are constrained to the tangent plane of the equi-observable hypersurface f −1 ({z}) at the point φ z (T ).
As in the previous sections, we use the model SDE ( 10) with z = 3 and ε = 0.5 to illustrate these findings.To do this, we compare the PDF of X ε t at different times t, when conditioning on f (X ε T ) = z, as obtained via sampling, to the Gaussian approximation N (φ z (t),εC z (t,t)) that we evaluate using the instanton as well as eigenvalues and eigenfunctions of A z that were computed previously.We use instanton-based importance sampling (Ebener et al. 2019) to generate 10 5 trajectories of (10) that satisfy f (X ε T ) = z up to a given precision f ((X ε T − φ z (T ))/ √ ε) < 0.05; the corresponding code, which again uses Euler steps with an integrating factor and a step size of ∆t = 5 ⋅ 10 −4 , can be found in the GitHub repository (Schorlepp, Tong, Grafke and Stadler 2023).Essentially, instead of using ( 9) directly, we shift the system by the instanton (cf.Tong et al. (2021) for a visualization and further analysis), solve and reweight the samples by The results are shown in figure 5, and we observe good agreement between the sampled conditioned distributions at times t ∈ {0.05,0.25,0.5,0.75,0.95}and the corresponding theoretical small-noise Gaussian approximations.In particular, the deformation of the fluctuation PDF along the instanton trajectory (φ z (t)) t∈[0,T ] is captured by the Gaussian approximation.It is not surprising that the Gaussian approximation works well for the parameters ε,z and T used here, since the probability P ε F (z) in section 1.1 as approximated by the Laplace method also matched the direct sampling estimate.

COMPUTATIONAL EXAMPLES
We now apply the numerical methods introduced in the previous section to two high-dimensional examples involving SPDEs: In section 4.1, we consider the Korteweg-De Vries equation in one spatial dimension, subject to spatially smooth Gaussian noise, and compute precise estimates for the probability to observe large wave heights at one instance in space and time.We compare our asymptotically sharp estimates to direct sampling, and also explicitly compare the two different prefactor computation strategies.Then, we focus on the stochastically forced three-dimensional incompressible Navier-Stokes equations in section 4.2.This is a much higher-dimensional problem, and we demonstrate that the eigenvalue-based prefactor computation indeed remains applicable in practice for this example.Note that both SPDE examples in this section have periodic boundary conditions in space, but this is not a restriction of the method and has merely been chosen for convenience.

Stochastic Korteweg-De Vries equation
To illustrate the instanton and prefactor computation, we study the Korteweg-De Vries (KdV) equation subject to large-scale smooth Gaussian noise.The KdV equation can be considered as a model for shallow water waves, so the problem we are interested in is to estimate the probability of observing large wave amplitudes.Since this is the first PDE example we study and the general theory in the previous sections has only been developed for ODEs, we explicitly state the instanton equations, second order adjoint equations and Riccati equation.We consider a field with constants ν = κ = 4 ⋅ 10 −2 and white-in-time, centered and stationary Gaussian forcing We choose χk = δ |k|,1 /(2π) as the spatial correlation function of the noise η in Fourier space, with ˆdenoting the spatial Fourier transform.Concretely, η(x,t) is then given by η(x,t) = π −1/2 ( Ḃ1 (t)sin(x) + Ḃ2 (t)cos(x)), where B 1 ,B 2 are independent standard one-dimensional Brownian motions.Hence, the forcing only acts on a single large scale Fourier mode, and excitations of all other modes are due to the nonlinearity of the SPDE.As our observable, we choose the wave height at the origin  10) with z = 3 and ε = 0.5 using instanton-based importance sampling (Ebener et al. 2019).We visualize the transition tube information obtained from the eigenvalues and eigenfunctions of the projected second variation operator.The upper left subfigure shows the histogram of the full data set for all times.The remaining subfigures show histograms of the transition paths at specific times t.The black lines, as a comparison, are the level sets of the normal PDF with covariance εCz(t,t), found by evaluating (37) numerically, and mean φz(t).Note that the deformation of the distribution of X ε t , conditioned on f (X ε T ) = z, is captured quite well using the quadratic, sampling-free approximation.
and we want to quantify the tail probability P ε F (z) = P[ f (u ε (⋅,T )) ≥ z] for different z > 0. Note that the effective dimension of the system when formulated in terms of the noise for our choice of noise correlation is small, and we have rankσ = 2 ≪ n = n x for typical spatial resolutions.Unless otherwise specified, we use n x = 1024 for all numerical results in this section, as well as n t = 4000 equidistant points in time, and we expect the prefactor computation in terms of eigenvalues of A z to be more efficient in this example, even though the Riccati approach still remains feasible.
We use a pseudo-spectral code and explicit second order Runge-Kutta steps in time with an integrating factor for the linear terms.The final-time constraint is treated with the augmented Lagrangian method.Denoting the state space instanton by u z with adjoint variable p z and Lagrange multiplier λ z , the first-order necessary conditions at the minimizers read Here, * denotes spatial convolution, which appears due to the stationarity of the forcing.As a starting point, we compute instantons for a range of equidistantly spaced observable values z ∈ [0,30].Knowledge of the instanton for different z gives us access to the rate function I F of the observable, which is shown on the left in figure 6.
In the table in figure 8, we show for fixed z how the value of I F (z) converges when increasing the spatio-temporal resolution, and in particular that the number of optimization steps needed to find the instanton is robust under changes of the numerical resolution, indicating scalability of the instanton computation.The numerical details for these instanton computations are as follows (cf.Schorlepp et al. (2022)): Initial control p ≡ 0 and initial Lagrange multiplier λ = 0; precise target observable value z = 8.39125; 6 logarithmically spaced penalty steps from 1 to 300 for augmented Lagrangian method; optimization is terminated upon reduction of gradient norm by 10 6 ; same (presumably) global minimizer was found for each resolution; discretize-thenoptimize; L-BFGS solver with 4 updates stored; Armijo line search.
Two comments on the instanton computations for this example are in order: Firstly, the observable rate function is non-convex for some z in the interval [1.5,5] (not visible  42) with height observable (44), as obtained from numerical instanton and prefactor computations.Note that the prefactor depends strongly, almost exponentially, on the observable value z in this example.Right: Comparison of LDT estimate (7) for different noise levels ε ∈ {0.1, 1, 10} to direct sampling for the SPDE (42).For each ε, we computed 4 ⋅ 10 4 samples of f (u ε (⋅, T )) to estimate the tail probabilities for various z.The shaded regions are 95% Wilson score intervals (Brown et al. 2001) for the sampling estimate of the tail probabilities.The solid lines show the asymptotically sharp estimate (7) without adjustable parameters.In comparison to this, the dashed lines show just the leading order LDT term exp {−I F (z)/ε} with a constant prefactor (chosen such that the curve matches (7) for large z), which shows that the prefactor C F is absolutely necessary to get useful results in this example at ε > 0.1 and can be understood in regard to the left column of the figure.Results use nx = 1024, nt = 4000 for the instanton computations, pseudo-spectral code with integrating factor, L-BFGS optimization with penalty term for observable; 80 eigenvalues with largest absolute value for Fredholm determinant; stochastic Heun steps with size ∆t = 10 −3 for direct sampling.in the figure).This poses a problem for the dual problem solved at fixed λ without penalty, but is not an issue for the penalty or augmented Lagrangian strategy that we used.Furthermore, this means that the Riccati prefactor computation is not directly applicable in this region, but the Fredholm expression remains valid.Secondly, since it is a priori unclear whether the minimization problem for the instanton has a unique solution (the target functional is quadratic, but the constraint is nonlinear), we started multiple optimization runs for the same z at different random initial conditions.In the KdV system, we found multiple subdominant minima that consist of multiple large wave crests (as opposed to just one for the dominant one, as shown in the top left of figure 9 for one z), but only took the (presumably) global minimizer for subsequent estimates.
To complete the asymptotic estimate of the wave height probability via (7), we further need the prefactor C F (z) for all z, which we compute by finding the dominant eigenvalues of A z as before.We specify the input and output of the linear operator A z only in terms of the two real Fourier modes of the noise that are relevant for this, to remove the memory cost of the eigenvalue solver.The second order adjoint equations (23) for noise fluctuations δ η∶[0,2π] × [0,1] → R for the KdV equation read with A z δ η = χ 1/2 * δ p.In our implementation, we supply the second variation operator with the two real Fourier coefficients (Re δ η 1 (t i )) i=0,...,n t and (Im δ η 1 (t i )) i=0,...,n t , assemble the full fluctuation vector δ η from it, and return χ 1/2 * δ p in the same format after solving (46).As the KdV solutions fit into memory, checkpointing, as discussed in section 2.4, is not necessary.In figure 7, we show the convergence of the determinant det(Id−A z ) for some z's based on the found eigenvalues, thereby demonstrating that a handful of eigenvalues suffices for an accurate approximation of the prefactor.The number of necessary eigenvalues increases only weakly with the observable value z in this example.In addition, figure 8 shows the effect of varying the spatio-temporal resolution (n x ,n t ) on the determinant det(Id−A z ) for one particular observable value of z = 8.4 at a fixed number of computed eigenvalues.We see that as long as the physical problem is resolved, the eigenvalue spectrum does not change much with the resolution, and the determinant converges when increasing the spatio-temporal resolution.This indicates that our methods are scalable, i.e., their cost does not increase with the temporal (and also spatial) discretization beyond the increased cost of the PDE solution.This is a crucial property of the eigenvalue-based prefactor computation and is in contrast with the Riccati approach.
The result for the prefactor C F as a function of z is shown on the bottom left of figure 6.Note that the vertical axis is scaled logarithmically, i.e. the prefactor strongly depends on the observable value.The importance of the prefactor is further confirmed by the comparison of the complete (dots: positive eigenvalues, crosses: negative eigenvalues).
z ) for different m as an approximation for the Fredholm determinant det(Id −Az).We see that the eigenvalues rapidly decay to zero for all z.Similarly, the cumulative product in the inset converges quickly, and the determinant is in fact well-approximated by less than 10 eigenvalues for all z.
asymptotic estimate (7) to the results of direct Monte Carlo simulations on the right in figure 6.For three values of ε ∈ {0.1,1,10}, we performed 4 ⋅ 10 4 respective simulations of the stochastic KdV equation (42) to estimate the tail probability P ε F (z) without approximations.Using both the rate function and prefactor, excellent agreement with the Monte Carlo simulations is obtained.In contrast to this, only using the leading order LDT term exp{−I F (z)/ε} with a constant prefactor leads to a much worse agreement with simulations, and in fact only works reasonably for ε = 0.1.Note also that one can see from these comparisons that the actual effective smallness parameter for the asymptotic expression (7) to be valid is ε/h(z) for some monotonically increasing function h, meaning that the estimate is also valid for large ε as long as suitably large z → ∞ are considered.In this sense, the estimate is truly an extreme event probability estimate, but we chose to work in terms of the formal parameter ε to have an explicit and general scaling parameter, in contrast to the example-specific function h(z).For works on large deviation principles directly in z → ∞, see e.g.Dematteis et al. (2019), Tong et al. (2021) In addition to the probability estimate itself, the instanton, eigenvalues and eigenfunctions of η z also carry physical information about the system, as discussed in general in section 3. Figure 9 shows the instanton u z , i.e. the most likely field realization to reach a large wave height of z = 8.4, and the dominant space-time fluctuations δ u    z ) for the different resolutions (nx, nt ) as an approximation for the Fredholm determinant det(Id −Az), which is seen to converge with increasing resolution.Note that only for the lowest resolution, the eigenvalue spectrum shows noticeable deviations from the results at (nx, nt ) = (1024, 4000).The latter resolution has been used for all other numerical results on the KdV equation in this paper.and reads using the same pseudospectral code and explicit second order Runge-Kutta steps with integrating factor.The result for the prefactor agrees with the one obtained using the Fredholm determinant expression, with C F (z = 8.4) ≈ 1.0793 ⋅ 10 −2 using the eigenvalues and C F (z = 8.4) ≈ 1.0794 ⋅ 10 −2  42) and height observable (44), and dominant 5 normalized state variable eigenfunctions δ u (i) z of the projected second variation operator Az.Due to the KdV nonlinearity and linear wave dispersion, the large-scale forcing input is transformed into a large wave with dominant peak at t = T , x = 0 for the instanton uz, i.e. the most likely field realization to obtain a large wave height z = 8.4 at t = T , x = 0.The strongest fluctuations around the instanton resemble the instanton itself, but are necessarily centered around 0 with final-time height δ u(0, T ) = 0 at the origin.Note that only two eigenvalues are larger than 0.1 in modulus, reflecting the small effective dimension of the system in the noise variable, and that δ u (4) z contains already higher modes in time.
from the Riccati approach with For this particular observable value, the Riccati equation could be integrated without numerical problems, but we encountered a removable singularity for larger observable values.The final-time covariance of the conditioned Gaussian fluctuations around the instanton, as predicted using either the Riccati solution (38) or the eigenfunctions and eigenvalues (37), indeed coincides for both approaches and is highly oscillatory (top row, center and right in figure 10).Denoting the eigenvalues and normalized eigenfunctions of the finaltime covariance operator C z (T,T ) by ν with Z i independent and identically standard normally distributed.All in all, this example demonstrates the practical relevance and ease of applicability of the asymptotically sharp LDT estimate including the prefactor in a nonlinear, one-dimensional SPDE.

Stochastically forced incompressible three-dimensional Navier-Stokes equations
As a challenging, high-dimensional example, we consider the estimation of the probability of a high strain event in the stochastically forced incompressible three-dimensional Navier-Stokes equations.Our main goal here is to demonstrate that in addition to instantons for this problem, which were computed by Schorlepp et al. (2022), it is also numerically feasible to compute the leading order prefactor using the Fredholm determinant approach (14).Our setup hence follows the one treated by Schorlepp et al. (2022).A comprehensive analysis of the problem, including the behavior of the prefactor in the vicinity of the critical points of the dynamical phase transitions observed in this example, is beyond the scope of this paper.For other works on instantons and large deviations for the three-dimensional stochastic Navier-Stokes equations, see Apolinário et al. (2022), Falkovich et al. (1996), Grafke, Grauer and Schäfer (2015), Moriconi (2004).We consider a velocity field u ε ∶[0,l = 2π] 3 × [0,T = 1] → R 3 with periodic boundary conditions   in space that satisfies Here, P denotes the pressure which is determined through the divergence constraint.The forcing η is centered Gaussian, large-scale in space, white in time, and solenoidal with covariance where a Mexican hat correlation function with correlation length 1 is used.Note that this corresponds to the situation rankσ ≪ 3n 3 x of section 2.4, where only a small number of degrees of freedom is forced due to the Fourier transform χ decaying exponentially.As our observable, we consider the strain f (u) = ∂ 3 u 3 (x = 0) at the origin.Denoting the Leray projection onto the divergence-free part of a vector field by P, the instanton equations for (u z , p z ,λ z ) are given by ] .
With the instantons computed, we are able to evaluate the application of the second variation operator A z to noise fluctuation vectors δ η∶[0,2π] 3 ×[0,1] → R 3 by solving the second order adjoint equations We focus on z = −25 here, where the unique instanton solution does not break rotational symmetry (Schorlepp et al. 2022).Numerically, we use a pseudo-spectral GPU code with a spatial resolution n x = n y = n z = 128, a temporal resolution of n t = 512, a nonuniform grid in time with smaller time steps close to T = 1, and second order explicit Runge-Kutta steps with an integrating factor for the diffusion term.We truncated χ in Fourier space by setting it to 0 for all k where | χk | < 10 −14 , leading to ∥k∥ ≤ 9 and an effective real spatial dimension, independently of n x , of approximately rankσ ≈ 2 ⋅ (2 ⋅ 9) 3 = 11664 for the noise (by taking a cube instead of sphere for the Fourier coefficients of the noise vectors that are stored, and noting that χk projects onto k ⊥ ).The evaluation of the second order adjoint equations is then possible with only a few GB of VRAM for this resolution when exploiting double checkpointing and low rank storage as described in section 2.4.We computed the 600 largest eigenvalues of operator A z , again realized as a scipy.sparse.linalg.LinearOperator, by using scipy.sparse.linalg.eigsas before.We transfer the data to the GPU to evaluate the second variation applied to δ η by solving (54) with PyCUDA (Klöckner et al. 2012), and transfer back χ 1/2 * δ p to the CPU afterwards.Computing 600 eigenvalues this way needs about 1200 operator evaluations, or about 30 hours on a modern workstation with Intel Xeon Gold 6342 CPUs at 2.80 GHz and an NVIDIA A100 80GB GPU.The main limitation for computing more eigenvalues is that the eigenvalue solver used stores all matrix vector products in RAM.This could be overcome by storing some of them on a hard disk, or using different algorithms that can be parallelized over multiple nodes such as randomized SVD (Maulik and Mengaldo 2021).
The results for the eigenvalues of A z are shown in figure 11.We see that the absolute value of the eigenvalues decays such that the product ∏ z ) converges as m increases, but that even more than 600 eigenvalues would be needed for a more accurate result.For smaller observable values z, faster convergence is expected.Also, the spectrum of A z shows a large number of doubly degenerate eigenvalues, which appear whenever the eigenfunctions break the axial symmetry of the instanton.This feature of the spectrum clearly depends on the domain and spatial boundary conditions that were chosen here.From the instanton computation, we obtain I F (z) ≈ 1900.7 for the rate function, and from the 600 eigenvalues of A z we estimate C F (z) ≈ 4.9⋅10 −3 .With this, we can estimate that e.g. for ε = 250, the probability to observe a strain event with ∂ 3 u 3 (x = 0,T = 1) ≤ −25 is approximately 1.5 ⋅ 10 −5 , which matches the sampling estimate of P 250 F (−25) ∈ [1.3 ⋅ 10 −5 ,1.7 ⋅ 10 −5 ] at 95% asymptotic confidence, as obtained from 10 4 direct numerical simulations of (50) (data set from Schorlepp et al. (2022)).For smaller ε, the event becomes more rare, and it quickly becomes unfeasible to estimate its probability via direct sampling, whereas the quadratic estimate using the rate function and prefactor can be computed for any ε and is known to become more precise as the event becomes more difficult to observe in direct simulations.In addition to these probability estimates, we can also analyze the dominant Gaussian fluctuations around the instanton now and easily sample high strain events within the Gaussian approximation. Figure 12 shows the instanton u z at final time, i.e. an axially symmetric pair of counterrotating vortex rings, as well as the dominant eigenfunctions of C z (T,T ), corresponding to the fluctuation modes that are most easily observed at final time in conditioned direct numerical simulations.Note that the Riccati equation ( 27) would be a PDE for a six-dimensional matrixvalued field Q z (x 1 ,x 2 ,x 3 ,y 1 ,y 2 ,y 3 ,t) here without obvious z ) for different m as an approximation for the Fredholm determinant det(Id −Az).We see in the main figure that the eigenvalues often appear in pairs, which happens whenever the eigenfunctions break the rotational symmetry of the problem, such that, due to the periodic box, there are two linearly independent eigenfunctions for the same eigenvalue.The inset shows that det(Id −Az) is approximately 11 in this example, but that even more eigenvalues would be needed to get an accurate result.sparsity properties.Solvers for such a problem are quite expensive, if feasible at all, and also not easy to scale to higher spatial resolutions, whereas this is possible for the dominant eigenvalue approach.

SUMMARY AND OUTLOOK
In this paper, we have presented an asymptotically sharp, sampling-free probability estimation method for extreme events of stochastic processes described by additive-noise SDEs and SPDEs.The method can be regarded as a pathspace SORM approximation.We have introduced and compared two different conceptual and numerical strategies to evaluate the pre-exponential factor appearing in these estimates, either through dominant eigenvalues of the second variation, corresponding to the standard formulation of precise Laplace asymptotics and SORM, or through the solution of matrix Riccati differential equations, which is possible for precise large deviations of continuous-time Markov processes.Highlighting the scalability of the first approach, we have shown that leading-order prefactors can be computed in practice even for very high-dimensional SDEs, and explicitly tested our methods in two SPDE examples.In all examples, the approximations showed good agreement with direct Monte Carlo simulations or importance sampling.We hope that the methods assembled in this paper are useful whenever sample path large deviation theory is used to obtain probability estimates in real-world examples.
There are multiple possible extensions of the methods presented in this paper.More general classes of SDEs and SPDEs could possibly be treated numerically within the eigenvalue-based approach, most notably SDEs with multiplicative Gaussian noise, but also SDEs driven by Levy noise or singular SPDEs.Furthermore, one could try to generalize the approach to include any additive Gaussian noise that is colored in time instead of white.This would potentially lead to further dimensional reduction for the instanton and prefactor computation for examples with a slowly decaying temporal noise correlation.It would also be interesting to apply the eigenvalue-based prefactor computation strategy to metastable non-gradient SDEs.Regarding the numerical applicability of the Riccati method in case of high-dimensional systems with low-rank forcing, there is an alternative formulation of the prefactor in terms of a backward-in-time Riccati equation (Grafke et al. 2021), which could be better suited for controlled low-rank ap-proximations.In general, improvements of the quadratic approximation used throughout this paper via loop expansions, resummation techniques or non-perturbative methods from theoretical physics could be investigated.In this regard, it would be desirable to obtain simple criteria that indicate whether the SORM approximation considered in this paper can be expected to be accurate for given ε and z.Finally, one could use the instanton and additional prefactor information for efficient importance sampling of extreme events for S(P)DEs.
It is a common strategy in large deviation theory to first study expectations of the type E[exp{ 1 ε F(φ ε )}] for a family of random variables φ ε satisfying a large deviation principle, and a real-valued function F. Only later will these results be transformed onto probabilities or other probabilistic quantities.We directly use the results of Ben Arous (1988) to conclude that the asymptotic behavior of the moment-generating function (MGF) A ε F ∶R → [0,∞], A ε F (λ ) = E[exp{ λ ε f (X ε T )}] of the observable f (X ε T ) for the additive-noise SDE (9) as ε ↓ 0 is given by with prefactor Here, I * F denotes the Legendre transform of the rate function I F , det is a Fredholm determinant, the second variation operator is trace class, and η λ is short for η z λ at the Legendre dual z λ of λ via I ′ F (z λ ) = λ .Note that for multiplicative noise, the result would be different, which can already be seen in the simple example of a one-dimensional geometric Brownian motion and f (x) = 1 2 log 2 x.Furthermore, Ben Arous (1988) also assumes that the vector field b, the observable f , and their respective derivatives are bounded.A remark by Deuschel et al. (2014) shows how one could relax this assumption via localization.
Evaluating the inverse Laplace transform from the MGF (A6) to the probability density function via a saddlepoint approximation, as well as a further integration to get the tail probability via a Laplace approximation yields the desired estimate with leading-order prefactor From the first-order necessary condition and λ z = I ′ F (z), we get via differentiation as claimed.The last equality is easy to see for finitedimensional matrices: For A ∈ R N×N invertible and a unit vector e ∈ R N , the adjugate is adj(A) = detA ⋅ A −1 , and applying e from the left and right yields detA⟨e,A −1 e⟩ N = ⟨e,adj(A)e⟩ N .The right-hand side is the (e,e) cofactor of A, which is equal to the determinant of the (N − 1) × (N − 1) matrix pr e ⊥ Apr e ⊥ with pr e ⊥ ∶R N → e ⊥ denoting the orthogonal projection.For the present infinite-dimensional case, an analogue of this relation can be verified using the series definition of the Fredholm determinant and adjugate as originally introduced by Fredholm himself (Fredholm 1903, McKean 2011).
3. From Fredholm determinants to zeta-regularized functional determinants In this section, we motivate (26) using purely formal manipulations, and only consider linear observables f for simplicity.For rigorous results on the relation between Fredholm determinants and zeta-regularized determinants for related classes of operators, see e.g.Forman (1987), Hartmann andLesch (2022).We start with the expression (A9) for the prefactor C F (z) in terms of the full second variation determinant without projection operators.According to the adjoint formulation of section 2.2, we write the second variation δ 2 (λ F)/δ η 2 as the composition of three linear operators Here, the critical step is in the second line where the operators are moved out of the Fredholm determinant to get a fraction of two zeta-regularized determinants, which is true for finite-dimensional matrices but non-trivial for general operators.We see that the boundary conditions of all appearing operators are which is the correct special case of the general boundary conditions FIG. 2. Visualization of five different sample paths (light orange)and the mean of 100 such paths (orange with black outline) of the model SDE (10) that satisfy f (X(T ),Y (T )) ≥ z with z = 3 (red set) and ε = 0.5.Using Euler-Maruyama steps with an integrating factor with step size ∆t = 5 ⋅ 10 −4 , we repeatedly simulated (10) until 100 such rare trajectories were found.The dashed blue line is the state variable instanton trajectory φz, solution of (19) with the optimal ηz as forcing.As in Figure1, the gray lines are field lines of the drift vector field b.
FIG. 3. Result of numerically computing 200 eigenvalues µ (i) z FIG.5.Results of numerically computing 10 5 transition paths from x = 0 to the target set f −1 ({z}) for the model SDE (10) with z = 3 and ε = 0.5 using instanton-based importance sampling(Ebener et al. 2019).We visualize the transition tube information obtained from the eigenvalues and eigenfunctions of the projected second variation operator.The upper left subfigure shows the histogram of the full data set for all times.The remaining subfigures show histograms of the transition paths at specific times t.The black lines, as a comparison, are the level sets of the normal PDF with covariance εCz(t,t), found by evaluating (37) numerically, and mean φz(t).Note that the deformation of the distribution of X ε t , conditioned on f (X ε T ) = z, is captured quite well using the quadratic, sampling-free approximation.
FIG.6.Left column: Rate function I F (top) and leading order prefactor C F (bottom) for the KdV equation (42) with height observable (44), as obtained from numerical instanton and prefactor computations.Note that the prefactor depends strongly, almost exponentially, on the observable value z in this example.Right: Comparison of LDT estimate (7) for different noise levels ε ∈ {0.1, 1, 10} to direct sampling for the SPDE (42).For each ε, we computed 4 ⋅ 10 4 samples of f (u ε (⋅, T )) to estimate the tail probabilities for various z.The shaded regions are 95% Wilson score intervals(Brown et al. 2001) for the sampling estimate of the tail probabilities.The solid lines show the asymptotically sharp estimate (7) without adjustable parameters.In comparison to this, the dashed lines show just the leading order LDT term exp {−I F (z)/ε} with a constant prefactor (chosen such that the curve matches (7) for large z), which shows that the prefactor C F is absolutely necessary to get useful results in this example at ε > 0.1 and can be understood in regard to the left column of the figure.Results use nx = 1024, nt = 4000 for the instanton computations, pseudo-spectral code with integrating factor, L-BFGS optimization with penalty term for observable; 80 eigenvalues with largest absolute value for Fredholm determinant; stochastic Heun steps with size ∆t = 10 −3 for direct sampling.

z
around it.We further computed the Gaussian fluctuations around the instanton for z = 8.4 at the final instance t = T in figure 10.First of all, we also solved the forward Riccati equation (27), which is a PDE for FIG. 8. Performance of instanton and prefactor computations forKdV problem with z = 8.4 for different spatio-temporal resolutions (nx, nt ) ∈ {(32, 125), .. ., (1024, 4000)}.The table shows the number of optimization iterations required to compute the instanton, and the value of the objective I F (z).The number of iterations does not increase with the resolution (nx, nt ).The bottom figure shows 80 eigenvalues µ (i) z with largest absolute value of Az.The main figure shows the absolute value of the eigenvalues µ FIG. 9. Example instanton field uz in space and time for z = 8.4 (top left) for the KdV equation (42) and height observable (44), and dominant 5 normalized state variable eigenfunctions δ u in the bottom left of figure 10 quickly decay.Using the eigenvalues and eigenfunctions, realizations of u ε (⋅,T ) when conditioning on u ε (0,T ) = z = 8.4 can now easily be sampled within the Gaussian approximation as FIG.10.Information on the conditioned final time Gaussian fluctuations around the KdV instanton for z = 8.4, calculated from the quantities used to evaluate the prefactor C F (z). Top, left: Riccati solution Qz(⋅, ⋅, T = 1) at final time.Top, center: Projection of the Riccati solution, such that the constraint δ u(0, T ) = 0 is satisfied.This way, the final time covariance Cz(T, T ) as given in (38) is obtained.Top, right: The same final time covariance Cz(T, T ) constructed from the eigenvalues and eigenfunctions of Az instead as in (37).The result is visually indistinguishable from the Riccati computations.Bottom, left: Eigenvalues ν (i) z (T ) of the covariance Cz(T, T ).We see that the eigenvalues quickly decay to zero, and less than 10 fluctuation modes are in fact relevant.Bottom, center: Eigenfunctions δ v (i) z for the 4 dominant eigenvalues ν (i) z (T ), i ∈ {1, 2, 3, 4}, which all necessarily satisfy δ v (i) z (x = 0) = 0. Bottom, right: Instanton uz(⋅, T ) at final time (dashed line), and variance of conditioned Gaussian fluctuations around it for ε = 0.1 (shaded area).
FIG. 11.Result of numerically computing 600 eigenvalues µ (i) z with largest absolute value of Az for the three-dimensional Navier-Stokes equations (50) with strain z = ∂ 3 u 3 (x = 0, T ) = −25, where the instanton is a rotationally symmetric pair of vortex rings.Main figure: absolute value of the eigenvalues µ