GitHub

DataDrivenSparse

DataDrivenSparse provides a universal framework to infer systems of equations using sparse regression. Assume the system:

\[y_{i} = f(x_{i}, p, t_i, u_{i})\]

We might then be able to express the unknown function $f$ as a linear combination of basis elements $\varphi_i : \mathbb R^{n_x} \times \mathbb R^{n_p} \times \mathbb R \times \mathbb R^{n_u} \mapsto \mathbb R$ .

\[y_i = \sum_{j=1}^k \xi_k ~ \varphi_k\left(x_i, p, t_i, u_i \right)\]

And simply solve the least squares problem

\[\Xi' = \min_{\Xi} \lVert Y - \Xi \varPhi \rVert_2^2\]

In the simplest case, we could use a Taylor expansion. However, if we want interpretable results, we need a key ingredient: sparsity! So, instead we aim to solve the problem

\[\Xi' = \min_{\Xi} \lVert\Xi \rVert_0 \\ \text{s.t.} \qquad \Xi \varPhi = Y\]

In its original version or via sufficient relaxation of the $L_0$ norm.

Similarly, implicit problems of the form

\[f(y_i, x_i, p, t_i, u_i) = 0\]

can be solved using an ImplicitOptimizer. Similar to the formulation above, we try to solve the corresponding optimization problem

\[\Xi' = \min_{\Xi} \lVert\Xi \rVert_0 \\ \text{s.t.} \qquad \Xi \varPhi_y = 0\]

Where the matrix of evaluated basis elements $\varPhi_y \in \mathbb R^{\lvert \varphi \rvert} \times \mathbb R^{m}$ now may also contain basis functions which are dependent on the target variables $y \in \mathbb R^{n_y}$.

Tuning parameters for sparse regression

The algorithms used by DataDrivenSparse are sensitive to the tuning of the hyperparameters! These are problem and coefficient specific, e.g., depend on the data and the unknown equations. While the examples used here are designed to work well, the used settings are not guaranteed to lead to success on other problems. Users who want to explore the space of possible hyperparameters further might be interested in using Hyperopt.jl.

Algorithms

The abstract algorithm and proximal operator entries below are developer interfaces for extending DataDrivenSparse. Application code should generally use the concrete algorithms and operators.

DataDrivenSparse.AbstractSparseRegressionAlgorithm — Type
AbstractSparseRegressionAlgorithm

Abstract interface for sparse-regression algorithms used by DataDrivenSparse.

Concrete subtypes are callable algorithm objects that solve a matrix regression problem for each target row. They are used directly as solve(prob, basis, alg) optimizers and through SparseLinearSolver.

Interface

A subtype Alg <: AbstractSparseRegressionAlgorithm must provide:

  • get_thresholds(alg): returns the scalar threshold or iterable threshold schedule evaluated by SparseLinearSolver.
  • alg(X, Y; options = DataDrivenCommonOptions(), kwargs...): returns (coefficients, optimal_thresholds, optimal_iterations) for feature matrix X and target matrix Y.

Algorithms that use SparseLinearSolver should also implement:

  • init_cache(alg, A::AbstractMatrix, B::AbstractMatrix): constructs the per-target cache for the design matrix A and target row matrix B.
  • step!(cache, threshold): performs one thresholded solver update.

Arguments

  • X::AbstractMatrix: feature or basis-evaluation matrix.
  • Y::AbstractMatrix: target matrix whose rows are fit independently.

Keywords

  • options::DataDrivenCommonOptions: convergence tolerances, iteration limits, and verbosity used by iterative sparse-regression solvers.
  • kwargs...: algorithm-specific options. Generic callers should not depend on any keyword that is not documented by the concrete algorithm.

Returns

The generic sparse-regression call returns a tuple (coefficients, optimal_thresholds, optimal_iterations). coefficients is a matrix with one row per target in Y; the other entries record the selected threshold and iteration count for each target row.

Examples

using DataDrivenSparse

X = [1.0 2.0 3.0; 1.0 4.0 9.0]
Y = [2.0 4.0 6.0]
alg = STLSQ([0.1])

coefficients, thresholds, iterations = alg(X, Y)
source
DataDrivenSparse.AbstractSparseRegressionCache — Type
AbstractSparseRegressionCache

Developer interface for the mutable cache used by SparseLinearSolver. This type is intended for packages implementing a new AbstractSparseRegressionAlgorithm, not for constructing user results directly.

Interface

A cache subtype must provide mutable fields Ã, B̃, X, X_prev, and active_set. Ã is the feature matrix, B̃ is the target vector or matrix, X is the current coefficient array, X_prev is the previous iterate, and active_set has the same shape as X. The generic solver calls step!(cache, λ), copies a winning cache with _set!, and checks convergence with the abstol and reltol values from DataDrivenCommonOptions.

The cache must also support the StatsAPI methods coef, rss, dof, and nobs; the default methods supplied here use the fields above. A custom cache should preserve the coefficient shape and return a numeric residual from rss.

source
DataDrivenSparse.get_thresholds — Function
get_thresholds(alg::AbstractSparseRegressionAlgorithm)

Return the scalar threshold or ordered threshold schedule explored by SparseLinearSolver. A custom algorithm must return either a scalar or an iterable that supports minimum and iteration.

source
DataDrivenSparse.get_relaxation — Function
get_relaxation(alg::AbstractSparseRegressionAlgorithm)

Return an optional relaxation parameter used by an algorithm. The default is nothing; algorithms that expose relaxation should specialize this method and document how it changes their update rule.

source
DataDrivenSparse.step! — Function
step!(cache, lambda)

Perform one thresholded update of a sparse-regression cache in place. The implementation must update cache.X, cache.X_prev, and cache.active_set consistently and return the cache or nothing.

source
DataDrivenSparse.STLSQ — Type
struct STLSQ{T<:Union{Number, AbstractVector}, R<:Number} <: DataDrivenSparse.AbstractSparseRegressionAlgorithm

STLSQ is taken from the original paper on SINDY and implements a sequentially thresholded least squares iteration. λ is the threshold of the iteration. It is based upon this Matlab implementation. It solves the following problem

\[\argmin_{x} \frac{1}{2} \| Ax-b\|_2 + \rho \|x\|_2\]

with the additional constraint

\[\lvert x_i \rvert > \lambda\]

If the parameter ρ > 0, ridge regression will be performed using the normal equations of the corresponding regression problem.

Arguments

  • threshold: positive scalar or iterable of positive thresholds to explore.
  • rho: nonnegative ridge-regression coefficient.

Returns

Return an algorithm object callable as alg(X, Y; options), producing coefficient matrices, selected thresholds, and iteration counts.

Fields

  • thresholds: Sparsity threshold

  • rho: Ridge regression parameter

Example

opt = STLSQ()
opt = STLSQ(1e-1)
opt = STLSQ(1e-1, 1.0) # Set rho to 1.0
opt = STLSQ(Float32[1e-2; 1e-1])

Note

This was formally STRRidge and has been renamed.

source
DataDrivenSparse.ADMM — Type
mutable struct ADMM{T, R<:Number} <: DataDrivenSparse.AbstractSparseRegressionAlgorithm

ADMM is an implementation of Lasso using the alternating direction methods of multipliers, and loosely based on this implementation. It solves the following problem

\[\argmin_{x} \frac{1}{2} \| Ax-b\|_2 + \lambda \|x\|_1\]

Fields

  • thresholds: Sparsity threshold parameter

  • rho: Augmented Lagrangian parameter

Arguments

  • threshold: positive scalar or iterable of positive sparsity thresholds.
  • ρ: positive augmented-Lagrangian parameter.

Returns

Return an algorithm object callable as alg(X, Y; options), producing coefficient matrices, selected thresholds, and iteration counts.

Example

opt = ADMM()
opt = ADMM(1e-1, 2.0)
source
DataDrivenSparse.SR3 — Type
mutable struct SR3{T, V, P<:DataDrivenSparse.AbstractProximalOperator} <: DataDrivenSparse.AbstractSparseRegressionAlgorithm

SR3 is an optimizer framework introduced by Zheng et al., 2018 and used within Champion et al., 2019. SR3 contains a sparsification parameter λ, a relaxation ν. It solves the following problem

\[\argmin_{x, w} \frac{1}{2} \| Ax-b\|_2 + \lambda R(w) + \frac{\nu}{2}\|x-w\|_2\]

Where R is a proximal operator, and the result is given by w.

Arguments

  • threshold: positive sparsity threshold or threshold schedule.
  • nu: positive relaxation parameter.
  • R: AbstractProximalOperator used for the relaxed update.

Returns

Return an algorithm object callable as alg(X, Y; options), producing coefficient matrices, selected thresholds, and iteration counts.

Fields

  • thresholds: Sparsity threshold

  • nu: Relaxation parameter

  • proximal: Proximal operator

Example

opt = SR3()
opt = SR3(1e-2)
opt = SR3(1e-3, 1.0)
opt = SR3(1e-3, 1.0, SoftThreshold())

Note

Opposed to the original formulation, we use nu as a relaxation parameter, as given in Champion et al., 2019. In the standard case of hard thresholding the sparsity is interpreted as λ = threshold^2 / 2, otherwise λ = threshold.

source
DataDrivenSparse.WyNDA — Type
struct WyNDA{T<:Number, C, IC} <: DataDrivenSparse.AbstractSparseRegressionAlgorithm

WyNDA implements the Wide-Array of Nonlinear Dynamics Approximation update as an online recursive least-squares estimator with exponential forgetting. Given a matrix of evaluated basis functions X and targets Y, it updates the coefficient matrix sample-by-sample so that Y ≈ coefficients * X.

The forgetting factor λ controls how quickly older samples are discounted. Values close to one recover a batch least-squares-like fit, while smaller values adapt faster to parameter drift.

Arguments

  • λ: forgetting factor satisfying 0 < λ <= 1.

Keywords

  • initial_covariance: positive scalar or square matrix used for the initial inverse covariance.
  • initial_coefficients: optional initial coefficient vector or target-by-feature matrix.

Returns

Return an online algorithm object callable as alg(X, Y; options), producing the coefficient matrix, the forgetting factor, and the number of observations.

Fields

  • λ: Exponential forgetting factor.

  • initial_covariance: Initial inverse covariance scale or matrix.

  • initial_coefficients: Optional initial coefficient vector or matrix.

Example

opt = WyNDA()
opt = WyNDA(0.998)
opt = WyNDA(0.998; initial_covariance = 1.0e4)
source
DataDrivenSparse.ImplicitOptimizer — Type
mutable struct ImplicitOptimizer{T<:DataDrivenSparse.AbstractSparseRegressionAlgorithm} <: DataDrivenSparse.AbstractSparseRegressionAlgorithm

Optimizer for finding a sparse implicit relationship via alternating the left-hand side of the problem and solving the explicit problem, as introduced here.

\[\argmin_{x} \|x\|_0 ~s.t.~Ax= 0\]

Arguments

Keywords

  • options::DataDrivenCommonOptions: convergence and selection settings.
  • necessary_idx: Boolean mask identifying coefficients that must participate in each candidate implicit relation.

Returns

Return an implicit sparse-regression algorithm object. Calling it returns the best cache, threshold, and iteration count for each candidate left-hand side.

Fields

  • optimizer: Explicit Optimizer

Example

ImplicitOptimizer(STLSQ())
ImplicitOptimizer(0.1f0, ADMM)
source
DataDrivenSparse.SparseLinearSolver — Type
struct SparseLinearSolver{A<:DataDrivenSparse.AbstractSparseRegressionAlgorithm, T<:Number}

Sparse regression solver that applies an AbstractSparseRegressionAlgorithm to one or more target variables.

Arguments

  • algorithm::AbstractSparseRegressionAlgorithm: algorithm used for each target.

Keywords

  • options::DataDrivenCommonOptions: tolerances, iteration limits, selector, and progress settings copied into the solver.

Returns

Return a solver object callable as solver(X, Y), where X has features in rows and Y has target variables in rows. The call returns one cache, selected threshold, and iteration count per target.

Fields

  • algorithm::DataDrivenSparse.AbstractSparseRegressionAlgorithm: Sparse-regression algorithm applied to each target.

  • abstol::Number: Absolute convergence tolerance for cache updates.

  • reltol::Number: Relative convergence tolerance for cache updates.

  • maxiters::Int64: Maximum number of iterations over all thresholds.

  • verbose::Bool: Whether progress information is printed.

  • progress::Bool: Whether the underlying algorithm reports progress.

  • selector::Function: Function used to select the best cache.

source

Proximal Operators

Custom proximal operators should subtype AbstractProximalOperator and implement the documented callable and active-set methods.

DataDrivenSparse.AbstractProximalOperator — Type
AbstractProximalOperator

Developer interface for thresholding operators used by sparse-regression algorithms such as SR3.

Interface

A subtype must implement operator(x, λ) as a callable object that updates x in place, operator(y, x, λ) as an out-of-place-buffer form, and active_set!(mask, operator, x, λ) to identify the nonzero coefficients. The two callable forms must preserve the shape and element type of the coefficient array. Concrete operators may store additional thresholds, but those fields and their defaults must be documented. The in-place form returns the modified x; the buffer form writes y and returns it. active_set! returns the modified mask.

source
DataDrivenSparse.active_set! — Function
active_set!(mask, operator, x, lambda)

Update the Boolean active-set mask for a sparse-regression proximal operator.

This is a developer extension point for AbstractProximalOperator. mask and x must have the same shape, and an active entry indicates that the corresponding coefficient survives thresholding.

Arguments

  • mask: Boolean array with the same shape as x.
  • operator::AbstractProximalOperator: thresholding operator.
  • x: coefficient array inspected by the operator.
  • lambda: nonnegative threshold parameter.

Returns

Return the modified mask.

source
DataDrivenSparse.SoftThreshold — Type
struct SoftThreshold <: DataDrivenSparse.AbstractProximalOperator

Proximal operator, which implements the soft thresholding operator.

sign(x) * max(abs(x) - λ, 0)

Arguments

  • x: coefficient array updated in place, or the source array for the buffered form.
  • λ: nonnegative threshold.

Returns

Return the modified array. The three-argument form writes the result into its first array argument.

See by Zheng et al., 2018.

source
DataDrivenSparse.HardThreshold — Type
struct HardThreshold <: DataDrivenSparse.AbstractProximalOperator

Proximal operator, which implements the hard thresholding operator.

abs(x) > sqrt(2*λ) ? x : 0

Arguments

  • x: coefficient array updated in place, or the source array for the buffered form.
  • λ: nonnegative threshold.

Returns

Return the modified array. The three-argument form writes the result into its first array argument.

See by Zheng et al., 2018.

source
DataDrivenSparse.ClippedAbsoluteDeviation — Type
struct ClippedAbsoluteDeviation{T} <: DataDrivenSparse.AbstractProximalOperator

Proximal operator, which implements the (smoothly) clipped absolute deviation operator.

abs(x) > ρ ? x : sign(x) * max(abs(x) - λ, 0)

Where ρ = 5λ per default.

Arguments

  • x: coefficient array updated in place, or the source array for the buffered form.
  • λ: nonnegative soft-threshold parameter.

Fields

  • ρ: optional hard cutoff; NaN selects the default 5λ cutoff.

Returns

Return the modified array. The three-argument form writes the result into its first array argument.

Fields

  • ρ: Upper threshold

Example

opt = ClippedAbsoluteDeviation()
opt = ClippedAbsoluteDeviation(1e-1)

See by Zheng et al., 2018.

source