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 [8] and [9]. An overview of numerical methods applying LDT is given in Grafke and Vanden-Eijnden [10].

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.

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.

To summarize, the following methods are available:

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

Method	Use when	Requirements (this package)	Not suitable when
MAM	You want a minimum action path for a specified travel time $T$ (FW/OM action minimization).	`CoupledSDEs` with additive, invertible, autonomous noise (constant covariance); discretized path uses equispaced time.	Multiplicative/state-dependent noise, degenerate/non-invertible noise, non-autonomous noise; when the transition time is unknown and you want a time-reparameterization-invariant formulation (prefer gMAM/sgMAM).
gMAM	You want a time-reparameterization-invariant minimum action path (no explicit optimization over $T$).	Same as MAM for `CoupledSDEs`: additive, invertible, autonomous noise (constant covariance).	Multiplicative/state-dependent, degenerate, or non-autonomous noise; if you need the Onsager–Machlup functional (only FW has a geometric formulation).
sgMAM	You want a Hamiltonian/simple gMAM formulation that can be efficient in practice.	Same as MAM and diagonal noise covariance (implementation restriction); assumes additive noise.	Non-diagonal covariance; multiplicative/state-dependent or degenerate noise; models where the required derivatives/Jacobian are not available or are too expensive.
String method	You want a minimum energy path / heteroclinic orbit driven by the deterministic drift (typical use: gradient systems).	Deterministic drift field (works for `ContinuousTimeDynamicalSystem`; does not rely on an SDE noise model).	In non-gradient systems if you need the most probable noise-induced transition path (string gives the deterministic heteroclinic orbit, which generally differs from the instanton).

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 [11].
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. [12].
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. [13].
Multiplicative / state-dependent noise: extensions of geometric action minimization to degenerate or multiplicative noise are discussed in Grafke et al. [14]; the current sgMAM implementation assumes additive noise.

Minimum action method (MAM)

Minimization of the specified action functional using the optimization algorithm of Optimization.jl. See also E et al. [15].

CriticalTransitions.min_action_method — Function

min_action_method(sys::ContinuousTimeDynamicalSystem, x_i, x_f, T::Real; 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

min_action_method(sys::ContinuousTimeDynamicalSystem, init::Matrix, T::Real; kwargs...).

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

Keyword arguments

functional = "FW": type of action functional to minimize. Defaults to fw_action, alternative: "OM" for om_action
points = 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"
optimizer = Optimisers.Adam(): minimization algorithm from Optimization.jl
ad_type = Optimization.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

min_action_method(sys::ContinuousTimeDynamicalSystem, init::Matrix, T::Real; 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$.

source

Geometric minimum action method (gMAM)

Minimization of the geometric action following Heymann and Vanden-Eijnden [16] and Heymann and Vanden-Eijnden [17]. 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.geometric_min_action_method — Function

geometric_min_action_method(
    sys::ContinuousTimeDynamicalSystem,
    init::Matrix;
    maxiters,
    abstol,
    reltol,
    optimizer,
    ad_type,
    stepsize,
    verbose,
    show_progress
) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Nothing, Nothing, Nothing} where {_A, _B<:Real, _C}

Runs the geometric Minimum Action Method (gMAM) to find the minimum action path (instanton) from an initial condition init, given a system sys and total arc length arclength.

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.

For more information see the main method, geometric_min_action_method(sys::CoupledSDEs, x_i, x_f, arclength::Real; kwargs...).

source

Simple geometric minimum action method (sgMAM)

Simplified minimization of the geometric action following Grafke et al. [18]. 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 assuming an autonomous system with additive Gaussian noise.

CriticalTransitions.simple_geometric_min_action_method — Function

simple_geometric_min_action_method(
    sys::ExtendedPhaseSpace,
    x_initial::Array{T, 2};
    stepsize,
    maxiters,
    show_progress,
    verbose,
    abstol,
    reltol
) -> MinimumActionPath{_A, _B, _C, Nothing, Nothing, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}} where {_A, _B<:Real, _C}

Performs the simplified geometric Minimal Action Method (sgMAM) on the given system sys. Our implementation is only valid for additive noise.

This method computes the optimal path in the phase space of a Hamiltonian system that minimizes the Freidlin–Wentzell action. The Hamiltonian functions H_x and H_p define the system's dynamics in a doubled phase. The initial state x_initial is evolved iteratively using constrained gradient descent with step size parameter stepsize over a specified number of iterations. The method can display a progress meter and will stop early if the absolute tolerance abstol or relative tolerance reltol is achieved.

The function returns a MinimumActionPath containing the final path, the action value, the Lagrange multipliers (.λ), the momentum (.generalized_momentum), and the state derivatives (.path_velocity). The implementation is based on the work of Grafke et al. (2019).

Keyword arguments

stepsize::Real=1e-1: step size for gradient descent. Default: 0.1
maxiters::Int=1000: maximum number of iterations before the algorithm stops
show_progress::Bool=false: if true, display a progress bar
verbose::Bool=false: if true, print additional output
abstol::Real=NaN: absolute tolerance for early stopping based on action change
reltol::Real=NaN: relative tolerance for early stopping based on action change

source

CriticalTransitions.ExtendedPhaseSpace — Type

A structure representing a extanded phase space system where your dissipative vector field is embedded in a doubled dimensional phase space. Given old phase space coordinates x of a vector field f(x), we can define the canonical momenta p, such that the new phase space coordinates are (x, p). The dynamics in this extended phase space are then governed by the Hamtiltonian system:

$H = p^2 + x \dot f(x)$

Hence, this system operates in an extended phase space where the Hamiltonian is assumed to be quadratic in the extended momentum.

The struct ExtendedPhaseSpace holds the Hamilton's functions H_x and H_p.

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:

ExtendedPhaseSpace(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 ExtendedPhaseSpace(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 = ExtendedPhaseSpace{false,2}(H_x, H_p)  # hardcoded analytic H_x/H_p

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

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

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

`MinimumActionPath`

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

CriticalTransitions.MinimumActionPath — Type

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

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

Calculates the Freidlin-Wentzell action of a given path with time points time in a drift field specified by the deterministic dynamics f = dynamic_rule(sys) and (normalized) noise covariance matrix covariance_matrix(sys).

The path must be a (D x N) matrix, where D is the dimensionality of the system sys and N is the number of path points. The time array must have length N.

Returns a single number, which is the value of the action functional

$S_T[\phi_t] = \frac{1}{2} \int_0^T || \dot \phi_t - f(\phi_t) ||^2_Q \text{d}t$

where $\phi_t$ denotes the path in state space, $b$ the drift field, and $T$ the total time of the path. The subscript $Q$ refers to the generalized norm $||a||_Q^2 := \langle a, Q^{-1} b \rangle$ (see anorm). Here $Q$ is the noise covariance matrix normalized by $D/L_1(Q)$, with $L_1$ being the L1 matrix norm.

source

Geometric Freidlin-Wentzell action

CriticalTransitions.geometric_action — Function

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

Calculates the geometric action of a given path with specified arclength for the drift field specified by the deterministic dynamics f = dynamic_rule(sys) and (normalized) noise covariance matrix covariance_matrix(sys).

For a given path $\varphi$, the geometric action $\bar S$ corresponds to the minimum of the Freidlin-Wentzell action $S_T(\varphi)$ over all travel times $T>0$, where $\varphi$ denotes the path's parameterization in physical time (see fw_action). It is given by the integral

$\bar S[\varphi] = \int_0^L \left( ||\varphi'||_Q \, ||f(\varphi)||_Q - \langle \varphi', \, f(\varphi) \rangle_Q \right) \, \text{d}s$

where $s$ is the arclength coordinate, $L$ the arclength, $f$ the drift field, and the subscript $Q$ refers to the generalized dot product $\langle a, b \rangle_Q := a^{\top} \cdot Q^{-1} b$ (see anorm). Here $Q$ is the noise covariance matrix normalized by $D/L_1(Q)$, with $L_1$ being the L1 matrix norm.

Returns the value of the geometric action $\bar S$.

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)

Calculates the Onsager-Machlup action of a given path with time points time for the drift field f = dynamic_rule(sys) at given noise_strength.

The path must be a (D x N) matrix, where D is the dimensionality of the system sys and N is the number of path points. The time array must have length N.

Returns a single number, which is the value of the action functional

$S^{\sigma}_T[\phi_t] = \frac{1}{2} \int_0^T \left( || \dot \phi_t - f(\phi_t) ||^2_Q + \sigma^2 \nabla \cdot f \right) \, \text{d} t$

where $\phi_t$ denotes the path in state space, $b$ the drift field, $T$ the total time of the path, and $\sigma$ the noise strength. The subscript $Q$ refers to the generalized norm $||a||_Q^2 := \langle a, Q^{-1} b \rangle$ (see anorm). Here $Q$ is the noise covariance matrix normalized by $D/L_1(Q)$, with $L_1$ being the L1 matrix norm.

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::Matrix,
    time,
    functional;
    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

String method

The string method is a technique for finding transition paths between two states in a dynamical system E et al. [19]. 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, ExtendedPhaseSpace},
    x_initial::Matrix;
    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. [19].

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. [19].

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

References

[19]

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