Scalability of the second-order reliability method for stochastic differential equations with multiplicative noise

TL;DR

This work generalizes the second-order reliability method to stochastic differential equations with multiplicative noise by deriving a renormalized, infinite-dimensional prefactor involving a Carleman–Fredholm determinant and a trace term. The authors establish a scalable, matrix-free computational framework enabled by forward-mode differentiation and operator-based determinants, enabling tail-probability estimates in high-dimensional SDEs and SPDEs. They validate the approach on a low-dimensional predator–prey system and a high-dimensional advection–diffusion SPDE, showing close agreement with Monte Carlo tails and robust convergence under increasing temporal resolution. The method, implemented in a public JAX-based codebase, opens practical avenues for accurate rare-event analysis in physics, engineering, and finance where multiplicative noise plays a central role.

Abstract

We show how to efficiently compute asymptotically sharp estimates of extreme event probabilities in stochastic differential equations (SDEs) with small multiplicative Brownian noise. The underlying approximation is known as sharp large deviation theory or precise Laplace asymptotics in mathematics, the second-order reliability method (SORM) in reliability engineering, and the instanton or optimal fluctuation method with 1-loop corrections in physics. It is based on approximating the tail probability in question with the most probable realization of the stochastic process, and local perturbations around this realization. We first recall and contextualize the relevant classical theoretical result on precise Laplace asymptotics of diffusion processes [Ben Arous (1988), Stochastics, 25(3), 125-153], and then show how to compute the involved infinite-dimensional quantities - operator traces and Carleman-Fredholm determinants - numerically in a way that is scalable with respect to the time discretization and remains feasible in high spatial dimensions. Using tools from automatic differentiation, we achieve a straightforward black-box numerical computation of the SORM estimates in JAX. The method is illustrated in examples of SDEs and stochastic partial differential equations, including a two-dimensional random advection-diffusion model of a passive scalar. We thereby demonstrate that it is possible to obtain efficient and accurate SORM estimates for very high-dimensional problems, as long as the infinite-dimensional structure of the problem is correctly taken into account. Our JAX implementation of the method is made publicly available.

Figures (9)

• Figure 1: A numerical example -- the probability estimation of an extreme prey concentration event $f(X_T^\varepsilon, Y_T^\varepsilon) = X_T^\varepsilon \geq z = 0.5$ in the stochastic Lotka--Volterra model with multiplicative noise \ref{['eq:predator-prey-SDE']} -- demonstrating that "naively" trying to compute the prefactor $C(z)$ in \ref{['eq:p-eps-asymp-form']} using a discretized version of \ref{['eq:tail-prob-prefac-sde-additive']} fails to be scalable with respect to the time discretization dimension $n_t$. Top left: state space instanton trajectories $\phi_z$ for different $n_t$ (dashed) found from numerical optimization, and a few sample paths of the SDE \ref{['eq:predator-prey-SDE']} at noise strength $\varepsilon = 0.01$ that reach the extreme event set at final time $T = 10$. The grey field lines visualize the drift vector field of the SDE \ref{['eq:predator-prey-SDE']}. Top right: table that summarizes the numerical results for the observable rate function $I(z)$, Lagrange multiplier $\lambda_z$, number of L-BFGS iterations in the instanton optimization, and prefactor $C(z)$ from \ref{['eq:prefac-wrong-discretize']} using all eigenvalues of the discrete projected Hessian $\text{pr}_{\eta_z^\perp} A_{\lambda_z} \text{pr}_{\eta_z^\perp}$. Bottom left: Full spectrum of all eigenvalues $\mu_z^{(i)}$ of the discretized projected Hessian at different resolutions. Positive eigenvalues are drawn as dots, negative eigenvalues as crosses. Bottom right: Prefactor $C(z)$ when approximating \ref{['eq:prefac-wrong-discretize']} using only the $m$ leading eigenvalues. No discretization-independent convergence is observed, and the number of necessary eigenvalues increases with $n_t$. The dashed lines, for references, are obtained by discretizing the correct continuum limit from theorem \ref{['thm:prefac']} instead of \ref{['eq:tail-prob-prefac-sde-additive']}, as explained later in section \ref{['sec:num-methods']}.
• Figure 2: Geometric Brownian motion toy example, where we want to estimate the MGF $J^\varepsilon(\lambda)$ of geometric Brownian motion \ref{['eq:geom-bm']} at $\lambda = -1$, and see that again the prefactor calculation according to the MGF equivalent of \ref{['eq:prefac-wrong-discretize']} is not scalable here. The other parameters are $\beta = T = 1$, such that the second variation in continuous time $A_\lambda$ has exactly one nonzero eigenvalue $2 \lambda T = -2$. Left: Eigenvalue spectrum of the discrete Hessians at different resolutions $n_t$, with the crosses indicating negative eigenvalues. Importantly, in addition to the predicted eigenvalue $\approx -2$, all other eigenvalues of the discrete Hessian are small but nonzero here. Right: Approximation of $R_\lambda$ through $m$ leading eigenvalues as $\left[\prod_{i = 1}^m \left(1 - \mu_\lambda^{(i)} \right) \right]^{-1/2}$. We see that this is not scalable, since all of the small but nonzero eigenvalues of the discrete Hessian are needed to "make up for the missing exponential term" compared to the exact result \ref{['eq:mgf-geom-bm-ex']} for the MGF prefactor $R_\lambda$. Dashed lines plotted according to the numerical evaluation of the correct prefactor expression \ref{['eq:mgf-prefac']} for the MGF, which immediately gives the correct result after one eigenvalue.
• Figure 3: Flow chart that visualizes how the second variation $A_\lambda$ of the noise-to-event map $F$ acts on $\delta \eta \in L^2\left([0,T], \mathbb{R}^n \right)$ for additive noise according to \ref{['eq:second-var-additive']} and \ref{['eq:second-order-adj-eq-additive']} and multiplicative noise as in \ref{['eq:second-var-mult']} and \ref{['eq:second-order-adjoint-mult-noise']}. Here, each integral symbol $\int \mathrm{d} t$ above an arrow corresponds to solving an ODE with input given by the previous box, and arrows without integral symbols correspond to directly plugging the entry of the previous box into an equation to calculate the next box. The dashed red arrows are only present for multiplicative noise.
• Figure 4: Rate function $I$ and prefactor $C$ (left column) as well as probability $P^\varepsilon(z)$ (right) of observing a large concentration of prey $f(x(T), y()T)) = x(T)$ in the Lotka--Volterra model (\ref{['eq:predator-prey-SDE']}), exceeding $z$ at time $t=T$ when starting at the fixed point $(x_0, y_0)$ of the system. The right subfigure shows a comparison between the sharp large deviation estimate \ref{['eq:tail-asymp']} from theorem \ref{['thm:prefac']} (solid lines) with prefactor according to \ref{['eq:tail-prefac-projected']} against Monte Carlo estimates (data points, with $99\%$ Wilson score intervals as shaded area) with $9 \cdot 10^8$ samples for each noise strength $\varepsilon \in \{0.001, 0.004, 0.01\}$. Temporal resolution $n_t = 1000$ for sampling and all asymptotic estimates, and all other parameters are as described in section \ref{['sec:predprey']}. The dashed lines are $\text{const} \cdot \exp \left\{-I(z) / \varepsilon \right\}$ with $\text{const}$ chosen in such a way that the curves match at large $z$. This is the best we can do without knowledge of the prefactor, and is not an a priori estimate without Monte Carlo data. The result is not too different from the full prediction including the prefactor $C(z)$ in this particular example, since $C(z)$ does not depend strongly on $z$ for large $z$ in this example.
• Figure 5: Numerical results for the leading $M=200$ eigenvalues of the spectra of the operators that are needed to evaluate the tail probability estimate \ref{['eq:tail-asymp']} from theorem \ref{['thm:prefac']} for the predator-prey model \ref{['eq:predator-prey-SDE']} for an observable value of $f(x,y) = z = 1.0$ at different time resolutions $n_t$. In all spectra, dots correspond to positive eigenvalues and crosses to negative eigenvalues. Left: Spectrum of the Hessian $A_{\lambda_z}$ without projection operators. Note that $I"(z = 1) < 0$ in this example, as evident from figure \ref{['fig:pred-prey-prob']}, and there is hence indeed a single eigenvalue $\mu_z^{(2)} > 1$ in the spectrum, and $\det_2 \left(\text{Id} - A_{\lambda_z}\right) < 0$. Center: Spectrum of the projected Hessian $\text{pr}_{\eta_z^\perp}A_{\lambda_z} \text{pr}_{\eta_z^\perp}$, as it actually enters into the prefactor formula \ref{['eq:tail-prefac-projected']}. The projection operators change the spectrum and in particular remove the eigenvalue $>1$, and hence the CF determinant shown in the inset is positive now. Right: Subtracting the HS part $\tilde{A}_{\lambda_z}$ from the full Hessian leads to a trace-class operator $\text{pr}_{\eta_z^\perp} \left( A_{\lambda_z} - \tilde{A}_{\lambda_z} \right) \text{pr}_{\eta_z^\perp}$ with eigenvalue decay $\propto i^{-2}$ instead of $\propto i^{-1}$, and the trace is indeed seen to converge in the inset. See table \ref{['tab:pred-prey-z-1.0']} for further details.

Theorems & Definitions (23)

• Theorem 2.1
• Remark 2.1
• Theorem D.1
• Definition D.1
• Definition D.2
• Definition D.3
• Definition D.4
• Example D.1
• Remark D.1
• Theorem D.2