A unified sparse optimization framework to learn parsimonious physics-informed models from data

TL;DR

This paper presents a unified sparse optimization framework for learning interpretable, physically consistent dynamical models from data, capable of generalizing beyond observed regimes and handling real-world challenges.

Contribution

It introduces a flexible, non-convex sparse optimization method that learns parsimonious, interpretable models with physical constraints from diverse datasets.

Findings

Successfully discovers simple dynamical models in example systems.

Handles outliers, parametric dependencies, and physical constraints.

Models generalize well to new parameter regimes.

Abstract

Machine learning (ML) is redefining what is possible in data-intensive fields of science and engineering. However, applying ML to problems in the physical sciences comes with a unique set of challenges: scientists want physically interpretable models that can (i) generalize to predict previously unobserved behaviors, (ii) provide effective forecasting predictions (extrapolation), and (iii) be certifiable. Autonomous systems will necessarily interact with changing and uncertain environments, motivating the need for models that can accurately extrapolate based on physical principles (e.g. Newton's universal second law for classical mechanics, $F=ma$). Standard ML approaches have shown impressive performance for predicting dynamics in an interpolatory regime, but the resulting models often lack interpretability and fail to generalize. We introduce a unified sparse optimization framework that learns governing dynamical systems models from data, selecting relevant terms in the dynamics from a library of possible functions. The resulting models are parsimonious, have physical interpretations, and can generalize to new parameter regimes. Our framework allows the use of non-convex sparsity promoting regularization functions and can be adapted to address key challenges in scientific problems and data sets, including outliers, parametric dependencies, and physical constraints. We show that the approach discovers parsimonious dynamical models on several example systems. This flexible approach can be tailored to the unique challenges associated with a wide range of applications and data sets, providing a powerful ML-based framework for learning governing models for physical systems from data.

Figures

Equations

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{x}(t)=\mathbf{f}(\mathbf{x}(t))$$

$$\dot{\mathbf{X}}=\mathbf{\Theta }(\mathbf{X})\mathbf{\Xi }.$$

$$\underset{\mathbf{\Xi }}{\min }\frac{1}{2}\Vert \dot{\mathbf{X}}-\mathbf{\Theta }(\mathbf{X})\mathbf{\Xi }{\Vert }^2+\lambda R(\mathbf{\Xi })$$

$$\underset{\mathbf{\Xi },\mathbf{W}}{\min }\frac{1}{2}\Vert \dot{\mathbf{X}}-\mathbf{\Theta }(\mathbf{X})\mathbf{\Xi }{\Vert }^2+\lambda R(\mathbf{W})+\frac{1}{2\nu }\Vert \mathbf{\Xi }-\mathbf{W}{\Vert }^2.$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$${\dot{x}}_3$$

$$\underset{\mathbf{\Xi },\mathbf{v}}{\min }\sum _{i=1}^m\frac{1}{2}{v}_i\Vert (\dot{\mathbf{X}}-\mathbf{\Theta }(\mathbf{X})\mathbf{\Xi }{)}_i{\Vert }^2+\lambda R(\mathbf{\Xi })$$

$$\text{mbox}s.t.0\le v_i\le 1,1^T\mathbf{v}=h,$$

$$\underset{\mathbf{\Xi },\mathbf{W},\mathbf{v}\in {\Delta }_h}{\min }$$

$$+\lambda R(\mathbf{W})+\frac{1}{2\nu }\Vert \mathbf{\Xi }-\mathbf{W}{\Vert }^2.$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$${\dot{x}}_3$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{X}$$

$$\dot{Y}$$

$$V(x,y)=-\frac{\omega }{2}(x^2+y^2)+\frac{\alpha }{4}(x^2+y^2{)}^2.$$

$$\dot{X}$$

$$\dot{Y}$$

$${\xi }_{1,3}$$

$${\xi }_{1,5}$$

$$2{\xi }_{1,6}$$

$${\xi }_{1,8}$$

$${\xi }_{1,9}$$

$$3{\xi }_{1,10}$$

$$\underset{\mathbf{\Xi }}{\min }\frac{1}{2}\Vert (\dot{\mathbf{X}}-\mathbf{\Theta }(\mathbf{X})\mathbf{\Xi }){\Vert }^2+\lambda R(\mathbf{\Xi })$$

$$\text{mbox}s.t.\mathbf{C}\mathbf{\xi }=\mathbf{d}.$$

$$\underset{\mathbf{\Xi },\mathbf{W}}{\min }\frac{1}{2}\Vert (\dot{\mathbf{X}}-\mathbf{\Theta }(\mathbf{X})\mathbf{\Xi }){\Vert }^2+\lambda R(\mathbf{W})+\frac{1}{2\nu }\Vert \mathbf{\Xi }-\mathbf{W}{\Vert }^2$$

$$\text{mbox}s.t.\mathbf{C}\mathbf{\xi }=\mathbf{d},\mathbf{\xi }=\mathbf{\Xi }(:).$$

$$\text{mbox}vec(\mathbf{\Theta }(\mathbf{X})\mathbf{\Xi })=(\mathbf{I}\otimes \mathbf{\Theta }(\mathbf{X}))\mathbf{\xi }:={\mathbf{\Theta }}_{\mathbf{X}}^{\otimes }\mathbf{\xi }$$

$$\left[\begin{array}{cc}{{\mathbf{\Theta }}_{\mathbf{X}}^{\otimes }}^T{\mathbf{\Theta }}_{\mathbf{X}}^{\otimes }+\frac{1}{\nu }\mathbf{I}& {\mathbf{C}}^T\\ \mathbf{C}& 0\end{array}\right]\left[\begin{array}{c}\mathbf{\xi }\\ \mathbf{\varphi }\end{array}\right]=\left[\begin{array}{c}{{\mathbf{\Theta }}_{\mathbf{X}}^{\otimes }}^T\dot{\mathbf{X}}+\frac{1}{\nu }{\mathbf{W}}^{k-1}\\ \mathbf{d}\end{array}\right]$$

$${\left({{\mathbf{\Theta }}_{\mathbf{X}}^{\otimes }}^T{\mathbf{\Theta }}_{\mathbf{X}}^{\otimes }+\frac{1}{\nu }\mathbf{I}\right)}^{-1}={\left(\mathbf{\Theta }(\mathbf{X}{)}^T\mathbf{\Theta }(\mathbf{X})+\frac{1}{\nu }\mathbf{I}\right)}^{-1}\otimes \mathbf{I}$$

$${\left(\mathbf{\Theta }(\mathbf{X}{)}^T\mathbf{\Theta }(\mathbf{X})+\frac{1}{\nu }\mathbf{I}\right)}^{-1}$$

$$\mathbf{\varphi }={\left(\mathbf{C}{\left(\mathbf{\Theta }(\mathbf{X}{)}^T\mathbf{\Theta }(\mathbf{X})+\frac{1}{\nu }\mathbf{I}\right)}^{-1}\otimes \mathbf{I}{\mathbf{C}}^T\right)}^{-1}\text{mbox}RHS$$

$$\text{mbox}RHS=\mathbf{d}-\mathbf{C}{\left({{\mathbf{\Theta }}_{\mathbf{X}}^{\otimes }}^T{\mathbf{\Theta }}_{\mathbf{X}}^{\otimes }+\frac{1}{\nu }\mathbf{I}\right)}^{-1}({{\mathbf{\Theta }}_{\mathbf{X}}^{\otimes }}^T\dot{\mathbf{X}}+\frac{1}{\nu }{\mathbf{W}}^{k-1}).$$

$$\mathbf{\xi }={\left({{\mathbf{\Theta }}_{\mathbf{X}}^{\otimes }}^T{\mathbf{\Theta }}_{\mathbf{X}}^{\otimes }+\frac{1}{\nu }\mathbf{I}\right)}^{-1}({{\mathbf{\Theta }}_{\mathbf{X}}^{\otimes }}^T\dot{\mathbf{X}}+\frac{1}{\nu }{\mathbf{W}}^{k-1}-{\mathbf{C}}^T\mathbf{\varphi })$$

$$\underset{\mathbf{\Xi },\mathbf{\alpha }}{\min }\frac{1}{2}\Vert \dot{\mathbf{X}}-\mathbf{\Theta }(\mathbf{X},\mathbf{\alpha })\mathbf{\Xi }{\Vert }^2+\lambda R(\mathbf{\Xi }).$$

Full text

††∗ Corresponding author: (kpchamp@uw.edu).

1 Introduction

With abundant data being generated across scientific fields, researchers are increasingly turning to machine learning (ML) methods to aid scientific inquiry. In addition to standard techniques in clustering and classification, ML is now being used to discover models that characterize and predict the behavior of physical systems. Unlike many applications in ML, interpretation, generalization and extrapolation are the primary objectives for engineering and science, hence we must identify parsimonious models that have the fewest terms required to describe the dynamics. This is in contrast to neural networks (NNs), which are defined by exceedingly large parametrizations which typically lack interpretability or generalizability. A breakthrough approach in model discovery used symbolic regression to learn the form of governing equations from data [4, 32]. Sparse identification of nonlinear dynamics (SINDy) [5] is a related approach that uses sparse regression to find the fewest terms in a library of candidate functions required to model the dynamics. Because this approach is based on a sparsity-promoting linear regression, it is possible to incorporate partial knowledge of the physics, such as symmetries, constraints, and conservation laws (e.g., conservation of mass, momentum, and energy) [17]. In this work, we develop a unified sparse optimization framework for dynamical system discovery that enables one to simultaneously discover models, trim corrupt training data, enforce known physics, and identify parametric dependency in the equations.

Although studied for decades in dynamical systems [10, 22], the universal approximation capabilities of NNs [13, 11] have generated a resurgence of approaches to model time-series data [20, 37, 38, 41, 34, 18, 26, 27, 2, 24, 6]. NNs can also learn coordinate transformations that simplify the dynamical system representation [20, 38, 41, 34, 18, 23, 6] (e.g. Koopman representations [15, 21]). However, NNs generally struggle with extrapolation, are difficult to interpret, and cannot readily enforce known or partially known physics.

In contrast, the SINDy algorithm [5, 6] has been shown to produce interpretable and generalizable dynamical systems models from limited data. SINDy has been applied broadly to identify models for optical systems [33], fluid flows [17], chemical reaction dynamics [12], plasma convection [7], structural modeling [16], and for model predictive control [14]. It is also possible to extend SINDy to identify partial differential equations [29, 31] and dynamical systems with parametric dependencies [30], to trim corrupt data [36], and to incorporate partially known physics and constraints [17]. Because the approach is fundamentally based on a sparsity-regularized regression, there is an opportunity to unify these innovations via the sparse relaxed regularized regression (SR3) algorithm [42], resulting in a unified sparse model discovery framework.

1.1 Basic problem formulation

SINDy [5] enables the discovery of nonlinear dynamical systems models from data. Assume we have data from a dynamical system

[TABLE]

where $\mathbf{x}(t)\!\in\!\mathbb{R}^{n}$ is the state of the system at time $t$. We want to find the terms in $\mathbf{f}$ given the assumption that $\mathbf{f}$ has only a few active terms: it is sparse in the space of all possible functions of $\mathbf{x}(t)$. Given snapshot data $\mathbf{X}=\begin{bmatrix}\mathbf{x}_{1}&\mathbf{x}_{2}&\cdots&\mathbf{x}_{m}\end{bmatrix}^{T}$ and $\mathbf{\dot{X}}=\begin{bmatrix}\mathbf{\dot{x}}_{1}&\mathbf{\dot{x}}_{2}&\cdots&\mathbf{\dot{x}}_{m}\end{bmatrix}^{T}$, we build a library of candidate functions $\boldsymbol{\Theta}(\mathbf{X})=[\boldsymbol{\theta}_{1}(\mathbf{X})\cdots\boldsymbol{\theta}_{p}(\mathbf{X})]$. We then seek a solution of

[TABLE]

where $\boldsymbol{\Xi}\!=\!(\boldsymbol{\xi}_{1}\ \boldsymbol{\xi}_{2}\ \cdots\ \boldsymbol{\xi}_{n})$ are sparse coefficient (loading) vectors. A natural optimization is given by

[TABLE]

where $R(\cdot)$ is a regularizer that promotes sparsity. When $R$ is convex, a range of well-known algorithms for (2) are available. The standard approach is to choose $R$ to be the sparsity-promoting $\ell_{1}$ norm, which is the convex relaxation of the $\ell_{0}$ norm. In this case, SINDy is solved via LASSO [35]. In practice, LASSO does not perform well at coefficient selection (see Section 2.1 for details). In the context of dynamics discovery we would like to use non-convex $R$, specifically the $\ell_{0}$ norm.

The standard SINDy algorithm performs sequential thresholded least squares (STLSQ): given a parameter $\eta$ that specifies the minimum magnitude for a coefficient in $\boldsymbol{\Xi}$, perform a least squares fit and then zero out all coefficients with magnitude below the threshold. This process of fitting and thresholding is performed until convergence. While this method works very well, it is customized to the least squares formulation and does not readily accommodate extensions including incorporation of additional constraints, robust formulations, or nonlinear parameter estimation. A number of extensions to SINDy have been developed and required adaptations to the optimization algorithm [36, 29, 31, 17, 14].

2 Formulation and approach

We extend the optimization formulation (2) to include additional structure, robustness to outliers, and nonlinear parameter estimation using the sparse relaxed regularized regression (SR3) approach that uses relaxation and partial minimization [42]. SR3 for (2) introduces the auxiliary variable $\mathbf{W}$ and relaxes the optimization to

[TABLE]

We can solve (3) using the alternating update rule in Algorithm 1. This requires only least squares solves and prox operators [42]. The resulting solution approximates the original problem (2) as $\nu\downarrow 0$. When $R$ is taken to be the $\ell_{0}$ penalty, the prox operator is hard thresholding, and Algorithm 1 is similar, but not equivalent, to thresholded least squares, and performs similarly in practice. However, unlike thresholded least squares, the SR3 approach easily generalizes to new problems and features.

2.1 Performance of SR3 for SINDy

SR3 for SINDy provides an optimization framework that both (i) enables the identification of truly sparse models and (ii) can be adapted to include additional features. We first compare SR3 to both STLSQ and the LASSO algorithm. While STLSQ works well for identifying sparse models that capture the behavior of a system, it is a standalone method without a true optimization cost function, meaning the algorithm must be reformulated to work with other adaptations to the SINDy problem [19]. LASSO provides a standard optimization approach but does not successfully identify sparse models. Even with noiseless (clean) data, LASSO models for SINDy typically have many coefficients that are small in magnitude but nonzero. Obtaining a sparse set of coefficients is key for interpretability. SR3 works with nonconvex regularization functions such as the $\ell_{0}$ norm, enabling the identification of truly sparse models.

In Fig. 2 we compare these algorithms using data from the canonical chaotic Lorenz system:

[TABLE]

We simulate the system from 20 initial conditions and fit a SINDy model with polynomials up to order 3 using the following optimization approaches: STLSQ with threshold $0.1$, SR3 with $\ell_{0}$ regularization, LASSO with a regularization weight of $0.1$, and LASSO with a regularization weight of $50$. For each model we analyze both the sparsity pattern of the coefficient matrix, and simulations of the resulting model on test trajectories. As shown in Fig. 2, STLSQ and SR3 yield the same correct sparsity pattern. In simulation, both track a Lorenz test trajectory for several trips around the attractor before eventually diverging. The eventual deviation is expected due to the chaotic nature of the Lorenz system, as a slight difference in coefficient values or initial conditions can lead to vastly different trajectories (although the trajectories continue to fill in the Lorenz attractor). These models also track the behavior well for a trajectory that starts off the attractor. The LASSO models both have many terms that are small in magnitude but still nonzero. The LASSO model with a low sparsity parameter has similar performance to the STLSQ and SR3 models, but is not a truly sparse model. As the regularization penalty is increased, rather than removing the unimportant terms in the dynamics the method removes many of the true coefficients in the Lorenz model. The LASSO model with a high sparsity parameter has a very poor fit for the dynamics.

3 Simultaneous Sparse Inference and Data Trimming

Many real world data sets contain corrupted data and/or outliers, which is problematic for model identification methods. For SINDy, outliers can be especially problematic, as derivative computations are corrupted. Many data modeling methods have been adapted to deal with corrupted data, resulting in “robust” versions of the methods (such as robust PCA). The SR3 algorithm for SINDy can be adapted to incorporate trimming of outliers, providing a robust optimization algorithm for SINDy. Starting with least trimmed squares [28], extended formulations that simultaneously fit models and trim outliers are widely used in statistical learning. Trimming has proven particularly useful in the high-dimensional setting when used with the LASSO approach and its extensions [39, 40].

The high-dimensional trimming extension applied to (2) takes the form

[TABLE]

where $h$ is an estimate of the number of ‘inliers’ out of the potential $m$ rows of the system. The set $\Delta_{h}:=\{\mathbf{v}:0\leq v_{i}\leq 1,\quad{\bf 1}^{T}\mathbf{v}=h\}$ is known as the capped simplex. Current algorithms for (5), such as those of [40], rely on LASSO formulations and thus have significant limitations (see previous section). Here, we use the SR3 strategy (3) to extend to the trimmed SINDy problem (5):

[TABLE]

We then use the alternating Algorithm 2 to solve the problem. The step size $\beta$ is completely up to the user, as discussed in the convergence theory (see Appendix S3). The trimming algorithm requires specifying how many samples should be trimmed, which can be chosen by estimating the level of corruption in the data. Estimating derivatives using central differences, for instance, makes derivative estimates on either side of the original corrupted sample corrupt as well, meaning that three times as many samples as were originally corrupted will be bad. Thus trimming will need to be more than the initial estimate of how many samples were corrupted. Trimming ultimately can help identify and remove points with bad derivative estimates, leading to a better SINDy model fit.

3.1 Example: Lorenz

We demonstrate the use of SINDy SR3 for trimming outliers on data from the Lorenz system (4). We randomly select a subset of samples to corrupt, adding a high level of noise to these samples to create outliers. We apply the SINDy SR3 algorithm with trimming to simultaneously remove the corrupted samples and fit a SINDy model. Fig. 3 shows the results of trimming on a dataset with 10% of the samples corrupted. The valid data points are shown in gray and the corrupt data points are highlighted in red. As derivatives are calculated directly from the data using central differences, this results in closer to 30% corruption (as derivative estimates on either side of each corrupt sample will also be corrupted). We find that the algorithm converges on the correct solution more often when a higher level of trimming is specified: in other words, it is better to remove some clean data along with all of the outliers than to risk leaving some outliers in the data set. Accordingly, we set our algorithm to trim 40% of the data. Despite the large fraction of corrupted samples, the method is consistently able to identify the Lorenz model (or a model with only 1-2 extra coefficients) from the remaining data in repeated simulations. As initial conditions starting off the attractor are included in the data, some samples of the off-attractor behavior are also trimmed. In contrast, if standard SINDy SR3 (without trimming) is applied to the data set with corrupted samples, the correct model is not identified.

3.2 Example: Rossler

The Rossler system

[TABLE]

exhibits chaotic behavior characterized by regular orbits around an attractor in the $x_{1},x_{2}$ plane combined with occasional excursions into the $x_{3}$ plane. The Rossler attractor is plotted in Fig. 3 with 1% of samples corrupted (highlighted in red). If standard SINDy SR3 is applied to this data, the system is not correctly identified. However with the algorithm set to trim 5% of the data, the system is consistently identified (or nearly identified, with only 1-2 incorrect coefficients) in repeated trials.

Much of the Rossler attractor lies in the $x_{1},x_{2}$ plane, which means that the excursions into the $x_{3}$ dimension can be seen as outliers in the dynamics and are more likely to be flagged as potential outliers. However, if the algorithm is run to convergence most of these points are typically recognized as part of the true dynamics and not removed.

4 Incorporating physical constraints

Physical systems are often subject to constraints based on first principles, such as conservation laws (e.g. conservation of mass or momentum). When performing model discovery, we may require a model that satisfies such a set of desired physical properties. In the case of SINDy, such requirements may manifest as constraints on the form of the coefficient matrix. This constrained approach has already been implemented with the standard thresholded least squares algorithms for SINDy [17], and it can also be straightforwardly implemented in the SR3 approach.

4.1 2D Duffing

Hamiltonian systems and gradient systems are commonly studied in physics and dynamical systems. In these systems, the equations are gradients of a Hamiltonian or potential function. This imposes a constraint on the form of the coefficient matrix, as the individual equations must be partial derivatives of the same function. We demonstrate the constrained optimization on a Hamiltonian system: the Duffing oscillator in two spatial dimensions. The full system in four variables is

[TABLE]

where $x,y$ represent spatial position and $X,Y$ represent momentum. The potential function for the 2D Duffing oscillator is

[TABLE]

Because the variables $X,Y$ are fixed to represent momentum, we apply SINDy only to fit the second two equations. If $V(x,y)$ is assumed to be a polynomial, the gradient functions will also be polynomials and a specific structure is imposed on the form of the coefficient matrix, with an interdependence between the two columns. In particular, if the equations are written as

[TABLE]

the constraints that must be imposed are

[TABLE]

By vectorizing the coefficient matrix $\boldsymbol{\xi}=\boldsymbol{\Xi}(:)$, this set of constraints can be written as a linear constraint $\mathbf{C}\boldsymbol{\xi}=\mathbf{d}$. The form of the coefficient matrix in the 2D Duffing example is illustrated in Figure 4.

The problem of interest is now given by

[TABLE]

Following the relaxation strategy used in the previous sections, we obtain the modified problem

[TABLE]

To solve this problem, we need to adapt the relax and split strategy to partially minimize the problem in $\boldsymbol{\Xi}$, which is a quadratic with affine equality constraints. This problem has a closed form solution, just as a quadratic problem, through solving an augmented system. To form this system, we first use the vec-kron identity:

[TABLE]

We can form a linear system corresponding to optimality conditions for (10) with respect to $\boldsymbol{\xi}$:

[TABLE]

Conceptually, (11) is just a large linear system. However, its dimension is very high, and so in the remainder of this section we explain how to use the structure of the matrices to solve it efficiently. The most important tools include the Kronecker product and associated identities. We use the fact that

[TABLE]

We can therefore compute the explicit inverse (in the lower dimensional space)

[TABLE]

and then use it to solve (11). Specifically, we can reduce (11) to a small problem (same dimension as number of constraints):

[TABLE]

where the right hand side is given by

[TABLE]

Once $\boldsymbol{\phi}$ is computed, we use (11) to solve for $\boldsymbol{\xi}$:

[TABLE]

These steps are summarized in Algorithm 3.

We apply both the standard (unconstrained) and constrained SINDy SR3 algorithms to data from the 2D Duffing system with varying levels of noise added. For very low noise levels, both systems correctly identify a coefficient matrix that satisfies the provided constraints. As the noise level increases, the coefficients identified by the unconstrained algorithm drift from satisfying the gradient constraint, while the coefficients identified by the constrained algorithm keep the desired form. In this example it is important to note that the constrained and unconstrained models have comparable performance in reproducing the system behavior via simulation. In other examples, constraints may provide additional benefits such as helping to identify a stable system [17]. Although the constrained SR3 algorithm imposes the constraints on the full coefficients $\boldsymbol{\Xi}$ and not the auxiliary coefficients $\mathbf{W}$, in this example the auxiliary coefficients still satisfy the constraints, producing a model of the desired form.

5 Parameterized library functions

In standard examples of SINDy, the library is chosen to contain polynomials, which make a natural basis for many models in the physical sciences. However, many systems of interest may include more complicated terms in the dynamics, such as exponentials or trigonometric functions, that include parameters that contribute nonlinearly to the fitting problem. In addition to parameterized basis functions, systems may be subject to parameterized external forcing: for example, periodic forcing where the exact frequency of the forcing is unknown. SINDy with unknown parameters is given by

[TABLE]

This is a regularized nonlinear least squares problem. The SR3 approach makes it possible to devise an efficient algorithm for this problem as well. The relaxed formulation is given by

[TABLE]

We solve (13) using Algorithm 4. The $\boldsymbol{\alpha}$ variable is updated using a true Newton step, where the gradient and Hessian are computed using algorithmic differentiation.

The joint optimization for the parameterized case yields a nonconvex problem with potentially many local minima depending on the initial choice of the parameter(s) $\boldsymbol{\alpha}$. This makes it essential to assess the fit of the discovered model through model selection criteria. While the best choice of model may be clear, this means parameterized SINDy works best for models with only a small number of parameters in the library, as scanning through different initializations scales combinatorially with added parameters.

5.1 Lorenz with parameterized forcing

We consider (4) with $x_{1}$ forced by a parameterized hyperbolic tangent function $\tanh(\alpha_{1}t-\alpha_{2})$. The parameters $\alpha_{1},\alpha_{2}$ determine the steepness and location of the sigmoidal curve in the forcing function. We simulate the system with forcing parameters $\alpha_{1}=0.8,\alpha_{2}=3$. Fig. 5 shows the results of fitting the SINDy model with and without the parameterized forcing term in the library. In the case without forcing, the equation for $x_{1}$ is loaded up with several active terms in an attempt to properly fit the dynamics. The model is not able to reproduce the correct system behavior through simulation. In the case with forcing, we start with an initial guess of $\alpha_{1}=5,\alpha_{2}=10$ and perform the joint optimization to fit both the parameters and the coefficient matrix. The algorithm correctly identifies the forcing and finds the correct coefficient matrix. The resulting system matches the true dynamics for several trips around the attractor.

6 Discussion

Machine learning for model discovery in physics, biology and engineering is of growing importance for characterizing complex systems for the purpose of control and technological applications. Critical for the design and implementation in new and emerging technologies is the ability to interpret and generalize the discovered models, thus requiring that parsimonious models be discovered which are minimally parametrized. Moreover, model discovery architectures must be able to incorporate the effects of constraints, provide robust models, and/or give accurate nonlinear parameter estimates. We here propose the SINDy-SR3 method which integrates a sparse regression framework for parsimonious model discovery with a unified optimization algorithm capable of incorporating many of the critical features necessary for real-life applications. This includes features for handling corrupt data, imposing constraints, and learning parametrizations, all tasks that the basic SINDy algorithm is unable to handle. We demonstrate its accuracy and efficiency on a number of example problems, showing that SINDy-SR3 is a viable framework for the engineering sciences.

All code and data can be found at github.com/kpchamp/SINDySR3

Acknowledgments

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1256082. SLB acknowledges support from the Army Research Office (W911NF-19-1-0045). JNK acknowledges support from the Air Force Office of Scientific Research (FA9550-17-1-0329). Aleksandr Aravkin is supported by the Washington Research Foundation Data Science Professorship.

S1 Choice of parameters for SR3

The SR3 algorithm requires the specification of two parameters, $\nu$ and $\lambda$. The parameter $\nu$ controls how closely the relaxed coefficient matrix $\mathbf{W}$ matches $\boldsymbol{\Xi}$: small values of $\nu$ encourage $\mathbf{W}$ to be a close match for $\boldsymbol{\Xi}$, whereas larger values will allow $\mathbf{W}$ to be farther from $\boldsymbol{\Xi}$.

The parameter $\lambda$ determines the strength of the regularization. If the regularization function is the $\ell_{0}$ norm, the parameter $\lambda$ can be chosen to correspond to the coefficient threshold used in the sequentially thresholded least squares algorithm (which determines the lowest magnitude value in the coefficient matrix). This is because the prox function for the $\ell_{0}$ norm will threshold out coefficients below a value determined by $\nu$ and $\lambda$. In particular, if the desired coefficient threshold is $\eta$, we can take

[TABLE]

and the prox update will threshold out values below $\eta$. In the examples shown here, we determine $\lambda$ in this manner based on the desired values for $\nu,\eta$. If the desired coefficient threshold is known (which is the case for the examples studied here, but may not be the case for unknown systems), this gives us a single parameter to adjust: $\nu$. With $\lambda$ defined in this manner, decreasing $\nu$ provides more weight to the regularization, whereas increasing $\nu$ provides more weight to the least squares model fit.

S2 Simulation details

S2.1 Performance of SR3 for SINDy

We illustrate a comparison of three algorithms for SINDy using data from the canonical example of the chaotic Lorenz system (4). To generate training data, we simulate the system from $t=0$ to 10 with a time step of $\Delta t=0.005$ for 20 initial conditions sampled from a random uniform distribution in a box around the attractor ($x\in[-36,36],y\in[-48,48],z\in[-16,66]$). This results in a data set with $40\times 10^{4}$ samples. We add random Gaussian noise with a standard deviation of $10^{-2}$ and compute the derivatives of the data using the central difference method. The SINDy library matrix $\boldsymbol{\Theta}(\mathbf{X})$ is constructed using polynomial terms through order 3.

We find the SINDy model coefficient matrix using the following optimization approaches: sequentially thresholded least squares (STLSQ) with threshold $0.1$, SR3 with $\ell_{0}$ regularization, LASSO with a regularization weight of $0.1$, and LASSO with a regularization weight of $50$. The STLSQ algorithm is performed by doing 10 iterations of the following procedure: (1) perform a least squares fitting on remaining coefficients, (2) remove all coefficients with magnitude less than $0.1$. The LASSO models are fit using the scikit-learn package [25]. LASSO models are fit without an intercept. For SR3 we initialize the coefficient matrix using least squares and we use parameters $\nu=1$ and $\lambda=0.005$ (which corresponds to a coefficient threshold of $0.1$, see Appendix S1). For all methods, an unbiasing procedure is performed at the end of the fitting: after the respective optimization algorithm is applied to select the correct terms in the dynamics, a least squares fit is performed on just these terms to find the coefficient values.

For each of the four resulting models we analyze (1) the sparsity pattern of the coefficient matrix and (2) the simulation of the resulting dynamical systems model. We compare the sparsity pattern of the coefficient matrix against the true sparsity pattern for the Lorenz system: SR3 and STLSQ identify the correct sparsity pattern, where as the LASSO models do not. For all models, we simulate the identified system on test trajectories using randomly selected initial conditions from two distributions: near attractor (using same distribution as the training set data) and off attractor (using initial conditions selected in a box outside of the training set data, farther from the attractor). In each case, 50 initial conditions are sampled and the simulations are run from $t=0$ to $5$ for each initial condition. For the near attractor data, the $R^{2}$ scores for the STLSQ, SR3, and LASSO model with low regularization are $0.82$, $0.82$, and $0.83$, respectively. For the off attractor data, the $R^{2}$ scores are $0.84$, $0.84$, and $0.82$. In both cases, the LASSO model with heavy regularization had very poor simulation results and had a strongly negative $R^{2}$ score. For all models, simulations for initial conditions are $(-8,7,27)$ (on attractor) and $(0.01,0.01,60)$ (off attractor) are shown in Figure 2 in the main text.

S2.2 Data trimming: Lorenz

We demonstrate the use of the SR3-trimming algorithm 2 on data from the Lorenz system (4). We simulate the system over the same time as in Appendix S2.1 from 5 randomly sampled initial conditions. This results in a data set with $10^{4}$ samples. We add Gaussian noise to the data with standard deviation $10^{-3}$. We then randomly choose 10% of the samples to corrupt ($1000$ total samples). For each state variable of each corrupted sample, noise chosen from a random uniform distribution over $[-50,50]$ is added. Derivatives are calculated from the data using central difference after the corruption is applied. We then apply the SR3-trimming algorithm, specifying that around 40% of the data points will be trimmed. We use SR3 parameters $\nu=1$ and $\lambda=0.005$ (corresponding to a coefficient threshold of $0.1$), and the step size is taken to be the default value $\beta=1$. With repeated testing we find that the algorithm is consistently able to correctly remove the outliers from the data set and identify the Lorenz system.

S2.3 Data trimming: Rossler

As an additional example, we test the trimming algorithm on data from the Rossler system (7). We generate sample data from 5 randomly sampled initial conditions around the portion of the attractor in the $x_{1},x_{2}$ plane, simulating trajectories from $t=0$ to 50 with a time step of $\Delta t=0.01$. Our data set consists of 25000 samples. We add Gaussian noise with a standard deviation of $10^{-3}$ and add outliers to 1% of the data in the same manner as in Appendix S2.2, with the noise level chosen from a random uniform distribution over $[-100,100]$. Derivatives are calculated from the corrupted data using central difference. We apply the SR3-trimming algorithm with parameters $\nu=1$ and $\lambda=1.25\times 10^{-3}$ (corresponding to a coefficient threshold of $0.05$) and default step size $\beta=1$. On repeated trials we find that if we trim 5% of the data, the system is correctly identified in most cases (or only 1 or 2 coefficients are misidentified).

S2.4 Physical constraints: 2D Duffing

We demonstrate the use of SINDy SR3 with constraints using the example of a 2D Duffing oscillator. We simulate the system (8) from time $t=0$ to 10 with a time step of $0.01$ for 20 initial conditions. The initial conditions are chosen from a uniform distribution with $x,y,X,Y\in[-\pi,\pi]$. We use constrained SINDy SR3 to discover the latter two equations in (8), imposing the constraints laid out in Section 4. We use a polynomial library with terms up through order 3 and parameters $\nu=1$ and $\lambda=0.005$ (corresponding to a coefficient threshold of $0.1$).

To assess the performance of the constrained algorithm, we apply noise at various levels and compare the performance of the constrained algorithm with the unconstrained algorithm. Gaussian noise is applied to the data at standard deviations of 0, 0.01, 0.05, and 0.1. With no noise and noise strength 0.01, both the constrained and unconstrained algorithms discover a model that satisfies the gradient constraint. At stronger noise strengths, both approaches identify the correct terms in the dynamics but only the models found via the constrained approach satisfy the gradient constraint. We also compare performance by simulating the identified systems on test trajectories from 50 randomly chosen initial conditions in the same distribution as the training set initial conditions. While the $R^{2}$ scores decrease slightly as noise increases, the constrained and unconstrained models have nearly identical scores at all noise levels.

S2.5 Parameterized library functions: Lorenz with parameterized forcing

To demonstrate the use of SR3 for SINDy with parameter estimation, we look at an example of the Lorenz system (4) forced by a parameterized hyperbolic tangent function $\tanh(\alpha_{1}t-\alpha_{2})$. The full set of equations for the system is

[TABLE]

The parameters $\alpha_{1},\alpha_{2}$ determine the steepness and location of the sigmoidal curve in the forcing function. We simulate the system as in Appendix S2.1 for a single initial condition $(8,-7,27)$ with forcing parameters $\alpha_{1}=0.8,\alpha_{2}=3$. We add Gaussian noise of standard deviation $10^{-3}$ and compute the derivatives via central difference.

We apply Algorithm 2 to perform a joint discovery of both the coefficients $\mathbf{W}$ and forcing parameters $\boldsymbol{\alpha}$. We use parameters $\nu=0.1$ and $\lambda=0.05$ (corresponding to coefficient threshold $0.1$). $\mathbf{W}$ is initialized using least squares, and as an initial guess for $\boldsymbol{\alpha}$ we use $\boldsymbol{\alpha}^{0}=(5,10)$. The algorithm discovers the correct parameters $\boldsymbol{\alpha}$ as well as the correct sparsity pattern in the coefficient matrix. We simulate the system and see that the discovered system tracks the behavior for several trips around the attractor. Results are shown in Figure 5.

For comparison, we apply the SR3 algorithm for SINDy with no forcing term in the library, using the same SR3 parameters as in the forcing case. The resulting model has many active terms in the equation for $\dot{x}_{1}$, as it attempts to capture the forcing behavior with polynomials of $x_{1},x_{2},x_{3}$. This model does not perform well in simulation, even from the same initial condition used in the training set. Figure 5 shows the coefficient matrix and model simulation for the discovered system.

S3 Convergence results

Here we state convergence results for Algorithm 1 and Algorithm 2. These algorithms fall under the framework of two classical methods, proximal gradient descent and the proximal alternating linearized minimization algorithm (PALM) [3]. While we demonstrate the use of Algorithm 4 on two example problems, this algorithm is much harder algorithm to analyze due to the complication from the Newton’s step. We leave obtaining theoretical guarantees of Algorithm 4 as future work.

S3.1 Convergence of Algorithm 1

Using the variable projection framework [9], we partially optimize out $\boldsymbol{\Xi}$ and then treat Algorithm 1 as the classical proximal gradient method on $\mathbf{W}$. The convergence result for Algorithm 1 is provided in [42, Theorem 2] and is restated here:

Theorem 1

Define the value function as,

[TABLE]

When $p$ is bounded below, we know that the iterators from Algorithm 1 satisfy,

[TABLE]

where $g^{k}\in\partial p(\mathbf{W}^{k})$ and $p^{*}=\min_{\mathbf{W}}p(\mathbf{W})$.

We obtain a sub-linear convergence rate for all prox-bounded regularizers $R$.

S3.2 Convergence of Algorithm 2

Following the same idea provided by the variable projection framework, the iterations from Algorithm 2 are equivalent with an alternating proximal gradient step between $\mathbf{W}$ and $\mathbf{v}$. This is the PALM algorithm, which is thoroughly analyzed in the context of trimming in [1] and [8]. We restate the convergence result here:

Theorem 2

Consider the value function,

[TABLE]

And we know that the iterators $(W^{k},v^{k})$ converge to the stationary point of $p$, with the rate,

[TABLE]

Algorithm 2 also requires the specification of a step size $\beta$ for the proximal gradient step for $\mathbf{v}$. Because the objective is linear with respect to $\mathbf{v}$, the step size will not influence the convergence result in the above theorem. However, because the objective is non-convex, $\beta$ will have an impact on where the solution lands. In this work we use a default step size of $\beta=1$ for all examples.

S3.3 Convergence of Algorithm 3

Algorithm 3 can be analyzed by bootstrapping on the analysis of Algorithm 1. As long as the constraint

[TABLE]

is feasible, it defines an affine subspace $\mathcal{S}$, which we can reparameterize in terms of a new variable

[TABLE]

Let $\mathbf{S}$ be the basis for the nullspace of $\mathbf{C}$. We can now rewrite the value function for (10) in terms of $\mathbf{S}$, parametrized by a new variable $\boldsymbol{\zeta}$:

[TABLE]

Since the solution for $(\boldsymbol{\phi},\boldsymbol{\xi})$ for the suproblem (10) is unique, Algorithm 3 is equivalent to optimizing this value function.

The convergence result for Algorithm 1 is therefore also provided by [42, Theorem 2], where the quadratic term to be used is

[TABLE]