Getting Started

Installation

using Pkg
Pkg.add(url = "https://github.com/JuliaRCM/SparseIdentification.jl")

The package requires Julia 1.11 or later.

The shape of a SINDy run

Every identification has the same four parts, whichever method you use:

  1. Data — states, and either their derivatives or their successors one time step on.
  2. A basis — the library of candidate functions to search.
  3. A method — SINDy or HamiltonianSINDy, carrying the sparsification threshold.
  4. A vector field — the identified model, callable and integrable.

Classical SINDy

The two-dimensional damped linear oscillator is the illustrative example of the original SINDy paper, and a good first check because the answer is known exactly:

\[\dot x = -0.1x + 2y, \qquad \dot y = -2x - 0.1y .\]

using SparseIdentification

# the true system
A = [-0.1  2.0
     -2.0 -0.1]

# sample the state space and evaluate the true vector field on it
x = randn(2, 500)
ẋ = A * x

# search polynomials up to degree 3
basis = CompoundBasis(polyorder = 3, trigonometric = 0)

# identify
result = identify(TrainingData(x, ẋ), SINDy(basis; λ = 0.05))

The coefficient matrix maps library terms to state components. Rows 2 and 3 hold the linear terms, so that block should be A transposed:

parameters(result)[2:3, :]
2×2 Matrix{Float64}:
 -0.1  -2.0
  2.0  -0.1

and everything else should be exactly zero — not merely small:

Ξ = parameters(result)
all(iszero, Ξ[1, :]) && all(iszero, Ξ[4:end, :])
true

Note that the data here is sampled from the state space, not taken along a single trajectory. Both work; uniform sampling over a region covers the library's domain more evenly and needs no integration, which is why the thesis uses it throughout.

Identifying a Hamiltonian

For a Hamiltonian system, identify the Hamiltonian instead. Take the harmonic oscillator, $H = \tfrac{1}{2}(q^2 + p^2)$, whose flow is a rotation in phase space:

using SparseIdentification

Δt = 0.01
R  = [ cos(Δt) sin(Δt)
      -sin(Δt) cos(Δt)]          # the exact flow map over one step

x = [randn(2) for _ in 1:60]     # states
y = [R * xⱼ for xⱼ in x]         # the same states one step later

method = HamiltonianSINDy(λ = 0.05, integrator_timestep = Δt, polyorder = 2)
result = identify(TrajectoryData(x, y, Δt), method)
vf     = HamiltonianSINDyVectorField(result)

The identified field should reproduce $\dot z = (p, -q)$ at points it never saw:

dz = zeros(2)
z  = [0.7, -0.3]
vf(dz, z)
dz            # should be ≈ [-0.3, -0.7]
2-element Vector{Float64}:
 -0.30000249043390187
 -0.7000058581073575

Because the field is a symplectic gradient by construction, this model conserves some Hamiltonian exactly, whether or not the coefficients are right.

Data layout

Two container types, distinguished by what they mean rather than by shape:

SparseIdentification.TrainingData — Type
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)
source
SparseIdentification.TrajectoryData — Type
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)
source

Both accept either a matrix whose columns are snapshots, or a vector of state vectors. Both check at construction that the pieces describe the same number of snapshots of the same dimension — a shape error here otherwise surfaces much later as a wrong answer rather than an exception.

Noise belongs to the data

The estimator does not add noise. If you want to test robustness, add it when you generate the data:

using Random
Random.seed!(1234)

A  = [-0.1 2.0; -2.0 -0.1]
x  = randn(2, 500)
ẋ  = A * x
η  = 0.05
ẋ_noisy = ẋ .+ η .* randn(size(ẋ))       # noise added here, deliberately and reproducibly

basis = CompoundBasis(polyorder = 3, trigonometric = 0)
Ξ = parameters(identify(TrainingData(x, ẋ_noisy), SINDy(basis; λ = 0.05)))
Ξ[2:3, :]
2×2 Matrix{Float64}:
 -0.0997889  -1.99879
  2.0021     -0.0983713

Earlier versions of this package injected noise inside the fitting routine, which made results irreproducible and could not be switched off. Two identical calls now give identical answers.