diff --git a/docs/images/museum/sindy/C2D.png b/docs/images/museum/sindy/C2D.png new file mode 100644 index 00000000..e1dc2395 Binary files /dev/null and b/docs/images/museum/sindy/C2D.png differ diff --git a/docs/images/museum/sindy/L2D.png b/docs/images/museum/sindy/L2D.png new file mode 100644 index 00000000..6fb55a8b Binary files /dev/null and b/docs/images/museum/sindy/L2D.png differ diff --git a/docs/images/museum/sindy/L3D.png b/docs/images/museum/sindy/L3D.png new file mode 100644 index 00000000..a2efd832 Binary files /dev/null and b/docs/images/museum/sindy/L3D.png differ diff --git a/docs/images/museum/sindy/LSQ.png b/docs/images/museum/sindy/LSQ.png new file mode 100644 index 00000000..2ee1636d Binary files /dev/null and b/docs/images/museum/sindy/LSQ.png differ diff --git a/docs/images/museum/sindy/Lorenz3D.png b/docs/images/museum/sindy/Lorenz3D.png new file mode 100644 index 00000000..84d2b4dc Binary files /dev/null and b/docs/images/museum/sindy/Lorenz3D.png differ diff --git a/docs/images/museum/sindy/Masking.png b/docs/images/museum/sindy/Masking.png new file mode 100644 index 00000000..bb132a3c Binary files /dev/null and b/docs/images/museum/sindy/Masking.png differ diff --git a/docs/images/museum/sindy/O3D.png b/docs/images/museum/sindy/O3D.png new file mode 100644 index 00000000..65057679 Binary files /dev/null and b/docs/images/museum/sindy/O3D.png differ diff --git a/docs/images/museum/sindy/P1.png b/docs/images/museum/sindy/P1.png new file mode 100644 index 00000000..a3a55cc1 Binary files /dev/null and b/docs/images/museum/sindy/P1.png differ diff --git a/docs/images/museum/sindy/P2.png b/docs/images/museum/sindy/P2.png new file mode 100644 index 00000000..91737928 Binary files /dev/null and b/docs/images/museum/sindy/P2.png differ diff --git a/docs/images/museum/sindy/SR.png b/docs/images/museum/sindy/SR.png new file mode 100644 index 00000000..1fc84667 Binary files /dev/null and b/docs/images/museum/sindy/SR.png differ diff --git a/docs/images/museum/sindy/SRin.png b/docs/images/museum/sindy/SRin.png new file mode 100644 index 00000000..977edf7b Binary files /dev/null and b/docs/images/museum/sindy/SRin.png differ diff --git a/docs/images/museum/sindy/STLSQ.png b/docs/images/museum/sindy/STLSQ.png new file mode 100644 index 00000000..0103d240 Binary files /dev/null and b/docs/images/museum/sindy/STLSQ.png differ diff --git a/docs/images/museum/sindy/Theta.png b/docs/images/museum/sindy/Theta.png new file mode 100644 index 00000000..1a6e1ec5 Binary files /dev/null and b/docs/images/museum/sindy/Theta.png differ diff --git a/docs/images/museum/sindy/Thresholding.png b/docs/images/museum/sindy/Thresholding.png new file mode 100644 index 00000000..b038446e Binary files /dev/null and b/docs/images/museum/sindy/Thresholding.png differ diff --git a/docs/images/museum/sindy/X_.png b/docs/images/museum/sindy/X_.png new file mode 100644 index 00000000..868f66ec Binary files /dev/null and b/docs/images/museum/sindy/X_.png differ diff --git a/docs/images/museum/sindy/X_dX_.png b/docs/images/museum/sindy/X_dX_.png new file mode 100644 index 00000000..715735c1 Binary files /dev/null and b/docs/images/museum/sindy/X_dX_.png differ diff --git a/docs/images/museum/sindy/Xtheta.png b/docs/images/museum/sindy/Xtheta.png new file mode 100644 index 00000000..98f147d2 Binary files /dev/null and b/docs/images/museum/sindy/Xtheta.png differ diff --git a/docs/images/museum/sindy/dX_.png b/docs/images/museum/sindy/dX_.png new file mode 100644 index 00000000..c1c5cb9c Binary files /dev/null and b/docs/images/museum/sindy/dX_.png differ diff --git a/docs/images/museum/sindy/d_xyz.png b/docs/images/museum/sindy/d_xyz.png new file mode 100644 index 00000000..5bd04836 Binary files /dev/null and b/docs/images/museum/sindy/d_xyz.png differ diff --git a/docs/images/museum/sindy/dx.png b/docs/images/museum/sindy/dx.png new file mode 100644 index 00000000..8991df0a Binary files /dev/null and b/docs/images/museum/sindy/dx.png differ diff --git a/docs/images/museum/sindy/dy.png b/docs/images/museum/sindy/dy.png new file mode 100644 index 00000000..2716a128 Binary files /dev/null and b/docs/images/museum/sindy/dy.png differ diff --git a/docs/images/museum/sindy/dz.png b/docs/images/museum/sindy/dz.png new file mode 100644 index 00000000..e9e65576 Binary files /dev/null and b/docs/images/museum/sindy/dz.png differ diff --git a/docs/images/museum/sindy/flow.png b/docs/images/museum/sindy/flow.png new file mode 100644 index 00000000..b32dcd20 Binary files /dev/null and b/docs/images/museum/sindy/flow.png differ diff --git a/docs/images/museum/sindy/flow_SR.png b/docs/images/museum/sindy/flow_SR.png new file mode 100644 index 00000000..ed9375dd Binary files /dev/null and b/docs/images/museum/sindy/flow_SR.png differ diff --git a/docs/images/museum/sindy/flow_full.jpg b/docs/images/museum/sindy/flow_full.jpg new file mode 100644 index 00000000..234b97c0 Binary files /dev/null and b/docs/images/museum/sindy/flow_full.jpg differ diff --git a/docs/images/museum/sindy/iterin.png b/docs/images/museum/sindy/iterin.png new file mode 100644 index 00000000..e182ef5d Binary files /dev/null and b/docs/images/museum/sindy/iterin.png differ diff --git a/docs/images/museum/sindy/xdx.png b/docs/images/museum/sindy/xdx.png new file mode 100644 index 00000000..3eef7e40 Binary files /dev/null and b/docs/images/museum/sindy/xdx.png differ diff --git a/docs/museum/sindy.md b/docs/museum/sindy.md index a00af2a6..c07ff654 100644 --- a/docs/museum/sindy.md +++ b/docs/museum/sindy.md @@ -1,82 +1,546 @@ -# Sparse Identification of Non-linear Dynamical Systems (SINDy)[1] + + -In this section, we teach, create, simulate, and visualize SINDy model implemented in NGC-Learn library. After going through this demonstration, you will: +# Sparse Identification of Non-linear Dynamical Systems (SINDy)[1] -1. Learn how to build a SINDy model of time-series dataset, generated using Ordinary Differential Equations (ODE) of known dynamical systems used in [1]. -2. Learn how to build polynomial libraries of given dataset with arbitrary order. -3. Learn how to solve the sparse regression problem by iteratively performing the least squares (LSQ) method followed by thresholding-- Sequential Thresholding Least Square (STLSQ)-- for the given model. +In this section, we teach, create, simulate, and visualize the Sparse Identification of Non-linear Dynamical Systems (SINDy) [1] model implemention in NGC-Learn library (JAX). After going through this demonstration, you will: +1. Learn how to discover the differential equation of a dynamical system using SINDy algorithm only by the system's stapshots. +2. Learn how to build polynomial libraries with arbitrary order out of the dataset. +3. Learn how to solve the sparse regression problem in 2 ways + - Iteratively finding the coefficient matrix by gradient descent. + - Iteratively performing the least squares (LSQ) method followed by thresholding-- Sequential Thresholding Least Square (STLSQ) for the given model. + + The model **code** for this exhibit can be found [here](https://github.com/NACLab/ngc-museum/exhibits/sindy/sindy.py). ## SINDy -SINDy is a data-driven algorithm that discovers the differential equation governing the dynamical systems. It uses symbolic regression to identify differential equation of the system and it solves sparse regression over the pre-defined library of candidate terms. It takes time series gathered dataset of the system and it gives you its describing differential equation. +SINDy is a data-driven algorithm that discovers the governing behavior of a dynamical system in terms of symbolic differential equation. It solves the sparse regression problem over the coefficients of pre-defined library that includes $p$ candidate predictors. It tries to find sparse model that only uses $s$ predictors out of $p$ where $s - -

+### SINDy Dynamics +If $\mathbf{X}$ is a system that only depends on variable $t$, a very small change in the independant variable ($dt$) can cause changing the system by $dX$ amount. +```math +d\mathbf{X} = \mathbf{Ẋ}(t)~dt +``` +SINDy models the derivative[^1] (a linear operation) as a linear transformations with: +[^1]: Derivative is a linear operation that acts on dt and gives the differential that is the linearization approximation of the taylor series of the function. +```math +\frac{d\mathbf{X}(t)}{dt} = \mathbf{Ẋ}(t) = \mathbf{f}(\mathbf{X}(t)) +``` +SINDy assumes thatt this linear operation, $\mathbf{f}(\mathbf{X}(t))$ is a matrix multiplication that linearly combines the relevant predictors to describe the system's equation. +```math +\mathbf{f}(\mathbf{X}(t)) = \mathbf{\Theta}(\mathbf{X})~\mathbf{W} +``` -### Inputs -> Time: $ts = [t_0, t_1, \dots, T]$ -> State matrix: $\mathbf{X}_{(m \times n)}$ (t measurements of n variables) -### Inputs +Given a group of candidate functions in the library $\mathbf{\Theta}(\mathbf{X})$, the coefficient $\mathbf{W}$ of choose the library terms is **sparse**. In other words, there are only a few functions that exist in the system's differential equation. Given these assumptions, SINDy solves the sparse regression problem to find the $\mathbf{W}$ that maps the library selected terms to each feature of the system. SINDy imposes parsimony constraints over symbolic regression (i.e., genetic programming) to describe a dynamical system's behavior by as few terms as possible. In order to select a sparse set of the given features, it adds the LASSO regularizarion (i.e., L1 norm) to the regression problem and solves the sparse regression or solves the regression problem by STLSQ. Here we desceibe STLSQ in third step of the SINDy dynamics. - **Time** -* $ts = [t_0,~t_1, \dots,~T]$ - - **State matrix** -* $\mathbf{X}(t)_{(m \times n)} = [x(t),~~y(t),~~z(t)]$ +SINDy's dynamics can be presented in 3 main phases according to the figure 1. +------------------------------------------------------------------------------------------ + +

+ +**Figure 1:** **Flow of three phases in SINDy.** **Phase-1)** Data collection: capturing system's states that are changing in time and making the state vector. **Phase-2A)** Library formation: manually creating the library of candidate predictors that could appear in the model. **Phase-2B)** Derivative computation: using the data collected in phase 1 and compute its derivative with respect to time. **Phase-3)** Solving the sparse regression problem. +

+ +------------------------------------------------------------------------------------------ + + + + + + + + + + +
+ +## Phase 1: Collecting Dataset → $\mathbf{X}_{(m \times n)}$ +This phase involves gathering the raw data points representing the system's states across time; In this example, capturing the x, y, and z coordinates of the system's states in this. Here, $m$ represents the number of data points (number of the spanshots/time length) and $n$ is the system's dimensions. + +

+ Dataset collection showing x, y, z coordinates + +

+
+ + + + + + + + + + + + + + + + + + + + + + -### Data Matrix +
+ +## Phase 2: Processing + +

+ +

+
+ +### 2.A: Making Library → $\mathbf{\Theta}_{(m \times p)}$ +In this step, using the dataset collected in step 1, given the pre-defined function terms, we construct the dictionary of candidate predictors for system's differential equations. These functions form the columns of our library matrix $\mathbf{\Theta}(\mathbf{X})$ and $p$ is the number of candidate predictors. To identify the dynamical structure of the system this library of candidate functions appear in the regression problem to propose the model's structure that later the coefficient matrix will give weight to them according to the problem setup. Assuming sparse models for the system, by sparsification (LASSO or thresholding weigths) decide which structure best describe the system's behavior using predictors. Given a set of time-series measurements of a dynamical system state variables ($\mathbf{X}_{(m \times n)}$) we construct: +Library of Candidate Functions: $\Theta(\mathbf{X}) = [\mathbf{1} \quad \mathbf{X} \quad \mathbf{X}^2 \quad \mathbf{X}^3 \quad \sin(\mathbf{X}) \quad \cos(\mathbf{X}) \quad ...]$ + +

+ +

+
+ +### 2.B: Compute State Derivatives → $\mathbf{Ẋ}_{(m \times n)}$ +Given a set of time-series measurements of a dynamical system state variables $\mathbf{X}_{(m \times n)}$ we construct the derivative matrix: $\mathbf{Ẋ}_{(m \times n)}$ (computed numerically) +In this step, using the dataset collected in step 1, we calculating the time derivatives of each state variable with respect to time. In this example, we compute ẋ, ẏ, and ż to capture how the system evolves over time. + +

+ +

+
+ -Given a set of time-series measurements of a dynamical system state variables ($\mathbf{X}_{(m \times n)}$) we construct: -Derivative matrix: $\dot{\mathbf{X}}_{(m \times n)}$ (computed numerically) + -Library of Candidate Functions: $\Theta(\mathbf{X}) = [\mathbf{1} \quad \mathbf{X} \quad \mathbf{X}^2 \quad \mathbf{X}^3 \quad \sin(\mathbf{X}) \quad \cos(\mathbf{X}) \quad ...]$ + + + + ------------------- + -

- - -

+ + + +
+ +## Phase 3: Solving Sparse Regression Problem → $\mathbf{W_s}_{(p \times n)}$ +Solving the Sparse Regression problem (SR) can be done with various method such as Lasso, STLSQ, Elastic Net, and many others. Here we describe STLSQ to solve the SR problem according to the SINDy method. + +

+ Dataset collection showing x, y, z coordincates +

+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + +### Solving Sparse Regression by Sequential Thresholding Least Square (STLSQ) +

- - - + + +**Figure 1:** **Flow of three phases in SINDy.** **Phase-1)** Data collection: capturing system's states that are changing in time and making the state vector. **Phase-2A)** Library formation: manually creating the library of candidate predictors that could appear in the model. **Phase-2B)** Derivative computation: using the data collected in phase 1 and compute its derivative with respect to time. **Phase-3)** Solving the sparse regression problem with STLSQ.

+------------------------------------------------------------------------------------------ +
+ +### Sequential Thresholding Least Square (STLSQ) +
+

+ State derivatives visualization +

+
+#### 3.A: Least Square method (LSQ) → $\mathbf{W}$ +Finds library coefficients by solving the following regression problem $\mathbf{Ẋ} = \mathbf{\Theta}\mathbf{W}$ analytically $\mathbf{W} = (\mathbf{\Theta}^T \mathbf{\Theta})^{-1} \mathbf{\Theta}^T \mathbf{Ẋ}$ + +

+ State derivatives visualization +

+
+ +#### 3.B: Thresholding → $\mathbf{W_s}$ +Sparsifies $\mathbf{W}$ by keeping only some terms in $\mathbf{W}$ that corresponds to the effective terms in the library. + +

+ State derivatives visualization +

+
+ +#### 3.C: Masking → $\mathbf{\Theta_s}$ +Sparsifies $\mathbf{\Theta}$ by keeping only the corresponding terms in $\mathbf{W}$ that are kept. + +

+ State derivatives visualization +

+
+ +#### 3.D: Repeat A → B → C until convergence +Solving LSQ with the sparse matrix $\mathbf{\Theta_s}$ and $\mathbf{W_s}$ and find the new $\mathbf{W}$ and repreat steps B and C everytime. + +

+ State derivatives visualization +

+
+ + + + + + + + + + + + +## Code + + + +```python + + + +import numpy as np +import jax.numpy as jnp +from ngclearn.utils.feature_dictionaries.polynomialLibrary import PolynomialLibrary +from ngclearn.utils.diffeq.ode_solver import solve_ode +from ngclearn.utils.diffeq.odes import lorenz, linear_2D + +jnp.set_printoptions(suppress=True, precision=5) + + +## system's ode function +dfx = lorenz + +x0 = jnp.array([-8, 7, 27], dtype=jnp.float32) ## initial state + +t0 = 0. ## starting time +dt = 1e-2 ## time steps +T = 2000 ## #of steps + +deg = 2 ## polynomial library degree +include_bias = False ## if include bias in making poly library +threshold = 0.02 ## sparaity threshold +max_iter=10 + +## Phase 1: Collecting Dataset (solving ode) +ts, X = solve_ode('rk4', t0, x0, T=T, dfx=dfx, dt=dt, params=None, sols_only=True) + +## Phase 2.A: Making Library (polynomial library) +lib_creator = PolynomialLibrary(poly_order=deg, include_bias=include_bias) +feature_lib, feature_names = lib_creator.fit([X[:, i] for i in range(X.shape[1])]) + +## Phase 2.B: Compute State Derivatives +dX = jnp.array(np.gradient(X, ts.ravel(), axis=0)) + +########## Solving Sparse Regression (for each dimension) ########## +#~~~~~~~~~~~~ By Seqyential Thresholding Least Square ~~~~~~~~~~~~~ + +for dim in range(dX.shape[1]): + ## 3.A: 'Initial' Least Square + coef = jnp.linalg.lstsq(feature_lib, dX[:, dim][:, None], rcond=None)[0] + + for i in range(max_iter): + coef_pre = jnp.array(coef) + coef_zero = jnp.zeros_like(coef) + + ## 3.B: thresholding + res_idx = jnp.where(jnp.abs(coef) >= threshold, + True, + False) + ## 3.C: masking + res_mask = jnp.any(res_idx, axis=1) ## residual mask + res_lib = feature_lib[:, res_mask] ## residual predictors + + ## 3.A: Least Square + coef_new = jnp.linalg.lstsq(res_lib, dX[:, dim][:, None], + rcond=None + )[0] ## least square + + coef = coef_zero.at[res_mask].set(coef_new) + + ## 3.B: 'Final' thresholding + coeff = jnp.where(jnp.abs(coef) >= threshold, coef, 0.) + + print(f"coefficients for dimension {dim+1}: \n", coeff.T) + + + +``` + + + + + +## Results + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ Model + + Results +
+ + ## Oscillator + +True model's equation \ +$\mathbf{ẋ} = \mu_1\mathbf{x} + \sigma \mathbf{xy}$ \ +$\mathbf{ẏ} = \mu_2\mathbf{y} + (\omega + \alpha \mathbf{y} + \beta \mathbf{z})\mathbf{z} - \sigma \mathbf{x}^2$ \ +$\mathbf{ż} = \mu_2\mathbf{z} - (\omega + \alpha \mathbf{y} + \beta \mathbf{z})\mathbf{y}$ + +```python +--- SINDy results ---- +ẋ = 0.050 𝑥 + 1.100 𝑥𝑦 +ẏ = 2.999 𝓏 -4.999 𝓏^2 + -0.010 𝑦 -1.998 𝑦𝓏 -1.100 𝑥^2 +ż = -0.010 𝓏 -3.000 𝑦 + + 5.000 𝑦𝓏 + 1.999 𝑦^2 + + [1, 𝓏, 𝓏^2, 𝑦, 𝑦𝓏, 𝑦^2, 𝑥, 𝑥𝓏, 𝑥𝑦, 𝑥^2] +[[ 0. 0. 0. 0. 0. 0. 0.049 0. 1.099 0.] + [ 0. 2.99 -4.99 -0.010 -1.99 0. 0. 0. 0. -1.099] + [ 0. -0.009 0. -2.99 4.99 1.99 0. 0. 0. 0.]] +``` + + +

+ + +

+
+ + ## Lorenz + +True model's equation \ +$\mathbf{ẋ} = 10(\mathbf{y} - \mathbf{x})$ \ +$\mathbf{ẏ} = \mathbf{x}(28 - \mathbf{z}) - \mathbf{y}$ \ +$\mathbf{ż} = \mathbf{xy} - \frac{8}{3}~\mathbf{z}$ + + + +```python +--- SINDy results ---- +ẋ = 9.969 𝑦 -9.966 𝑥 +ẏ = -0.972 𝑦 + 27.833 𝑥 -0.995 𝑥𝓏 +ż = -2.657 𝓏 + 0.997 𝑥𝑦 + + [𝓏, 𝓏^2, 𝑦, 𝑦𝓏, 𝑦^2, 𝑥, 𝑥𝓏, 𝑥𝑦, 𝑥^2] +[[ 0. 0. 9.968 0. 0. -9.965 0. 0. 0.] + [ 0. 0. -0.971 0. 0. 27.832 -0.995 0. 0.] + [-2.656 0. 0. 0. 0. 0. 0. 0.996 0.]] +``` + + +

+ + +

+
+ + ## Linear-2D + +True model's equation \ +$\mathbf{ẋ} = -0.1\mathbf{x} + 2.0\mathbf{y}$ \ +$\mathbf{ẏ} = -2.0\mathbf{x} - 0.1\mathbf{y}$ + +```python +--- SINDy results ---- +ẋ = 2.000 𝑦 -0.100 𝑥 +ẏ = -0.100 𝑦 -2.000 𝑥 + +$[𝑦, 𝑦^2, 𝑥, 𝑥𝑦, 𝑥^2]$ +[[ 1.999 0. -0.100 0. 0.] + [-0.099 0. -1.999 0. 0.]] +``` + + + +

+ + +

+
+ + ## Linear-3D + +True model's equation \ +$\mathbf{ẋ} = -0.1\mathbf{x} + 2\mathbf{y}$ \ +$\mathbf{ẏ} = -2\mathbf{x} - 0.1\mathbf{y}$ \ +$\mathbf{ż} = -0.3\mathbf{z}$ + +```python +--- SINDy results ---- +ẋ = 2.000 𝑦 -0.100 𝑥 +ẏ = -0.100 𝑦 -2.000 𝑥 +ż = -0.300 𝓏 + +[1, 𝓏, 𝓏^2, 𝑦, 𝑦.𝓏, 𝑦^2, 𝑥, 𝑥𝓏, 𝑥.𝑦, 𝑥^2] +[[ 0. 0. 1.999 0. 0. -0.100 0. 0. 0.] + [ 0. 0. -0.100 0. 0. -1.999 0. 0. 0.] + [-0.299 0. 0. 0. 0. 0. 0. 0. 0.]] +``` + + +

+ + +

+
+ + ## Cubic-2D + +True model's equation \ +$\mathbf{ẋ} = -0.1\mathbf{x}^3 + 2.0\mathbf{y}^3$ \ +$\mathbf{ẏ} = -2.0\mathbf{x}^3 - 0.1\mathbf{y}^3$ + +```python +--- SINDy results ---- +ẋ = 1.999 𝑦^3 -0.100 𝑥^3 +ẏ = -0.099 𝑦^3 -1.998 𝑥^3 + +[𝑦, 𝑦^2, 𝑦^3, 𝑥, 𝑥𝑦, 𝑥𝑦^2, 𝑥^2, 𝑥^2𝑦, 𝑥^3] +[[ 0. 0. 1.99 0. 0. 0. 0. 0. -0.100] + [ 0. 0. -0.099 0. 0. 0. 0. 0. -1.99]] +``` + + +

+ + +

+
+ + + + + + + + + + ## References [1] Brunton SL, Proctor JL, Kutz JN. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc Natl Acad Sci U S A. 2016 Apr 12;113(15):3932-7. doi: 10.1073/pnas.1517384113. Epub 2016 Mar 28. PMID: 27035946; PMCID: PMC4839439. -