# Scalable Sparse Regression for Model Discovery: The Fast Lane to Insight

Matthew Golden<EMAIL_ADDRESS>School of Physics, Georgia Institute
of Technology

###### Abstract

There exist endless examples of dynamical systems with vast available data and
unsatisfying mathematical descriptions. Sparse regression applied to symbolic
libraries has quickly emerged as a powerful tool for learning governing
equations directly from data; these learned equations balance quantitative
accuracy with qualitative simplicity and human interpretability. Here, I
present a general purpose, model agnostic sparse regression algorithm that
extends a recently proposed exhaustive search leveraging iterative Singular
Value Decompositions (SVD). This accelerated scheme, Scalable Pruning for
Rapid Identification of Null vecTors (SPRINT), uses bisection with analytic
bounds to quickly identify optimal rank-1 modifications to null vectors. It is
intended to maintain sensitivity to small coefficients and be of reasonable
computational cost for large symbolic libraries. A calculation that would take
the age of the universe with an exhaustive search but can be achieved in a day
with SPRINT.

††preprint: APS/123-QED

##  I Introduction

There has been a recent explosion in the application of model discovery
algorithms to learn symbolic governing equations directly from temporal and
spatiotemporal data. These methods rely on some form of symbolic regression
[1, 2, 3, 4] to fit data with a set of allowed symbolic operations (i.e.,
prescribed nonlinear functions or linear summation). Sparse regression applied
to a symbolic library has been particularly successful at recovering dynamics
[3, 5, 6, 4, 7, 8, 9, 10]. In these methods, the dynamics of a system
$\dot{\bf z}$ are assumed to be a sparse linear combination of nonlinear
functions $\mathcal{L}=\\{{\bf f}_{n}\\}.$ A sparse coefficient vector ${\bf
c}$ balances model simplicity with quantitative accuracy:

$\dot{\bf z}=\sum_{n}c_{n}{\bf f}_{n}({\bf z}).$ (1)

Ideally, ${\bf c}$ should be chosen to be maximally descriptive while
containing no irrelevant information. $\mathcal{L}$ is typically made up of
polynomials or simple trigonometric functions. This is not a particularly
restrictive assumption: all data is bounded and all smooth functions on
bounded subsets of $\mathbb{R}^{n}$ are approximated arbitrarily well by
polynomials [11].

Model discovery has been extended to partial differential equations (PDEs)
using spatiotemporal data [6, 12, 8, 9, 10]. In these libraries, spatial
gradients also enter the library. If there is a known symmetry of the system
(rotational, Galilean, etc.) then $\mathcal{L}$ should contain terms that
transform identically to $\dot{\bf z}$ under the action of the symmetry group
[12, 13, 10]. Such symmetry considerations significantly reduce the size of
the library and improve the robustness of sparse regression.

A prevalent sparse regression paradigm for learning the dynamics is the Sparse
Identification of Nonlinear Dynamics (SINDy) algorithm. In this approach, the
vector is identified numerically by constructing a system of equations ${\bf
G}{\bf c}={\bf b}$ using empirical data. The $n$th column of the feature
matrix ${\bf G}$ corresponds to a set of observations of ${\bf f}_{n}$ and
${\bf b}$ corresponds to observations of $\dot{\bf z}$.

$\displaystyle{\bf G}=\begin{pmatrix}|&|&|\\\ {\bf f}_{1}&{\bf f}_{2}&{\bf
f}_{3}&\cdots\\\ |&|&|\end{pmatrix},$ $\displaystyle{\bf
b}=\begin{pmatrix}|\\\ \dot{\bf z}\\\ |\end{pmatrix}.$ (2)

In early model discovery efforts, ${\bf G}$ was constructed via pointwise
evaluation of these library terms [3, 6]. Pointwise sampling is susceptible to
numerical error or experimental noise, especially for terms containing
derivatives. A more robust evaluation is in the weak formulation in which
terms are evaluated as integrals [5, 7, 14, 9, 10]. While the vector ${\bf b}$
usually corresponds to the time derivative $\dot{\bf z}$, other physically
motivated terms can be used instead [9]. Once this linear system is
constructed, a sparse approximate solution ${\bf c}$ is generated using
sequential thresholding: the linear system ${\bf G}{\bf c}={\bf b}$ is solved
iteratively with coefficients set to zero if they are less than a threshold
$|c_{n}|<\epsilon$, for some small $\epsilon$ [3, 15]. In practice, sequential
thresholding needs surprisingly few iterations to converge.

This paradigm of sparse regression exactly minimizes a loss $\|{\bf G}{\bf
c}-{\bf b}\|_{2}$ at each stage and then makes sparsity promoting
modifications (thresholding) to arrive at a sparse coefficient vector ${\bf
c}.$ This is in constrast to sparse regression methods like LASSO and elastic
net, which modify the loss to contain sparsity-promoting penalties $\|{\bf
G}{\bf c}-{\bf b}\|_{2}+\lambda_{1}\|{\bf c}\|_{1}+\lambda_{2}\|{\bf c}\|_{2}$
[16, 17]. The latter approaches inevitably introduce hyper-parameters that
must be tuned to obtain the desired output. Furthermore, the alignment of the
$\ell_{1}$ sparsity-promoting penalty with the true goal of sparse regression
is questionable. Why should the magnitude of coefficients matter? An
$\ell_{0}$ penalty is certainly more desirable [18], although this problem is
NP hard [19] and tuning of a hyperparameter is still required.

Sparse regression methods face three generic difficulties: i) in principle the
dynamics of a system may not be described by equations of the form (1) [9,
10]; ii) the symbolic expression for ${\bf b}$ is not unique if there is a
degeneracy in the library; and iii) thresholding can discard small
coefficients like dissipation that are important for understanding the
dynamics [7]. The first two problems are addressed in part by implicit-SINDy
and SINDy-PI [20, 4] where each column of ${\bf G}$ is solved for
independently.

A robust Singular Value Decomposition (SVD) based alternative was introduced
during the development of the Sparse Physics-Informed Discovery of Empirical
Relations (SPIDER) algorithm [12, 7, 10]. The sparse regression algorithm is
an exhaustive search for optimal, rank-one modifications to a coefficient
vector ${\bf c}$. This exhaustive search does not assume any right hand side
${\bf b}$, does not threshold coefficients, and can find multiple relations
thereby finding any degeneracy in the library. Instead of directly solving a
system of equations, it sequentially minimizes the homogeneous residual

$r({\bf c})\equiv\frac{1}{\sqrt{m}}\frac{\|{\bf G}{\bf c}\|_{2}}{\|{\bf
c}\|_{2}}$ (3)

for sparse ${\bf c},$ where $m$ is the number of rows in ${\bf G}.$ All of
these features can be seen directly in the example explained in Ref [7]: in
simulated 3D turbulence the exhaustive search was able to correctly identify
the pressure-Poisson equation, the energy equation, the incompressibility
condition, and the Navier-Stokes equation as well as boundary conditions. The
discovered equations include two spatial constraints rather than dynamical
relations, and recovered three of these relations from the same scalar library
displaying robustness to degeneracy. Furthermore, a viscous term with
coefficient $5\times 10^{-5}$ was identified accurately using data from two
dynamically distinct regions of the flow. In another application, the
exhaustive search with a modified residual found the governing equations for
an experimental active nematic suspension [10]. This application displayed
these same strengths: both dynamic relations and spatial constraints were
identified. Six libraries were searched to obtain nine PDEs, which logically
reduced to a minimal set of three governing equations.

While the exhaustive search has attractive model agnosticism and the ability
to indentify accurate small coefficients, it has displayed
$O(|\mathcal{L}|^{4})$ scaling with library size. Calculations for
$|\mathcal{L}|>10^{4}$ become prohibitively expensive. Here, I propose an
accelerated variant of this search algorithm: Scalable Pruning for Rapid
Identification of Null vecTors (SPRINT). This acceleration hinges on a simple
correspondence: adding or removing the coefficient $c_{k}$ from ${\bf c}$ is
equivalent to a performing a rank-1 update to the feature matrix ${\bf
G}\leftarrow{\bf G}\pm{\bf g}_{k}\otimes{\bf e}_{k}$. (Here ${\bf g}_{k}$ is
the $k$th column of ${\bf G}$ and ${\bf e}_{k}$ is a unit vector.) Finding the
modified singular values after a rank-1 update can be done efficiently with
bisection [21]. SPRINT comes in two flavors: SPRINT$\pm$. In SPRINT–, terms
are removed from ${\bf c}$ greedily to minimize $\|{\bf G}{\bf c}\|_{2}$ at
each iteration. In SPRINT+, one starts with a sparse guess and adds terms in
the same greedy fashion. It is unnecessary to add terms until the full library
is used, so in practice SPRINT+ is halted when a maximum number of terms
$\|{\bf c}\|_{0}$ is reached. In both cases, bisection is used to quickly
determine the optimal rank-1 modification to the coefficient vector ${\bf c}$.
SPRINT- and SPRINT+ empirically scale like $O(|\mathcal{L}|^{3.38})$ and
$O(|\mathcal{L}|^{1.65})$, respectively.

##  II The Library Catastrophe

Despite their success, library based methods have an unfortunate scaling with
allowed complexity of library terms. This scaling is significantly reduced by
symmetry considerations [7, 10], but I will ignore these to illustrate the
potentially catastrophic scaling of library size with complexity. Consider an
ideal Magnetohydrodynamic (MHD) system with eight fields: a velocity field
$(u_{x},u_{y},u_{z})$, a magnetic field $(B_{x},B_{y},B_{z})$, a density
$\rho$, and a pressure $P$. In this example, I will estimate the computational
requirement for model discovery in such a system. These fields can vary in all
four spacetime directions $(t,x,y,z)$, so one can take four different partial
derivatives $(\partial_{t},\partial_{x},\partial_{y},\partial_{z}).$ These
fields and derivatives form a symbolic alphabet $\mathcal{A}$:

$\mathcal{A}=\\{u_{x},u_{y},u_{z},B_{x},B_{y},B_{z},\rho,P,\partial_{t},\partial_{x},\partial_{y},\partial_{z}\\}.$
(4)

Each element of this alphabet will be a letter. A word is a combination of
letters, and its length is the total number of letters used in a word. For
example, $\rho\partial_{x}\rho$ is a 3-letter word. Since the number of
symbols is $|A|=12$, there can be at most $12^{n}$ $n$-letter words. If
$\mathcal{L}_{n}$ is the library of words length $n$ and less, then
$|\mathcal{L}_{n}|$ is bounded by a geometric sum of $|\mathcal{A}|^{n}$.

$|\mathcal{L}_{n}|\leq|\mathcal{A}|\frac{|\mathcal{A}|^{n}-1}{|\mathcal{A}|-1}\sim|\mathcal{A}|^{n}.$
(5)

In reality, both physicality conditions and permutation symmetries
(commutativity of multiplication and partial derivatives) makes
$|\mathcal{L}_{n}|$ smaller than this bound. The exact magnitude of the
library can be found by counting words that meet the following conditions:

  1. 1.

The final letter of a word cannot be a derivative: a derivative must act on
something. If the final letter is a field, then this word can be evaluated.

  2. 2.

The order that fields appear in a word should coincide with the ordering
$(u_{x},u_{y},u_{z},B_{x},B_{y},B_{z},\rho,P)$ such that $u_{y}u_{x}\wedge
u_{x}u_{y}$ are redundant.

  3. 3.

Partial derivatives acting on an individual field should coincide with the
canonical ordering $(\partial_{t},\partial_{x},\partial_{y},\partial_{z})$
such that $\partial_{x}\partial_{t}u_{x}\wedge\partial_{t}\partial_{x}u_{x}$
are redundant.

  4. 4.

Products of the same field with distinct derivatives, i.e.,
$\partial_{y}\rho\partial_{t}\partial_{x}\rho$ should be ordered in a
canonical way. Providing a canonical ordering is beyond the scope of this
paper, but any ordering on $\mathbb{Z}^{n}$ is sufficient.

The true scaling of $\mathcal{L}_{n}$ with $n$ for the MHD system discussed
here is shown in Figure 1. This $\mathcal{L}_{n}$ was computed exhaustively by
constructing every combination of $n$ (or fewer) symbols from $\mathcal{A}$
and keeping only those that satisfy the four physicality and uniqueness
conditions. By the modest size $n=5,$ one has a library with O$(10^{4})$
terms. Enforcing continuous symmetry covariance, such as rotational symmetry
[7], can reduce the library size substantially. However, large libraries can
be desirable if symmetry-breaking is suspected in the data or if no symmetry
is known a priori. Such symmetry breaking terms can appear with very small
coefficients, which highlights the need for a sparse regression method without
thresholding.

Figure 1: The size of a library $|\mathcal{L}_{n}|$ using 3+1D ideal MHD
variables. A 12-letter symbolic alphabet is made up of eight physical fields
and four partial derivatives. The library size is trivially bounded by
equation (5) in red. Permutation symmetries make many words redundant but this
does not prevent the catastrophic scaling of $|\mathcal{L}_{n}|$. A complete
library allowing only up to 4 letters per word will require thousands of
library terms.

##  III The SPRINT Algorithm

In the existing approach, the exhaustive search algorithm finds increasingly
sparse coefficient vectors ${\bf c}$ that minimize the residual $r({\bf c})$
defined in equation (3). Since $r({\bf c})$ is invariant under the rescaling
of ${\bf c}\rightarrow\lambda{\bf c}$, we assume ${\bf c}$ is a unit vector
without loss of generality. The residual is exactly minimized by performing
the Singular Value Decomposition (SVD) ${\bf G}={\bf U}{\bm{\Sigma}}{\bf
V}^{T}$, where ${\bm{\Sigma}}$ is a diagonal matrix and ${\bf U}$ and ${\bf
V}$ are rotation matrices. The diagonal elements $\sigma_{i}$ of
${\bm{\Sigma}}$ are the singular values of ${\bf G}$ and the columns of ${\bf
U}$ and ${\bf V}$ are the left and right singular vectors, respectively.
Efficient and stable algorithms for computing the SVD have been developed and
standardized [22, 23]. The exact minima of $r({\bf c})$ is given by the last
column of ${\bf V}$; it is the right singular vector associated with the
smallest singular value $\sigma_{\textrm{min}}$ since $r({\bf
c})=\sigma_{\textrm{min}}$. This vector will be dense in general, and we are
interested in sparse approximate minima. Rather than modify the residual to be
sparsity promoting, we manually remove terms one at a time such that $r({\bf
c})$ stays as low as possible at each iteration. Let ${\bf c}_{(k)}$ denote
the coefficient vector with $k$ nonzero elements. These form approximate
Pareto-optimal sets of coefficients.

This sequence of vectors are in nested sparse subspaces of the full library.
Once a coefficient is set to zero, it is never allowed to become nonzero. The
residuals $r_{k}\equiv r({\bf c}_{(k)})$ monotonically increase as $k$
decreases. One decides which term to eliminate from ${\bf c}_{(k)}$ by
recomputing the SVD with each column removed ${\bf G}$. Removing the $n$th
column of ${\bf G}$ before the SVD computation is equivalent to setting
$c_{n}=0.$ The greedy modification is selected such that $r_{k}$ increases the
least. This means that to make a sparsification step $c_{(k)}\rightarrow
c_{(k-1)}$ an SVD must be computed $k$ times. This elimination process is
repeated until a single term remains. The resulting optimization curve $r_{k}$
shows the greedy trade-off between sparsity and quantitative accuracy in the
coefficient vector ${\bf c}$.

A selection rule for the final sparse relation is ${\bf c}_{k}$, where
$r_{k}>\gamma r_{k+1}$ for some $\gamma>1$. The choice $\gamma=1.25$ was
empirically determined [7]. This is the only hyperparameter of the method, but
it does not change the underlying optimization curve $r_{k}$ or the
coefficients ${\bf c}_{(k)}$. Adjusting the hyperparameter $\gamma$ has
effectively no computational cost. Once a sparse relation is identified, the
term with the largest contribution $\|c_{n}{\bf g}_{n}\|_{2}$ can be
permanently removed, where ${\bf g}_{n}$ is the $n$th column of ${\bf G}$.
This reduced library library can be searched for further relations. Searching
is halted when $\sigma_{min}$ of the reduced library suggests further
relations do not exist.

For large libraries, recomputing the SVD for every possible rank-one
modification is expensive. Most computational time is spent at large $k$ when
${\bf c}_{k}$ is dense. Most library terms are not needed in this stage and it
seems wasteful to determine the optimal way to discard thousands of unhelpful
terms. There are two major considerations that speed up this method.

  1. 1.

All one cares about is the smallest singular value after a rank-1 modification
${\bf G}$. Calculation of the full SVD can be avoided.

  2. 2.

One does not need to start from scratch in the SVD evaluation. Helpful
theoretical results for updating the SVD after a rank-one modification are
presented in Ref [21]. This update does require full SVD information but can
be used to compute individual modified singular values quickly, which I
describe below.

Suppose that the SVD of ${\bf G}$ has already been computed ${\bf G}={\bf
U}{\bm{\Sigma}}{\bf V}^{T}$. Consider a rank-one modification ${\bf
G}^{\prime}={\bf G}\pm{\bf g}_{k}\otimes{\bf e}_{k}$ such that the column
${\bf g}$ is removed or added to the matrix. The matrix determinant lemma
implies that the new singular values of ${\bf G}^{\prime}$ are the roots of
the simple secular function $f_{\pm}(\sigma)$ [21] with

$\displaystyle
f_{-}(\sigma)=1-\frac{1}{\alpha^{2}}\sum_{j=1}^{n}\frac{w_{j}^{2}}{\sigma_{j}^{2}-\sigma^{2}},$
(6) $\displaystyle
f_{+}(\sigma)=1+\frac{1}{\alpha^{2}}\sum_{j=1}^{n}\frac{w_{j}}{\sigma_{j}^{2}-\sigma^{2}}-\frac{1-\|{\bf
w}\|^{2}_{2}}{\alpha^{2}\sigma^{2}},$ (7)

where $\alpha\equiv 1/\|{\bf g}\|_{2}$, ${\bf w}=\alpha{\bf U}^{T}{\bf g}$,
$n$ is the rank of ${\bf G}$ and $\sigma_{j}$ are the singular values of ${\bf
G}.$ $f_{+}(\sigma)$ and $f_{-}(\sigma)$ will be used for removing and adding
a column, respectively. An example of $f_{-}(\sigma)$ for a low-dimensional
matrix is shown in Figure 2. Note the poles at the previous singular values
are very useful: they imply the new singular values are uniquely placed
between the old ones. The new smallest singular value
$\sigma^{\prime}_{\textrm{min}}\in(\sigma_{n},\sigma_{n-1})$, when removing a
column, and $\sigma^{\prime}_{\textrm{min}}\in(0,\sigma_{n})$, when adding a
column.

Figure 2:  An evaluation of the secular function $f_{-}(\sigma)$ for an
example matrix ${\bf G}$ in $\mathbb{R}^{3\times 3}$. The last column of the
matrix is set to be removed. (a) $f_{-}(\sigma)$ has poles at the singular
values with definitive signs: the new singular values $\sigma^{\prime}_{i}$
will be between the previous singular values
$\sigma_{i}>\sigma^{\prime}_{i}>\sigma_{i+1}$ when there is no degeneracy in
${\bf G}$. The two new singular values are the roots of $f(\sigma)$
represented by red points. (b) The regularized function
$\tilde{f}_{-}(\sigma)$, which is now finite at the smallest two singular
values. Note that $\tilde{f}_{-}(\sigma)$ is not monotonic. The blue points
represent the bounds on the singular value, and the values of $\tilde{f}_{-}$
at these points will always be of opposite sign. (c) The numerical convergence
of the bisection. $\sigma_{(k)}$ is the $k$th iteration of bisection.

The functions $f_{\pm}(\sigma)$ have poles at these respective bounds but they
are trivially removable. Indeed, regularized secular functions
$\tilde{f}_{\pm}$ are better suited for numerical root-finding, i.e.,

$\displaystyle\tilde{f}_{-}(\sigma)=\alpha^{2}\frac{(\sigma^{2}-\sigma_{n}^{2})(\sigma^{2}-\sigma_{n-1}^{2})}{\sigma_{n}^{2}-\sigma_{n-1}^{2}}f_{-}(\sigma),$
(8)
$\displaystyle\tilde{f}_{+}(\sigma)=\alpha^{2}\frac{(\sigma^{2}-\sigma_{n}^{2})\sigma^{2}}{\sigma_{n}^{2}}f_{+}(\sigma).$
(9)

Let $l$ and $u$ be the lower and upper bounds of
$\sigma^{\prime}_{\textrm{min}}$. One can show that $f_{\pm}(l)>0$ and
$f_{\pm}(u)<0$. By the intermediate value theorem, a root $l<\sigma<u$ exists.
Bisection can be used to converge a root of these functions. One constructs a
guess $\sigma_{\textrm{guess}}=(u+l)/2$ and evaluates $f_{\pm}$ at this guess.
If the sign of the function is positive (negative), then replace the lower
(upper) bound. Since bisection halves the candidate interval at each
iteration, one can expect the error in the root at iteration $k$ to be
$\delta\equiv|\sigma^{\prime}_{\textrm{min}}-\sigma_{\textrm{guess}}|\propto
2^{-k}$. The value of the function will exhibit this same scaling:
$f(\sigma_{\textrm{guess}})\approx\delta
f^{\prime\pm}(\sigma_{\textrm{min}})\propto 2^{-k}$ if the derivative is
nonvanishing at the root. This predicted scaling can be seen in Figure 2(c).
In this example, $f_{\pm}$ will be zero to machine precision after $O(50)$
function evaluations.

SPRINT$\pm$ uses the bisection method to determine the best column
modification to a matrix. After the best modification is found, the economy
SVD of the modified matrix is recalculated. In SPRINT$-$, one begins with a
full feature matrix and removes columns one at a time such that
$\sigma_{\textrm{min}}$ increases the least at each stage. In SPRINT$+$, an
initially small feature matrix has columns added such that the
$\sigma_{\textrm{min}}$ decreases the most at each stage. SPRINT$+$ is not as
model agnostic as SPRINT$-$ since the resulting optimization curve is heavily
determined by the initial guess of sparse ${\bf c}$. One should use a dominant
balance or physical domain knowledge to guess an initial relation that is
refined by adding terms sequentially. While constructing this initial guess
requires more thought, one does not need to add terms to the relation until
the full feature matrix is recovered. One halts SPRINT$+$ once a satisfactory
number of terms have been added to the relation. Indeed, most computation time
in the model agnostic SPRINT$-$ is spent when the feature matrix has most of
its columns, which is undesirable. SPRINT$+$ avoids this regime altogether by
halting at a maximal model size.

Lastly, neither the residual $r_{k}$ nor the coefficient vectors ${\bf c}_{k}$
are affected by applying a rotation ${\bf Q}$ to the left of ${\bf
G}\rightarrow{\bf Q}{\bf G}$. This rotational symmetry can be used to
drastically reduce the size of ${\bf G}$ for regression. If ${\bf G}$ is a
large over-determined matrix, one can compute the reduced QR factorization
such that ${\bf G}={\bf Q}{\bf R}$ [24]. ${\bf R}$ is an upper-triangular
square matrix and, using it in place of ${\bf G}$, eliminates any strong
scaling with the number of observations in the feature matrix. This can be
treated as a preprocessing step.

##  IV Quantifying Algorithm Efficiency

To quantify the efficientcy of SPRINT$\pm$, I now explore a particular example
of model discovery using the Kuramoto–Sivashinsky (KS) equation, which is a
classical example of a chaotic partial differential equation. The KS equation
has been a canonical testing grounds for model discovery studies [14, 6, 25].
The chaotic nature of the KS equation is appealing because it provides rich
timeseries data. The purpose of SPRINT is to scale to large library sizes and
to learn small but important coefficients.

In principle, such small modifications are inevitable when solving a PDE on a
numerical mesh. Grid effects manifest as modifications to the governing PDE
with coefficients proportional to grid spacing. To demonstrate that sparse
regression can identify real small terms, I instead use a high accuracy
integrator and modify the KS equation to contain several higher order
nonlinearities with very small coefficients

$\displaystyle\partial_{t}u+u\partial_{x}u+\partial_{x}^{2}u+\partial_{x}^{4}u-\varepsilon\sum_{k=3}^{6}\partial_{x}u^{k}=0.$
(10)

The small coefficient is taken to be $\varepsilon=10^{-6}$ so that the
solution closely resembles the unaltered KS trajectory. This choice is made so
that the target for sparse regression is unambiguous yet it retains the flavor
of identifying grid effects.

I solve equation (10) on a spatial domain of length $L=22$ with periodic
boundary conditions for a temporal length of $T=400$. The geometry of the KS
state space has been well-studied at this length and the Lyupanov exponent is
$\sim 1/23$ [26, 27]. I use 128 gridpoints in space and 5120 in time. I
approximate spatial derivatives psuedospectrally with Fast Fourier Transforms.
In order to ensure that the temporal discretization error is minimal, I have
used a sixth order Gauss-Legendre Runge-Kutta scheme [28]. This scheme is
fully implicit, so the Runge-Kutta constraints are determined by iterating a
Newton-Krylov scheme until their $\ell_{2}$ error is below $10^{-9}$ [29].
Achieving this accuracy requires O(10) Newton iterations per timestep. The
initial condition is $u_{0}(x)=\cos(3\tilde{x})-\sin(\tilde{x})/2$ where
$\tilde{x}$ is the nondimensionalized position $\tilde{x}\in[0,2\pi]$. I chose
this initial condition to not have translational symmetry, although it does
possess shift-reflection symmetry $u_{0}(\tilde{x}+\pi)=-u_{0}(\tilde{x})$.
The nonlinear modifications of equation (10) do not respect this symmetry, so
it should be violated slowly in time. The resulting dynamic field $u(x,t)$ can
be seen in Figure 3.

Figure 3:  A numerical solution $u(x,t)$ to equation (10). The domain size
$L=22$ is chosen to make the evolution chaotic. The red rectangles show the
sizes and random distribution of weak-form subdomains used for the weak
evaluation of library terms.

I now use sparse regression to recover the modified KS equation. However, a
completely general library will not be helpful for this demonstration. This is
because, if equation (10) is true, so are derivatives of this equation. This
means sparse regression should be able to recover linearly independent
relations for $\partial_{t}^{2}u,$ $\partial_{t}^{3}u$, $\cdots$, or some
linear combination of these. Further, simple multiplication of equation (10)
by an arbitrary term $\phi$ produces true equations. To fix the library and
only produce a single dynamical relation, I only include the time derivative
once as $\partial_{t}u.$ I then allow for the rest of the library to contain
words made from the symbolic alphabet $\mathcal{A}=\\{u,\partial_{x}\\}$ up to
length 10.

$\mathcal{L}=\\{\partial_{t}u\\}\cup\textrm{words}_{10}(\mathcal{A})$ (11)

The resulting library contains $139$ terms. Homogeneous regression on this
library will either produce a dynamic closure or a spatial constraint
capturing the geometry of the attractor: that is either a closure
$\partial_{t}u=f(u,\partial_{x}^{n}u)$ or a spatial PDE
$f(u,\partial_{x}^{n}u)=0$ that describes the observed data despite its
dynamic nature. The numerical solution $u(x,t)$ is perturbed by additive
uncorrelated Gaussian noise $\varepsilon N(\mu,\sigma)$ with amplitude
$\varepsilon=10^{-7}$, mean $\mu=0$, and a standard deviation of
$\sigma=0.218$.

The observations of library terms are taken to be standard weak-formulation
integrals over rectangular spacetime subdomains $\Omega_{m}$ [12, 14, 7, 10].
Each element of ${\bf G}$ corresponds to an appropriately scaled integral over
the data with a weight function $\phi$ chosen so that integration by parts is
convenient. Any smooth weight function can be used, including data-driven
weights [10], but I have chosen envelope polynomials

$\displaystyle
G_{mn}\equiv\frac{1}{V_{m}S_{n}}\int_{\Omega_{m}}d\bar{x}d\bar{t}\,\,\phi
f_{n},$ (12) $\displaystyle\phi(x,t)=(1-\bar{x}^{2})^{8}(1-\bar{t}^{2})^{8}.$
(13)

Here, $\bar{x}$ and $\bar{t}$ are rescaled coordinates that lie in the
canonical interval $[-1,1]$. The nondimensionalization factors $V_{m}$ and
$S_{n}$ ensure that the elements of ${\bf G}$ take on reasonable values for
all library terms [7]. $V_{m}$ is the integral of the weight function.

$V_{m}=\int_{\Omega_{m}}d\bar{x}d\bar{t}|\phi|.$ (14)

$S_{n}$ are combinations of characteristic scales and statistics of the data.
Specifically, $S_{n}$ is a product of scales defined such that the scale of a
undifferentiated field $u$ is its mean $\mu_{u}$ and the scale of a
differentiated field $u$ is proportional to its standard deviation
$\sigma_{u}$, i.e.,

$\displaystyle S[u]=\mu_{u},$ $\displaystyle
S[\partial_{x}^{n}u]=\frac{\sigma_{u}}{L_{u}^{n}}.$ (15)

$L_{u}$ is a length scale associated with the variation of $u$. When a library
term is a product of such subterms, the overall scale is the product of
subscales. For reasonably smooth data, this rescaling should produce elements
$G_{mn}$ that are O(1). I chose the physical size of the subdomains
$\Omega_{m}$ to be the order of the length and time scales of characteristic
process in the KS dynamics: 1024 subdomains $\Omega_{m}$ chosen from a uniform
distribution (see Figure 3).

Once I compute the matrix ${\bf G}$, I apply both an exhaustive search and
SPRINT$\pm$ to produce optimization curves. The resulting curves are shown in
Figure 4 displaying the monotonic trade-off between quantitative accuracy and
model complexity.

Figure 4:  Optimization curve showing the greedy residual $r_{k}$ as a
function of nonzero terms in ${\bf c}$. The starting subspace for SPRINT+ was
the $k=1$ choice of $\partial_{x}^{2}u$ to coincide with the exhaustive
search. (a) Optimization curve for the library including the time derivative
$\partial_{t}u$. A dynamic closure is found as demonstrated by the elbow at
$k=8$. (b) The optimization curve when the time derivative is not included in
the library. Instead of a dynamic closure, an approximate spatial constraint
is found on the solution manifold at $k=13$.

Figure 5: (a) Empirical walltime scaling for various sparse regression
techniques. Power law exponents were empirically determined by performing a
least squares fit to $\log(t)=\alpha\log(n)+\beta$ for $\alpha$ and $\beta$.
Each point shows the walltime of an algorithm applied to a real square matrix
with elements drawn from a uniform distribution over $[-1,1]$. The reported
walltime is the mean of eight independent trials. Four algorithms are
considered: the exhaustive search, SPRINT–, SPRINT+, and an implementation of
implicit-SINDy available at https://github.com/dynamicslab/SINDy-PI [4]. The
exhaustive search algorithm displays the $n^{\sim 4}$ scaling, while the
SPRINT– algorithm has approximately $n^{\sim 3}$ scaling. These algorithms
start with similar cost for small libraries, but quickly depart as the library
size $n$ grows beyond O(100). (b) The extrapolated walltimes of sparse
regression for the library sizes of Figure 1. Calculations that would take the
exhaustive search the order of a Hubble time can be accomplished in one day
with SPRINT+.

Figure 4 shows that the exhaustive search and SPRINT+ produce the same
optimization curves, while SPRINT– discards useful information too early. All
three methods find an elbow in the residual at $k\approx 8$. SPRINT+ and the
exhaustive search exactly identify the full modified equation (10) at $k=8$ as
the primary elbow. A numerical instability in SPRINT– produces an artificial
elbow at k=10. All three methods recover the secondary elbow $k=4$
corresponding to the unmodified KS equation.

The SPRINT+ method is demonstrably capable of capturing terms with small
coefficients and doing so accurately. The identified coefficients of equation
(10) are $c_{3}=0.999\times 10^{-6}$, $c_{4}=0.993\times 10^{-6}$,
$c_{5}=1.001\times 10^{-6}$, and $c_{6}=1.001\times 10^{-6}$. The ratio of
residuals after including the four modifications is $r_{4}/r_{8}=264$.
Capturing these produces a substantial change in the accuracy of the relation.

An interesting situation arises if the time derivative is discarded from
$\mathcal{L}$ so that only terms algebraic in $u$ and its spatial derivative
remain. The library can then contain effective spatial descriptions of the
state-space of this modified KS system. An elbow for this library can be seen
in the right panel of Figure 4 at $k=13$. The resulting constraint is

$\displaystyle
C(u,\partial_{x}^{n}u)=c_{1}u^{2}\partial_{x}^{3}u+c_{2}u\partial_{x}^{1}u\partial_{x}^{2}u+c_{3}u\partial_{x}^{4}u$
$\displaystyle+c_{4}(\partial_{x}^{1}u)^{3}+c_{5}\partial_{x}^{1}u\partial_{x}^{3}u+c_{6}(\partial_{x}^{2}u)^{2}+c_{7}\partial_{x}^{5}u$
$\displaystyle+c_{8}u\partial_{x}^{6}u+c_{9}\partial_{x}^{1}u\partial_{x}^{5}u+c_{10}\partial_{x}^{2}u\partial_{x}^{4}u+c_{11}(\partial_{x}^{3}u)^{2}$
$\displaystyle+c_{12}\partial_{x}^{7}u+c_{13}\partial_{x}^{9}u=0.$ (16)

This learned constraint respects the reflection symmetry
$u(x,t)\rightarrow-u(-x,t)$ of the KS equation and the initial data despite
this symmetry not being enforced in the library. The modified equation (10)
does not possess this reflection symmetry exactly, but the effect of symmetry
breaking is well below the accuracy of relation (16). This highlights the
versatility of homogeneous regression.

##  V Conclusion

Bisection greatly accelerates the calculations of optimal rank-one
modifications of systems of equations. These rank-one modifications include
the addition and removal of coefficients in sparse regression. The two
algorithms SPRINT– and SPRINT+ presented here utilize bisection to find sparse
approximate null vectors. For well conditioned matrices, SPRINT– should
exactly reproduce the exhaustive search [7], although a numerical instability
can lead to early removal of useful coefficients as seen in Figure 4. There
can be numerical difficulties in accurately evaluating the secular equation
(9) when $\sigma_{\textrm{min}}$ becomes exceptionally close to zero. If
enough noise is added to the numerical solution $u$, then SPRINT– is capable
of reproducing the same optimization curve as the exhaustive search. This
instability is not observed in the additive variant SPRINT+.

The only complication is that SPRINT+ needs a sparse initial guess to add
terms to. This is not a problem for debugging if the suspected governing
equations are known a priori. For data with unknown governing equations, few-
term dominant balances can be cheaply computed with a combinatoric search.
SPRINT+ will generally not reproduce the same curve as the exhaustive search
or SPRINT– because the optimization curve depends on the initial sparse
coefficient vector. However, we see excellent performance of SPRINT+ in
reconstructing the desired nonlinear dynamics with small coefficients and
without enforcing the inclusion of the time derivative.

SPRINT contains no hyper-parameters other than some model selection criterion
such as $r_{k-1}>\gamma r_{k}$ as was done in [7]. A different choice of such
hyper-parameters does not require the costly aspects of the algorithm to be
repeated. This is because the primary output is the entire optimization curve
(and associated coefficients) rather than a single sparse model as in
sequential thresholding or LASSO. Th ultimate choice of sparse model is up to
the modeler, which lends this method to being highly interactive. There can be
multiple reasonable cutoffs depending on acceptable accuracy of a relation. In
Figure 4 the four term KS equation provides a reasonable description of the
data with a residual of $r<10^{-4}.$

The benefit of searching for a null vector ${\bf G}{\bf c}=0$ rather than
fitting a right hand side ${\bf G}{\bf c}={\bf b}$ is the model agnosticism.
This modification makes model discovery general purpose while not sacrificing
the ability to find dynamical equations. Furthermore, the null vector approach
can find single term equations like $\nabla\cdot{\bf B}=0$ [12], which the
inhomogeneous method is incapable of doing. Note that implicit-SINDy [20] was
used as a benchmark rather than SINDy-PI since SPRINT$\pm$ is investigated
here as a serial algorithm. SPRINT$\pm$ can be further accelerated by
parallelizing the SVD [30, 31] and the embarrassingly parallel bisection
calculations.

The additive variant SPRINT+ offers a staggering computational speedup
compared to SPRINT– and the exhaustive search. SPRINT+ is a promising
candidate for large-library sparse regression. It displays sensitivity to
small coefficients, efficient scaling with library size, and weak model
agnosticism. The resulting optimization curve provides more insight than a
single sparse solution.

Since the symbolic form of modifications are unknown and potentially symmetry-
breaking, it is important for sparse regression to scale well with library
size. The accelerated method SPRINT presented here is an order of magnitude
faster than the exhaustive search at the O(100) library. Figure 5(a)
investigates the empirical scaling of several sparse regression methods on
uniform random matrices with elements in $[-1,1]$. These walltimes were
measured with MATLAB’s benchmarking tools. All methods display power law
scaling with the library size with distinct exponents. Figure 5 displays the
extropolated walltimes of various methods when considering the motivating MHD
library. For small word sizes, the exhaustive search, SPRINT–, and an
implementation of implicit-SINDy become prohibatively expensive. The
exhaustive search in particular requires more than the Hubble time for a word
size of eight, while SPRINT+ reproduces the same effective result in a day.

All code for reproducing these results is available at
https:/github.com/mgolden30/SPIDER. A minimal implementation of SPRINT$\pm$ in
both MATLAB and Python is available at
https:/github.com/mgolden30/FastSparseRegression.

## Acknowledgments

I am thankful to Dimitrios Psaltis and Daniel Gurevich for helpful
conversations.

## References

  * Bongard and Lipson [2007] J. Bongard and H. Lipson, Automated reverse engineering of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 104, 9943 (2007).
  * Schmidt and Lipson [2009] M. Schmidt and H. Lipson, Distilling free-form natural laws from experimental data, science 324, 81 (2009).
  * Brunton  _et al._ [2016] S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the national academy of sciences 113, 3932 (2016).
  * Kaheman  _et al._ [2020] K. Kaheman, J. N. Kutz, and S. L. Brunton, Sindy-pi: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics, Proceedings of the Royal Society A 476, 20200279 (2020).
  * Schaeffer and McCalla [2017] H. Schaeffer and S. G. McCalla, Sparse model selection via integral terms, Physical Review E 96, 023302 (2017).
  * Rudy  _et al._ [2017] S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Data-driven discovery of partial differential equations, Science advances 3, e1602614 (2017).
  * Gurevich  _et al._ [2021] D. R. Gurevich, M. R. Golden, P. A. Reinbold, and R. O. Grigoriev, Learning fluid physics from highly turbulent data using sparse physics-informed discovery of empirical relations (spider), arXiv preprint arXiv:2105.00048 (2021).
  * Alves and Fiuza [2022] E. P. Alves and F. Fiuza, Data-driven discovery of reduced plasma physics models from fully kinetic simulations, Physical Review Research 4, 033192 (2022).
  * Joshi  _et al._ [2022] C. Joshi, S. Ray, L. M. Lemma, M. Varghese, G. Sharp, Z. Dogic, A. Baskaran, and M. F. Hagan, Data-driven discovery of active nematic hydrodynamics, Physical review letters 129, 258001 (2022).
  * Golden  _et al._ [2023] M. Golden, R. O. Grigoriev, J. Nambisan, and A. Fernandez-Nieves, Physically informed data-driven modeling of active nematics, Science Advances 9, eabq6120 (2023).
  * Stone [1948] M. H. Stone, The generalized weierstrass approximation theorem, Mathematics Magazine 21, 237 (1948).
  * Gurevich  _et al._ [2019] D. R. Gurevich, P. A. Reinbold, and R. O. Grigoriev, Robust and optimal sparse regression for nonlinear pde models, Chaos: An Interdisciplinary Journal of Nonlinear Science 29 (2019).
  * Reinbold  _et al._ [2021] P. A. Reinbold, L. M. Kageorge, M. F. Schatz, and R. O. Grigoriev, Robust learning from noisy, incomplete, high-dimensional experimental data via physically constrained symbolic regression, Nature communications 12, 3219 (2021).
  * Messenger and Bortz [2021] D. A. Messenger and D. M. Bortz, Weak sindy for partial differential equations, Journal of Computational Physics 443, 110525 (2021).
  * Zhang and Schaeffer [2019] L. Zhang and H. Schaeffer, On the convergence of the sindy algorithm, Multiscale Modeling & Simulation 17, 948 (2019).
  * Santosa and Symes [1986] F. Santosa and W. W. Symes, Linear inversion of band-limited reflection seismograms, SIAM journal on scientific and statistical computing 7, 1307 (1986).
  * Zou and Hastie [2005] H. Zou and T. Hastie, Regularization and variable selection via the elastic net, Journal of the Royal Statistical Society Series B: Statistical Methodology 67, 301 (2005).
  * Huang  _et al._ [2018] J. Huang, Y. Jiao, Y. Liu, and X. Lu, A constructive approach to $l\\_0$ penalized regression, Journal of Machine Learning Research 19, 1 (2018).
  * Natarajan [1995] B. K. Natarajan, Sparse approximate solutions to linear systems, SIAM journal on computing 24, 227 (1995).
  * Mangan  _et al._ [2016] N. M. Mangan, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Inferring biological networks by sparse identification of nonlinear dynamics, IEEE Transactions on Molecular, Biological and Multi-Scale Communications 2, 52 (2016).
  * Bunch and Nielsen [1978] J. R. Bunch and C. P. Nielsen, Updating the singular value decomposition, Numerische Mathematik 31, 111 (1978).
  * Stewart [1993] G. W. Stewart, On the early history of the singular value decomposition, SIAM review 35, 551 (1993).
  * Anderson  _et al._ [1999] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen,  _LAPACK Users’ Guide_ , Vol. 9 (SIAM, 1999).
  * Trefethen and Bau III [1997] L. N. Trefethen and D. Bau III,  _Numerical Linear Algebra_ , Vol. 50 (SIAM, 1997).
  * Köster  _et al._ [2023] F. Köster, D. Patel, A. Wikner, L. Jaurigue, and K. Lüdge, Data-informed reservoir computing for efficient time-series prediction, Chaos: An Interdisciplinary Journal of Nonlinear Science 33 (2023).
  * Cvitanović  _et al._ [2010] P. Cvitanović, R. L. Davidchack, and E. Siminos, On the state space geometry of the kuramoto–sivashinsky flow in a periodic domain, SIAM Journal on Applied Dynamical Systems 9, 1 (2010).
  * Edson  _et al._ [2019] R. A. Edson, J. E. Bunder, T. W. Mattner, and A. J. Roberts, Lyapunov exponents of the kuramoto–sivashinsky pde, The ANZIAM Journal 61, 270 (2019).
  * Butcher [1964] J. C. Butcher, Implicit runge-kutta processes, Mathematics of computation 18, 50 (1964).
  * Zhigunov and Grigoriev [2023] D. Zhigunov and R. O. Grigoriev, Exact coherent structures in fully developed two-dimensional turbulence, Journal of Fluid Mechanics 970, A18 (2023).
  * Andez  _et al._ [2008] V. H. Andez, J. E. R. AN, and A. É. T. AS, A robust and efficient parallel svd solver based on restarted lanczos bidiagonalization, Electronic Transactions on Numerical Analysis 31, 68 (2008).
  * Jessup and Sorensen [1994] E. R. Jessup and D. C. Sorensen, A parallel algorithm for computing the singular value decomposition of a matrix, Siam Journal on Matrix Analysis and Applications 15, 530 (1994).