Extracting Sparse High-Dimensional Dynamics from Limited Data
Hayden Schaeffer Giang Tran Rachel Ward
math.OC 34F05, 37H99, 65P99, 65L09, 65L99, 37N30


[
    [
    December 30, 2023
=====================





Extracting governing equations from dynamic data is an essential task in model selection and parameter estimation. The form of the governing equation is rarely known a priori; however, based on the sparsity-of-effect principle one may assume that the number of candidate functions needed to represent the dynamics is very small. In this work, we leverage the sparse structure of the governing equations along with recent results from random sampling theory to develop methods for selecting dynamical systems from under-sampled data. In particular, we detail three sampling strategies that lead to the exact recovery of first-order dynamical systems when we are given fewer samples than unknowns.  The first method makes no assumptions on the behavior of the data, and requires a certain number of random initial samples. The second method utilizes the structure of the governing equation to limit the number of random initializations needed.  The third method leverages chaotic behavior in the data to construct a nearly deterministic sampling strategy. Using results from compressive sensing, we show that the strategies lead to exact recovery, which is stable to the sparse structure of the governing equations and robust to noise in the estimation of the velocity. Computational results validate each of the sampling strategies and highlight potential applications.
 







§ INTRODUCTION

Since the scientific revolution in the seventeenth century, scientists have endeavored to extract increasingly sophisticated physical models from experimental and observational data. This process is usually done manually in the sense that a scientist, with a strong expertise in the field, must examine the data to uncover meaningful physical laws. Parsing through data by hand is often time-consuming, expensive, and infeasible – especially when the dimension of the system grows. However, with advances in computing and machine learning algorithms, automated discovery of physical models and mathematical equations directly from data is becoming increasingly possible. Some applications of data-based modeling include (but are certainly not limited to) weather predictions, controls for fluid flows <cit.>, construction of climate trend models <cit.>, and disease control models <cit.>.


Over the last decade, data-driven methods have seen strong growth due to the abundance of data and the development of sophisticated analytical tools. In this work, a computational method for identifying high-dimensional differential equations (which we will refer to as the `model') from under-sampled dynamical data is developed. We will restrict our attention to quadratic models, which contain a large subset of potential physical laws and applications; for example, Lorenz-like systems used in atmospheric science, Fisher's equation and related reaction-diffusion systems, Vlasov-Poisson equations from plasma physics, and fluid dynamics models like the Navier-Stokes equations. The governing equations for such systems are often moreover sparse quadratic polynomials, in that the dynamics depend on a small number of variables and second-order interactions.  Indeed, this allows us to recover the governing equations from under-sampled data: we only need to observe a number of samples of the data roughly on the order of the sparsity level of the system. Our method learns the governing models via a sparse optimization problem over a large set of potential candidate functions. It will be shown the our method selects the correct governing dynamical system (under certain conditions) even when the size of the candidate set far exceeds the size of the data set. In this way, these methods are one of the first applications of compressive sensing to model extraction for dynamical systems.



Recently, regression based methods have been developed and applied to model selection and parameter estimation of dynamic data. The authors of <cit.> first proposed the use of regression to select physical laws from synthetic and experimental data. In particular, a symbolic regression method was developed that compared computed derivatives of the data to analytic derivatives of trial functions, while controlling for the total number of trial functions selected (as to avoid overfitting).  In <cit.>, a sparsity-promoting method was proposed for extracting dynamical systems by comparing the computed velocity to a large set of potential trial functions. A sequential thresholded least-squares algorithm was used to fit a (redundant) set of trial functions (typically in the form of monomials) to the velocity.  In <cit.> sparsity-based methods were proposed for learning nonlinear partial differential equations from spatio-temporal data. In <cit.>, a non-convex sparse optimization method was developed to identify the underlying dynamical system from noisy data sets using an integrated trial set. In <cit.>, the authors utilized the Akaike information criterion to rank different sparse solutions using the method from <cit.>, allowing for the automated selection of different models when varying the method's free parameter. Using a sparse regression method along with the Takens' embedding theorem, a data-driven method was proposed in <cit.> to decompose chaotic systems into intermittently forced linear systems. A sparse convex optimization method was proposed for joint outlier detection and model selection in <cit.>, when the observed data is locally corrupted by high variance noise. Unlike the previous work in the literature, which was mostly empirical, the authors of <cit.> proved that the separation between the outliers and the `clean' data can be exactly recovered. The use of L^1 minimization has also appeared in various data-driven and scientific computing applications, for example <cit.>. Other data-driven methods for learning data structure and approximating dynamics include: the proper orthogonal decomposition 
<cit.>, Koopman representations
<cit.>, diffusion maps
<cit.>, and dynamic mode decomposition
<cit.>. 


These data-based regression methods often use the ℓ^0 or ℓ^1 penalty to promote sparsity in the learned models (i.e. to select a small number of active candidates from the large set of potential trial functions).  Soft-thresholding (related to the ℓ^1 penalty) for sparse recovery and denoising was first proposed in <cit.>. The ℓ^1 regularized least squares problem (referred to as the least absolute shrinkage and selection operator or LASSO) was introduced in <cit.> to reconstruct a sparse vector from linear observations. Conditions under which ℓ^1 penalized problems admit sparse solutions are detailed in <cit.> and have led to many applications in imaging and signal processing.  Refined conditions for the setting of function interpolation and approximation were developed in <cit.>.  


In this work, we develop a model selection and parameter estimation method for learning quadratic high-dimensional differential equations from under-sampled data. Using results from compressive sensing and sampling theory, we show that given a certain sampling of the initial data, a convex optimization problem can recover the coefficients and the governing equations exactly even when the data is under-sampled. We detail three sampling strategies, depending on the level of knowledge of the data or the governing equation. In particular, one can decrease the number of samples needed by adding assumptions to the system. It is important to note that the methodology presented in this work can be extended to higher-order polynomial systems and, more generally, governing equations which are sparse with respect to any given bounded orthogonal system. Also, this formulation benefits from well-known numerical methods for solving ℓ^1 penalized problems. [Our code is available on our github page: 
 <https://github.com/GiangTTran/ExtractingSparseHighDimensionalDynamicsFromLimitedData>.]





§ PROBLEM STATEMENT

Consider the dynamical variable x(t)∈ℝ^n (n≫ 1) governed by the equation ẋ=f(x) with initial data x(t_0)=x_0.  Assume that f(x) is a quadratic vector-valued equation in x, which can be written component-wise as:

    ẋ_1    = f_1(x_1,…, x_n)
    ẋ_2    = f_2(x_1,…, x_n)
        ⋮
    ẋ_n   = f_n(x_1,…, x_n).

The goal is to learn f_1, …, f_n, given x and ẋ. For a given initial condition x_0, assume that we can obtain a sequence of discrete measurements { x(t_1), x(t_2), …, x(t_m-1)} through either simulations or observations; however, the function f(x) is unknown.  The total number of samples along a given trajectory (including the initialization), denoted as m, can be small and thus we will refer to these measurements as a single burst. The key here is that by using a small number of bursts corresponding to random initial data, we are able to provide conditions on the recovery of the underlying dynamics. Although we will show the construction and results for quadratic systems, it is important to note that it is possible to generalized to high-order systems.

We construct the under-sampled measurements as follows. Let k be the index for a given burst; that is, if we sample the initial data x(t_0;k) from some random distribution, we can observe the k^th burst, i.e {x(t_0;k), x(t_1;k), x(t_2;k), …, x(t_m-1;k)}. The corresponding velocity along the burst, denoted by {ẋ(t_0;k), ẋ(t_1;k), ẋ(t_2;k), …, ẋ(t_m-1;k)}, is either observed or computed to some level of accuracy. The collection of trial functions corresponding to the k^th burst is denoted as:


    A^(k) =[                        1               x_1(t_0;k)               x_2(t_0;k)                        ⋯               x_n(t_0;k)             x^2_1(t_0;k)     x_1(t_0;k)x_2(t_0;k)                        ⋯             x^2_n(t_0;k);                        1               x_1(t_1;k)               x_2(t_1;k)                        ⋯               x_n(t_1;k)             x^2_1(t_1;k)     x_1(t_1;k)x_2(t_1;k)                        ⋯             x^2_n(t_1;k);                        1               x_1(t_2;k)               x_2(t_2;k)                        ⋯               x_n(t_2;k)             x^2_1(t_2;k)    x_1(t_2;k) x_2(t_2;k)                        ⋯             x^2_n(t_2;k);                                                                                                   ⋯                                                                           ;                                                                                                   ⋯                                                                           ;                                                                                                   ⋯                                                                           ;                        1             x_1(t_m-1;k)             x_2(t_m-1;k)                        ⋯             x_n(t_m-1;k)           x^2_1(t_m-1;k) x_1(t_m-1;k)x_2(t_m-1;k)                        ⋯           x^2_n(t_m-1;k);                          ]

    
and the velocity matrix is denoted as:

    V^(k)=[   ẋ_1(t_0;k)   ẋ_2(t_0;k)            ⋯   ẋ_n(t_0;k);   ẋ_1(t_1;k)   ẋ_2(t_1;k)            ⋯   ẋ_n(t_1;k);   ẋ_1(t_2;k)   ẋ_2(t_2;k)            ⋯   ẋ_n(t_2;k);                                      ⋯             ; ẋ_1(t_m-1;k) ẋ_2(t_m-1;k)            ⋯ ẋ_n(t_m-1;k) ] .


Let c_j be the vector of coefficients corresponding to the j^th governing equation, 1≤ j ≤ n, and define the coefficient matrix:

    C= [   |   |       |; c_1 c_2   ⋮ c_n;   |   |       | ].

Using the k^th burst data leads to the subproblem: find c such that V^(k) = A^(k) C. Next, we combine the data over all bursts k from 1,…, K  by simply concatenating the burst arrays vertically: 

    V = [ V^(1); V^(2);     |; V^(K) ]and
    A = [ A^(1); A^(2);     |; A^(K) ]

and thus the full problem is to find C such that V = A C. Let N:=n^2+3n+22   be the number of monomials of n variables up to degree two, then the size of matrix A is mK × N with mK<N therefore the linear inverse problem is ill-posed. Since the coefficient matrix C is sparse, the inversion could be regularized by solving the non-convex optimization problem:

    min_C  ||C||_0   subject to  AC-V≤σ,

where ||· ||_0 is the ℓ^0 penalty, i.e. the number of nonzero terms and the norm ||· || is the maximum of the ℓ^2 norm of each row. In this way, the inversion is component-wise separable.  The parameter σ>0 controls for the error between the computed derivative and the true velocity. For any fixed value of σ, the general ℓ^0 penalized problem is NP hard <cit.>; and thus we use the convex relaxation known as the ℓ^1 basis pursuit method:

    (M-BP_σ):  min_C  ||C||_1   subject to  AC-V≤σ.

Note that if V and A are given with high accuracy then we can solve the following:

    min_C  ||C||_1   subject to  AC=V.


The procedure and optimization problem above will be referred to as the monomial basis pursuit M-BP, since the trial set contains monomials up to degree two. We can also repeat the process using the tensorized quadratic Legendre polynomials, which corresponds to changing Equation (<ref>) to: 

    A_L^(k) =[                         1            √(3)x_1(t_0;k)                         ⋯            √(3)x_n(t_0;k)      √(5)3x^2_1(t_0;k)-12    3x_1(t_0;k) x_2(t_0;k)                         ⋯;                         1            √(3)x_1(t_1;k)                         ⋯            √(3)x_n(t_1;k)      √(5)3x^2_1(t_1;k)-12     3x_1(t_1;k)x_2(t_1;k)                         ⋯;                         1            √(3)x_1(t_2;k)                         ⋯            √(3)x_n(t_2;k)      √(5)3x^2_1(t_2;k)-12     3x_1(t_2;k)x_2(t_2;k)                         ⋯;                                                                             ⋯                                                                              ;                         1          √(3)x_1(t_m-1;k)                         ⋯          √(3)x_n(t_m-1;k)   √(5) 3x^2_1(t_m-1;k)-12 3x_1(t_m-1;k)x_2(t_m-1;k)                         ⋯;                           ]

          
The resulting inverse problem will be referred to as the Legendre basis pursuit L-BP, and we write it concisely as:

    min_C_L ||C_L||_1   subject to  A_L C_L=V.

 Note that sparsity of the solution of L-BP can be slightly larger than M-BP, since the `pure' quadratic terms now include a constant term.  In particular, if a component of the governing equation is s-sparse with respect to A, then it is at most (s+1)-sparse with respect to A_L. A noise-robust extension of this is given by:

    (L-BP_σ):  min_C_L  ||C_L||_1   subject to  A_LC_L-V≤σ.





§ RECONSTRUCTION GUARANTEE

Compressive sensing theory provides reconstruction guarantees for sparse solutions to ill-posed linear inverse problems via the solution to ℓ_1 optimization problems like L-BP – let us recall the basics of this theory.  Consider a general system of linear equations, y = A x.  If A is underdetermined, or has fewer rows than columns, it is in general impossible to determine x given only y and A.   However, if we know a priori that x ∈ℝ^N is s-sparse, or only has s ≪N non-zero entries, the locations of which are unknown, and if the underdetermined matrix A is suitably incoherent, it is possible to recover such an x from only M = O(s logN) measurements as the unique solution to the ℓ_1-minimization problem x = min_z  z _1  subject to A z = y.    Moreover, if there is noise on the measurements y = Ax + η, then x is well-approximated by the solution x^#∈min_z  z _1  subject to  A z - y _2 ≤σ for σ appropriately chosen. Thus, by leveraging sparsity, we can potentially overcome the curse of dimensionality in the number of measurements we need to take, which often renders high-dimensional inverse problems intractable.  

In our setting, the matrix of coefficients C (equivalently, each of the n columns c_1, c_2, …, c_n) is s-sparse, and so we can ask: under what conditions on the measurement  matrix A (or the Legendre-transformed A_L) can we apply compressive sensing results and conclude that C is exactly recovered as the solution to M-BP_σ (or L-BP_σ)?  As it turns out, if each of the K bursts is initialized at a uniformly random point in [-1,1]^n, then the K× N measurement matrix A_L corresponding to only the K initialization measurements, has precisely the incoherence properties that provide optimal compressive sensing results.  Thus, if each of the coefficient vectors c_k is s-sparse, then we only need to measure on the order of K ∼ s log N bursts to recover the coefficients exactly.  Each burst need only be measured long enough to get an accurate approximation to the initial velocity. Informally, the size of the burst should be large enough to get a stable approximation to the velocity via a finite difference approximation.

The remarkable fact that sparse vectors can be exactly recovered from vastly underdetermined linear systems of equations cannot be possible for just any underdetermined system of equations Ax = y – for example, the M × N matrix  A consisting of the first M rows of the N × N identity matrix maps all sparse vectors whose support does not intersect the first M coordinates to its null space, rendering it impossible to distinguish them. Indeed, the matrix A must have the incoherence property which implies that its null space only intersects the set of sparse vectors trivially, and is sufficiently “well-conditioned" over the set of sparse vectors to permit stability to noise.   The Legendre-transformed matrix A_L with burst length m = 0 and with initializations x(t_0;1), …, x(t_0;K) taken as independent and identically distributed as uniform random variables [-1,1]^n satisfies these requirements; in particular, it has two key properties which permit theoretical results on sparse recovery: 


  * Its rows ω_1, ω_2, …, ω_K are independent realizations of an isotropic random variable ω; that is, its covariance matrix is the identity matrix 𝔼[ωω^*] =  1_N × N,

  * Its rows are uniformly bounded: max_1 ≤ i ≤ N, 1 ≤ j ≤ K | A_L(i,j) |^2 ≤ 9.

Using these properties, we may apply Theorem 1.2 from <cit.> (see also Theorem 12.22 in <cit.>) to arrive at the reconstruction guarantee.



Consider a dynamical system ẋ=f(x) to be recovered from snapshots 
x(t_0; k), … x(t_m-1; k) and corresponding velocities ẋ(t_0; k), …ẋ(t_m-1;k), k=1,2,…, K.   Assume that each component f(x) = (f_1(x), f_2(x), …, f_n(x)) is a quadratic vector-valued equation in x, and that each of the f_k has at most s out of N polynomial coefficients non-zero.  Suppose that the initialization x(t_0;k)
for each of the k = 1,2,… K bursts is chosen independently at random from the uniform distribution over [-1,1]^n.  
Suppose the number of bursts satisfies

    K    ≥ 9 c_* s log(N) log( ε^-1)

where c_* is a universal constant.
Then with probability 1-ε, any particular component of the n governing equations in the system ẋ=f(x) is recovered exactly via the polynomial coefficients as the unique solution to (L-BP).



More generally, under the same conditions as above and with the same probability, if the measured gradient terms are only approximate,
 
    ẋ(t_0;k) = ẋ(t_0; k) + τ_0;k, … ,  ẋ(t_m-1;k) = ẋ(t_m-1; k) + τ_m-1;k,
 
such that √(1/K∑_k=1^K |τ_0;k|^2)≤η, then considering (L-BP_σ) with matrix A_L consisting of only the initial burst data and with σ = √(K)η, any particular coefficient vector c_k is approximated by a minimizer c_k^# of  (L-BP_σ) according to

    c_k^# - c_k _2  ≤ c^* √(s)η

where c^* is a universal constant.

The proof of Theorem  <ref> is deferred to the Appendix.  We note that by the union bound (Boole's inequality), if we observe K' ≥ 9 c_* s log(N) log( nε^-1) bursts, then with probability 1 - ε, the error bounds are guaranteed to hold uniformly for each of the n governing equations. Also, the arguments require a uniform bound on the elements of the matrix A_L restricted to the initial data, which can be chosen to be uniform random variables on [-1,1]^n. In general, trajectories of polynomial governing systems can become unbound in finite-time, which would make a uniform bound on A_L not possible. However, with the burst approach, one should expect the trajectory to exist for a short time. 


  This sparse recovery result can be generalized in several ways; namely,


  * The initializations should be independent and identically distributed random variables, but it is not crucial that they be  uniform random variables on [-1,1]^n.  Using results from <cit.>, the initializations could instead be taken to be i.i.d. according to the Chebyshev measure on [-1,1]^n, any measure interpolating between the Chebyshev measure and the uniform measure. 

  * The cubic case is an important structure that can arise from special symmetries. Although it is not directly considered here, our approach can be applied to extract cubic governing equations from data. Theorem <ref> holds if we modify Equation (<ref>) to:

    K    ≥ 27 c_* s log(N) log( ε^-1)

 since one can show that the matrix A_L restricted to uniformly random measurements satisfies the bound max_1 ≤ i ≤ N, 1 ≤ j ≤ K | A_L(i,j) |^2 ≤ 27. However, in the cubic case the number of unknowns N can be quite large since N = n +3  3 compared to the quadratic case where  N = n +2  2. Therefore, the number of basis elements N will be extremely large for high-dimensional cubic systems, rendering the ℓ_1 reconstruction algorithm quite costly without taking into account additional assumptions on the sparsity structure. Constructing a computational efficient approach under reasonable structural conditions will be explored in a future work.
 
  * We could have derived reconstruction guarantees for higher-order polynomial systems.  As with the cubic case, the only difference in the theoretical result is that the uniform bound max_1 ≤ i ≤ N, 1 ≤ j ≤ K | A_L(i,j) |^2 will grow with the maximal polynomial degree; the constant 9 in Equation (<ref>) will increase accordingly.    

  * We are not even restricted to consider governing equations which are sparse with respect to the polynomial basis; one could consider a different orthonormal basis such as e.g., sines and cosines which are uniformly bounded, and take the initializations of the bursts as i.i.d. random variables according to the orthogonalization measure (or a measure which is “close to" the orthogonalization measure) of that basis.
 






§ ERGODICITY AND THE NUMBER OF BURSTS


Theorem <ref> says that any dynamical system described as a system of sparse quadratic ordinary differential equations can be recovered exactly from a small number of randomly initialized bursts.  In particular, the number of bursts need only scale with the sparsity level of the system, and only logarithmically with the ambient dimension.  This result is true for any such dynamical system, independent of the behavior of the trajectories along the bursts.  In this sense, it is a “worst-case"  result.  In many situations, the data x(t_0), x(t_1), …, x(t_m-1) along a single burst behaves chaotically, mimicking the behavior of a random sequence, and in such cases the number of bursts actually required to achieve exact recover should be far smaller – as our numerical evidence shows, even a single burst often suffices.


While theoretical results concerning the behavior of high-dimensional dynamical systems have remained elusive, recent large-scale simulation studies, such as <cit.>  demonstrate that high-dimensional dynamical systems described by polynomial systems of equations often exhibit chaotic behavior; in fact, such behavior becomes more and more “probable" as the dimension of the system increases.   If we are in the regime of chaotic behavior, and if we measure snapshots of the system at the time scale of the chaotic dynamics, then for a smooth function F, the sequence F(x(t_0)), F(x(t_1)), …, F(x(t_m-1)) will satisfy a deterministic form of the law of large numbers – the so-called Birkhoff Ergodic Theorem – which says that the time average of the sequence will converge to the space average of F(x) with respect to an underlying invariant measure <cit.>.  For strongly ergodic systems, a stronger Central Limit Theorem holds, and the convergence rate can be quantified.  In particular, the rate of convergence is faster if the correlation between successive observed values F(x_k) and F(x_k+q) (for index q>0) is smaller, starting with the works <cit.>.  There is a rich theory on the decay of correlations for certain classes of low-dimensional dynamical systems, beyond the scope of this article. Morally speaking, if we measure a single burst which has the property that a smooth functional applied to the sequence exhibits a fast decay of correlation, then the resulting matrix A_L will have rows which are only weakly correlated, and thus, should almost fit the theoretical requirements for a compressive sensing result of the form of Theorem <ref>.  It is worth noting that a more recent paper <cit.> generalizes the results from  <cit.> to accommodate matrices with only weakly-correlated rows.  However, even so, the existing compressive sensing theory does not extend to this situation exactly as the sparse signal to be recovered is not independent of the measurement matrix.   Exploring a more precise relationship between the level of ergodicity in the system to recover and its effect on reducing the required number of bursts remains an intriguing direction for future study. 



§ SAMPLING STRATEGIES AND COMPUTATIONAL RESULTS


In this section, we present three sampling strategies that will lead to exact recovery (with high probability). In particular, Strategy 1 uses the construction from Section <ref>, Strategy 2 is a reduction of Strategy 1 when one can limit the number of unknowns to an ℓ-sized neighbor (i.e. when f_j only depends on  x_i, for i∈ [j-ℓ-1/2, j+ℓ-1/2] for odd ℓ), and lastly Strategy 3 shows that chaotic trajectories greatly reduce the number of random initializations. It is worth noting that Strategy 1 makes no assumptions on the trajectories x(t) or the support of the monomial basis for f(x), Strategy 2 incorporates information on the support of the monomial basis for f(x) but makes no assumptions on the trajectories x(t), and Strategy 3 uses information about the trajectories x(t). In practice, we expect that some combination of these three approaches could be optimal when dealing with various types of data. Note that although we will show examples with quadratic governing equations,  the number of unknowns, N, is non-trivial and grows quickly with dimension.


In all cases, we have under-sampled the data and it is easy to check that standard methods, like the least-squares algorithm, do not produce meaningful results. To validate our approaches, we will apply these strategies to the Lorenz 96 system and a quadratic reaction-diffusion equation. The Lorenz 96 system, introduced in <cit.> as an atmosphere model, contains n>3 variables x_1,⋯, x_n and satisfies

    _k = -x_k-2  x_k-1 + x_k-1  x_k+1 - x_k + F,   k=1,⋯, n.

The constant F is independent of k and x_0=x_n and x_n+1=x_1.  In both the standard monomial basis and in the Legendre monomial basis, the system is 4-sparse (component-wise). Unless otherwise stated, the number of variables is fixed to be n = 50 and the constant is set to F=8, where chaotic behavior is expected, see  <cit.>.


The second system we consider is known as the Fisher's equation, a quadratic reaction-diffusion equation, with reaction term F(x)=x-x^2. This equation has applications ranging from population dynamics to combustion physics. Its finite difference discretization of the Fisher's equation with n nodes is given by:

    _k = x_k+1-2x_k + x_k-1 + γ (x_k-x_k^2),   k=1,…, n,

where the coefficients and time-scale have been re-scaled by the grid-spacing and we impose the periodicity condition: x_0=x_n and x_n+1=x_1. Reaction-diffusion systems have multi-scale phenomena and can produce traveling waves, for more details see, for example, <cit.>. In our examples, the solutions remain bounded since the initial data is sampled in the unit box. Note that when transformed to the Legendre basis, each component of Equation (<ref>) is 5-sparse.


Throughout this section, we denote K to be the number of initializations, m is the size of each burst. For each burst, the time-derivative is approximated using central differencing except at the first and last time-step where we use forward and backward differencing (respectively). The dynamical data is generated by solving the ODEs numerically using Runge-Kutta 45 with a finer tolerance than the sampled time-step dt. To solve the L-BP problem we use the spgl1 algorithm <cit.> for all examples here. Other algorithms that could be successfully applied include: primal-dual <cit.>, Douglas-Rachford <cit.>, SpaRSA <cit.>, or the convex optimization package cvx <cit.>.





 §.§ Strategy 1: K ∼ c  s log(N) Random Initializations
 

A direct consequence of Theorem <ref> is that one can recover the governing equation with K random initializations as long as K≥ c  s log(N). We first determine the number of uniformly random initializations based on the theorem and collect a burst of size m generated from the system starting with each initialization. The data matrix is concatenated vertically over each burst. The parameters used are specified in each experiment. 

   



First, we validate the recovery results by varying the number of initialization K and measuring the probability that the governing system can be recovered. To test this we use the 10th component  of the Lorenz 96 system (see Equation (<ref>)):

    dx_10dt = -x_8  x_9 + x_9  x_11 - x_10 + 8.

 For each K=5, …, Nm=1326m, we repeat the simulation 100 times and record the number of successes to calculate the probability p, see Figure <ref>. The maximum K is chosen so that we do not oversample the system. From Figure <ref>, we can see that for dt = 0.001, K=80 is needed to achieve 90% probability of success for both m = 5 and m = 10. For a larger time-step, dt = 0.01, K=130 is needed to achieve 90% probability of success for m = 5. This is expected since the larger time-steps yield less accurate time-derivative approximations. Note that the recovery degrades when the burst size (in terms of m and dt) increases due to the propagation of error between time-steps. When the data is less correlated in time, the degradation may not occur, see Strategy 3. It is worth noting that although the theoretical universal constant c_* in Equation (<ref>) has a large upper bound, it is much smaller in practice.

 We also verify the recovery results for Fisher's equation for various γ. Table <ref> displays the values of the coefficients (as well as their support) when applying the proposed approach to the first component of Equation <ref>. The initial data is sampled from the uniform distribution over [0,1]^n and the full data matrix is transformed to [-1,1]^n. We set n=200 (N= 20301), the number of random samples to 159, and the size of the burst is 5 – this corresponds to the bound provided in Theorem <ref>. By varying the model parameter γ, we show that the method selects the correct basis terms and accurately approximates the parameter values (within a few significant digits). It is worth highlighting that as γ decreases, the relative scale between the maximum and minimum coefficient increases; however, the recovery results remain relatively stable.

 


Note that once we have identified the support S_L with respect to the Legendre dictionary matrix, it is possible to “correct” the results from Table <ref> by solving:

    min_CA_M|_S C|_S - V_2

where A_M|_S is the monomial dictionary A_M restricted to the corresponding support set S.  This debiasing step, when applied to the examples from Table <ref>, produces exact results (up to the fourth significant digit).



 §.§ Strategy 2: K < c  s log(N) Random Initializations with Localization
 


For many ODEs, especially those related to finite dimensional approximations of local PDEs, it is safe to assume that each variable x_j only relates to its ℓ-neighbors x_i, i∈ [j-ℓ-1/2, j+ℓ-1/2] (ℓ odd). In particular, the assumption is that the governing equation satisfies:

    ẋ_j = f_j(x_j-ℓ-1/2,…, x_ j+ℓ-1/2),

for all j. With this additional assumption, the number of unknowns in the dictionary matrix is reduced, and thus we can decrease the number of initial conditions needed to guarantee exact recovery. In particular, with the same conditions in Theorem <ref>, if additionally it is known that the sparse support is restricted to an ℓ-sized neighborhood, then if we sample the initial data K times, with K∼ c  s log(ℓ)log(ε^-1), then with probability 1-ε, the system ẋ=f(x) is recovered exactly by the unique solution to (L-BP).


Incorporating this additional condition simply amounts to downsampling the column space of the dictionary A. In Table <ref>, we consider the n=1000 dimensional Fisher's equation, which yields 501501 unknowns.
By varying the size of the neighborhood, ℓ, we calculate the minimum number of random initial samples needed for exact recovery. To validate the theoretical scaling, the ratio between the samples needed versus the log of the neighborhood size is shown to be nearly constant. Since the sampling rate is related to the log of the number of unknowns, we see substantial gains between using all terms versus restricting to a neighborhood, but further refinement is not needed. This highlights the benefit of restricting the optimization to an ℓ-sized neighborhood, even when ℓ can only be estimated. One key consequence from this strategy is that the sampling rate is independent of the dimension of the ODE system n, allowing this approach to be applied to very large systems.






 §.§ Strategy 3: Chaotic Systems with K Small and m Large
 

In Strategies 1 and 2, we made no assumption on the behavior of the trajectories x(t).  For Strategy 3, we will assume that the data exhibits chaotic behavior, or uncorrelated long-time behavior. Here we show that if this is the case, we can reduce the number of random initial data needed. In particular, we will show an example of the extreme case where K=1, since including additional trajectories with random initial data will only improve the recovery.  Consider collecting data along one trajectory: 
{ x(t_0), x(t_1),  …, x(t_m-1)} with time-steps dt that are  large enough so that x(t_i) and x(t_j) are sufficiently uncorrelated for i<j. The velocity along the trajectory {ẋ(t_0), ẋ(t_1), …, ẋ(t_m-1)} is either observed directly (possibly with some error) or calculated by using a fine time-step (smaller than dt). In either case, the data is under-sampled, i.e., m<N. In essence, taking large enough time-steps of a chaotic system mathematically resembles a random “re-sampling" of the data, thus fitting in with Strategies 1 and 2.

We test this strategy on the first component of the Lorenz 96 system, Equation (<ref>), with n=50. In Table <ref>, we show that with 500 measurements and dt = 1.0, the solution of L-BP identifies the correct terms and approximates the coefficients within the expected error. This shows that it is possible to use the randomness of the data along one trajectory to learn the governing equation. Adding more trajectories with initial data sampled i.i.d. from the uniform distribution, while keeping the total size of the data fixed, only helps the recovery processes. 











 §.§ Comparisons

For comparison, we apply Strategy 1, the standard least-square algorithm, and the linear regression method with sequential thresholding proposed in <cit.> on both the Lorenz 96 and Fisher's equation. In Figure <ref>, the coefficients extracted using the L-BP method (left), the least-square algorithm (middle), and the sequential thresholding algorithm (right) for the 35th component of the Lorenz 96 equation with n=50, dt = 0.001. The threshold parameter for the least-square algorithm and the sequential thresholding algorithm is set to λ=0.05. In fact, the solutions of the least-square and the linear regression are the same for this problem. This is the case for all reasonable λ>0.

 





Next, we compare the recovery of the coefficients from the first component of the Fisher's Equation (<ref>). The model parameter is set to γ=0.1 and dimension is set to n=100, the number of random samples is set to 138, and the size of the burst is set to 5. The least-squares solution is 653-sparse and has coefficients on the order of 10^4. Applying the sequential thresholding algorithm proposed in <cit.> with threshold parameter λ∈ [5, 5000] results in a similar solution. The sparsity matches that of the least-square solution, with s∈[ 579, 653] and has coefficients on the order of 10^4. Increasing λ to achieve a better sparsity level will yield the trivial solution (i.e. all zeros). This is likely an effect of their algorithm's dependence on the least-square solution. It is important to note that the algorithm proposed in <cit.> was not intended for the case of under-sampling.





 §.§ Measurement Noise



We consider the effects of noisy state-space measurements on the reconstruction of the sparse coefficient vector. Let X be the data matrix:


    X =[   x_1(t_0;1)   x_2(t_0;1)            ⋯   x_n(t_0;1);   x_1(t_1;1)   x_2(t_1;1)            ⋯   x_n(t_1;1);                                      ⋯             ; x_1(t_m-1;1) x_2(t_m-1;1)            ⋯ x_n(t_m-1;1);                                      ⋯             ;   x_1(t_0;k)   x_2(t_0;k)            ⋯   x_n(t_0;k);                                      ⋯             ; x_1(t_m-1;k) x_2(t_m-1;k)            ⋯ x_n(t_m-1;k);                                      ⋯             ; x_1(t_m-1;K) x_2(t_m-1;K)            ⋯ x_n(t_m-1;K) ]


and let Y = X + η  be the matrix of noisy measurements, where η is random Gaussian noise. The noise ratio is defined as:

    Noise Ratio =  X-Y_2/ X_2× 100%

and the relative ℓ^2 error is define as:

    relative ℓ^2 error =  c-c_true_2/ c_true_2× 100%,

where c is the computed/learned vector of coefficients and c_true is the true coefficients. 
We generate the data from Equation (<ref>) using F=8 and dimension equal to 50. The parameters are set to dt=0.001,  K=200, and m=3. By adding noise directly to the state-space, we can measure the recovery when the measurements are corrupted. The state variable is corrupted by random Gaussian noise before the derivatives are calculated, thus the value of V and the matrix A will be inaccurate. This makes the problem challenging when the noise is large. The results are summarized in Table <ref>. As we vary the noise level, we measure the relative ℓ^2 error and the recovery of the support set.  For noise under 5%, the method is stable with respect to the noise. This is consistent with Theorem <ref>. After 5% noise, the ℓ^2 error jumps and after 6%, we cannot reliably recover the support set. In particular, since we know that the true sparsity of Equation (<ref>) is 4,  we can check if the largest four learned coefficients (in magnitude) coincide with the correct support set. After about 6%, the largest four coefficients do not represent the correct support set. This is likely due to the large inaccuracies in V, which is not stable to the noise, and inaccuracies in A, which scales nonlinearly with the noise.




  
 


§ CONCLUSION AND DISCUSSION

Extracting dynamical systems remains a difficult task with many open areas of research. Recent work in sparse model selection for dynamical systems has focused on the overdetermined case, where regression must be controlled so as to prevent overfitting. In this work, we utilized the  fundamental idea from compressive sensing to develop several sampling strategies for extracting governing equations from high-dimensional dynamic data. In all cases, the number of measurements is less than the number of unknowns. The main differences between these strategies is the degree of prior knowledge about the data or the governing equations.  If no assumptions on the evolution can be made, then randomizing the initial data is sufficient. If some assumptions on the governing equations are provided, such as locality induced by discretizing a local PDE, then the number of random initial data can be reduced further, to be nearly independent of the dimension of the dataset.  In the third case, if the data is chaotic (or shows a low temporal correlation), then we may reduce the number of initial samples to a fixed number. Using results from compressive sensing, the first two strategies are shown to hold; however, more theory is needed to verify the last case. In several of our experiments,  we have shown that the three strategies are robust to various factors as well as highlighted the benefits of this approach over existing methods. The effective combination of reconstruction guarantees from compressive sensing with sparse learning for dynamical systems presented here opens a wide range of applications. We are currently investigating the use of group sparsity and random sampling approaches for learning dynamic models from multiple data sources. In future work, we also would like to extend the current framework to other bounded orthonormal bases and to learn the dynamics from noisy data.





§ ACKNOWLEDGMENTS

The authors would like to thank Amit Singer and Scott McCalla for helpful discussions that improved this manuscript. H.S. acknowledges the support of AFOSR, FA9550-17-1-0125. R.W. and G.T. acknowledge the support of NSF CAREER grant #1255631. 



§ APPENDIX


For purposes of being as self-contained as possible, we first recall some background on the theory of sparse recovery in bounded orthonormal systems via random sampling.  We refer the reader to the text <cit.> for more details.




 §.§ Random Sampling in Bounded Orthonormal Systems


Let D⊂ℝ^n be endowed with a probability measure μ.  Suppose that {ϕ_1, ϕ_2, …, ϕ_d } (d ≤ n) is a (possibly complex-valued) orthonormal system  on D:

    ∫_ Dϕ_j(t) ϕ_k(t) dμ(t) = δ_j,k   =  {[        0 if j ≠ k;        1   if j=k;          ].

We call {ϕ_1, ϕ_2, …, ϕ_d } a bounded orthonormal system with constant B ≥ 1 if moreover 

    ϕ_j _∞ := sup_t ∈ D | ϕ_j(t) | ≤ B   for all j ∈ 1,2,…, d.

Suppose that t_1, t_2, …, t_m ∈ D are sampling points which are drawn i.i.d. according to the orthogonalization measure μ, and consider the sampling matrix 
A ∈ℂ^m × n with entries

    A_ℓ,k = ϕ_k(t_ℓ),   ℓ∈ [m], k ∈ [n].

With high probability, a random matrix formed as such permits stable “inversion" of the (possibly highly underdetermined) system y= Ax if x is sufficiently sparse, and if “inversion" is carried out through e.g. solving an ℓ_1-minimization problem.  The following is a restatement of Theorem 12.22 in <cit.>, which is a restatement of a result from <cit.>.



Let x ∈ℂ^n and let A ∈ℂ^m × n to be the random sampling matrix associated to a BOS with constant B ≥ 1.  For y = Ax + e with e _2 ≤η√(m) for some η≥ 0, let x^# be a solution to

    min_z ∈ℂ^n z _1   subject to  Az - y _2 ≤η√(m).

If 

    m ≥ C B^2 s log(n) log(ε^-1),

then with probability at least 1 - ε, the reconstruction error satisfies 

    x - x^#_2 ≤ C_1 σ_s(x)_1 + C_2 √(s)η

where σ_s(x)_1 = inf_u: u  is s-sparse x - u _1, and the constants C, C_1, C_2 > 0 are universal.

We apply this result to prove our reconstruction guarantee.



 §.§ Proof of Theorem <ref>

Recall that A_L is the K × N matrix consisting of those rows corresponding to the K initializations.  The tensor product of univariate Legendre polynomials, normalized as in  the construction of A_L, forms a bounded orthonormal system with respect to the uniform measure over [-1,1]^n; precisely, dμ = 1/2dx.  The Legendre polynomials up to degree 2 are uniformly bounded in magnitude over the domain [-1,1]^n by 3, as realized by the terms 3 x_j x_k at x_j, x_k = ± 1. Thus, A_L satisfies the requirements of Proposition <ref> with K = 3.   Moreover, by assumption, each column in the coefficient matrix C associated to the underlying dynamical system consists of at most s nonzero terms.  Thus, after transforming to the Legendre basis, each column in the the Legendre coefficient matrix C_L has at most s' = s+1 nonzero terms. Thus, we apply Proposition <ref> with these parameters to any particular one of the n ℓ_1-minimization problems. 

The recovery guarantee in the noiseless case (η = 0) moreover extends to the optimization problem L-BP over the full set of measurement constraints, not just those constraints corresponding to the burst initializations, since the domain corresponding to the full measurement set is strictly included in the domain corresponding to the reduced measurement set, and the unique minimizer over the reduced measurement set belongs to this subdomain. 
   



10

ABW17
Ben Adcock, Simone Brugiapaglia, and Clayton G. Webster.
Polynomial approximation of high-dimensional functions via compressed sensing.
arXiv preprint  arXiv:1703.06987, 2017.


anderson1992infectious
Roy M. Anderson and Robert M. May.
Infectious diseases of humans: dynamics and control, Oxford university press, 1992.



  berkooz1993proper
Gal Berkooz, Philip Holmes, and John L. Lumley. 
The proper orthogonal decomposition in the analysis of turbulent flows.
Annual Review of Fluid Mechanics, 25(1):539–575, 1993.

Birkhoff31
George Birkhoff.
Proof of the ergodic theorem.
Proceedings of the National Academy of Sciences, 17(12): 656–660, 1931.

bongard2007
Josh Bongard and Hod Lipson.
Automated reverse engineering of nonlinear dynamical systems.
Proceedings of the National Academy of Sciences, 104(24):9943–9948, 2007.

Bowen75
Rufus Bowen.
Equilibrium states and the ergodic theory of Anosov diffeomorphisms. 
Lect. Notes Math. 470.
Springer-Verlag, Berlin, 1975.

bright2013compressive
Ido Bright, Guang Lin, and J. Nathan Kutz.
 Compressive sensing based machine learning strategy for
  characterizing the flow around a cylinder with limited pressure measurements.
Physics of Fluids, 25(12):127102, 2013.


brunton2016
Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz.
Discovering governing equations from data by sparse identification of nonlinear dynamical systems.
Proceedings of the National Academy of Sciences, 113(15):3932–3937, 2016.

brunton2017chaos
Steven L. Brunton, Bingni W. Brunton, Joshua L. Proctor, Eurika Kaiser, and J. Nathan Kutz. 
Chaos as an intermittently forced linear system.
Nature Communications, 8, 2017.


caflisch2013pdes
Russel E Caflisch, Stanley J Osher, Hayden Schaeffer, and Giang Tran.
PDEs with compressed solutions.
Communications in Mathematical Sciences, 13(8):2155–2176, 2015.


candes2006
Emmanuel J. Candès, Justin K. Romberg, and Terence Tao.
 Stable signal recovery from incomplete and inaccurate measurements.
 Communications on Pure and Applied Mathematics, 59(8):1207–1223, 2006.



candes2007
Emmanuel Candès and Justin Romberg.
 Sparsity and incoherence in compressive sampling.
 Inverse Problems, 23(3):969–985, 2007.
 
 CP11
Emmanuel Candès and Yaniv Plan.
A probabilistic and RIPless theory of compressed sensing.
IEEE Transactions on Information Theory, 57(11):7235–7254, 2011.


chambolle2011first
Antonin Chambolle and Thomas Pock. 
  A first-order primal-dual algorithm for convex problems with applications to imaging.
  Journal of Mathematical Imaging and Vision, 40(1):120–145, 2011.

chen2001atomic
 Scott Shaobing Chen, David L. Donoho, and Michael A. Saunders. 
 Atomic decomposition by basis pursuit.
 SIAM Review 43(1):129–159, 2001.

coifman2005geometric
Ronald R. Coifman, Stéphane Lafon, Ann B. Lee, Mauro Maggioni, Boaz Nadler, Frederick Warner, and Steven W. Zucker. 
Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps.
Proceedings of the National Academy of Sciences of the United States of America, 102(21):7426–7431, 2005.


coifman2006diffusion
Ronald R. Coifman and Stéphane Lafon. 
Diffusion maps.
Applied and Computational Harmonic Analysis, 21(1):5–30, 2006.


donoho1995noising
David L. Donoho.
 De-noising by soft-thresholding.
 IEEE Transactions on Information Theory, 41(3):613–627, 1995.

  donoho2006most
  David L. Donoho.
 For most large underdetermined systems of linear equations the minimal ℓ1-norm solution is also the sparsest solution.
 Communications on Pure and Applied Mathematics, 59(6):797–829, 2006.

eckstein1992douglas
Jonathan Eckstein and Dimitri P. Bertsekas.
On the Douglas–Rachford splitting method and the proximal point algorithm for maximal monotone operators.
Mathematical Programming, 55(1):293–318, 1992.



FR11
Simon Foucart and Holger Rauhut.
A mathematical introduction to compressive sensing.
Springer, 2011.

giannakis2015data
Dimitrios Giannakis. 
Data-driven spectral decomposition and forecasting of ergodic dynamical systems.
arXiv preprint arXiv:1507.02338, 2015.

grant2008cvx
Michael Grant, Stephen Boyd, and Yinyu Ye.
 CVX: Matlab software for disciplined convex programming, 2008.

 
 holmes2012turbulence
 Philip Holmes, John L. Lumley, Gahl Berkooz, and Clarence W. Rowley. 
  Turbulence, coherent structures, dynamical systems and symmetry. 
 Cambridge University Press, 2012.

hou2015sparse
 Thomas Y. Hou, Qin Li, and Hayden Schaeffer.
Sparse+ low-energy decomposition for viscous conservation laws.
Journal of Computational Physics, 288:150–166, 2015.

 chaos2
 Iaroslav Ispolatov,  Vaibhav Madhok, Sebastian Allende, and Michael Doebeli.
 Chaos in high-dimensional dissipative dynamical systems.
 Scientific Reports, 5, 2015.
 
keeling2011modeling
 Matt J. Keeling and Pejman Rohani. M
  Modeling infectious diseases in humans and animals, Princeton University Press, 2011.
 
 
  RIPless2
Richard  Kueng and David Gross.
RIPless compressed sensing from anisotropic measurements.
Linear Algebra and its Applications, 441: 110–123, 2014.


lorenz1996predictability
Edward N. Lorenz.
Predictability: A problem partly solved.
Proc. Seminar on Predictability, 1(1), 1996.

mackey2014
 Alan Mackey, Hayden Schaeffer, and Stanley Osher. 
 On the compressive spectral method.
 Multiscale Modeling & Simulation, SIAM, 12(4):1800–1827, 2014.
 
 
mangan2017model
Niall M. Mangan, J. Nathan Kutz, Steven L. Brunton, and Joshua L. Proctor.
Model selection for dynamical systems via sparse regression and information criteria.
arXiv preprint arXiv:1701.01773, 2017.

majda2009normal
Andrew J. Majda,  Christian Franzke, and Daan Crommelin. 
Normal forms for reduced stochastic climate models.
Proceedings of the National Academy of Sciences 106(10): 3649-3653, 2009.


mckean1975application
Henry P. McKean.
Application of brownian motion to the equation of Kolmogorov?Petrovskii?Piskunov.
Communications on Pure and Applied Mathematics 28(3):  323-331, 1975.


  chaos1
  Z. Musielak and D. Musielak. 
  High-dimensional chaos in dissipative and driven dynamical systems.
  International Journal of Bifurcation and Chaos, 19(9): 2823–2869, 2009.

 nadler2006diffusion
Boaz Nadler, Stéphane Lafon, Ronald R. Coifman, and Ioannis G. Kevrekidis. 
Diffusion maps, spectral clustering and reaction coordinates of dynamical systems.
Applied and Computational Harmonic Analysis, 21(1):113–127, 2006.

nadler2006diffusion2
Boaz Nadler, Stéphane Lafon, Ioannis Kevrekidis, and Ronald R. Coifman. 
Diffusion maps, spectral clustering and eigenfunctions of Fokker-Planck operators.
Advances in Neural Information Processing Systems. 955–962, 2006.

ozolinvs2013compressedCPW
Vidvuds Ozolinš, Rongjie Lai, Russel Caflisch, and Stanley Osher.
Compressed plane waves-compactly supported multiresolution basis for
  the laplace operator.
Proceedings of the National Academy of Sciences, 2013.

proctor2016dynamic
Joshua L. Proctor, Steven L. Brunton, and J. Nathan Kutz. 
Dynamic mode decomposition with control.
SIAM Journal on Applied Dynamical Systems 15(1): 142-161, 2016.


PHD16
Ji Peng, Jerrad Hampton, and Alireza Doostan.
On polynomial chaos expansion via gradient-enhanced ℓ_1-minimization. 
Journal of Computational Physics, 310: 440-458, 2016.
 

R10
Holger Rauhut.
Compressive sensing and structured random matrices.
Theoretical Foundations and Numerical Methods for Sparse Recovery, ed. by M. Fornasier. 
Radon Series on Computational and Applied Mathematics, 9:1–92, 2010.
 
   RWLegendre
  Holger Rauhut and Rachel Ward.
  Sparse Legendre Expansions via ℓ_1 minimization.
  Journal of Approximation Theory,  164(5):517–533, 2012.
  
  RWWeighted
  Holger Rauhut and Rachel Ward.
  Interpolation via weighted ℓ_1 minimization.
  Applied and Computational Harmonic Analysis, 40(2): 321-351, 2016. 

 RV08
 Mark Rudelson and Roman Vershynin.
 On sparse reconstruction from Fourier and Gaussian measurements.
Communications on Pure and Applied Mathematics, 61:1025–1045, 2008.

 
  rudy2017data
Samuel H. Rudy, Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. 
Data-driven discovery of partial differential equations.
 Science Advances, 3(4): e1602614, 2017.
 
Ruelle68
David Ruelle.
Statistical mechanics of a one-dimensional lattice gas. 
Commun. Math. Phys., 9:267–278. 1968.

Ruelle76
David Ruelle.
A measure associated with Axiom A attractors. 
Amer. J. Math., 98:619–654.  1976.

Sinai72
Yakov Sinai.
Gibbs measures in ergodic theory. 
Russ. Math. Surveys 27:21–69. 1972.

schaeffer2013sparse
Hayden Schaeffer, Russel Caflisch, Cory D. Hauck, and Stanley Osher.
Sparse dynamics for partial differential equations.
Proceedings of the National Academy of Sciences,
  110(17):6634–6639, 2013.


 
 schaeffer2017learning
Hayden Schaeffer.
Learning partial differential equations via data discovery and sparse optimization.
Proc. R. Soc. A, 473(2197):20160446, 2017.

schaeffer2017sms
Hayden Schaeffer and Scott G. McCalla.
Sparse model selection via integral terms.
To appear in Phys. Rev. E, 2017.

schmid2010dynamic
Peter J. Schmid. 
Dynamic mode decomposition of numerical and experimental data.
Journal of Fluid Mechanics 656: 5–28, 2010.

schmid2012decomposition
Peter J. Schmid, Daniele Violato, and Fulvio Scarano. 
Decomposition of time-resolved tomographic PIV.
Experiments in Fluids 52(6):1567–1579, 2012.

  
  schmidt2009
Michael Schmidt and Hod Lipson.
Distilling free-form natural laws from experimental data.
Science, 324(5923):81–85, 2009.

 


tibshirani1996
Robert Tibshirani. 
Regression shrinkage and selection via the lasso. 
Journal of the Royal Statistical Society. Series B (Methodological), 267–288, 1996.


tran2015penalty
Giang Tran, Hayden Schaeffer, William M. Feldman, and Stanley J. Osher.
An L^1 penalty method for general obstacle problems.
SIAM Journal on Applied Mathematics, 75(4):1424–1444, 2015.


tran2016exact
Giang Tran and Rachel Ward.
Exact recovery of chaotic systems from highly corrupted data.
Multiscale Modeling & Simulation, SIAM, 15(3):1108–1129, 2017.



van2008probing
Ewout Van Den Berg and Michael P. Friedlander.
Probing the Pareto frontier for basis pursuit solutions.
SIAM Journal on Scientific Computing, 31(2):890–912, 2008.
  

van2003front
Wim Van Saarloos.
Front propagation into unstable states.
Physics reports, 386(2-6): 29-222, 2003.



 








wright2009sparse
Stephen J. Wright, Robert D. Nowak, and Mário AT Figueiredo.
  Sparse reconstruction by separable approximation.
  IEEE Transactions on Signal Processing, 57(7):2479–2493, 2009.
  
  


 
 



williams2014data
Matthew O. Williams, Ioannis G. Kevrekidis, and Clarence W. Rowley. A data-driven approximation of the Koopman operator: extending dynamic mode decomposition. arXiv preprint arXiv:1408.4408, 2014.

  
  



