Sparse identification of non-linear dynamics

Sparse identification of nonlinear dynamics (SINDy) is a data-driven algorithm for obtaining dynamical systems from data. Given a series of snapshots of a dynamical system and its corresponding time derivatives, SINDy performs a sparsity-promoting regression (such as LASSO and spare Bayesian inference) on a library of nonlinear candidate functions of the snapshots against the derivatives to find the governing equations. This procedure relies on the assumption that most physical systems only have a few dominant terms which dictate the dynamics, given an appropriately selected coordinate system and quality training data. It has been applied to identify the dynamics of fluids, based on proper orthogonal decomposition, as well as other complex dynamical systems, such as biological networks.

Mathematical Overview

First, consider a dynamical system of the form

$$\dot{\mathbf{\text{x}}}=\frac{d}{dt}\mathbf{\text{x}}(t)=\mathbf{\text{f}}(\mathbf{\text{x}}(t)),$$

where $\mathbf{\text{x}}(t)\in {ℝ}^n$ is a state vector (snapshot) of the system at time $t$ and the function $\mathbf{\text{f}}(\mathbf{\text{x}}(t))$ defines the equations of motion and constraints of the system. The time derivative may be either prescribed or numerically approximated from the snapshots.

With $\mathbf{\text{x}}$ and $\dot{\mathbf{\text{x}}}$ sampled at $m$ equidistant points in time ($t_1,t_2,\cdots ,t_m$), these can be arranged into matrices of the form

$$𝐗\mathbf{=}\left[\begin{array}{c}{𝐱}^{𝐓}\mathbf{(}{𝐭}_{\mathrm{𝟏}}\mathbf{)}\\ {𝐱}^{𝐓}\mathbf{(}{𝐭}_{\mathrm{𝟐}}\mathbf{)}\\ \mathbf{⋮}\\ {𝐱}^{𝐓}\mathbf{(}{𝐭}_{𝐦}\mathbf{)}\end{array}\right]\mathbf{=}\left[\begin{array}{cccc}{𝐱}_{\mathrm{𝟏}}\mathbf{(}{𝐭}_{\mathrm{𝟏}}\mathbf{)}& {𝐱}_{\mathrm{𝟐}}\mathbf{(}{𝐭}_{\mathrm{𝟏}}\mathbf{)}& \mathbf{\cdots }& {𝐱}_{𝐧}\mathbf{(}{𝐭}_{\mathrm{𝟏}}\mathbf{)}\\ {𝐱}_{\mathrm{𝟏}}\mathbf{(}{𝐭}_{\mathrm{𝟐}}\mathbf{)}& {𝐱}_{\mathrm{𝟐}}\mathbf{(}{𝐭}_{\mathrm{𝟐}}\mathbf{)}& \mathbf{\cdots }& {𝐱}_{𝐧}\mathbf{(}{𝐭}_{\mathrm{𝟐}}\mathbf{)}\\ \mathbf{⋮}& \mathbf{⋮}& \mathbf{\ddots }& \mathbf{⋮}\\ {𝐱}_{\mathrm{𝟏}}\mathbf{(}{𝐭}_{𝐦}\mathbf{)}& {𝐱}_{\mathrm{𝟐}}\mathbf{(}{𝐭}_{𝐦}\mathbf{)}& \mathbf{\cdots }& {𝐱}_{𝐧}\mathbf{(}{𝐭}_{𝐦}\mathbf{)}\end{array}\right],$$

and similarly for $\dot{\mathbf{\text{X}}}$.

Next, a library $\mathrm{\Theta }\mathbf{(}\mathbf{\text{X}}\mathbf{)}$ of nonlinear candidate functions of the columns of $\mathbf{\text{X}}$ is constructed, which may be constant, polynomial, or more exotic functions (like trigonometric and rational terms, and so on):

$$\mathrm{\Theta }\mathbf{(}𝐗\mathbf{)}\mathbf{=}\left[\begin{array}{cccccccc}|& |& |& |& & |& |& \\ \mathrm{𝟏}& 𝐗& {𝐗}^{\mathrm{𝟐}}& {𝐗}^{\mathrm{𝟑}}& \mathbf{\cdots }& \sin \mathbf{(}𝐗\mathbf{)}& \cos \mathbf{(}𝐗\mathbf{)}& \mathbf{\cdots }\\ |& |& |& |& & |& |& \end{array}\right]$$

The number of possible model structures from this library is combinatorically high. $\mathbf{\text{f}}(\mathbf{\text{x}}(t))$ is then substituted by $\mathrm{\Theta }\mathbf{(}\mathbf{\text{X}}\mathbf{)}$ and a vector of coefficients $\mathrm{\Xi }\mathbf{=}\left[{\mathit{\xi }}_{\mathrm{𝟏}}{\mathit{\xi }}_{\mathrm{𝟐}}\mathbf{\cdots }{\mathit{\xi }}_{𝐧}\right]$ determining the active terms in $\mathbf{\text{f}}(\mathbf{\text{x}}(t))$:

$$\dot{𝐗}=\mathrm{\Theta }\mathbf{(}𝐗\mathbf{)}\mathrm{\Xi }$$

Because only a few terms are expected to be active at each point in time, an assumption is made that $\mathbf{\text{f}}(\mathbf{\text{x}}(t))$ admits a sparse representation in $\mathrm{\Theta }\mathbf{(}\mathbf{\text{X}}\mathbf{)}$. This then becomes an optimization problem in finding a sparse $\mathrm{\Xi }$ which optimally embeds $\dot{\mathbf{\text{X}}}$. In other words, a parsimonious model is obtained by performing least squares regression on the system Template:EquationRef with sparsity-promoting ($L_1$) regularization

$${\mathit{\xi }}_{𝐤}\mathbf{=}\underset{\mathit{\xi }{\text{'}}_{𝐤}}{\arg \min }{\left||{\dot{𝐗}}_k-\mathrm{\Theta }\mathbf{(}𝐗\mathbf{)}\mathit{\xi }{\text{'}}_{𝐤}|\right|}_{\mathrm{𝟐}}\mathbf{+}\mathit{\lambda }{\left|\left|\mathit{\xi }{\text{'}}_{𝐤}\right|\right|}_{\mathrm{𝟏}},$$

where $\lambda$ is a regularization parameter. Finally, the sparse set of ${\mathit{\xi }}_{𝐤}$ can be used to reconstruct the dynamical system:

$${\dot{x}}_k=\mathrm{\Theta }\mathbf{(}𝐱\mathbf{)}{\mathit{\xi }}_{𝐤}$$