Library

HODEProblem(result::HamiltonianSINDyResult, timespan, timestep, q₀, p₀; kwargs...)

The identified Hamiltonian system as a GeometricEquations.HODEProblem, ready for integrate with a symplectic integrator.

The problem carries the identified Hamiltonian, so the energy behaviour of the integrated solution can be checked directly with GeometricSolutions.compute_invariant_error.

Examples

result = identify(TrajectoryData(x, y, Δt), HamiltonianSINDy(λ = 0.05, polyorder = 2))
prob   = HODEProblem(result, (0.0, 10.0), 0.01, q₀, p₀)

using GeometricIntegrators
sol = integrate(prob, ImplicitMidpoint())

ODEProblem(result::SINDyResult, timespan, timestep, ics...; kwargs...)

The identified system as a GeometricEquations.ODEProblem, ready for integrate.

Examples

result = identify(TrainingData(x, ẋ), SINDy(CompoundBasis(polyorder = 3)))
prob   = ODEProblem(result, (0.0, 10.0), 0.01, x₀)

using GeometricIntegrators
sol = integrate(prob, ExplicitMidpoint())

AbstractBasis

Supertype of the candidate-function libraries a sparse regression selects from.

A basis is defined symbolically by basis_functions, which returns its candidate functions as Symbolics expressions in the state. Everything else follows from that one definition: evaluate compiles the expressions into a fast numerical evaluator, and the Hamiltonian methods differentiate them to build J∇φₖ. There is deliberately no second, numeric definition that could drift from the symbolic one.

Concrete bases are PolynomialBasis, TrigonometricBasis, ExponentialBasis, LogarithmicBasis, RationalBasis and CompoundBasis.

BasisArguments

Supertype of the argument selections a univariate basis can be applied to: StateComponents and Differences.

CompoundBasis(bases...)
CompoundBasis(; polyorder = 5, trigonometric = 0)

A basis assembled from several others, evaluated in the order given.

The keyword form builds the common case: the constant and all monomials up to polyorder, optionally followed by trigonometric terms up to wavenumber trigonometric.

Bases are combined with ⊕ as well:

julia> basis = CompoundBasis(polyorder = 2) ⊕ ExponentialBasis(; rates = (-1.0,));

julia> nterms(basis, 2)
8

Differences(indices; consecutive = false)

Apply the basis functions to differences of the state components named by indices.

With consecutive = true only neighbouring differences z[i₊₁] - z[i] are formed, which is what a lattice with nearest-neighbour interaction needs. Otherwise all pairs z[i] - z[j] with i > j are formed, which is what an all-to-all interaction needs.

Examples

A Toda lattice of four particles interacts through consecutive differences of the positions, the first four of eight phase-space components:

julia> args = Differences(1:4; consecutive = true);

julia> basis = ExponentialBasis(args; rates = (-1.0,));

julia> nterms(basis, 8)
3

ExponentialBasis(args = StateComponents(); rates = (1.0,))

The basis functions exp(α u) for each rate α and each argument u.

The Toda lattice needs exp(-(qₙ₊₁ - qₙ)), so both a negative rate and a Differences argument selection:

julia> basis = ExponentialBasis(Differences(1:3; consecutive = true); rates = (-1.0,));

julia> nterms(basis, 6)
2

HamiltonianFunctions

The compiled functions of a parametrised Hamiltonian H(z; a), in the form GeometricEquations expects.

Fields

• H: the Hamiltonian, callable as H(t, q, p, params) with params.a the coefficients
• v: q̇ = ∂H/∂p, callable as v(v, t, q, p, params)
• f: ṗ = -∂H/∂q, callable as f(f, t, q, p, params)
• ż: the combined field J∇H, callable as ż(out, z, a) — the form the regression uses
• d: the number of degrees of freedom
• nparam: the number of coefficients

Splitting v and f is what lets an identified system become a HODEProblem; the combined ż is kept because the fitting loop evaluates the whole field at once.

HamiltonianSINDy(; λ, polyorder, trigonometric, integrator_timestep, nloops, picard_iterations)

Sparse identification of a Hamiltonian system (Khan 2023).

A scalar Hamiltonian is parametrised over a library of candidate functions, H(z; a) = Σₖ aₖ φₖ(z), and the vector field is obtained as ż = J ∇H(z). The identified dynamics is therefore Hamiltonian by construction rather than by penalty: every candidate the fit considers is a symplectic gradient field, whatever the coefficients happen to be.

The coefficients are fitted by matching the flow map — minimising Σⱼ ‖Φ_Δt(zⱼ; a) − zⱼ₊₁‖² where Φ_Δt is an implicit-midpoint step. This is nonlinear in a and needs an optimiser.

Apply it to a TrajectoryData problem with identify.

The basis may be given directly, which is what the systems the thesis treats require — a Toda lattice needs exp of a difference of positions, a point vortex log of one:

HamiltonianSINDy(hamiltonian_basis(polyorder = 2) ⊕
                 ExponentialBasis(Differences(1:4; consecutive = true); rates = (-1.0,)))

The keyword form HamiltonianSINDy(; polyorder, trigonometric) builds the polynomial and trigonometric library and is equivalent to passing hamiltonian_basis.

HamiltonianSINDyResult

The outcome of identify with HamiltonianSINDy.

Reach the coefficients with parameters and the compiled Hamiltonian and its symplectic gradient with functions. Convert to a GeometricEquations.HODEProblem to integrate the identified system with a symplectic integrator.

HamiltonianSINDyVectorField(result)

The identified Hamiltonian vector field J∇H, callable as f(dz, z).

IdentificationProblem

Supertype of the data an identification method is applied to.

An identification problem plays the role that a GeometricEquations problem plays for an integrator: it is the thing a method is applied to. identify is the verb, so that

identify(problem, method)

reads as GeometricIntegrators' integrate(problem, method) does.

Two concrete problems, distinguished by what the data means rather than by its shape: TrainingData for matching a vector field, TrajectoryData for matching a flow map.

JuliaLeastSquare()

Solve the regression Θ x = ẋ by the ordinary least-squares solution Θ \ ẋ.

This is the right solver whenever the model is linear in the coefficients, which is the case for every method in this package that fits a fixed library of candidate functions.

LogarithmicBasis(args = StateComponents())

The basis functions log(abs(u)) over the arguments args.

A point-vortex Hamiltonian is built from log|qᵢ - qⱼ|, so this is normally paired with Differences. abs is applied inside the logarithm so the basis is defined on both signs of the argument; it is singular where the argument vanishes, which for a difference means two coordinates coinciding.

OptimizerSolver(; algorithm = BFGS(), linesearch = Backtracking())

Minimise a nonlinear least-squares loss with GeometricOptimizers.

Only needed where the model is not linear in its coefficients — in this package that is the flow-map formulation of HamiltonianSINDy, where the coefficients enter through an integrator. Prefer JuliaLeastSquare wherever the problem is linear.

PolynomialBasis(p)

All monomials of degree exactly p, each one once.

Variables may repeat within a monomial, so degree 2 in two variables is z₁², z₁z₂, z₂². Degree 0 is the constant. For a state of dimension d there are binomial(d + p - 1, p) terms of degree p.

RationalBasis(args = StateComponents(); powers = (1,))

The basis functions abs(u)^-k for each k in powers, over the arguments args.

An N-body gravitational Hamiltonian is built from 1/|qᵢ - qⱼ|, so this is normally paired with Differences. abs is applied inside the power, as it is in LogarithmicBasis, so that the basis functions depend on the distance between two coordinates and not on their order. Singular where the argument vanishes, which for a difference means two coordinates coinciding.

SINDy(basis; λ = 0.05, nloops = 10)

Sparse Identification of Nonlinear Dynamics (Brunton, Proctor & Kutz, PNAS 2016).

Fits ẋ = Θ(x) Ξ over the fixed library basis and sparsifies Ξ by sequentially thresholded least squares: threshold every coefficient below λ to zero, refit on the surviving support, and repeat until the support stops changing or nloops is reached.

λ is a hard threshold on coefficient magnitude applied after each least-squares fit — it is not an ℓ¹ penalty and does not appear in the objective being minimised.

Apply it to a TrainingData problem with identify.

Examples

julia> A = [-0.1 2.0; -2.0 -0.1];

julia> x = randn(2, 500);

julia> result = identify(TrainingData(x, A * x), SINDy(CompoundBasis(polyorder = 3)));

julia> isapprox(parameters(result)[2:3, :], A'; atol = 1e-10)
true

SINDyResult

The outcome of identify with SINDy.

Carries the coefficient matrix and the method that produced it. Reach the coefficients with parameters and the basis they refer to with basis; convert to a GeometricEquations.ODEProblem to integrate the identified system.

SINDyVectorField(result)

The identified vector field, callable as f(dy, y, params, t).

SparsificationMethod

Supertype of the identification methods, e.g. SINDy and HamiltonianSINDy.

A method carries the basis it searches and the parameters of its sparsification. It is applied to an IdentificationProblem with identify, mirroring the way a GeometricIntegrators method is applied to a problem with integrate.

Traits

SparsificationMethod extends GeometricBase.AbstractMethod, and the ecosystem's trait functions are answered where they mean something for a regression method:

• issymplectic and isenergypreserving report whether the identified model is symplectic and energy-preserving by construction. This is the substantive question for a structure-preserving method, and the reason HamiltonianSINDy exists.
• isexplicit and isimplicit report whether the regression is a direct solve or goes through an implicit integrator.
• name, description and reference identify the method and its source.

issymmetric, isstifflyaccurate and order describe a Runge–Kutta tableau and have no meaning for a regression, so they are left as missing. GeometricBase.isAbstractMethod therefore returns false for these methods, which is correct: it is a conformance check for integrator methods.

StateComponents()

Apply the basis functions to each state component separately: f(z₁), f(z₂), ….

TrainingData(x, ẋ)
TrainingData(solution::GeometricSolution)

States and the corresponding time derivatives, for methods that match a vector field.

x and ẋ are either matrices whose columns are snapshots, or vectors of state vectors. Both must describe the same number of snapshots of the same dimension; that is checked at construction, because the alternative is a shape error surfacing much later as a wrong answer.

Examples

julia> data = TrainingData(randn(2, 20), randn(2, 20));

julia> nsamples(data), statedimension(data)
(20, 2)

TrainingData(solution::GeometricSolution)

Training data from an integrated solution, pairing each stored state with the vector field evaluated there.

This closes the loop: integrate a known problem, identify it from the solution, and compare.

TrajectoryData(x, y, Δt)
TrajectoryData(solution::GeometricSolution)

Consecutive states, for methods that match a flow map: y[j] is the state one step of size Δt after x[j].

Use this where the derivatives are not available and only sampled trajectories are — matching the flow map avoids differentiating the data, at the cost of a nonlinear regression.

Examples

julia> data = TrajectoryData([randn(2) for _ in 1:8], [randn(2) for _ in 1:8], 0.01);

julia> nsamples(data), statedimension(data), timestep(data)
(8, 2, 0.01)

TrajectoryData(solution::GeometricSolution)

Consecutive states from an integrated solution, for flow-map matching.

Pairs each stored state with its successor, taking Δt from the solution's own time step.

TrigonometricBasis(n, args = StateComponents())

The basis functions sin(k u) and cos(k u) for 1 ≤ k ≤ n, over the arguments args.

basis(method::SparsificationMethod)

The library of candidate functions the method searches.

nsamples(problem)

The number of snapshots in an IdentificationProblem.

b₁ ⊕ b₂

Concatenate two bases into a CompoundBasis.

_ispartitioned(solution)

Whether solution stores a conjugate momentum alongside the position.

GeometricSolutions defines hasproperty on the type of the data series, so this is answered at compile time rather than by probing the object.

basis_functions(basis, z)

The candidate functions of basis as symbolic expressions in the state vector z.

This is the definition of a basis; evaluate is generated from it.

calculate_nparams(d, polyorder, trig_wave_num)

The number of coefficients a Hamiltonian basis carries for d degrees of freedom.

The phase space has 2d dimensions, so the polynomial part counts the monomials of degree up to polyorder in 2d variables. The trigonometric part, when trig_wave_num > 0, adds sin and cos at each wave number for each variable.

degreesoffreedom(result)

The number of degrees of freedom d; the phase space has 2d dimensions.

evaluate(data, basis)

Evaluate every candidate function of basis on every snapshot of data.

data is a matrix whose columns are snapshots, a vector of state vectors, or a single state vector. The result Θ has one row per snapshot and one column per candidate function.

Examples

julia> Θ = evaluate([1.0 2.0; 3.0 4.0], CompoundBasis(polyorder = 1));

julia> size(Θ)      # 2 snapshots × (1 constant + 2 linear) terms
(2, 3)

hamilGrad_func_builder(d, polyorder, trig_wave_num)

The symplectic gradient J∇H of the parametrised Hamiltonian, callable as out = f(out, z, a).

A thin wrapper over hamiltonian_functions for the regression, which needs only the combined field.

hamiltonian_basis(; polyorder = 3, trigonometric = 0)

The polynomial (and optionally trigonometric) library used for a Hamiltonian ansatz.

The constant is omitted: it contributes nothing to ∇H, so its coefficient is unidentifiable and its column of the regression matrix is identically zero.

hamiltonian_functions(basis, d)
hamiltonian_functions(d, polyorder, trig_wave_num)

Build the symbolic Hamiltonian H(z; a) = Σₖ aₖ φₖ(z) over d degrees of freedom — so 2d phase-space variables — and compile it, together with its symplectic gradient, into HamiltonianFunctions.

Constant terms are dropped, since they are unidentifiable; see strip_constants.

hamiltonian_poly(z, order, inds...)

All monomials of degree exactly order in the variables z, each one once.

A variable may repeat within a monomial, so degree 2 in two variables gives z₁², z₁z₂, z₂².

Shared by PolynomialBasis and by the symbolic Hamiltonian construction, so that the two cannot disagree about which terms exist or in what order.

identify(problem::IdentificationProblem, method::SparsificationMethod; kwargs...)

Identify the governing equations of problem with method.

This mirrors GeometricIntegrators' integrate(problem, method): the problem carries the data, the method carries the basis and the sparsification parameters, and the result carries the identified coefficients.

Which methods accept which problems is not interchangeable, because the two formulations need different data:

Applying a method to a problem it does not accept raises an ArgumentError naming both, rather than a MethodError from several layers down.

Examples

julia> A = [-0.1 2.0; -2.0 -0.1];

julia> x = randn(2, 400);

julia> result = identify(TrainingData(x, A * x), SINDy(CompoundBasis(polyorder = 3); λ = 0.05));

julia> nterms(result)
4

minimize(loss, x₀, solver::OptimizerSolver)

Minimise loss starting from x₀, returning the minimiser. x₀ is not modified.

nloops(method::SparsificationMethod)

The cap on thresholding passes.

nterms(basis, d)

The number of candidate functions the basis contributes for a state of dimension d.

nterms(result)

The number of Hamiltonian basis terms retained after sparsification.

nterms(result)

The number of library terms retained after sparsification.

sparsify(method::HamiltonianSINDy, hfuns, problem::TrajectoryData, solver; verbose = false)

Sequentially thresholded regression of the Hamiltonian coefficients against the flow map.

sparsify(method::SINDy, Θ, ẋ, solver)

Sequentially thresholded least squares.

Returns the coefficient matrix Ξ with size(Ξ) == (size(Θ, 2), size(ẋ, 1)).

sparsity_threshold(method::SparsificationMethod)

The sparsification threshold: coefficients below it are set to zero after each fit.

statedimension(problem)

The dimension of a single state.

strip_constants(basis, z)

The basis functions of basis whose gradient is not identically zero.

A term with vanishing gradient — a constant — cannot be identified from ż = J∇H: it contributes an all-zero column, which makes the linear formulation singular and wastes a parameter in the nonlinear one. Filtering on the gradient rather than on the type of the term catches every such case, whatever basis it came from.