Large deviation theory

In the weak-noise limit, the probability of a stochastic system following a given path obeys a large deviation principle (LDP): transitions between metastable states concentrate, with overwhelming probability, on a single most probable path called the instanton or maximum likelihood path. The instanton minimizes a rate functional (the action), and its minimal value is the quasipotential, which sets the exponential transition rate and the most likely exit location.

This page follows that chain:

The action: the rate functional being minimized.
Finding the instanton: minimizing the action to obtain the most probable path.
The quasipotential: the minimized action as a field, and the transition rate.

For the underlying theory see Freidlin and Wentzell (1998) and (Börner, 2025); for a survey of numerical methods see Grafke and Vanden-Eijnden (2019).

System class

The methods on this page 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. The convention relating $\sigma$, $\mathbf{Q}$ and noise_strength is described under Noise strength and the trace convention.

The action

The action (or rate functional) assigns a cost to every candidate path; the instanton is its minimizer and the minimal cost is the leading-order $-\log$ probability of the transition. Three forms are available:

Freidlin-Wentzell action, fw_action: the standard rate functional for a path traversed in a fixed transition time $T$.
Geometric Freidlin-Wentzell action, geometric_action: a time-reparameterization-invariant form, used when $T$ is not fixed.
Onsager-Machlup action, om_action: the finite-noise functional that adds the $\mathcal{O}(\sigma^2)$ drift-divergence correction; it reduces to the Freidlin-Wentzell action as $\sigma \to 0$. Additive noise only.

All three evaluate the integrand in the canonical covariance metric fixed by the trace convention. A practical consequence: fw_action and geometric_action are independent of the noise_strength chosen at construction, whereas om_action takes $\sigma$ as an explicit argument because its correction scales with $\sigma^2$; pass the $\sigma$ at which you want $-\log P[\phi]$ evaluated.

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

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

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

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

Finding the instanton

The instanton is the minimizer of the action. Several methods have been proposed to compute it; in the weak-noise limit, the minimum action path coincides with 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. In gradient systems, minimum action paths between attractors coincide with heteroclinic orbits, which can be computed via the 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.

Method	Use when
MAM (`minimize_action`)	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 (`minimize_geometric_action`)	$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.
Multiple shooting (`MultipleShooting`)	Endpoints are hyperbolic fixed points and you want a high-accuracy boundary-value solution; best warm-started from a gMAM/sgMAM solve.
String method (`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.

Action minimization as an optimal control problem

All of the action minimizers 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.

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 and sgMAM. The diffusion tensor $a(x)$ is classified once when the sgMAM cache is built and the resulting inner loop dispatches at compile time on the cache type (see Developer notes & internals); the Hamiltonian formulation is given in Grafke et al. (2017). 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 / sgMAM)

gMAM minimizes the geometric action following Heymann and Vanden-Eijnden (2008), Heymann and Vanden-Eijnden (2008). It reformulates MAM to avoid the double optimization over both the action and the transition time, 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.

The simple gMAM (sgMAM) following Grafke et al. (2017) 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 performs a constrained gradient descent on the Hamiltonian system, supporting autonomous diffusions with additive, diagonal-multiplicative, or general-matrix multiplicative noise.

The gradient step size is set by the GeometricGradient optimizer passed to minimize_geometric_action; for difficult problems the probe-based AdaptiveGeometricGradient variant is more robust (see the adaptive step-size example).

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. Rank-deficient $a(x)$ is rejected when the per-path cache is built. The compile-time classification of $a(x)$ is described under Developer notes & internals.

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

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 (internal helper _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 (internal helper 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. The 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.

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)

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 BVP is segmented into nshoots shooting segments stitched by Newton-iterated continuity; the solver internals are described under Developer notes & internals.

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 for this. 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

String method

The string method finds transition paths between two states by representing the path as a "string" of points connecting the states and evolving it to minimize the drift along the path E et al. (2007). 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 (where detailed balance is broken) the heteroclinic orbit differs from the transition path, but it correctly captures the deterministic dynamics from the saddle point onward (the "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

`MinimumActionPath`

gMAM/sgMAM and multiple shooting return their output as a MinimumActionPath:

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

The quasipotential

The quasipotential $U_A(\mathbf{x})$ is the minimal action required to reach a state $\mathbf{x}$ from an attractor $\mathbf{x}_A$. It generalizes the potential of gradient systems to non-gradient drift and sets both the exponential transition rate ($\sim e^{-U_A/\sigma^2}$) and the most likely exit location. There are two ways to obtain it:

Along a single instanton: minimize the action between two states with one of the minimizers above and read off the action value.
As a field: the Ordered Line Integral Method (OLIM) (Dahiya and Cameron, 2018) computes the entire quasipotential field on a Cartesian grid in a single Dijkstra-style sweep, given a stable attractor $\mathbf{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.
regularization::Real = default_regularization(grid): for a system with a single noiseless coordinate (rank-1 diffusion, e.g. momentum-only noise in a second-order Langevin or van der Pol oscillator) the trace-normalized diffusion is singular and has no metric. This amount is added to that coordinate's diagonal (a structural vanishing-viscosity perturbation that leaves the noisy block exact) to make it invertible; the degenerate Freidlin-Wentzell quasipotential is recovered as it tends to 0. The default scales like 1 / minimum(grid.nbox) (about 0.04 at 80 cells), so it shrinks under refinement; smaller values are more accurate but stiffer. The escape sheet and saddle barrier carry negligible bias, the non-escape region carries the regularization bias. Ignored for full-rank systems; two or more noiseless coordinates (sub-Riemannian) or a non-coordinate-aligned null space are rejected.
verbose::Bool = false
show_progress::Bool = true