The problem

A large chunk of applied mathematics aims to understand dynamical systems and to characterize their behavior.

But what happens if you have a physical system for which you don’t know the governing equations?

You might suspect that there is a set of differential equations pulling the strings, but maybe you are unable or unwilling to derive them analytically.

If you can get your hands on measurements of the system, then you may be able to use a data-driven approach.

A solution

One recently developed method for nonlinear model discovery is the Sparse Identification of Nonlinear Dynamical systems (SINDy).

How it works

The main idea is that the right-hand sides of many dynamical systems of interest do not include many terms, implying that they are sparse with respect to an appropriately chosen basis.

Suppose we have measurements $x(t)$ of a physical system over time. SINDy seeks to represent the time evolution of $x(t)$ in terms of a nonlinear function $f$, defining a dynamical system model for $x(t)$: $$\frac{d}{dt}x(t)=f(x(t)).$$

The vector $x(t)=[x_1(t),x_2(t),\dots x_n(t){]}^{\mathrm{\top }}$ gives the state of the physical system at time $t$. The function $f(x(t))$ constrains how the system evolves in time.

In many problems of interest, $f$ turns out to be sparse in the set of possible basis functions. For example, the function

$$\frac{d}{dt}x=f(x)=\left(\begin{array}{c}{f}_1(x)\\ f_2(x)\end{array}\right)=\left(\begin{array}{c}1-x_1+3x_1{x}_2\\ x_2^2-5x_1^3\end{array}\right)$$

is sparse with respect to the set of polynomials in two variables in the sense that just a few of the infinite number of polynomial terms are present on the right-hand side of the equation. So if we pick an appropriate basis ${\theta }_1,{\theta }_2,\dots ,{\theta }_k$ then it is often possible to express each component of $f$ (note that $f$ is often vector-valued) as

$$f_i(x)\approx {\xi }_1{\theta }_1(x)+{\xi }_2{\theta }_2(x)+\dots {\xi }_k{\theta }_k(x)$$

with most coefficients ${\xi }_j$ being zero. SINDy leverages sparse regression to identify the active terms in the dynamics (the nonzero coefficients ${\xi }_j$).

Mathematical formulation

To apply SINDy one needs a set of measurement data collected at times $t_1,t_2,\dots ,t_n$, and the time derivatives of these measurements (either measured directly or numerically approximated). These data are aggregated into the matrices $X$ and $\dot{X}$, respectively

$$X=\left[\begin{array}{cccc}{x}_1(t_1)& x_2(t_1)& \dots & x_n(t_1)\\ x_1(t_2)& x_2(t_2)& \dots & x_n(t_2)\\ ⋮& ⋮& & ⋮\\ x_1(t_m)& x_2(t_m)& \dots & x_n(t_m)\end{array}\right],$$ $$\dot{X}=\left[\begin{array}{cccc}\dot{x_1}(t_1)& \dot{x_2}(t_1)& \dots & \dot{x_n}(t_1)\\ \dot{x_1}(t_2)& \dot{x_2}(t_2)& \dots & \dot{x_n}(t_2)\\ ⋮& ⋮& & ⋮\\ \dot{x_1}(t_m)& \dot{x_2}(t_m)& \dots & \dot{x_n}(t_m)\end{array}\right].$$

Next, one forms a library matrix $\mathrm{\Theta }(X)$ whose columns consist of a chosen set of basis functions applied to the data

$$\mathrm{\Theta }(X)=\left[\begin{array}{cccc}\mid & \mid & & \mid \\ {\theta }_1(X)& {\theta }_2(X)& \dots & {\theta }_{\ell }(X)\\ \mid & \mid & & \mid \end{array}\right].$$

We seek a set of sparse coefficient vectors (collected into a matrix)

$$\mathrm{\Xi }=\left[\begin{array}{cccc}\mid & \mid & & \mid \\ {\xi }_1& {\xi }_2& \dots & {\xi }_n\\ \mid & \mid & & \mid \end{array}\right].$$

The vector ${\xi }_i$ gives the coefficients for a linear combination of basis functions ${\theta }_1(x),{\theta }_2(x),\dots ,{\theta }_{\ell }(x)$ representing $f_i$, the $i$th component function of $f$, i.e. $f_i(x)\approx \mathrm{\Theta }\left(x^{\mathrm{\top }}\right){\xi }_i$ (note that the preceding equation should be interpreted as a symbolic expression).

Finally, we are ready to write down the approximation problem underlying SINDy:

$$\dot{X}\approx \mathrm{\Theta }(X)\mathrm{\Xi }.$$

A sparse regression algorithm is then employed to solve the related minimization problem

$$\underset{\mathrm{\Xi }}{min}\frac{1}{2}\Vert \dot{X}-\mathrm{\Theta }(X)\mathrm{\Xi }{\Vert }_F^2+R(\mathrm{\Xi }),$$ where $R$ is a sparsity-promoting regularization function such as $\Vert \cdot {\Vert }_1$.

The following image shows an overview of the SINDy method applied to the Lorenz equations.

Resources

• The original paper introducing SINDy: Discovering governing equations from data by sparse identification of nonlinear dynamical systems
• A recent paper of mine applying SINDy to a noisy real-world data set: Discovery of physics from data: universal laws and discrepancies
• PySINDy: a Python package for SINDy