Basis Libraries

The library of candidate functions is the single most consequential choice in a SINDy run. SINDy cannot discover a term you did not offer it, and every term you do offer costs conditioning and raises the chance of a spurious fit. This page covers what is available and how to choose.

Available bases

SparseIdentification.AbstractBasis — Type
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.

source
SparseIdentification.PolynomialBasis — Type
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.

source
SparseIdentification.ExponentialBasis — Type
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
source
SparseIdentification.LogarithmicBasis — Type
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.

source
SparseIdentification.RationalBasis — Type
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.

source
SparseIdentification.CompoundBasis — Type
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
source
SparseIdentification.basis_functions — Function
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.

source
SparseIdentification.evaluate — Function
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)
source

One definition, not two

A basis is defined symbolically, by basis_functions(basis, z), and everything else is derived from that single definition: evaluate compiles the expressions into a fast numerical evaluator (cached per basis and state dimension), and the Hamiltonian methods differentiate the same expressions to build J∇φₖ. There is deliberately no second, hand-written numeric path that could drift from the symbolic one.

using SparseIdentification, Symbolics

z = Symbolics.variables(:z, 1:2)
basis_functions(CompoundBasis(polyorder = 2), z)
6-element Vector{Symbolics.Num}:
     1
    z₁
    z₂
  z₁^2
 z₁*z₂
  z₂^2

Polynomials

CompoundBasis(polyorder = n) assembles the constant, then all monomials of degree 1 through $n$, without repetition. For $d$ degrees of freedom the count is $\binom{d+n}{n}$.

using SparseIdentification

basis = CompoundBasis(polyorder = 3, trigonometric = 0)
x = randn(2, 5)
Θ = evaluate(x, basis)
size(Θ)      # 5 snapshots × 10 terms: 1 + 2 + 3 + 4
(5, 10)

The column order is: the constant, then degree 1 in state order, then degree 2 as $x_1^2, x_1x_2, x_2^2$, and so on.

Θ[:, 1] ≈ ones(5), Θ[:, 2] ≈ x[1, :], Θ[:, 4] ≈ x[1, :] .^ 2
(true, true, true)

Growth is combinatorial. In four variables — a two-degree-of-freedom Hamiltonian system — a cubic library has 34 terms, a quintic one 125. The thesis routinely works with libraries of roughly 100 to 170 terms, which is a reasonable ceiling.

Trigonometric terms

TrigonometricBasis(n) supplies $\sin(k x_i)$ and $\cos(k x_i)$ for $1 \le k \le n$ on every component:

btrig = TrigonometricBasis(2)
Θt = evaluate(x, btrig)
size(Θt)     # 5 × 8:  2 wavenumbers × {sin, cos} × 2 components
(5, 8)

These are essential for pendulum-like systems, whose Hamiltonians carry $\cos q$ and which no polynomial library can represent. CompoundBasis(polyorder = p, trigonometric = n) combines both.

Choosing a library

Start from what you know about the system. The thesis's framing is worth adopting: supply enough to represent plausible dynamics, but treat the library as an expression of genuine prior belief rather than a fishing expedition. Its nonlinear-oscillator run started from 97 candidate functions and selected 4.

Watch the conditioning. Near-duplicate terms make the least-squares refit ill-conditioned. Some redundancy is unavoidable — $\sin$ and $\cos$ at different wavenumbers are not orthogonal on a finite sample — but adding both a high-degree polynomial and an exponential of the same variable invites trouble.

Sample where the library is distinguishable. Terms that differ only at large amplitude cannot be told apart from data clustered near the origin. Uniform sampling over a broad range, as used throughout the thesis, distinguishes them better than a single trajectory that lingers in one region.

Arguments: applying a function to differences

This is the piece that makes interacting systems expressible, and it is easy to miss. Applying $\exp$ to individual state components gives $e^{q_1}, e^{q_2}, \dots$, which is useless for a lattice: a Toda chain interacts through $e^{-(q_{n+1} - q_n)}$, an exponential of a difference. The same is true of a point vortex ($\log\lvert q_i - q_j\rvert$) and of the $N$-body problem ($1/\lvert q_i - q_j\rvert$).

Every univariate basis therefore takes an argument selection:

SparseIdentification.Differences — Type
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
source

Differences(indices) forms all pairs $z_i - z_j$ with $i > j$, which is what an all-to-all interaction needs; Differences(indices; consecutive = true) forms only neighbouring differences, which is what a nearest-neighbour lattice needs. The indices select which components take part — for a Hamiltonian system in $z = (q, p)$ the interaction is usually among the positions alone, i.e. the first half.

# a Toda chain of three particles: interaction among the positions only
basis_functions(ExponentialBasis(Differences(1:3; consecutive = true); rates = (-1.0,)),
                Symbolics.variables(:z, 1:6))
2-element Vector{Symbolics.Num}:
 exp(z₁ - z₂)
 exp(z₂ - z₃)
# a point vortex: log of every pairwise separation
basis_functions(LogarithmicBasis(Differences(1:3)), Symbolics.variables(:z, 1:6))
3-element Vector{Symbolics.Num}:
 log(abs(-z₁ + z₂))
 log(abs(-z₁ + z₃))
 log(abs(-z₂ + z₃))
# an N-body gravitational term
basis_functions(RationalBasis(Differences(1:3)), Symbolics.variables(:z, 1:6))
3-element Vector{Symbolics.Num}:
 1 / abs(-z₁ + z₂)
 1 / abs(-z₁ + z₃)
 1 / abs(-z₂ + z₃)

Combining bases

⊕ concatenates bases, so a library is assembled from the pieces a system actually needs rather than from one monolithic keyword:

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

nterms(basis, 4)
15
SparseIdentification.hamiltonian_basis — Function
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.

source
SparseIdentification.strip_constants — Function
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.

source

For a Hamiltonian ansatz the constant term is dropped, because it contributes an identically-zero column to $J\nabla H$ and so cannot be identified. hamiltonian_basis omits it, and strip_constants removes any term whose gradient vanishes — filtering on the gradient rather than on the type of the term catches every such case, whatever basis it came from.

What is still out of reach

The bases above take scalar arguments. A genuinely three-dimensional $N$-body problem needs $1/\lVert \mathbf{q}_i - \mathbf{q}_j \rVert$ — the norm of a difference of position vectors — which needs a block structure over components that Differences does not express. In one spatial dimension the rational basis above is exactly right; in three it is not.

Where a system's Hamiltonian is not in the span of any fixed library, an evolutionary symbolic-regression search is the better tool, since it composes operators rather than selecting from a list.