Large deviation theory

This section applies results of large deviation theory (LDT), particularly action minimization problems for computing most probable transition paths in stochastic dynamical systems driven by weak noise. For a description of the theory, see Freidlin and Wentzell (1998) and (Börner, 2025). An overview of numerical methods applying LDT is given in Grafke and Vanden-Eijnden (2019).

Info

The methods in this section apply to $D$-dimensional stochastic dynamical systems of the form

\[\text{d} \mathbf{x} = \mathbf{b} (\mathbf{x}) \text{d}t + \sigma \mathbf{\Sigma} \text{d}\mathbf{W}_t \,,\]

where the drift field $\mathbf{b}$ may be non-gradient but the noise term must consist of Gaussian noise ($\mathbf{W}_t$ is a $D$-dimensional vector of independent standard Wiener processes) and a constant covariance matrix $\mathbf{Q} = \mathbf{\Sigma}\mathbf{\Sigma}^\top$.

This is a special case of the broader class of noise types supported by CoupledSDEs.

Noise convention used by the action functionals

Starting point. A CoupledSDEs exposes its noise via the diffusion function $g$ and, for additive invertible noise, the noise rate covariance

\[\mathbf{\Sigma} \;=\; \texttt{covariance\_matrix(sys)} \;=\; g\,g^\top.\]

For an SDE built as $\mathrm{d}\mathbf{x} = \mathbf{b}\,\mathrm{d}t + \sigma\,\mathbf{\Sigma}_0\,\mathrm{d}\mathbf{W}_t$ with $\mathbf{\Sigma}_0\mathbf{\Sigma}_0^\top = \mathbf{Q}$, this gives $\mathbf{\Sigma} = \sigma^2\,\mathbf{Q}$. The Freidlin–Wentzell rate function only depends on the product $\sigma^2\mathbf{Q}$; the individual $(\sigma, \mathbf{Q})$ split is not determined by the SDE alone ($(\sigma, \mathbf{Q}) \leftrightarrow (c\sigma, \mathbf{Q}/c^2)$ is the same physical system).

Convention. The package picks a unique $(\sigma, \mathbf{Q})$ pair from $\mathbf{\Sigma}$ by adopting the trace convention $\mathrm{tr}(\mathbf{Q}) = D$:

\[\sigma_{\mathrm{eff}}^2 \;=\; \frac{\mathrm{tr}(\mathbf{\Sigma})}{D}, \qquad \mathbf{Q}_{\mathrm{can}} \;=\; \frac{D}{\mathrm{tr}(\mathbf{\Sigma})}\,\mathbf{\Sigma}.\]

The accessor noise_strength(sys) returns $\sigma_{\mathrm{eff}}$. The trace is invariant under orthogonal changes of basis and reduces to the per-direction average noise variance, making this the natural choice among rotation-invariant normalizations.

Action functionals use the canonical covariance. All action functionals in this section evaluate the FW integrand in the $\mathbf{Q}_{\mathrm{can}}^{-1}$ metric:

\[\texttt{fw\_action(sys, path, time)} \;=\; \Phi_{FW}^{\mathbf{Q}_{\mathrm{can}}}[\phi] \;=\; \tfrac{1}{2}\!\int (\dot\phi - \mathbf{b})^\top \mathbf{Q}_{\mathrm{can}}^{-1}(\dot\phi - \mathbf{b})\,\mathrm{d}t.\]

Two consequences:

The returned action is independent of the noise_strength keyword chosen at construction (both $\sigma$ and the magnitude of $\mathbf{Q}$ are absorbed into $\sigma_{\mathrm{eff}}$).
The returned action is invariant under orthogonal changes of basis, matching the coordinate-independence of the FW action as a geometric quantity.

Converting to the user's parameterization. If you built your system as $\mathrm{d}\mathbf{x} = \mathbf{b}\,\mathrm{d}t + \sigma_{\mathrm{user}}\,\mathbf{\Sigma}_0\,\mathrm{d}\mathbf{W}_t$ with $\mathbf{Q}_{\mathrm{user}} = \mathbf{\Sigma}_0\mathbf{\Sigma}_0^\top$ and want the FW action in your metric $\mathbf{Q}_{\mathrm{user}}^{-1}$, multiply the returned action by $\mathrm{tr}(\mathbf{Q}_{\mathrm{user}})/D$:

\[\Phi_{FW}^{\mathbf{Q}_{\mathrm{user}}}[\phi] \;=\; \frac{\mathrm{tr}(\mathbf{Q}_{\mathrm{user}})}{D}\cdot \texttt{fw\_action(sys, path, time)}.\]

The factor is $1$ whenever $\mathbf{Q}_{\mathrm{user}}$ is isotropic ($c\,\mathbf{I}$ for any $c>0$) or already trace-normalized ($\mathrm{tr}(\mathbf{Q}_{\mathrm{user}}) = D$), in which case fw_action returns $\Phi_{FW}^{\mathbf{Q}_{\mathrm{user}}}$ directly. All examples and tests in the package fall in this case.

fw_action vs om_action. fw_action and geometric_action are fully $\sigma$-independent under this convention. om_action takes noise_strength as an explicit positional argument because the Onsager–Machlup correction $(\sigma^2/2)\int \nabla\cdot \mathbf{b}\,\mathrm{d}t$ is a finite-noise correction that scales with $\sigma^2$; only the FW part of OM is $\sigma$-independent. Pass the $\sigma$ at which you want $-\log P[\phi]$ evaluated (typically the same noise_strength you used at construction). As $\sigma \to 0$ the OM correction vanishes and OM $\to$ FW, recovering the leading-order LDP rate function.

Action minimizers

Several methods have been proposed to calculate transition paths that minimize a given action functional. In the weak-noise limit, this minimum action path (or instanton) corresponds to the most probable transition path. While the minimum action method (MAM) is the most basic version, it is often beneficial to minimize the geometric action via a time-independent version called gMAM. The problem can also be cast in a Hamiltonian form, implemented as simple gMAM (sgMAM), which can have numerical advantages.

These methods apply to non-gradient systems driven by Gaussian noise. In gradient systems, minimum action paths between attractors coincide with heteroclinic orbits, which can be computed via the so-called string method.

Action minimization as an optimal control problem

All of the action minimizers listed below (MAM, gMAM, sgMAM) can equivalently be cast as optimal control problems: find a path $\mathbf{x}$ and a control $\mathbf{u}$ that minimize a Freidlin-Wentzell-type action

\[\int \| \mathbf{u}(t) - \mathbf{b}(\mathbf{x}(t)) \|^2 \, \text{d}t\]

subject to the path dynamics $\dot{\mathbf{x}}(t) = \mathbf{u}(t)$ and the boundary conditions $\mathbf{x}(0) = \mathbf{x}_A$, $\mathbf{x}(T) = \mathbf{x}_B$ (with $T$ either fixed, as in MAM, or eliminated by reparameterization, as in gMAM/sgMAM). In this form the problem can be solved with dedicated optimal-control packages such as OptimalControl.jl; see the control-toolbox.org MAM tutorial for a worked example on the Maier-Stein system.

To summarize, the following methods are available:

Minimum action method (MAM)
Geometric minimum action method (gMAM)
Simple geometric minimum action method (sgMAM)
String method

All three action-minimizers (MAM, gMAM, sgMAM) require only autonomous noise and reject rank-deficient diffusion; they all support additive, diagonal-multiplicative, and general-matrix multiplicative diffusion. The Onsager-Machlup functional (selectable in MAM via functional = "OM") is the one exception: it is implemented only for additive noise and throws otherwise. Pick by what you want out of the result:

Method	Use when
MAM	The transition time $T$ is fixed and you want the action at that $T$ (or you specifically need the Onsager-Machlup functional, which only has a time-parameterized form).
gMAM / sgMAM	$T$ is unknown and you want a time-reparameterization-invariant instanton. Prefer sgMAM (Hamiltonian picture) when an analytic `FreidlinWentzellHamiltonian(H_x, H_p)` is available or when AD-based `jacobian(ds)` evaluations are cheap; prefer gMAM for the closer-to-textbook geometric-action formulation.
String method	You want the deterministic heteroclinic orbit (typical use: gradient systems, where it coincides with the instanton). In non-gradient systems it generally differs from the most probable transition path.

Variants and extensions

The literature contains a number of extensions of MAM-type methods that may be relevant depending on the model class and numerical difficulties. These variants are not currently implemented in CriticalTransitions.jl, but serve as useful pointers:

tMAM / optimal linear time scaling: avoids explicit optimization over the transition time by introducing an optimal linear time scaling; can be combined with adaptivity in time discretization Wan (2015).
Adaptive MAM: uses a moving-mesh strategy to concentrate grid points in dynamically important portions of the path, improving efficiency and robustness Zhou et al. (2008).
Non-Gaussian (jump / Lévy) noise: for systems driven by jump noise, the rate function and path optimization problem differ from the Freidlin–Wentzell diffusive setting; see e.g. an optimal-control-based approach in Wei et al. (2023).
Multiplicative / state-dependent noise: state-dependent diagonal and full-matrix multiplicative noise is supported by both gMAM (via minimize_geometric_action(::CoupledSDEs, ...)) and sgMAM (via FreidlinWentzellHamiltonian(::CoupledSDEs)). The diffusion tensor a(x) is classified once when the sgMAM cache is built (constant vs state-dependent, diagonal vs coupled) and the resulting inner loop dispatches at compile time on the cache type; see Grafke et al. (2017) for the underlying Hamiltonian formulation. Rank-deficient (degenerate) noise is rejected at cache build; see Grafke et al. (2017) for the rank-deficient extension that would be needed to support it.

Minimum action method (MAM)

Minimization of the specified action functional using the optimization algorithm of Optimization.jl. See also E et al. (2004).

CriticalTransitions.minimize_action — Function

minimize_action(sys::ContinuousTimeDynamicalSystem, x_i, x_f, T::Real, optimizer=Optimisers.Adam(); kwargs...)

Minimizes an action functional to obtain a minimum action path (instanton) between an initial state x_i and final state x_f in phase space.

This algorithm uses the Optimization.jl package to minimize the specified action functional (either fw_action or om_action) for the system sys over paths connecting x_i to x_f in time T.

The path is initialized as a straight line between x_i and x_f, parameterized in time via N equidistant points and total time T. Thus, the time step between discretized path points is $\Delta t = T/N$. To set an initial path different from a straight line, see the multiple dispatch method

minimize_action(sys::ContinuousTimeDynamicalSystem, init::AbstractMatrix, T::Real, optimizer=Optimisers.Adam(); kwargs...).

Returns a MinimumActionPath object containing the optimized path and the action value.

The optional positional argument optimizer selects the Optimization.jl solver. It defaults to Optimisers.Adam().

Keyword arguments

functional = "FW": type of action functional to minimize. Defaults to fw_action, alternative: "OM" for om_action
npoints = 100: number of path points to use for the discretization of the path
noise_strength = nothing: noise strength for the action functional. Specify only if functional = "OM"
ad_type = OptimizationBase.AutoFiniteDiff(): type of automatic differentiation to use (see Optimization.jl)
maxiters = 100: maximum number of iterations before the algorithm stops
abstol::Real=NaN: absolute tolerance of action gradient to determine convergence
reltol::Real=NaN: relative tolerance of action gradient to determine convergence
verbose = false: whether to print Optimization information during the run
show_progress = false: whether to print a progress bar

source

minimize_action(sys::ContinuousTimeDynamicalSystem, init::AbstractMatrix, T::Real, optimizer=Optimisers.Adam(); kwargs...)

Minimizes the specified action functional to obtain a minimum action path (instanton) between fixed end points given a system sys and total path time T.

The initial path init must be a matrix of size (D, N), where D is the dimension of the system and N is the number of path points. The physical time of the path is specified by T, such that the time step between consecutive path points is $\Delta t = T/N$.

The optional positional argument optimizer selects the Optimization.jl solver. It defaults to Optimisers.Adam().

source

Geometric minimum action method (gMAM)

Minimization of the geometric action following Heymann and Vanden-Eijnden (2008), Heymann and Vanden-Eijnden (2008). gMAM reformulates MAM to avoid double optimization of both the action and the transition time. It achieves this by using a geometric action functional that is independent of the time parametrization of the path. This reparameterization invariance makes the method more robust and computationally efficient, particularly for multiscale systems.

CriticalTransitions.minimize_geometric_action — Function

minimize_geometric_action(sys::FreidlinWentzellHamiltonian, x_initial, optimizer = GeometricGradient(; stepsize = 1e3); kwargs...)

Runs the simplified geometric Minimum Action Method (sgMAM, (Grafke et al., 2017)) on the FreidlinWentzellHamiltonian sys, returning the instanton (minimizer of the Freidlin-Wentzell action) that connects the two endpoints of x_initial.

sgMAM works directly in the doubled $(x, p)$ phase space and minimizes the symplectic action $\int p \cdot \mathrm{d}x$ under the zero-energy constraint $H(x, p) = 0$. Here $p = a(x)^{-1}(\dot x - b(x))$ is the conjugate momentum, i.e. the (dimensionless) noise force that the Wiener process must supply to drive the trajectory over the barrier.

Arguments

sys::FreidlinWentzellHamiltonian: the Hamiltonian, typically obtained from a CoupledSDEs via FreidlinWentzellHamiltonian(ds).
x_initial::AbstractMatrix{T} (or StateSpaceSet): initial guess for the instanton, shape D × Nt with state points in columns. Endpoints x_initial[:, 1] and x_initial[:, end] are held fixed; the interior is reparameterized to uniform arclength on the first iteration.
optimizer: step-size control, either
- GeometricGradient (default; backtracking projected gradient), or
- AdaptiveGeometricGradient (multi-phase probe variant; more robust on underdamped systems at higher per-iteration cost).

Keyword arguments

maxiters::Int = 1000: maximum outer iterations.
abstol::Real = NaN, reltol::Real = NaN: convergence on the absolute / relative action change between accepted iterations. NaN disables that criterion.
verbose::Bool = false: log convergence and backtracking events.
show_progress::Bool = false: show a ProgressMeter bar.

Returns a MinimumActionPath

source

minimize_geometric_action(sys::CoupledSDEs, x_i, x_f, optimizer=GeometricGradient(); kwargs...)
minimize_geometric_action(sys::CoupledSDEs, init::AbstractMatrix, optimizer=GeometricGradient(); kwargs...)

Computes the minimizer of the geometric Freidlin-Wentzell action based on the geometric minimum action method (gMAM), using optimizers of Optimization.jl or the original formulation by Heymann and Vanden-Eijnden (2008) (projected gradient descent with optional backtracking).

The minimizer is computed for system sys over all paths from x_i to x_f. To set an initial path different from a straight line, see the multiple-dispatch method minimize_geometric_action(sys::CoupledSDEs, init::AbstractMatrix, optimizer; kwargs...).

Returns a MinimumActionPath.

Keyword arguments

npoints::Int = 100: number of discretization points (endpoint form only)
maxiters::Int = 100: maximum optimizer iterations
abstol::Real = NaN, reltol::Real = NaN: action-change convergence criteria
ad_type = OptimizationBase.AutoFiniteDiff() (Optimization.jl path only)
verbose::Bool = false, show_progress::Bool = true

source

Simple geometric minimum action method (sgMAM)

Simplified minimization of the geometric action following Grafke et al. (2017). The simple gMAM reduces the complexity of the original gMAM by requiring only first-order derivatives of the underlying Hamiltonian optimization formulation. This simplifies the numerical treatment and computational complexity.

The implementation below performs a constrained gradient descent on the Hamiltonian system, supporting autonomous diffusions with additive, diagonal-multiplicative, or general-matrix multiplicative noise.

Hamiltonian picture

The Freidlin-Wentzell rate function admits an equivalent Hamiltonian formulation with conjugate momentum $p$:

\[H(\varphi, p) = \langle b(\varphi), p \rangle + \tfrac{1}{2}\, \langle p,\, a(\varphi)\, p \rangle,\]

where $a(x) = \sigma(x)\sigma(x)^{\top}$. The action along an instanton (where $H \equiv 0$) reduces to $S[\varphi] = \int_0^T \langle p, \mathrm{d}\varphi\rangle$, the form used by minimize_geometric_action when the input is a FreidlinWentzellHamiltonian.

Compile-time dispatch on the shape of $a(x)$ is driven by the type of sys.a (a Base.Returns wrapper marks constant-in-x diffusion) together with the type of a(x_ref) (a LinearAlgebra.Diagonal marks diagonal coupling). Classification happens once when the per-path cache is built; the resulting cache type then selects the diagonal or coupled inner loops in update_p!, update_x!, and geometric_gradient_step!.

Rank-deficient $a(x)$ is rejected at cache build (a is probed at a reference state x_ref (the first path point) and at the $2D$ neighbors x_ref ± h·eₗ).

CriticalTransitions.FreidlinWentzellHamiltonian — Type

FreidlinWentzellHamiltonian{IIP, D, ...}

Freidlin-Wentzell Hamiltonian for the small-noise SDE $\mathrm{d}X_t = b(X_t)\,\mathrm{d}t + \sigma\,\Sigma(X_t)\,\mathrm{d}W_t$, with diffusion tensor $a(x) = \Sigma(x)\Sigma(x)^\top$:

\[H(x, p) \;=\; \langle b(x),\, p\rangle \;+\; \tfrac{1}{2}\,\langle p,\, a(x)\,p\rangle.\]

p is the conjugate momentum (Legendre dual of $\dot x$); Hamilton's equations are $\dot x = b(x) + a(x)\,p$ and $\dot p = -\partial_x b^\top p - \tfrac{1}{2}\langle p, \partial_x a\,p\rangle$. Freidlin-Wentzell instantons (action minimizers between invariant sets) live on the zero-energy shell $H \equiv 0$; the simplified geometric MAM ((Grafke et al., 2017)) minimizes the action directly in (x, p) space.

The stored a(x) is trace-normalized (see _trace_normalized_a) so that the action is invariant to the overall scale of noise_strength.

Fields

H_x, H_p: callables (x, p) -> ∂_x H, (x, p) -> ∂_p H, returning D × N matrices.
a: trace-normalized diffusion callable; a Base.Returns for constant noise.
x_ref: reference state used by Base.show; nothing if not supplied.

Constructors

FreidlinWentzellHamiltonian(ds::ContinuousTimeDynamicalSystem): builds H_x, H_p, a from ds. Central finite differences are used for $\partial_x a$ when a is state-dependent. Validation and diagonal-vs-coupled classification happen at cache build (build_sgmam_cache); rank-deficient a is rejected there.
FreidlinWentzellHamiltonian{IIP, D}(H_x, H_p; a = Returns(Diagonal(ones(D))), x_ref = nothing): for hand-rolled Hamiltonians; you are responsible for matching the convention above.

See (Freidlin and Wentzell, 1998) for the underlying theory.

source

Performance notes (sgMAM)

sgMAM repeatedly evaluates H_p(x, p) and H_x(x, p) along a discretized path. If this allocates, sgMAM will be slow.

Key reason for performance differences:

FreidlinWentzellHamiltonian(ds::CoupledSDEs) typically relies on jacobian(ds) (often automatic differentiation unless you provide an analytic Jacobian) and evaluates it pointwise along the path.
A hardcoded FreidlinWentzellHamiltonian(H_x, H_p) with analytic expressions operating on the full D×Nt path matrix usually allocates far less.

Benchmark pattern:

using CriticalTransitions
using BenchmarkTools

sys_fast = FreidlinWentzellHamiltonian{false,2}(H_x, H_p)  # hardcoded analytic H_x/H_p

ds = CoupledSDEs(KPO, zeros(2), ())
sys_generic = FreidlinWentzellHamiltonian(ds)              # uses jacobian(ds)

opt = GeometricGradient(; stepsize=0.5)
@btime minimize_geometric_action($sys_fast,    $x_initial, $opt; maxiters=100, show_progress=false)
@btime minimize_geometric_action($sys_generic, $x_initial, $opt; maxiters=100, show_progress=false)

Aside: the same “vectorized + allocation-free inner loop” principle also tends to make string_method faster when used with FreidlinWentzellHamiltonian.

Multiple shooting

The MultipleShooting optimizer treats the instanton as a boundary value problem on arclength-reparameterized Hamilton equations

\[\frac{\mathrm{d}\varphi}{\mathrm{d}s} = \alpha\,H_p(\varphi, p),\qquad \frac{\mathrm{d}p}{\mathrm{d}s} = -\alpha\,H_x(\varphi, p),\qquad \alpha = L\,/\,\|H_p\|\]

on s ∈ [0, 1], with the total path length L a Newton unknown so that ‖dφ/ds‖ ≡ L pointwise. The boundary states y_0, y_{N_seg} are parameterized by the unstable / stable eigenvectors of the Hamiltonian Jacobian M = [J A; 0 -Jᵀ] at each fixed-point endpoint. The BVP is segmented into nshoots shooting segments stitched by Newton-iterated continuity, solved by NonlinearSolveFirstOrder.NewtonRaphson with a sparse ForwardDiff Jacobian (block-bidiagonal sparsity contracted by SparseMatrixColorings.GreedyColoringAlgorithm) and LinearSolve.UMFPACKFactorization().

MultipleShooting dispatches on FreidlinWentzellHamiltonian, not CoupledSDEs, because the shooting method operates on the LDT-derived deterministic Hamiltonian system rather than the original stochastic system. Passing a CoupledSDEs raises an ArgumentError pointing at the explicit construction:

H = FreidlinWentzellHamiltonian(sys)
res = minimize_geometric_action(H, x_initial, MultipleShooting(; nshoots = 10))

Endpoint constraints. Both endpoints must be hyperbolic fixed points of the drift (‖b(x*)‖ ≈ 0, ∂b(x*) has no purely-imaginary eigenvalues). Non-fixed-point endpoints and non-hyperbolic fixed points are rejected with an ArgumentError. Use GeometricGradient (gMAM/sgMAM) in either regime.

No interior fixed-point crossings. The arclength reparameterization develops a singularity at any (φ, p) where H_p = 0, which on the H = 0 shell is exactly (x*, 0) for x* a drift fixed point. The BVP-integrated portion of the path therefore must not pass through such a fixed point. The classical 1D bistable transition xa = −1 → xb = +1 (which crosses the saddle x = 0) does not work as a single BVP; the user must split it into the noise-driven attractor → saddle leg and the deterministic saddle → attractor leg. Only the uphill attractor → saddle leg carries action; the downhill saddle → attractor leg follows the deterministic drift and contributes zero Freidlin-Wentzell action. The transition action is therefore the action of the single uphill leg (and, when several competing barrier legs exist, the maximum over them), not the sum of both attractor → saddle solves:

res_left  = minimize_geometric_action(H, x_init_left,  MultipleShooting())  # -1 → 0 (uphill)
res_right = minimize_geometric_action(H, x_init_right, MultipleShooting())  # +1 → 0 (uphill)
# -1 → +1 transition: uphill -1 → 0, then deterministic 0 → +1 (zero action).
S_transition = res_left.action
# Symmetric case: res_left.action ≈ res_right.action.

There is no runtime guard (cheap interior-fixed-point detection requires integrating, which the residual function does each Newton iteration anyway). If a user passes a through-saddle path, the BVP typically converges to a degenerate point with near-zero action; the symptom is res.action close to 0 when a nontrivial transition was expected.

Warm-starting. The BVP Newton iteration has a finite basin of attraction. Warm-starting from a GeometricGradient (gMAM/sgMAM) solve helps on nontrivial problems; x_initial follows the same D × N matrix convention as gMAM, so the swap is a one-line change.

The output MinimumActionPath carries the BVP-integrated path in path, the conjugate momentum in generalized_momentum, and the converged path length in λ.

CriticalTransitions.MultipleShooting — Type

struct MultipleShooting{S, J, T} <: CriticalTransitions.GMAMOptimizer

Multiple-shooting BVP optimizer for the Freidlin-Wentzell instanton on a FreidlinWentzellHamiltonian. Integrates the arclength-reparameterized Hamilton equations

\[\frac{\mathrm{d}\varphi}{\mathrm{d}s} = \alpha\,H_p,\qquad \frac{\mathrm{d}p}{\mathrm{d}s} = -\alpha\,H_x,\qquad \alpha = L / \|H_p\|\]

on s ∈ [0, 1] with path length L a Newton unknown. Boundary states are parameterized by the unstable / stable eigenvectors of the Hamiltonian Jacobian at each fixed-point endpoint. The heteroclinic must not cross a drift fixed point in its interior; through-saddle problems must be user-split into attractor → saddle legs and the actions summed.

Fields

nshoots::Int64: Number of shooting segments (≥ 2).
ode_solver::Any: ODE integrator passed to the per-segment solve.
nlsolve::Any: Override for the Newton solver; nothing builds a default sparse-AD NewtonRaphson.
abstol::Any: Newton absolute tolerance.
reltol::Any: Newton relative tolerance.
eps_lin::Any: Anchor magnitude: distance of y_0 from (xa, 0) in 2D-space.
invariant_tol::Any: Warning threshold on max_s |H(φ(s), p(s))| along the converged path.
maxiters::Int64: Maximum Newton iterations.

source

`MinimumActionPath`

(gMAM) and (sgMAM) return their output as a MinimumActionPath type:

CriticalTransitions.MinimumActionPath — Type

struct MinimumActionPath{D, T<:Real, V, Phis, Ahis, Lambda, GPV, PV}

The minimum action path between two points in a D-dimensional phase space.

Fields

path::StateSpaceSet{D, T, V} where {D, T<:Real, V}: The path matrix.
action::Real: The action value associated to the path.
path_history::Any: The history of action of the paths in the optimisation algorithm (optional).
action_history::Any: The history of action of the paths in the optimisation algorithm (optional).
λ::Any: The Lagrange multiplier parameter for the minimum action path (optional).
generalized_momentum::Any: The generalized momentum of the phase space variables (optional).
path_velocity::Any: The path velocity (optional).

Constructors

MinimumActionPath(
    path,
    action;
    path_history,
    action_history,
    λ,
    generalized_momentum,
    path_velocity
)

defined at /home/runner/work/CriticalTransitions.jl/CriticalTransitions.jl/src/largedeviations/MinimumActionPath.jl:29.

source

Action functionals

Freidlin-Wentzell action

CriticalTransitions.fw_action — Function

fw_action(sys::CoupledSDEs, path, time) -> Any

Freidlin-Wentzell action of path (a D × N matrix with D = dimension(sys)) with time points time (length N) for the drift b = dynamic_rule(sys):

\[S_T[\phi] \;=\; \tfrac{1}{2}\int_0^T \big\| \dot\phi - \mathbf{b}(\phi)\big\|_{\mathbf{Q}(\phi)}^2 \,\mathrm{d}t\]

where $\|a\|_\mathbf{Q}^2 = \langle a, \mathbf{Q}^{-1} a\rangle$ (see anorm), $T$ is the total time of the path, and $\mathbf{Q}(x)$ is the trace-normalized diffusion tensor a(x) / s with s = tr(a(u₀))/D pinned at current_state(ds); for state-dependent diffusions s depends on the chosen u₀.

source

Geometric Freidlin-Wentzell action

CriticalTransitions.geometric_action — Function

geometric_action(sys::CoupledSDEs, path) -> Any
geometric_action(sys::CoupledSDEs, path, arclength) -> Any

Geometric Freidlin-Wentzell action of path (a D × N matrix) with the given arclength, for the drift b = dynamic_rule(sys). Equals $\inf_T S_T[\varphi]$ of the time-parameterized fw_action on the instanton:

\[\bar S[\varphi] \;=\; \int_0^L \Big( \|\varphi'\|_\mathbf{Q}\,\|\mathbf{b}(\varphi)\|_\mathbf{Q} \;-\; \langle \varphi', \mathbf{b}(\varphi)\rangle_\mathbf{Q} \Big)\,\mathrm{d}s\]

where $s$ is the arclength coordinate, $L$ the arclength, and $\|a\|_\mathbf{Q}^2 = \langle a, b\rangle_\mathbf{Q} = a^\top \mathbf{Q}^{-1} b$ (see anorm). $\mathbf{Q}$ is the trace-normalized covariance_matrix(sys). As with fw_action, the returned value is independent of the noise_strength keyword and rotation-invariant.

source

geometric_action(b::Function, path, arclength=1.0; A=nothing)

Geometric Freidlin-Wentzell action for a drift field b and a discrete path.

This is a drift-only convenience overload that uses the same integrand as geometric_action(sys::CoupledSDEs, ...), but requires an explicit inverse covariance (A = Q^{-1}) if you want anything other than the identity metric.

If A === nothing, it defaults to the identity metric.

The path must be a (D × N) matrix with points in columns.

source

Onsager-Machlup action

CriticalTransitions.om_action — Function

om_action(sys::CoupledSDEs, path, time, noise_strength)

Onsager-Machlup action of path (a D × N matrix) with time points time (length N) for the drift b = dynamic_rule(sys), at the given noise_strength σ:

\[S^\sigma_T[\phi] \;=\; \tfrac{1}{2}\int_0^T \Big(\big\|\dot\phi - \mathbf{b}(\phi)\big\|_\mathbf{Q}^2 \;+\; \sigma^2 \,\nabla\cdot \mathbf{b}(\phi)\Big)\,\mathrm{d}t\]

where $\|a\|_\mathbf{Q}^2 = \langle a, \mathbf{Q}^{-1} a\rangle$ (see anorm), $T$ is the total time, and $\mathbf{Q}$ is the trace-normalized covariance_matrix(sys). The first term is exactly fw_action and is independent of the noise_strength keyword under the trace-normalize convention. The second term is the Onsager-Machlup finite-noise correction parameterized by the explicit noise_strength argument: pass the σ at which you want $-\log P[\phi]$ evaluated (typically the same noise_strength you used at construction). As $\sigma \to 0$, om_action → fw_action.

source

For convenience, a general action function is available where the type of functional is set as an argument:

CriticalTransitions.action — Function

action(
    sys::CoupledSDEs,
    path::AbstractMatrix,
    time,
    functional::AbstractString;
    noise_strength
) -> Any

Computes the action functional specified by functional for a given CoupledSDEs sys and path parameterized by time.

functional = "FW": Returns the Freidlin-Wentzell action (fw_action)
functional = "OM": Returns the Onsager-Machlup action (om_action)

source

Grid-based quasipotential (OLIM)

The Ordered Line Integral Method (Dahiya and Cameron, 2018) computes the entire Freidlin-Wentzell quasipotential field $U_A(x)$ on a Cartesian grid in a single Dijkstra-style sweep, given a stable attractor $x_A$. The dimension $D$ is taken from the CoupledSDEs type parameter; the method is most accurate for $D = 2, 3$.

using CriticalTransitions, StaticArrays
f(x, p, t) = SVector(x[1] - x[1]^3 - 5 * x[1] * x[2]^2,
                      -(1 + x[1]^2) * x[2])
sys  = CoupledSDEs(f, [-1.0, 0.0]; noise_strength = 0.3)
grid = CartesianGrid((-1.5, 1.5, 121), (-1.0, 1.0, 81))
qp   = quasipotential(sys, grid, [-1.0, 0.0])
qp.U[CartesianIndex(61, 41)]   # U at the saddle (0, 0)

See examples/quasipotential_maierstein.jl for a worked example with contour plot.

CriticalTransitions.quasipotential — Function

quasipotential(sys, grid, attractor; band_radius, near_source_layers,
               verbose, show_progress) -> QuasiPotential{D}

Compute the Freidlin-Wentzell quasipotential field U_A(x) from attractor using the Ordered Line Integral Method (Dahiya and Cameron 2018). The state dimension D is taken from sys::CoupledSDEs{IIP, D, I, P} and must match grid::CartesianGrid{D}. A warning is emitted for D > 4.

Keyword arguments

band_radius::Int = default_K(grid): accepted-band radius in grid cells.
near_source_layers::Int = 3: size of the analytic CARE seed box; 0 disables analytic seeding.
verbose::Bool = false
show_progress::Bool = true

String method

The string method is a technique for finding transition paths between two states in a dynamical system E et al. (2007). The method represents the path as a "string" of points that connects the states and evolves it to minimize the drift along the path. The resulting tangent path is parallel to the drift of the system, i.e., the string method computes the heteroclinic orbit. For non-gradient systems (detailed -balance is broken), the heteroclinic orbit differs from the transition path, it does correctly predict, it correctly captures the deterministic dynamics from the saddle point onward ("downhill" portion of the path).

CriticalTransitions.string_method — Function

string_method(
    sys::Union{Function, FreidlinWentzellHamiltonian},
    x_initial::AbstractMatrix;
    stepsize,
    integrator,
    maxiters,
    show_progress
) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}

Compute the string method for a given system using the method of E et al. (2007).

The string method is an iterative algorithm used to find minimum energy path (MEP) between two points in phase space. It works by evolving a discretized path (string) according to the system's drift while maintaining equal arc-length parametrization between points.

This implementation allows for computation between arbitrary points, not just stable fixed points.

References

Keyword arguments

stepsize::Real=1e-1: Step size for the evolution step. Default: 0.1
integrator=Euler(): SciML integrator algorithm (e.g. Euler(), Tsit5())
maxiters::Int=1000: Maximum number of iterations for path convergence
show_progress::Bool=false: Whether to display a progress meter during computation

Returns

A MinimumActionPath containing the converged path and action value.

For sys::CoupledSDEs the .action field is computed as the geometric Freidlin-Wentzell action via geometric_action. For drift-only systems (e.g. sys::Function) the same geometric action is computed assuming an identity noise covariance.

source

string_method(
    sys::ContinuousTimeDynamicalSystem,
    init;
    kwargs...
) -> MinimumActionPath{_A, _B, _B1, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _B1}

Compute the string method for a given system using the method of E et al. (2007).

This implementation allows for computation between arbitrary points, not just stable fixed points.

References

(E et al., 2007)

Keyword arguments

stepsize::Real=1e-1: Step size for the evolution step. Default: 0.1
integrator=Euler(): SciML integrator algorithm (e.g. Euler(), Tsit5())
maxiters::Int=1000: Maximum number of iterations for path convergence
show_progress::Bool=false: Whether to display a progress meter during computation

Returns

A MinimumActionPath containing the converged path and action value.

For sys::CoupledSDEs, the action is computed using the system's noise covariance.

source