to classical results from independent component analysis (ICA) (Hyvärinen et
al., 2004). Specifically, a linear SEM can be represented by,

$\bm{X}=\mathbf{B}\bm{X}+\bm{\varepsilon},$ (2.4)

where $\bm{X}$ is the vector of system variables, $\bm{\varepsilon}$ the noise
vector and $\mathbf{B}$ is the matrix of coefficients. Thus solving Equation
2.4 for $\bm{X}$ we obtain,

$\bm{X}=(\mathbf{I}-\mathbf{B})^{-1}\bm{\varepsilon}.$

The above equation is the an ICA problem, and its theory states that when
$\bm{\varepsilon}$ are non-Gaussian noises variables, the mixing matrix
$\mathbf{A}=(\mathbf{I}-\mathbf{B})^{-1}$ is identifiable in the large sample
limit (Hyvärinen et al., 2004). LiNGAM-ICA works by first exploiting the ICA
results to obtain the mixing matrix $\mathbf{A}$, and secondly by permuting
and normalising $\mathbf{A}$ to obtain the appropriate matrix $\mathbf{B}$ and
the corresponding causal DAG. The LiNGAM can also be solved directly without
using ICA (Shimizu et al., 2011). The direct-LiNGAM iteratively selects the
variable which is the most independent from the residual (in a linear
regression onto the remaining variables) and replaces the original data with
the residuals matrix. In the end, the procedure gives a causal order; and
ultimately, simple linear regressions (e.g. estimated with least squares) are
used to obtain the final estimation of the lower triangular $\mathbf{B}$
matrix.

VAR-LiNGAM (Hyvärinen et al., 2008; Zhang, 2008a; Hyvärinen et al., 2010) is
the extension of the classical LiNGAM model to time series data. It considers
a linear structural vector auto-regressive (VAR) model, with possibly acyclic
instantaneous relationships, and assumes non-Gaussian disturbances. In
particular, the process $\bm{X}_{t}$ is assumed to evolve following

$\bm{X}_{t}=\sum_{\tau=0}^{k}\mathbf{B}_{\tau}\bm{X}_{t-\tau}+\bm{\varepsilon}_{t}.$

Where $\mathbf{B}_{0}$ is the matrix of instantaneous relationships and its
sparsity pattern corresponds to a DAG while $\mathbf{B}_{\tau}$ for $\tau>0$
are the matrices of lagged relationships; moreover the non-Gaussian noise
vector $\bm{\varepsilon}_{t}$ has independent components and is assumed
independent over time. Hyvärinen et al. (2010) propose two methods to estimate
the coefficients of the model above: (1) a two-stage method which combines
least-square estimation of the autoregressive model and classical LiNGAM
estimation; (2) a method based on multichannel blind deconvolution. In Figure
2.2, we illustrate that if the assumptions are satisfied, VAR-LiNGAM is able
to discover and orient all causal links.

##### Structural equation models for time series, TiMiNo.

Peters et al. (2013) extended the classical structural equation modelling
framework to the time series setting by introducing time series models with
independent noise (TiMINo). In particular, a multi-variate time series
$X_{t}=(X_{t}^{i})_{i\in V}$ satisfies a TiMINo if there exists a $p>0$ and,
for every $i\in V$, there are subsets $\mathbf{PA}_{0}^{i}\subseteq
X^{V\setminus\\{i\\}}$ and $\mathbf{PA}_{k}^{i}\subseteq X^{V}$ for
$k=1,\ldots,p$ such that,

$X_{t}^{i}=f_{i}\left((\mathbf{PA}_{p}^{i})_{t-p},\ldots,(\mathbf{PA}_{1}^{i})_{t-1},(\mathbf{PA}_{0}^{i})_{t},N_{t}^{i}\right),$
(2.5)

where $N_{t}^{i}$ are assumed to be jointly independent over $i$ and time and
for each $i$, $N^{i}_{t}$ are identically distributed in time. Thus, the
assumed model is extremely general, instantaneous relationships are allowed
and the noise contribution is not restricted by allowing arbitrarily mixing
through the $f_{i}$ functions. Nevertheless, to prove the identifiability of
the full-time graph and the causal summary time graph, additional assumptions
have to be made. Specifically, in Peters et al. (2013) it is assumed that
either (i) Equations 2.5 follows an identifiable functional model class (e.g.
nonlinear functions with additive Gaussian noise or linear functions with
additive non-Gaussian noise (Peters et al., 2011)) or (ii) each function
$f_{i}$ exhibits a time structure, that is the union of the causal parents of
$X^{i}_{t}$ contains at least one $X_{t-k}^{i}$ and moreover the joint
distribution is faithful with respect to the full-time graph with the summary
time graph being acyclic. Under either assumption (i) or (ii), the complete
causal graph is proved to be identifiable in the large sample limit.
Practically, to estimate a TiMINo, methods for causal discovery for additive
noise models with i.i.d. data (Mooij et al., 2009) can be adapted to the time
series setting. As in the direct-LiNGAM, it iteratively selects the causal
order between the variables by fitting regression models and evaluating the
independence between the variables and the residuals. To fit the regression
models for $f_{i}$, various methods can be used, such as vector autoregressive
models (linear), generalised additive models or Gaussian processes. Moreover,
to test independence from the residuals, the HSIC (Gretton et al., 2007) can
be applied to all possible shifted time series up to the maximum lag order.

##### Invariant Causal Prediction.

Invariant causal prediction deals with the setting of independent and
identically distributed samples of the random vectors
$X=(X_{1},X_{2},\ldots,X_{p})^{\top}\in\mathbb{R}^{p}$,
$E=(E_{1},E_{2},\ldots,E_{q})^{\top}\in\mathbb{R}^{q}$, and $Y\in\mathbb{R}$.
For a variable or vector of interest $Y$, $E$ is a set of environment
variables that may be causes of $X$ but are not direct causes or effects
(direct or indirect) of $Y$. ICP assumes that $Y$ is generated _causally_ from
a subset $S^{*}\subseteq\\{1,\ldots,p\\}$ of the $p$ variables considered, so
that there is causal sufficiency, and $Y$ is generated from a Structural
Causal Model obeying:

$\displaystyle Y=g(X_{S^{*}},\epsilon),\phantom{bl}\epsilon\sim
F,\phantom{bl}\epsilon\rotatebox[origin={c}]{90.0}{$\models$}X_{S^{*}},$ (2.6)

where $g$ and $F$ are arbitrary functions and distributions, respectively. A
set $S\subseteq\\{1,\ldots,p\\}$ is a generic subset of the full set of
candidate causes. We refer to $S$ and $X_{S}$ interchangeably for brevity. The
task of ICP is to infer the set $S^{*}$ of direct causes of the variable of
interest $Y$.

The Invariance Causal Prediction (ICP) framework (Peters et al., 2016) is
based on the observation that $Y$ is independent of $E$ given $X_{S^{*}}$,
denoted $Y\rotatebox[origin={c}]{90.0}{$\models$}E|X_{S^{*}}$. Assuming we
have a set of candidate causes $X$ that includes $S^{*}$, the causal subset,
and an environment variable $E$ that we know does not directly cause $Y$ or is
an effect of $Y$. We can search for $S^{*}$ by applying a conditional
independence test $Y\rotatebox[origin={c}]{90.0}{$\models$}E|X_{S}$ on $Y$,
$E$ and subsets $S\subseteq\\{1,\ldots,p\\}$ of $X$. ICP then selects as
causal variables the intersection of all those subsets $S$ where the
corresponding conditional independence test is not rejected:

$\displaystyle\hat{S}^{*}=\bigcap\limits_{S:p_{S}>\alpha}S.$ (2.7)

Here $p_{S}$ is the $p$-value associated with the conditional independence
test of $Y\rotatebox[origin={c}]{90.0}{$\models$}E|X_{S}$, with the null
hypothesis corresponding to conditional independence. We do not reject
conditional independence at significance level $\alpha$ if $p_{S}>\alpha$.

One way to interpret the problem setting is that causal associations are more
robust than other associations. So ICP finds the causes of a variable of
interest by investigating which associations are invariant across
environments. Another interpretation is that the environment variables define
data generated under different interventions to the system (SCM).

In physical sciences, regional and temporal variables are good candidates
since these often describe changes in environments that alter the conditions
under which physical processes occur.

Peters et al. (2016) introduce an algorithm to implement ICP that assumes
linear relationships between causes $X$ and effects $Y$ and a categorical,
univariate variable $E$. In (Heinze-Deml et al., 2018a), more general
algorithms are presented that allow for nonlinear relationships between cause
and effects and for continuous and multivariate environment variables $E$.
Pfister et al. (2019) provide an ICP variant for time series data. The
proposed time-series ICP relies on the causal invariance assumption across
time points, thus removing the requirement of environment knowledge. This
setting is also partially robust to hidden confounders, similar to the
original ICP (Peters et al., 2016) framework, and in general, the ICP is
expected to be conservative with respect to violations of its assumptions
(Pfister et al., 2019).

##### Causal frameworks for continuous-time systems.

As we already reviewed, discrete-time causal systems fit directly as an
extension of i.i.d. and DAG framework. When considering discrete-time systems,
we can express the value of a variable at a time $t$ as a function of other
variables (and itself) observed at past instants, thus the complete causal
graph can be seen as a DAG extended (possibly infinitely) in time. Classical
methods for causal discovery in the i.i.d. setting can then be adapted to
discrete-time systems quite straightforwardly.

Conversely, considering continuous-time systems, raise the issue that a time-
extended DAG is not feasible, since the included variables would be
uncountable (Hansen and Sokol, 2014). Nevertheless, modelling continuous
dynamical system helps in dealing with non-uniform sampling and extrapolating
among different sampling frequencies, two major drawbacks of causal discovery
methods for discrete-time systems. The following are the major causal
frameworks available for continuous-time systems:

  1. 1.

Causal interpretation of ODE and SDE. Various efforts have been made to
describe causal systems with ODEs and SDEs. First, causal discourses around
ODE were used to obtain different justifications for the cyclic SEM (Hyttinen
et al., 2012; Mooij et al., 2013a; Bongers et al., 2018; Peters et al., 2022).
Rubenstein et al. (2018) described Dynamical Structural Causal Models (DSCM)
as extensions of SEM where each equation or assignment is a relationship
between a set of causal parent trajectories and an effect trajectory. Under
some stability conditions, such DSCM can be obtained from ODE systems. SDEs
have also been studied from a causal perspective (Hansen and Sokol, 2014;
Peters et al., 2022). The advantage of SDEs in modelling physical systems is
that they allow incorporating an inherent source of stochasticity, a common
assumption in numerous real-world systems (Abbati et al., 2019). Graphical
parameterisations of SDE equilibrium distributions leading to models allowing
for cycles have been investigated and different structure learning algorithms
have been proposed (Varando and Hansen, 2020).

  2. 2.

Dynamic Causal Models. Dynamic Causal Models (DCM) (Friston et al., 2003;
Marreiros et al., 2010; Stephan et al., 2010; Abbati et al., 2019; Friston et
al., 2013), in short, is a Bayesian framework for fitting and comparing causal
models for coupled dynamical systems. DCM was introduced and applied mostly in
Neuroscience, and particularly in the problem of estimating the connectivity
between brain regions from neuroimaging data, as discussed in detail in
Section 4.2.2. Recently, and it has been even employed to model the COVID-19
pandemic dynamics (Friston et al., 2020, 2022).

  3. 3.

Local independence graphs. Local independence is a notion of conditional
independence for stochastic processes (both discrete and continuous time ones)
which can be (in)dependent on each other pasts. In detail, for real-valued
stochastic processes $X_{t}=(X^{1}_{t},\ldots,X^{p}_{t})$ and
$A,B,C\subseteq\\{1,\ldots,p)$, we say that $X^{B}$ is locally independent of
$X^{A}$ given $X^{C}$ at time $t$ if the past of $X^{C}$ until time $t$
provides the same information, to predict
$E[X_{t}^{\beta}|\mathcal{F}_{t}^{A\cup C}]$ 111where $\mathcal{F}_{t}^{A\cup
C}$ is a right-continuous and complete filtration which represents the history
of the processes $X^{A\cup C}$ see e.g. (Mogensen and Hansen, 2020) for a
detailed description, as the past of $X^{A\cup C}$ until time $t$, for each
$\beta\in B$.

Didelez (Didelez, 2008, 2006) studied graphical representations of local
independence with directed graphs together with $\delta$-separation and proved
the equivalence of the pairwise and global Markov properties for multivariate
counting processes. Directed graphs and $\delta$-separation has been then
extended to mixed graphs and $\mu$-separation (Mogensen and Hansen, 2020) to
model partially unobserved systems. A constrained-based algorithm has been
proposed (Mogensen et al., 2018), which is proven to be sound and complete
under faithfulness assumption. Local independence graphs can be applied to
multivariate processes which are solutions of SDEs (such as the multivariate
Ornstein-Uhlenbeck process) or event and counting processes such as Hawkes
processes (Mogensen, 2022).

###  2.2 Challenges

In this section, we discuss several challenges for causal discovery that are
frequent in real-world applications. Following Runge et al. (2019a), we
distinguish between challenges related to the data-generating process itself
(see Sec. 2.2.1), challenges associated with the available data (see Sec.
2.2.2), and challenges of statistical or computational nature (see Sec.
2.2.3). Users should carefully consider the challenges they face in their
application and choose a suitable causal discovery method. More generally, how
to reason and deal with typical challenges in causal inference is further
discussed in a time series context in Runge et al. (2023).

####  2.2.1 Process challenges

“Discovering causal relations from observational data is impossible without
assumptions about the mechanisms and faces important challenges related to
data and statistical characteristics. The field will need to incorporate
domain knowledge and post-selection inference.”

Non-linearities pose challenges for causal discovery in both independence-
based and score-based methods. Non-linear conditional independence tests, such
as the GPDC test (Rasmussen and Williams, 2006) and tests based on conditional
mutual information (Runge, 2018), are computationally more expensive and tend
to have lower statistical power than linear tests. Non-linear functional
relationships in score-based methods require more complex score functions,
which can decrease finite-sample performance. However, non-linear functional
relationships can enhance identifiability in some asymmetry-based causal
discovery methods (Peters et al., 2017a). Overall, nonlinearities increase
model complexity and require careful consideration when selecting appropriate
causal discovery methods.

Most time series causal discovery methods assume data generated from a
causally stationary process with a unique equilibrium distribution. However,
many real-world processes are non-stationary. Recently there have been works
on devising causal discovery methods that first detect the variables afflicted
with non-stationary driving mechanisms and subsequently infer the entire
causal graph, including possibly proxy variables corresponding to the driving
force of non-stationarity (Mooij et al., 2020; Huang et al., 2020).
Furthermore, techniques exist to detect regime or context changes in non-
stationary data (Huang et al., 2015; Saggioro et al., 2020; Huang et al.,
2020), but it remains an an active area of research. Domain experts may be
able to identify the source of non-stationarity in time series data and
preprocess the data to remove it.

Time series data is common in physical sciences, and it has a distinctive
feature of auto-correlation. Many causal discovery algorithms are not designed
for time series data and show decreased performance when applied without
modification (Runge et al., 2019b; Runge, 2020). However, some methods are
specifically designed for time series data, reducing the detrimental effect of
auto-correlation. See the PCMCI algorithm above for a discussion on the
challenges of auto-correlated data. In many domains, space adds to time as
well. In principle, one could feed different spatial locations of the same
variable as distinct variables into causal discovery methods. However, this
naive approach ignores spatial correlations and quickly results in a high-
dimensional problem. Another workaround, employed for example in (Tibau et
al., 2022; Diego Bueso, 2020), is to perform dimension-reduction as a
preprocessing step and then perform causal discovery on the dimensionally-
reduced space. The development of causal discovery inherently designed for
spatio-temporal data is an active area of research; for example, see
(Christiansen et al., 2022). Dealing with variables whose dynamics operate on
different time scales, such as e.g. fast atmospheric and slow oceanic
processes, pose important challenges. Granger causality in the frequency
domain, e.g. (Bressler and Seth, 2011; Chicharro, 2011; Faes et al., 2012),
and a combination of wavelet analysis with transfer entropy (Lungarella et
al., 2007) are examples of approaches to deal with the time scales of causal
influences.

Finally, it is worth mentioning that many causal discovery methods assume the
data-generating SCM to be acyclic (see axis VII in Sec. 2.1.2). However, in
real-world applications, one can often not exclude the existence of feedback
that acts on time scales below the measurement interval. Such feedbacks make
the time series graph cyclic. In the non-temporal setting, one might often not
be able to exclude the existence of causal cycles. As discussed and referenced
above, there are causal discovery methods that can handle cyclic causal
graphs.

####  2.2.2 Data challenges

Various data challenges arise when tackling the problem of causal discovery in
practice. This is mainly due to the discrepancy between the assumed hypothesis
needed from each method or framework and the real-world data. One of the most
common assumptions of most causal discovery methods is causal sufficiency,
which is the hypothesis that all relevant variables are observed. Unobserved
variables are especially problematic when they are possible confounders
between system variables since omitting confounders from the causal discovery
could lead to learning spurious or wrong relationships. There are some
available methods which do not assume causal sufficiency, such as LPCMCI
(Gerhardus and Runge, 2020), SVAR-FCI (Malinsky and Spirtes, 2018) and GPS
(Claassen and Bucur, 2022) (see Sec. 2.1.3).

Missing data and selection bias are other common issues in real-world
applications, and there have been some efforts in developing causal discovery
methods which are resilient to these challenges (Strobl et al., 2018; Gain and
Shpitser, 2018; Tu et al., 2019; Versteeg et al., 2022). Not always,
especially in the physical realm, the data follow predictable and well-behaved
distributions. One such example is zero-inflated data, which is common, for
instance, in gene expression data, where single-cell expressions lack
detectable values of transcripts that appear abundant on bulk (thousands of
cells) gene expression experiments. Recent advances have developed graphical
models and causal discovery methods in such scenarios (McDavid et al., 2019;
Yu et al., 2020).

####  2.2.3 Statistical and computational challenges

The high dimensionality of data in physical systems, such as spatiotemporal
data, and small sample sizes are central statistical challenges for causal
discovery. On the other hand, large sample sizes raise issues of unaffordable
computational time, which can scale up to cubically for kernel methods
typically used for independence testing (Schölkopf and Smola, 2008; Gretton et
al., 2007). High-dimensional data leads to large conditioning sets in
particular algorithms, effectively reducing the sample size available to test
the hypothesis.

As noted in Sec. 2.2.1, non-linearity is a common characteristic of processes
in the physical sciences. For the case of independence-based or hybrid casual
discovery techniques, this calls for devising non-parametric tests of
independence, for instance, tests based on measures of conditional mutual
information (Runge, 2018), or on Gaussian process regression or other kernel-
based measures on independence (Gretton et al., 2007), or using quantile
regression (Petersen and Hansen, 2021) and copula-based methods (Bouezmarni et
al., 2012) (also applied to Granger causality), etc. The no-free-lunch theorem
of Shah and Peters (2020) states that no single conditional independence test
can have power against all alternatives. Here, the challenge lies in devising
and applying (a combination of) conditional independence tests that are the
most suited for a particular physical system.

The concept of post-selection inference (Berk et al., 2013) involves
performing statistical inference on a model that was selected based on data-
driven methods rather than being pre-selected. While there are some advances
in solving the post-selection inference problem for regression and causal
effect estimation, few solutions have been proposed for the inference after
causal discovery setting (Berk et al., 2013; Belloni et al., 2014; Rinaldo et
al., 2019). One possible solution is sample-splitting, but this is often
statistically inefficient. A recent development is the randomised version of
the greedy equivalence search (GES) algorithm, which allows for finite-sample
correction of classical confidence intervals (Gradu et al., 2022).

###  2.3 Opportunities for the physical sciences

The field of causal discovery from observational data is still in its infancy
but growing in methodologies, theoretical guarantees of performance, and
empirical evidence. Causal inference, in general, is a vast field that offers
alternatives and scientific opportunities that we review in what follows. Only
revisiting the whole body of empirical science based on association would take
a village, but the advances would pay off.

####  2.3.1 Causal hypothesis testing and targeted interventions

“Observational causal discovery offers revolutionary opportunities to test
hypotheses, evaluate the impact of interventions, attribute extreme events
with counterfactuals, and characterise complex systems by deriving causal
pathways and robust forecasting models.”

Scientists need a principled way to test different hypotheses against each
other. A causal hypothesis is a supposition or theory about how things
interact, specifically on whether one thing causes another. Causal studies aim
to confirm or reject any given causal hypothesis. The problem is that
hypotheses in the physical sciences are often presented as narratives giving a
chain of causal factors that lead to the studied phenomenon. Without a causal
vocabulary and analytical tools, it is often impossible to precisely state the
hypothesis, which leads to several competing hypotheses or, even worse, a
false hypothesis, which is accepted as true due to its compelling narrative
quality. Testing hypotheses have been conducted in myriad ways (Pearl, 2009c;
Peters et al., 2017a; Robins and Hernan, 2020).

Causal graphs, as graphical representations of assumed or learned causal
relations, provide a more principled way to talk about causal hypotheses.
Learned graphs imply causal links and pathways and provide evidence for
deciding between rivalling causal hypotheses, in Kretschmer et al. (2016), for
instance, regarding competing hypotheses of Arctic climate teleconnections.

The conclusiveness and interpretability of discovered causal graphs from
purely observational data sets on the often untestable validity of the
methods’ assumptions and the statistical complexity of the task. But
observational causal discovery can help inform more targeted subsequent
interventions, which are often too expensive to employ on a large scale
(Pearl, 2009c; Robins and Hernan, 2020). Incorporating interventions, if
performed meaningfully, could thus make the causal discovery process much more
efficient and robust (discovered DAGs not being confined to the Markov
equivalence class). Interestingly, interventions could be differentiable,
i.e., ‘learnable’ from data (Brouillard et al., 2020).

####  2.3.2 Cause-effect estimation

Causal discovery results in qualitative causal graphs, or often Markov
equivalence classes of graphs. But often, the target question is a
quantitative estimate of a causal effect of one variable $X$ on another
variable $Y$, as pioneered by Pearl (2009c). This topic is discussed in a time
series context in Runge et al. (2023). The quantity of interest then is the
(interventional) distribution of $Y$ given an intervention in $X$,
$p(Y=y\leavevmode\nobreak\ |\leavevmode\nobreak\ do(X=x))$. The fundamental
problem is that typically $p(Y=y\leavevmode\nobreak\ |\leavevmode\nobreak\
do(X=x))\neq p(Y=y\leavevmode\nobreak\ |\leavevmode\nobreak\ X=x)$.
Confounders, for example, can introduce a non-causal association between the
treatments and the outcome. Randomised experiments would be the gold standard
by eliminating the unwanted non-causal associations (Fisher, 1935; Imbens and
Rubin, 2015; Pearl, 2009c). The goal of causal effect estimation is to do so
without access to interventions by expressing $p(Y=y\leavevmode\nobreak\
|\leavevmode\nobreak\ do(X=x))$ as a function of the observational
distribution $p(\mathbf{x})$:

$\displaystyle p(Y=y\leavevmode\nobreak\ |\leavevmode\nobreak\ do(X=x))$
$\displaystyle=\text{function of}\leavevmode\nobreak\ p(\mathbf{x})\,.$ (2.8)

If such a re-expression is possible, one calls the causal effect identifiable
and obtains a causal estimand, which involves only the observational
distribution. The most well-known method for causal effect estimation from
data without parametric assumptions is _covariate adjustment_ (Pearl, 2009c),
which refers to de-confounding the causal relationship by adjusting for a set
of variables $\mathbf{Z}$. In the general case, the adjustment formula is

$\displaystyle p(y\leavevmode\nobreak\ |\leavevmode\nobreak\ do(X=x))=\int
p(y\leavevmode\nobreak\ |\leavevmode\nobreak\
x,\mathbf{z})p(\mathbf{z})d\mathbf{z}\,.$ (2.9)

Recent work has focused on finding statistically optimal adjustment sets
(Runge, 2021), i.e., for which the estimators have minimal variance. Using the
do-calculus (Pearl, 1995; Huang and Valtorta, 2006; Shpitser and Pearl, 2006,
2008), it is possible to determine whether a causal effect is, in principle,
identifiable from observational data or not. To this end, causal effect
estimation requires fully specified causal graphs from causal discovery (with
its inherent reliance on further assumptions) or domain expertise that can
qualitatively specify a causal graph. For example, it is known that
temperature influences ecosystem respiration, but one may want to quantify how
much when given a graph of other observed and unobserved confounding
variables. The graph then encodes assumptions about the absence and presence
of causal relations.

Different variants of causal effects can be defined based on the
interventional distribution $p(y\leavevmode\nobreak\ |\leavevmode\nobreak\
do(X=x))$, and an estimate then involves further parametric assumptions. For
example, in a linear model, the total causal effect on the expected value of
$Y$ when setting $X$ by intervention to $x^{\prime}$ as opposed to $x$ is
given by

$\Delta_{X\rightarrow Y}(x^{\prime},x)=\,\Delta x\cdot\beta_{X{\to}Y}\,,$
(2.10)

where $\Delta x=x^{\prime}-x$ and $\beta_{X{\to}Y}$ can be estimated as the
regression parameter of $X$ in the linear regression of $Y$ on
$X\cup\mathbf{Z}$.

####  2.3.3 Causal pathway analysis and mediation

Next to quantifying the overall causal effect of $X$ on $Y$, a relevant
follow-up question is often about the causal pathways: the mechanisms by which
this effect propagates. In complex systems, it is often interesting to analyse
how perturbations spread throughout the systems and through which subprocesses
perturbations are mediated (Runge et al., 2015b; Runge, 2015). Within the
structural causal model framework, mediation formally leads to counterfactual
quantities, see, for example, VanderWeele (2015) and briefly below in Sec.
2.3.8. But for linear models, the mediated causal effect (MCE) of $X$ on $Y$
that passes through a mediator $M$ (here $X,Y,M\in\mathbf{V}$) can be computed
by summing up the contributions along all paths passing through it:

$\displaystyle{\rm MCE}(X,Y|M)$
$\displaystyle=\sum_{\pi^{M}_{k}}\prod_{\lambda_{i\to
j}\in\pi^{M}_{k}}\beta_{i{\to}j}\,,$ (2.11)

where the summand iterates over causal paths $\pi^{M}_{k}$ from $X$ to $Y$
through $M$ and the product is over all links $\lambda_{i\to j}$ on each path.
The link coefficient $\beta_{i{\to}j}$ can be estimated as the regression
coefficient of $V^{j}_{t}$ in the linear regression of $V^{j}_{t}$ on the
parents of $V^{j}_{t}$. Mediation analysis can also answer the complementary
question: how strong is the direct effect of $X$ on $Y$.

####  2.3.4 Identifying causes and pathways leading to anomalies

Anomaly detection (Chandola et al., 2009) is an important problem in
engineering and the physical sciences. In Earth sciences, extreme events form
a subclass of anomalies and can be structured across different dimensions,
such as compound extremes (Zscheischler et al., 2020). While detecting
anomalies is an important problem, it does not answer the often relevant
question of what causes a particular anomaly or, more generally, what causes
the anomalous process. In engineering, business science, and healthcare, a
related problem is _root cause analysis_ (Andersen and Fagerhaug, 2006).
Causal discovery can address the problem of identifying causal drivers
(parents) or indirect mediating pathways and facilitate quantitative analyses
to analyse the contribution of different physical drivers in causing an
extreme.

####  2.3.5 Causal complex network analysis

Complex systems are often viewed as networks of interacting subprocesses, for
example, the human brain (Bullmore and Sporns, 2009), or the Earth system
(Donges et al., 2009a; Ludescher et al., 2021). Tools of network theory
(Boccaletti et al., 2006) have been used to analyse quantities such as the
information flow as it propagates through the system or the stability of
subprocesses (Gozolchiani et al., 2011). A common network measure is the node
degree, which quantifies the number of processes linked to a node. A more
involved measure is betweenness centrality, which quantifies the number of
shortest paths through a particular node. A crucial question is then to define
what these paths mean. In works where the networks are based on pairwise
correlation or mutual information (Donges et al., 2009a; Ludescher et al.,
2021), one may associate paths with a transfer of information.

However, there is a difference between information being transferred versus
perturbations propagating through the network. Here a question can be to
identify how critical individual subprocesses are in spreading and mediating
perturbations in such dynamic complex systems. The propagation of
perturbations, aka interventions, relates to a causal question requiring a
causal definition of network links able to distinguish direct from indirect
interactions.

In addition, the toolbox of classical network measures is not rich enough for
quantifying gateways and mediators of perturbations. Essentially, these
measures—with many originating from the social sciences (Freeman, 1977)—are
based on a different definition of links, for example, two persons knowing
each other, as opposed to dynamical interactions in a complex system. Hence,
here measures based on causal pathways on which perturbations propagate in a
complex system’s interaction network can be utilised, such as those studied in
Runge et al. (2015b); Runge (2015): Identifying the nodes of the causal graph
with the components of the complex system, the average causal effect can be
defined as a causal version of the out-degree or closeness centrality, which
quantifies by how much an individual component causes any of the remaining
components. This serves as a quantitative measure of how much a component is a
gateway of perturbations. On the other hand, the average causal susceptibility
measures how much a component is changed on average by a perturbation in any
of the remaining components as a causal version of the in-degree or in-
closeness. Finally, the average mediating causal effect measures how much of
the pairwise causal effects between any pair are mediated through a particular
variable, which can be seen as a causal version of betweenness centrality. In
Runge (2015), these measures are generalised in an information-theoretic
framework.

####  2.3.6 Causally robust forecasting models

Forecasting a time series from multivariate predictors constitutes another
problem where causal knowledge helps. Even considering the case of forecasting
inside the same distribution, that is, assuming a stationary distribution, it
can be proven information-theoretically that causal predictors maximise the
mutual information with the target variable and, by the Markov property, any
further predictors do not add further information. More formally (Brown et
al., 2012), the negative log-likelihood can be decomposed as follows

$\displaystyle\lim_{n\to\infty}-{l}$
$\displaystyle=\lim_{n\to\infty}-\frac{1}{n}\sum_{t=1}^{n}\log\widehat{p}\left(X_{t+1}\leavevmode\nobreak\
|\leavevmode\nobreak\ \mathcal{P};\theta\right)$ (2.12)
$\displaystyle=\mathbb{E}\left[\log\frac{p(X_{t+1}\leavevmode\nobreak\
|\leavevmode\nobreak\ \mathcal{P})}{\widehat{p}(X_{t+1}\leavevmode\nobreak\
|\leavevmode\nobreak\
\mathcal{P};\theta)}\right]\quad+\quad\underbrace{I\left(X_{t+1};\mathbf{X}^{-}_{t+1}\setminus\mathcal{P}\leavevmode\nobreak\
|\leavevmode\nobreak\ \mathcal{P}\right)}_{=0}\quad+\quad
H\left(X_{t+1}\leavevmode\nobreak\ |\leavevmode\nobreak\
\mathbf{X}^{-}_{t+1}\right)\,,$ (2.13)

where $n$ is the sample size, $\widehat{p}\left(X_{t+1}\leavevmode\nobreak\
|\leavevmode\nobreak\ \mathcal{P};\theta\right)$ the prediction model for
$X_{t+1}\in\mathbf{X}_{t+1}$ of the true underlying
$p\left(X_{t+1}\leavevmode\nobreak\ |\leavevmode\nobreak\ \mathcal{P}\right)$
given its causal parents $\mathcal{P}$ and model parameters $\theta$. As
shown, the log-likelihood decomposes into the model approximation error given
$\mathcal{P}$ (first term), the conditional mutual information between the
target and unselected variables $\mathbf{X}^{-}_{t+1}\setminus X_{t+1}$ given
$\mathcal{P}$ (second term, zero by the Markov condition), and the irreducible
entropy or uncertainty (last term). This is especially relevant for finding
optimal sets of predictors in the case where greedy selection strategies do
not work because the predictors cause the target variable synergistically, for
example, $X_{t+1}=Z^{1}_{t}\cdot Z^{2}_{t}+\eta_{t+1}$. As shown in Runge et
al. (2015a), an optimal subset selection can be better performed on the
smaller subset of causal predictors. In Kretschmer et al. (2017); Di Capua et
al. (2019), causal pre-selection was used in a climate context. Beyond
stationary distributions, Huang et al. (2019) address the task of causally-
informed forecasting under nonstationary environments through state-space
models.

####  2.3.7 Physical simulation model evaluation

Causal graphs and causal effects can be utilised to intercompare the output of
physical models and evaluate and validate them against observations at the
level of causal dependencies (Eyring et al., 2020, 2019; Nowack et al., 2020;
Pérez-Suay and Camps-Valls, 2019). One approach in this direction is to
compare the causal graph obtained from the observational data to those
obtained from simulated data. This procedure has been proposed in the climate
sciences to compare climate model simulations and observational data through
their corresponding causal graphs derived from PCMCI (Nowack et al., 2020),
cf. Section 4. The methodology could be adapted and applied to other physical
science problems where one typically has complex datasets and simulations to
confirm hypotheses.

####  2.3.8 Counterfactual causal attribution of extreme events

Counterfactual questions are not about the distributions of a target variable
due to possible (future) interventions, but about the distribution of a target
variable for an alternative past intervention, given that a particular outcome
was observed. Formally, just like interventional causal queries are
represented by interventional SCMs, counterfactual queries are represented in
counterfactual SCMs (Pearl, 2009c; Peters et al., 2017a). Given an SCM over
variables $\mathbf{V}$ and observations $\mathbf{v}$, in the _counterfactual
SCM_ , the noise distribution is updated such that the $\mathbf{V}=\mathbf{v}$
holds. Then the noise terms may not be independent anymore. Counterfactual
queries are then _do_ -statements in the counterfactual SCM. One example of a
counterfactual distribution query is $p(y^{\prime}_{x^{\prime}}|y_{x})$, which
specifies the probability of observing $Y=y^{\prime}\neq y$ under the
hypothetical past intervention $do(X=x^{\prime})$ when, in fact, $Y=y$ was
observed under the intervention $do(X=x)$. Such queries can be computed in
different ways (Pearl, 2009c; Correa et al., 2021) and generally require more
assumptions about the underlying structural causal model than causal effect
questions or causal discovery. Next to counterfactual distributional queries,
Halpern (2016) discusses causation and counterfactuals regarding single
events.

An example of a counterfactual question in climate is the causal attribution
of extreme events (Hannart et al., 2016). The above query
$p(y^{\prime}_{x^{\prime}}|y_{x})$ is one specific type of a counterfactual
question and is sometimes called “probability of necessity” (PN), which is
typically the quantity of interest in lawsuits. Extreme event attribution
requires to study of anthropogenic forcings compared to their absence, that
is, solely natural forcings or internal variability of the climate system. If
the probability of necessity is high enough, then a human-caused extreme event
is established.

####  2.3.9 Signal tracking for the discovery of proximal causes

Many phenomena, such as extreme events in complex systems, such as El Niño
events in the climate system and extreme volatility in the financial system,
are caused by an initial anomaly that triggers a travelling cascade of events
(Press, 1967). This phenomenon is often called the “butterfly effect”222The
term is attributed to Lorenz when he noted that a weather model failed to
reproduce the results of runs with the unrounded initial conditions. However,
the idea was earlier recognised by Poincaré and further formalised by Wiener.
The analogy became popular and originated the quantitative science of
characterising instability in complex systems undergoing nonlinear dynamics
and deterministic chaos., characterised by an anomaly in one part of a system
having extreme consequences in another space and time. Such cascaded events
are challenging to detect, predict, understand and characterise (Zscheischler
et al., 2018; Menzly et al., 2004), and has led to the development of the
field of the science surrounding the concept of predictability in complex
systems (Grunberg and Modigliani, 1954; Bialek et al., 2001; Boffetta et al.,
2002). A potential strategy for uncovering the cause of notable events is
causal discovery, for instance, by conducting a simulation which begins from
the start of the event in question and tracing the initiating signal back to
its source. However, this presents difficulty in determining the initiating
trigger. It is challenging to have a dialogue about recognising the drivers of
intense impacts because the amount of correlated drivers is usually much
greater than the number of causally pertinent drivers, which may only have a
substantial effect when combined (synergy) (Runge et al., 2019a). These kinds
of relations are hard to portray with a pairwise network.

####  2.3.10 Causal benchmarks, software and platforms

Method development and comparison require benchmark datasets with known causal
ground truth for validation. Ideally, such ground truth comes from expert
knowledge of real data or actual experiments that can also be used to falsify
causal relationships predicted from observational causal inference methods.
Unfortunately, in many fields, such as Earth system sciences, such datasets
exist only for expert-labelled causal relations among a few variables (e.g.,
some bivariate examples (Mooij et al., 2016)). A tractable approach is to
generate synthetic data with simple model systems that mimic properties and
challenges of data from the system under study but where the underlying ground
truth is known. These can then be used to study the performance of causal
discovery (and causal inference methods more generally) for different
challenges in realistic finite sample situations. From a practitioner’s
perspective, it is essential to determine which method is best suited for a
particular task with particular challenges and for a specific set of
assumptions. Synthetic data, adapted to the problem at hand, can be used to
choose the suitable method, including method parameters. A list of key methods
for causal discovery and the available software and platforms is given in
Table 2.2. An example is the SAVAR model (Tibau et al., 2022) that mimics
spatio-temporal features of climate data. The website causeme.net (Runge et
al., 2019a, 2020) aims to provide an open platform with synthetic models
mimicking real data challenges on which causal discovery methods can be
compared. Next to method comparison, the platform also calls for submissions
of actual and modelled data sets where the causal structure is known with high
confidence and was used on the Causality 4 Climate NeurIPS competition (Runge
et al., 2020). That competition sparked the investigation of a particular
property of synthetic data and models called _var-sortability_ , which led to
new insights in causal discovery methods (Reisach et al., 2021).

Table 2.2: Methods and open-source software for causal discovery. Method |  Software
---|---
Granger causality (GC) (Granger, 1969), kernel GC (Marinazzo et al., 2008a), explicit KGC (Diego Bueso, 2020) |  causal-learn, statsmodels, KGC, XKGC
CCM (Sugihara et al., 2012; Ye et al., 2015) |  rEDM
PC (Spirtes and Glymour, 1991; Spirtes et al., 2000), FCI Entner and Hoyer (2010) |  Tetrad, causal-learn, pcalg, Tigramite, MXM, bnlearn, dbnlearn, PyWhy
PCMCI (Runge et al., 2019b), PCMCI+ (Runge, 2020), LPCMCI (Gerhardus and Runge, 2020) |  Tigramite
DYNOTEARS Pamfil et al. (2020) |  Causalnex
TiMINo Peters et al. (2013) |  R script
VARLiNGAM (Hyvärinen et al., 2010) |  causal-learn, lingam, original R code
ICP Peters et al. (2016); Heinze-Deml et al. (2018b); Pfister et al. (2019) |  seqICP

###  2.4 Perspectives

This section reviewed causal discovery in the physical sciences, describing
the main methods, challenges, and opportunities for future research. We laid
out the fundamental elements of the causal discovery framework –SCMs, graphs,
and associated distributions–, and gave an overview of the methodological
concepts of learning qualitative causal graphs (Runge et al., 2019a, 2023). We
deliberated on commonplace difficulties encountered in the field, such as
determining and preprocessing causal variables, addressing non-stationarity,
contemporaneous causation and hidden confounding, and selecting parametric
models for nonlinear dependencies and non-Gaussian distributions. Section 4
will illustrate causal discovery methods in neurosciences and Earth sciences
case studies.

The body of causal inference has traditionally been embedded in several
communities, mainly statistics, social sciences, econometrics and health
sciences. Irreconcilable positions and long-standing discussions exist (Dowd,
2011). Pearl argues that it is essential to distinguish between causal and
statistical information, as they refer to two separate concepts333Statistical
information deals with the probability of certain variables being observed. In
contrast, causal information deals with hypothetical relationships in new
situations., and suggests that clear distinctions should be established in the
notation used, and each should be subject to different means of calculation
(Pearl, 2011). Arguably, nonparametric SCMs (as a natural generalisation of
those used by econometricians and social scientists in the 1950-1960s) have
developed the field of causality in new mathematical underpinnings: explicate
and enumerate causal assumptions, test implications, decide measurements and
experiments, recognise and generate equivalent models, recognise instrumental
variables, generalise structural equation models and solve the mediation and
external validity problems. These tools, methods and solutions help to
determine the accuracy and validity of causal claims in the analysis. The
machine learning community is approaching the field of causal discovery in
innovative ways by leveraging data, assumptions and models collectively. In
recent decades, mathematical foundations have been established to address
questions of causality in various scientific fields, mainly emerging for
statistics and machine learning (Pearl, 2009c; Spirtes et al., 2000; Peters et
al., 2017a). The causal machine learning (CausalML) field has recently
introduced (Kaddour et al., 2022) as an umbrella for machine learning methods
based on SCMs. It aims to advance the field in several directions: causal
supervised learning, causal generative modelling, causal explanations, causal
fairness, and causal reinforcement learning. Applications of the new methods
are vast and promise advances in computer vision, natural language processing,
and graph representation learning. Therefore, the field of causal discovery is
growing in methods, approaches and impactful applications. A unified agenda
for Causal Inference is built and deployed in the wild.

However, despite the significant advances in the last decades, many unresolved
philosophical and methodological issues remain for causal discovery from
observational data. All such challenges also create avenues of research. On
the one hand, we identified and discussed algorithmic and data challenges and
summarised possible ways to address them in §2.2. Indeed, we must develop more
effective methods for incorporating (uncertain) expert knowledge, determining
the spatio-temporal complexity of the underlying dynamic phenomena, and
creating more reliable and statistically efficient algorithms.

On the other hand, perhaps the most critical challenge is the theoretical
impossibility of causal discovery from purely observational data (Pearl,
2009a). There are, however, ways to tackle the challenge. For example,
specifying a causal DAG using domain knowledge can help mitigate the potential
inaccuracy of their assumptions of sufficiency and faithfulness. Another
possibility to learn about a certain equivalence class may consider
incorporating domain knowledge into structure learning algorithms by using
“allow lists” and “deny lists” to determine which edges should or should not
be included in a DAG or creating a Bayesian prior to assigning varying levels
of probability to certain causal relationships (Castelo and Siebes, 2000; Ness
et al., 2017). This is very much related to using inductive biases (such as
Occam’s razor) (Janzing, 2007) and causal invariances (such as parameter
modularity and independence of mechanism) (Hoyer et al., 2008) to learn
structure beyond likelihood-based scores and conditional independence
constraints. Finally, if data from natural experiments, such as $do(A=a)$, is
available, this intervention information can be incorporated into the
algorithm (Cooper and Yoo, 1999) to automate this reasoning process.
Incorporating the abundant domain knowledge within the causal discovery
routine can address the identifiability and faithfulness assumptions (very
much in line with the basis or sparsity priors used in equation discovery, cf.
§3). By joining forces, both can contribute to resolving pressing scientific
issues, ranging from process comprehension to evaluating and upgrading the
physics included in physics models.

Many problems in the physical sciences can be framed as causal questions. Yet
many researchers in economics and health services, and even many computer
scientists in machine learning, have been trained to be reluctant to use the
language of causality (Dowd, 2011; Pearl, 2011; Hernán, 2018; Schölkopf et
al., 2021). This is a cognitive barrier to resolve in the future. Besides, the
language barrier between the methodological and domain science communities is
a significant challenge in the causal discovery endeavour. Bridging this
divide by translating domain questions into actionable and precisely stated
causal inference tasks seems reasonable. Additionally, hesitance to adopt
causal inference can be attributed to the lack of suitable benchmarks to help
choose an appropriate method. A benchmarking platform (https://causeme.net)
was introduced that covers the causal discovery problem setting. It is
necessary to have more of these benchmark platforms and easily accessible
databases to facilitate better collaboration between the two communities. To
successfully address a causal inference problem, it is essential for the
domain scientist and computer scientist to work together; assumptions must be
discussed and formalised, data characteristics must be jointly analysed, and
conclusions must be assessed from both perspectives.

It is commonly believed that most research questions in science can be
interpreted as causal inference problems. Sentences such as “we find that $X$
[increases/ decreases/ lags/ leads/ affects/ drives/ impacts] $Y$” are often
found in papers, creating an impression of causality. Nevertheless, scientists
should be more explicit and transparent when making assumptions which lead to
causal conclusions. This is not only sensible but is also necessary when
analysing complex systems like the Earth, the Brain, or the Economy, since the
outcomes of such research may have significant economic, environmental, and
social implications. The field of Causal Inference provides a comprehensive
arsenal of tools grounded in rigorous mathematical principles and a vibrant
interdisciplinary milieu to confront the challenge at hand successfully.

##  Chapter 3 Learning physical laws from data

“For centuries, scientists observed Nature to extract simple laws and
equations to explain the world mechanisms, to anticipate and predict
behaviours and gain faith in their interventions/actions.”

Distilling mechanistic models of the world is how physicists have successfully
understood and explained the natural world. The prototypical process starts
with an experiment or observation. A mathematical model is hypothesised, which
can predict a new experiment’s outcome. The observations will either support
or falsify the hypothesis, leading to more experiments and refined models.

For many complex systems, we have poor models because certain interaction
terms are unknown. Another case is when we know some microscopic interaction
laws. Still, the emergent properties at a larger scale do not directly follow,
such that predictions at the larger scale need new coarse-scale interactions.
To cope with these situations and make sense of the sheer amount of data
produced by modern instruments, researchers have been looking into automating
the processes involved in model building and creating new insights. Many
motivations and inspirations have been adopted to guide the scientific method
development (see Fig. 3.1): e.g. a physical law or equation describing the
system should be compositional, thus following the “divide and conquer” rule,
should be as simple as possible (but not simpler) thus following Occam’s
razor, further developed and formalised by Solomonoff, eventually focus on the
interesting parts of the system and disregard the rest, create understanding,
generalise and unify different working theories and models, and the learned
representation (and its parameters) should be self-explainable, amenable and
intuitive.

Figure 3.1: Historical Inspirations/Motivations in Law/Equation Discovery.

This section reviews the state-of-the-art in equation discovery from data.
Unlike in traditional law discovery à la Kepler where trial-and-error was
dominant, modern statistics and machine learning techniques exploit the
regularities found in the data to discover plausible, simple and explainable
equations, to learn feature representations that describe (typically dynamic)
systems. We will first consider the explicit discovery of equations that
describe observed data, also called _Symbolic Regression_. Second, we look
into implicit discovery through dimensionality reduction techniques and
transfer operators. We finish the section by discussing the main challenges
and research opportunities.

###  3.1 Explicit equation discovery with symbolic regression

Symbolic Regression (SR) refers to a class of machine learning techniques that
aim to discover mathematical relationships and patterns in data. SR aims to
find a compact, human-readable mathematical expression that accurately
reflects the underlying relationships in the data. This approach is
particularly useful in cases where the relationships between variables are
complex or unknown, and traditional statistical methods may not provide
adequate explanations. SR operates under the assumption that the underlying
data-generating mechanism can be described by a sparse and algebraic input-
output relationship.

There is a major divide in the method to achieve that. One class of methods
use genetic algorithms or other discrete search methods to find mathematical
terms, typically represented by a graph of mathematical operations. Another
class of methods uses continuous space search methods and solves a relaxation
of the discrete problem of finding compact equations. A third and most recent
addition to the family of symbolic regression methods uses massive amounts of
synthetic data for pretraining a system that can quickly guess a suitable
expression at test time.

More formally, we are trying to solve the following optimisation problem:

$\displaystyle\operatorname*{arg\,min}_{f}\mathbb{E}_{(x,y)\sim\mathcal{D}}\|f(\bm{x})-\bm{y}\|^{2}+\lambda
C(f)$ (3.1)

where we attempt to find a low-complexity function (equation) $f$ that
best111We use the squared error here for simplicity, but other notions of
distance are possible. maps the inputs $\bm{x}\in\mathbb{R}^{n}$ to their
corresponding outputs $\bm{y}\in\mathbb{R}^{m}$ in the data distribution
$\mathcal{D}$, $C(f)$ refers to a measure of complexity of $f$, and $\lambda$
is the weighting factor. For instance, the complexity measure could be the
number of terms in the equation.

A central problem when performing symbolic regression is selecting an
appropriate weighting factor. More generically, the question is which level of
complexity is right. There is probably no definite answer to this question.
Instead, we consider the solution to a symbolic regression problem as a family
of Pareto-optimal solutions:

$\displaystyle
f^{(c)}=\operatorname*{arg\,min}_{f}\mathbb{E}_{(\bm{x},\bm{y})\sim\mathcal{D}}\|f(\mathbf{x})-\mathbf{y}\|^{2}\qquad\text{s.t.}\quad
C(f)=c$ (3.2)

where $f^{(c)}$ refers to the best fitting expression with complexity $c$. A
more complex expression will be able to fit the given data at least as well as
a less complex expression. Fig. 3.2 illustrates a typical Pareto curve along
the optimisation objectives: goodness of fit (error reduction) and function
complexity. In addition to the training error, which monotonically decreases
with function complexity, the illustration also shows a hypothetical test
error that shows the overfitting of too complex functions.

Figure 3.2: Illustration of a Pareto curve of solutions to an SR problem.

We will now look more closely into different methods that have been proposed
to solve the optimisation problem (Eq. 3.2) in practice.

####  3.1.1 Symbolic regression using discrete search methods

The problem of symbolic regression is, at its core, a search for suitable
functions $f$ in (3.2). Since those functions should have low complexity, it
is natural to attempt to perform a search for functions. The first attempt to
do that was proposed by Cramer (1985) by inventing Genetic Programming, which
got popularised and applied through Koza (Koza, 1990, 1994). The idea is
simple: search for computer programs to solve a particular problem by
iteratively creating many random programs and selecting the best fit, and
create a new pool of candidates by recombination and random modification. This
mimics the biological evolution process of nature to create the genetic
material of living organisms. Applied to symbolic regression: The functions
are represented as a graph of input variables, operators and basic algebraic
functions.

In the paper by Schmidt and Lipson (2009b), this approach was refined and
applied to the discovery of physical laws. As the method was indeed able to
discover Lagrangian and Hamiltonian formulations from data, it stimulated a
growing interest in symbolic regression and sparked the development of many
methods. These general search methods are also referred to as _evolutionary
algorithms_. The approach of Schmidt and Lipson (2009b) is illustrated in Fig.
3.3. The general method for symbolic regression was implemented in a tool
called _Eureqa_ (Praksova, 2011) that is now only available as an online
service (DataRobot Inc, 2023).

1 Observational time series $\\{x(t),y(t),z(t)\\}$ |  2 Partial derivatives $\Big{\\{}\dfrac{\Delta x}{\Delta y},\dfrac{\Delta z}{\Delta x},\dfrac{\Delta y}{\Delta z}\Big{\\}}$ |  3 Candidate symbolic functions |  4 Symbolic partial derivatives |  5 Compare predictive par.der. |  6 Accurate and sparsest solutions
---|---|---|---|---|---
|  |  $f=x\cos(y)$ $f=0.7\exp(y/z)$ $f=y^{2}-9.8\cos(x)$ … |  $\dfrac{\partial f}{\partial y}=y+\sin(x)\dfrac{\Delta x}{\Delta y}$ $\dfrac{\partial y}{\partial x}\bigg{|}_{f(x,y)}=\dfrac{\partial f}{\partial x}/\dfrac{\partial f}{\partial y}$ | $\dfrac{\Delta y}{\Delta y}\bigg{|}_{D_{i}}=\dfrac{\partial y}{\partial x}\bigg{|}_{f_{i}}$ |  $f=z+9.8\sin(x)$ $f=y^{2}-9.8\cos(x)$
Figure 3.3: Schematic view of the symbolic regression method for discovering
physical laws in (Schmidt and Lipson, 2009b). Starting from observational data
(1), partial derivatives are computed numerically for all pairs of variables
(2). A set of candidate symbolic functions $f$ is derived (3), whose symbolic
partial derivatives are computed (4) and compared to the predictive ones (5).
The process 3-5 is iterated until, finally, a small set of the most accurate
and simple equations is returned (6).

There are several publicly available and open-source implementations, such as
the PySR (Cranmer, 2020), gplearn (Stephens, 2022), Glyph (Quade et al., 2019)
and Operon (Burlacu et al., 2020). A more detailed overview of genetic
algorithm-based methods and their combination with gradient descent can be
found in Kommenda et al. (2020).

##### Feynman AI.

An approach that exploits physical knowledge such as units and makes
reasonable assumptions for equations in physics is Feynman AI (Udrescu and
Tegmark, 2020; Udrescu et al., 2020). The method augments the genetic
algorithm searching for expressions by enforcing fitting physical units,
decomposing the problem using symmetries and checking separability. To check
for symmetries, a neural network is trained on the data to allow accessing
whether the underlying function is symmetric. It is worth noting that a large
amount of data is used here. From a set of 100 equations taken from the
Feynman Lectures, the method was able to recover all of them whereas Eureqa
only solves 71.

##### Search with Deep Reinforcement Learning.

A method that uses Deep Reinforcement Learning to search for a suitable
solution to the symbolic regression problem is Deep Symbolic Regression (DSR)
(Petersen et al., 2021). The key idea is to treat the search for expressions
as an exploration problem in reinforcement learning (RL). The functional
expressions are represented as a sequence of tokens corresponding to a depth-
first graph traversal and are generated by a recurrent neural network.
Numerical constants are fitted using the BFGS optimiser. This generative
network is trained on the given dataset using RL to find a highly fitting
solution. An interesting contribution is a formulation of a risk-seeking
policy gradient that tries to optimise for the best-case scenario (a good
solution can be found) rather than the typical average case. The method was
able to solve 83% of the standard Nguyen-1 dataset.

####  3.1.2 Sparse linear regression and neural network approach

The symbolic regression problem can also be tackled via traditional regression
methods. In contrast to the search in a discrete set of functions, the search
is performed in a dense set, typically represented by a real-valued
parameterised function. So (3.1) is solved by choosing a large enough function
class described by $f_{w}$ with $\bm{w}\in\mathbb{R}^{p}$. The optimisation is
then performed over the space of parameter-values $w$. Linear regression is a
special case, where the function $f_{w}(\bm{x})=\bm{w}\cdot\bm{x}^{\top}$.

What about the complexity regularisation term $C(f_{w})$ in (3.1)? Ideally,
the term should count the number of non-zero parameters in $\bm{w}$, expressed
as $|\bm{w}|_{0}$ and referred to as $L_{0}$ norm. However, this term
jeopardises the efficient solution of the regression problem because it is
non-linear and non-differentiable. One practical alternative is to use the
$L_{1}$ norm instead, i.e. the sum of absolute values, which also leads to
sparse solutions (see Fig. 3.4). In the case of linear regression, this is
termed LASSO regression (Tibshirani, 1996).

Figure 3.4: $L_{1}$ regularisation typically leads to sparse solutions. The
lines show the isolines of quadratic loss (red) and $|w|$ (blue). Instead of
the data-only solution (red star), a sparse solution (green) is found.

The methods differ in the function class $f_{w}$, the regularisation term and
the optimisation method used.

##### SINDy: Sparse identification of dynamical systems.

In some cases, the class of building blocks that might occur as summands in
the analytical description of the data are known. Then a rather simple but
effective method can be employed that is called _sparse identification of
dynamical systems_ , SINDy for short. It was proposed in Brunton et al.
(2016b) to find differential equations of dynamical systems from observations.
For the general symbolic regression problem, the FFX method by McConaghy
(2011) was already earlier proposing the same idea. The input data is passed
through a predefined library of base functions and interaction terms. Then the
resulting high-dimensional representation is fit to the data using sparse
linear regression. All relevant terms keep a non-zero weight and constitute
the final expression.

Let us unpack this in more detail for a dynamical system of $n$ variables
described by the system of ordinary differential equations
$\frac{d}{dt}\bm{x}=\bm{g}(\bm{x})$, where $\bm{x}(t)\in\mathbb{R}^{n}$. Each
component of $\bm{g}$ can now be substituted by a linear combination of
library functions:

$\begin{array}[]{rclcl}\dfrac{d}{dt}x_{1}&=&g_{1}(x_{1},x_{2},\ldots,x_{n})&=&w_{11}l_{1}(x_{1},...,x_{n})+\ldots+w_{1m}l_{m}(x_{1},...,x_{n})\\\\[5.0pt]
\dfrac{d}{dt}x_{2}&=&g_{2}(x_{1},x_{2},\ldots,x_{n})&=&w_{21}l_{1}(x_{1},...,x_{n})+\ldots+w_{2m}l_{m}(x_{1},...,x_{n})\\\
&\vdots&\\\
\dfrac{d}{dt}x_{n}&=&g_{n}(x_{1},x_{2},\ldots,x_{n})&=&w_{n1}l_{1}(x_{1},...,x_{n})+\ldots+w_{nm}l_{m}(x_{1},...,x_{n}),\end{array}$

where $l_{1},\ldots,l_{m}$ is the _predefined finite library of candidate
functions_ and $w_{ij}$ are the scalar coefficients to be learned following
our objective (3.1). To obtain a sparse solution, the $L_{1}$ regularisation
(LASSO) can be used, i.e. $C(f)=|\bm{w}|_{1}$, as described above.
Alternatively, the (squared) $L_{2}$ norm of the weights $\|\bm{w}\|^{2}_{2}$
can be used, corresponding to classical ridge regression that permits a
closed-form solution. However, an iterative pruning of small weights must be
used to obtain a sparse solution.

Illustration of sparse identification of dynamical systems (SINDy) (Brunton et
al., 2016b) using a synthetic dataset of the well-known Lotka-Volterra system,
as shown in Adsuara et al. (2020). The Lotka-Volterra system models the
interaction between prey and its predator in ecology and is given by the
following equations: $\begin{array}[]{rcl}\dfrac{d}{dt}x_{1}&=&\alpha
x_{1}-\beta x_{1}x_{2}\\\\[4.5pt] \dfrac{d}{dt}x_{2}&=&-\gamma x_{2}+\delta
x_{1}x_{2}\end{array}$ with the coefficients, $\alpha$ and $\gamma$, being the
intrinsic growth/decrease rates of $x_{1}$ and $x_{2}$, and $\beta$ and
$\delta$ are cross terms taking into account the interaction between species.
In our particular case, we will set $\alpha={3}/{2}$, $\beta=1$, $\gamma=3$,
and $\delta={1}/{2}$. We show the results of the identification of these
parameters using SINDy in the table below for two levels of additive white
Gaussian noise of the signal-to-noise ratio of $5$ (high noise level) and $40$
dB (low noise level). As usual, the created data was split into train/test
data ($75\%,25\%$), respectively. The ODE coefficients are recovered
sufficiently well to achieve a high correlation coefficient $R$ but are
generally underestimated due to the sparsity regularisation. Learned
Coefficients True Library functions 40 dB 5 dB Coefficients
$\frac{d}{dt}x_{1}$ $\frac{d}{dt}x_{2}$ $\frac{d}{dt}x_{1}$
$\frac{d}{dt}x_{2}$ $\frac{d}{dt}x_{1}$ $\frac{d}{dt}x_{2}$ $x_{1}$ 1.3822 0
1.1404 0 1.5 0 $x_{2}$ 0 -2.9123 0 -2.7946 0 -3 $x_{1}x_{2}$ -0.9797 0.4849
-0.9520 0.4710 -1 0.5 $x_{2}^{3}$ 0 0 0 -0.0001 0 0 R $0.9999$ $0.8674$

For fitting dynamical systems, the temporal derivatives need to be computed.
Finite differences are often too sensitive to noise such that kernel
regression (aka Gaussian processes), which allow for explicit derivative
computation, are preferred (Camps-Valls et al., 2016; Johnson et al., 2018). A
more recent approach is to solve noise estimation and model identification in
one joint optimisation procedure Kaheman et al. (2022). Intuitively, for every
data point, the corruption by noise is estimated. Since this optimisation
problem is highly underdetermined, an additional constraint is used, namely
that when integrating the estimated dynamical system model a small error
should occur. This trick separates noise from the signal and leads to an
improved estimation quality.

##### Neural network approach: Equation Learner.

Enlarging the function class $f$ is possible using neural networks. Probably
the first work in this direction is the Equation Learner (EQL) introduced in
Martius and Lampert (2016) that uses a neural network with algebraic base
functions and a particular regularisation scheme to solve (3.1) and 3.2. The
function $f$ is represented by a neural network, modified only to contain
elementary operations that should appear in a potential solution. Figure
3.5(left) shows the architecture of the Equation Learner (EQL) in a simplified
form.

Example:
$y=\nicefrac{{1}}{{3}}\left((1+x_{2})\sin(\pi x_{1})+x_{2}x_{3}x_{4}\right)$

Figure 3.5: Equation Learning Architecture and example equation. Left:
Illustration of the EQL Architecture, reproduced from (Sahoo et al., 2018), a
feed-forward neural network with special _activation_ functions (sin,
multiplication etc.). Note that each unit type will occur many times. Right:
Example system with four inputs $x_{1,2,3,4}$ and one output $y$. Training is
only in the $[-1,1]^{4}$. EQL recovers the equation and extrapolates (Sahoo et
al., 2018).

The input variables are mapped with a dense layer to multiple instances of
trigonometric functions, identity, multiplication and division, but more base
functions, such as squares or exponentials, are possible. The resulting values
are again mapped with a matrix to another layer of elementary functions, and
so forth, until the last layer corresponds to the output (containing only
division operators in the picture).

The network is trained using stochastic gradient descent (e.g. Adam) on the
mean squared error loss and $L_{1}$ regularisation on the weights to induce
sparsity:

$\displaystyle\mathbb{E}_{(\bm{x},\bm{y})\sim\mathcal{D}}\|f_{W}(\bm{x})-\bm{y}\|^{2}+\lambda|\bm{W}|$
(3.3)

where $f_{W}$ denotes the neural network with parameters $\bm{W}$.

Toy data

Mauna Atmospheric CO2 concentration dataset

Figure 3.6: Illustration of uncertainty estimates using a mixture of Laplace
approximations of learned equations. Top row: toy example $y=0.8\cos
x-0.4+\epsilon$ where $\epsilon\sim\mathcal{N}(0,0.03^{2})$ with just 6
datapoints. Bottom row: Atmospheric CO2 concentration at Mauna Loa Observatory
(Observatory, 2020) (concentration vs. time, both in arbitrary units). The
left panels show the predictive distributions. The panel in the middle shows
individual local Laplace approximations with $2\sigma$ (shaded area) for the
toy data and a zoomed density for the Mauna dataset. The colour represents the
weight and aligns with the Pareto plots on the right side, showing RMSE over
the complexity of each equation. Reproduced from (Werner et al., 2022).

Note that the system needs to be differentiable for training, such that a pure
complexity term, such as $L_{0}$ regularisation that would count the number of
non-zero weights, does not work out of the box. Although methods have been
developed since then Louizos et al. (2018), the $L_{1}$ regularisation in
(3.3) is effective but creates an undesired trade-off between error and
sparsity. Something that we also encountered in the illustrating example when
applying SINDy. The EQL method introduces an additional regularisation phase
after converging with $L_{1}$ that clamps all $|W_{i}|<\epsilon\ll 1$ to zero
and optimises without regularisation. This yields a practical approach to
optimising for sparsity without trade-offs. In Kim et al. (2021) an
alternative to $L_{1}$ with $L_{0.5}$ was used.

The reader may wonder how the system is successfully trained with elementary
functions such as division or square root. Indeed, a naive application would
fail due to exploding values or gradients. In Sahoo et al. (2018) and Werner
et al. (2021), suitable parameterisations and training steps are proposed.
Choosing different $\lambda$ (Eq. 3.3) will create differently sparse
resulting networks. Each represents a particular symbolic expression
resembling the Pareto curve illustrated in Fig. 3.2. In Fig. 3.5 (right), a
synthetic example system is shown. The training data is only generated in the
$[-1,1]^{4}$ hypercube. The correct equation was discovered, and perfect
extrapolation is possible in this case, see Sahoo et al. (2018) for details.

Instead of manually selecting a particular solution, which might be a good
procedure when structural insights are to be obtained when investigating some
unknown phenomenon, one can also use several or all solutions along the Pareto
curve to estimate uncertainty about the predictions for extrapolation, as
proposed in Werner et al. (2022). Figure 3.6 illustrates this approach. For
each found equation, a Laplace approximation allows to approximate the
uncertainty due to parameter estimation errors and yields a Gaussian
posterior. Combining these using a weight based on the validation error and
the complexity yields the estimated density. Note how the uncertainty in
extrapolation shows clearly the structure of the discrete set of automatically
generated hypotheses.

####  3.1.3 Learning to solve symbolic regression

All methods so far treat every symbolic regression problem in isolation – the
search or optimisation algorithm was applied to a new dataset from scratch. We
are now looking into the idea of learning to solve a particular problem
quickly by using data from a whole class of symbolic regression instances.
Generally, the idea is to approximate the inverse mapping from data to a
suitable equation. The Dreamcoder paper by Ellis et al. (2021) showed the
first instantiation of this idea. Provided with the language of algebraic
expressions (arithmetic operations, variables, base functions) and a
_simulator_ to generate data for a particular equation instance, the method
learns a probabilistic mapping from data to equation terms and a library of
common equation building blocks. Given a particular instance of data, the
system can relatively quickly guess and verify suitable explaining equations.
Developing the idea of pretraining further and specialising it for symbolic
regression was done by the following method.

##### NeSymReS.

The approach in Biggio et al. (2021) is to use a high-capacity transformer
model pretrained to solve the symbolic regression problem. The method is
called _Neural Symbolic Regression that Scales_ (NeSymReS). The method uses a
large set of symbolic regression problems to approximate the inverse mapping
from data to equations. After this pretraining phase, given new data, the
inverse mapping can generate likely candidate equations. Intuitively, an
experienced data scientist might solve the problem similarly: looking at the
data and postulating a particular functional form that might explain it,
testing it, and potentially trying a different plausible hypothesis.

Figure 3.7: Overview of the NeSymReS method. Left: step with randomly
generated training data. Right: inference of candidate equations for unseen
data. Reproduced from (Biggio et al., 2021).

Let us look closer at the method. As visualised in Fig. 3.7, the core is a
transformer222The transformer architecture is the basic building block of many
large-scale machine learning systems, such as GPT-3 (Brown et al., 2020).
architecture (Vaswani et al., 2017) that can generate algebraic expressions
symbol-by-symbol given a set of data points.

The input is represented as a set of $(\bm{x},\bm{y})$ pairs (1024 in Biggio
et al. (2021)), which are processed through a set-encoder. The latter is
invariant to permutations of the data points. The output of the transformer
are tokens that correspond to the typical symbols of input variables, base
functions, operators and constant placeholders, which resemble the skeleton of
the predicted function. Importantly, the transformer does not have to guess
the right constants, just their location in the expression, as these constant
placeholders are fit to the data using non-linear optimisation (here BFGS).

Trained on millions of synthetically generated pairs of random expressions
with corresponding data, the transformer does a remarkable job guessing likely
equations. Importantly, the prediction is not deterministic but allows
sampling of possible functions. Thus, new potential solutions can be generated
and validated when new data is presented at test time until a sufficiently
good fit is found or the Pareto.

Figure 3.8: Performance of NeSymReS, DSR, classical SR (using gplearn
(Stephens, 2022)), and Gaussian processes on the AI Feynman and Nguyen
datasets (equations unseen during training). Reproduced from (Biggio et al.,
2021).

Figure 3.8 shows the accuracy of different SR methods for unseen equations
from the Feynman and Nguyen benchmarks. The performance is presented in
dependence on wall-clock time. NeSymReS is remarkably fast at finding a well-
fitting expression for the data.

As a downside, the method was only shown for three input variables, and it
remains to be seen how much it can be scaled in this respect. Also, the
dataset used to guess an equation at test time cannot be big (currently in the
order of 1000 data points) because the set transformer encoder cannot yet
handle larger sets well.

####  3.1.4 Comparison

As there are quite a number of methods, we aim to discuss their differences,
strengths and weaknesses by comparing them along a set of axes. We start with
using domain knowledge, as we seldom face a completely uninformed setting in
physics. We continue with aspects of the embedding, scaling, speed and
usability, summarised in Table 3.1.

##### Using domain knowledge.

A common form of domain knowledge is the base functions and their approximate
frequency of occurrence in describing the system under consideration. In
standard symbolic regression with genetic algorithms, the number and kind of
base functions are very flexible, and each term can have its individual
penalty in terms of complexity. More specific domain knowledge, such as
monotonicity, function image constraints and derivative constraints, can also
be considered, as presented in Kronberger et al. (2022). Another recent work
is presented in Cornelio et al. (2023) that allows to incorporate axiomatic
contraints.

In FFX/SINDy, the library of functions is the prime way to specify domain
knowledge. Relative preferences could be implemented by different
regularisation strengths.

For the EQL framework, the choice of base functions is a bit more complicated,
as the systems need to remain optimisable with gradient descent. In Werner et
al. (2021), a suitable relaxation for functions with divergences (in function
value or derivative) is proposed, and a way to specify preferences among base
functions is analysed. The control of the relative frequency of used terms is
possible but less direct than in genetic algorithm-based methods. In NeSymReS,
domain knowledge can be embedded by selecting/generating the training set with
appropriate synthetic problems, although this was not explicitly demonstrated.

##### Scaling.

Most SR methods are for small-scale problems with a few hundred to a few
thousand data points and low-dimensional problems, i.e. 1-10 input and output
variables. Classical search methods scale unfavourably with dimensionality as
the search space grows exponentially. That is why most SR methods are good at
finding relatively small and compact equations for low-dimensional systems but
fail for both high-dimensional systems or those where larger equations are the
most compact solution.

FFX/SINDy can handle large output dimensions easily, using some form of ridge
regression. However, it also suffers from large input dimensions as the
library bank becomes exponentially large unless a factorisation or other
simplifying structure is known, for instance, a strong locality assumption in
PDEs.

NeSymReS is also limited to several variables and small data sets. Although,
the limitation of the data set size can be lifted by sampling a smaller subset
of data points for guessing the skeletons and using all data for the parameter
tuning with BFGS.

EQL is the only scalable method, as gradient descent works on all dimensions
simultaneously, and machine learning methods are developed to scale. A larger
initial neural network should be used for larger systems, but no specific
adaptations are required.

Table 3.1: Comparison of symbolic regression methods. See Sec. 3.1.4 for more details. | embeddable | scaling | speed |  restriction |  domain knowledge
---|---|---|---|---|---
Genetic Programming (Schmidt and Lipson, 2009b) | ✗ | ✗ | slow |  for small systems |  base-functions, complexity of terms
AI Feynman (Udrescu and Tegmark, 2020) | ✗ | ✗ | slow |  for physical systems in canonical form |  physic: units, symmetries
DSR Petersen et al. (2021) | ✗ | ✗ | medium |  small input dim |  training domain
FFX (McConaghy, 2011), SINDy (Brunton et al., 2016b) | ✓ |  1✓333It scales well to large output sizes; for high-dimensional input strong structural assumptions are required. | blazing |  needs known library |  training domain
EQL (Sahoo et al., 2018) | ✓ | ✓ | slow |  base functions limited, sometimes less concise |  base-functions, complexities
NeSymReS (Biggio et al., 2021) | ✗ | ✗ | fast |  small input dim |  training set

##### Differentiability and Embeddability.

An interesting feature is whether the SR system can be used as a module in a
larger computational pipeline. For instance, if observations are available as
images and the causal variables have to be first extracted from the images
before they can be used for a concise SR prediction module. An example using
SINDy inside an autoencoder architecture (Champion et al., 2019) learns the
coordinate frame and dynamics equations at the same time, as detailed in Sec.
3.2.4 below. EQL is conceptually easy to embed into larger architectures, such
as deep networks, as it is end-to-end differentiable. An example is
discovering PDEs (Long et al., 2019), or learning an energy function given
observations of the dynamics in the context of density functional theory (Lin
et al., 2020), discussed in more detail in Sec. 4.5. The other methods are
more difficult to embed.

##### Speed.

As shown in Fig. 3.8, symbolic regression methods are best compared in
performance per compute-time, because search methods can, in principle, find
the global optimum given enough time (although this time might be longer than
the age of the universe), so just the “final” performance is difficult to
measure. SINDy is not in this comparison because it requires knowledge of the
base functions. However, it would be the fastest method, followed by NeSymReS,
DSR and classic GPs. EQL is likely the slowest method on small systems because
it requires a long training time. However, the time does not significantly
increase with system size and amount of data.

##### When to use which method?

For systems where we have a good idea about the occurring modules of the
functional form, FFX and SINDy are probably the method of choice. They are
simple and effective. For dynamical systems SINDy is the most specialized
method. The other SR methods are good if the functional building blocks are
unknown or nested structures are expected. Genetic Programming based methods
generally shine on small problem settings. Modern implementations can also fit
constants, but they result is sometimes complex nested structures. DSR and
NeSymReS are faster than standard SR methods and can yield potentially less
complex equations. However, less software is built for them, and they are less
easy to use. For high-dimensional data or when high fitting accuracy is
required, EQL might be the right choice. Also, when SR should be embedded,
only FFX, SINDy and EQL are practically usable.

###  3.2 Implicit equation discovery: dimensionality reduction and transfer
operators

This section reviews state-of-the-art methods for recovering implicit feature
representations of systems from data. We will review connections among methods
and emphasise the role and examples in the broad discipline of physics. Unlike
in the previous section, an explicit equation is not discovered but rather an
operator that encapsulates the system’s characteristics (typically spatio-
temporal dynamics). The field is tightly related to dimensionality reduction
and feature extraction in machine learning and signal processing (Arenas-
García et al., 2013), but also to transfer operators in functional analysis
(Khodkar et al., 2019).

####  3.2.1 Reduced-order models

In many domains, the goal is to study system dynamics from model simulations.
In this case, equations are encapsulated in the model itself, but large-scale,
high-fidelity nonlinear models can be challenging to simulate and require
significant computational power. In such cases, reduced order models (ROMs)
can simplify analysis and control design by trading off model accuracy for
computational complexity reduction. ROM can combine complex component-level
simulation models into system-level simulations used for control analysis and
design.

Two main classes of techniques for building ROMs: model-based and data-driven.
Model-based methods rely on a mathematical or physical understanding of the
underlying model and are designed for specific PDE-based models. In contrast,
data-driven methods use input-output data from the original high-fidelity
first-principles model to construct a ROM that accurately represents the
underlying system. Data-driven ROMs can be either static or dynamic models.
Static ROMs can be developed using techniques such as curve fitting and lookup
tables (LUT), while dynamic ROMs can be developed using deep learning
techniques such as LSTM, feedforward neural nets, and neural ODEs. The
obtained ROM ideally contains the essential physical mechanisms of the
original system while exhibiting simpler dynamics that can enhance
interpretability. In addition to the simpler physics, the ROM would be much
cheaper from the computational point of view than the numerical integration of
the governing equations, which can help obtain more efficient control and
optimisation techniques.

Developing a ROM typically requires two steps. The first step is to find a set
of coordinates where the original dynamics can be expressed in a compact form.
This is typically done in terms of modes, associated with coherent features in
the original system (Taira et al., 2017; Rowley and Dawson, 2017). In the
second step, it is necessary to find a set of differential equations governing
the temporal evolution of the amplitudes of the aforementioned modes. These
equations, which constitute a dynamical system, enable shedding light on the
physics of the (reduced) phenomenon under study. A widely-used approach to
perform the modal decomposition is the so-called proper-orthogonal
decomposition (POD) (Lumley, 1967), which is also known as principal component
analysis (PCA) in statistics and empirical orthogonal function (EOF) analysis
in meteorology (Arenas-García et al., 2013), and is closely connected with the
singular-value decomposition (SVD) (Brunton and Kutz, 2022).

Proper-orthogonal decomposition (POD), aka PCA or EOF  (Lumley, 1967; Arenas-
García et al., 2013) decomposes a dataset or a high-dimensional (spatio-
temporal) field into a set of orthogonal basis functions called modes or
eigenfunctions, which capture the dominant features of the data. The first few
modes explain most of the variability in the data, while the later modes
explain smaller and smaller amounts of variability. By truncating the number
of modes, one can obtain a low-dimensional representation of the data that
preserves the essential features of the original system. Given a set of data
snapshots
$\mathbf{X}=[\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{n}]\in\mathbb{R}^{m\times
n}$, where $m$ is the dimension of the data and $n$ is the number of
snapshots, we seek to decompose the data into a set of $r$ orthogonal modes
$\\{{\mathbf{u}_{i}}\\}_{i=1}^{r}$, such that
$\mathbf{X}\approx\mathbf{U}\bm{\Sigma}\mathbf{V}^{\top}$, where
$\mathbf{U}=[\mathbf{u}_{1},\mathbf{u}_{2},\ldots,\mathbf{u}_{r}]\in\mathbb{R}^{m\times
r}$, $\bm{\Sigma}\in\mathbb{R}^{r\times r}$ is a diagonal matrix containing
the singular values, and $\mathbf{V}^{\top}\in\mathbb{R}^{r\times n}$ is the
matrix of temporal coefficients. The modes $\\{\mathbf{u}_{i}\\}_{i=1}^{r}$
can be computed by performing a singular value decomposition (SVD) of the data
matrix
$\mathbf{X}=\mathbf{U}\bm{\Sigma}\mathbf{V}^{\top}=\sum_{i=1}^{r}\sigma_{i}\mathbf{u}_{i}\mathbf{v}_{i}^{\top}$
where $\sigma_{i}$ is the $i$-th singular value, and $\mathbf{v}_{i}$ is the
$i$-th right singular vector. The modes are then given by
$\mathbf{u}_{i}=\frac{1}{\sqrt{\sigma_{i}}}\mathbf{X}\mathbf{v}_{i}$.

The POD/PCA framework enables decomposing spatio-temporal data (e.g. flow
velocities, weather or climate variables which depend on the spatial
coordinates and time) into a set of spatial modes (which only depend on the
spatial coordinates), multiplied by their temporal coefficients (which define
their change of amplitude with time). POD/PCA ensures that components are
orthogonal and optimality with respect to the variance explained by a reduced
number of modes, cf. Fig. 3.9 for an example of PCA on a toy spatio-temporal
data flow. Other alternative multivariate methods, like partial least squares
(PLS) or canonical correlation analysis (CCA) seek projections that maximise
covariance or correlation, respectively (Arenas-García et al., 2013). Still,
all these projection methods are linear and thus cannot cope with nonlinear
spatio-temporal feature relations and complex dynamics. This can be addressed
with kernel machines (Arenas-García et al., 2013). Oblique and nonlinear
transformations can also be learned by embedding Varimax in Reproducing Kernel
Hilbert Spaces (RKHS) explicitly (Bueso et al., 2020; Diego Bueso, 2020).
Other ways to obtain non-linear transformations from observation space to ROM
space, e.g. using neural networks, will be discussed below.

####  3.2.2 Transfer operators for learning nonlinear dynamics

Transfer operators are related to the abovementioned methods and allow the
characterisation and modelling of complex dynamic systems. These operators’
eigenfunctions can decompose a system given by an ergodic Markov process into
fast and slow dynamics and identify modes of the stationary measure called
metastable sets.

The Koopman operator is a linear operator that describes the dynamics of a
system by lifting the state variables into an infinite-dimensional Hilbert
space. Thus, it enables us to effectively linearise complex temporal
trajectories and hence is a compelling approach in dynamical systems research
(Khodkar et al., 2019). Its application is expanding in both theoretical
(Khodkar et al., 2019; Brunton et al., 2017), and practical domains from
molecular dynamics and fluid dynamics, atmospheric sciences, and control
theory (Schmid, 2010; Brunton et al., 2016a; Klus et al., 2018).

The advantages of the Koopman operator are numerous. First, it is a powerful
tool for analysing and predicting the behaviour of a system over time. By
lifting the state variables into a higher dimensional space, the Koopman
operator can identify patterns in a system’s behaviour that may otherwise be
difficult to detect. This can be especially useful for uncovering hidden
dynamical structures in chaotic systems.

Second, the Koopman operator allows us to develop data-driven models of
dynamical systems. Using the operator’s eigenfunctions as basis functions, it
is possible to develop models of dynamical systems operating on these
summarised coordinates, without solving or even understanding the underlying
equations of motion. This makes the Koopman operator an attractive tool for
model-based control and optimisation. Third, the Koopman operator is useful to
discern key properties of highly nonlinear dynamical systems. In short, with
the expansion of original state variables to infinite dimensions, the operator
can uncover subtle nonlinear behaviour in a system that would otherwise be
difficult or impossible to detect (Khodkar et al., 2019).

Koopman operator (Khodkar et al., 2019) The Koopman operator is a linear
operator that describes the evolution of an observable function of a dynamical
system. Let $\mathcal{M}$ be a manifold of dimension $n$ and
$f:\mathcal{M}\rightarrow\mathbb{R}$ be a real-valued function. The Koopman
operator, denoted by $\mathcal{K}$, is defined as an infinite-dimensional
linear operator that acts on the space of observable functions $f$ such that
for any $f\in L^{2}(\mathcal{M})$, ${\mathcal{K}}f({\bf x})=f(T({\bf x})),$
where $T:\mathcal{M}\rightarrow\mathcal{M}$ is the evolution operator that
maps each point ${\bf x}\in\mathcal{M}$ to its next iterate in time. Now, let
$f:\mathcal{M}\rightarrow\mathbb{R}$ be a bounded, measurable observable and
$\mathcal{K}$ be the Koopman operator associated with a dynamical system.
Then, there exists a sequence of eigenfunctions
$\psi_{j}:\mathcal{M}\rightarrow\mathbb{C}$ and a corresponding sequence of
eigenvalues $\lambda_{j}\in\mathbb{C}$ such that ${\mathcal{K}}\psi_{j}({\bf
x})=\lambda_{j}\psi_{j}({\bf x}).$ The Koopman operator preserves the linear
structure of the space of observables, provides a linear representation of
nonlinear dynamics, and its eigenfunctions provide a useful basis for
approximating dynamical systems.

There are, however, some drawbacks associated with the Koopman operator.
First, it is intrinsically an infinite-dimensional operator, and although
there are efficient finite-dimensional approximations available, the accurate
computation of their eigenvalues and eigenfunctions can be computationally
expensive (Kaiser et al., 2020). This can be a problem for real-time
applications, such as model-based control (Kaiser et al., 2021a). Second,
related to this drawback, the operator’s eigenfunctions are often difficult to
interpret, hampering the capacity of the learned representations to explain
the underlying system’s dynamics. Finally, the Koopman operator assumes that a
locally-linear behaviour can represent nonlinearities sufficiently accurately.
Thus, the Koopman operator is an important theoretical and applied research
tool for understanding and predicting the behaviour of complex dynamical
systems. Its data-driven approach to model-based control and optimisation has
opened up new possibilities for real-time applications (Kaiser et al., 2021a).
Moreover, its ability to uncover subtle nonlinear behaviour in chaotic systems
has made it invaluable for studying chaotic dynamical systems (Takeishi et
al., 2017).

More specifically, given a space in which the dynamics is linear, a successful
approach to approximate transfer operators (such as the Koopman operator) from
the data is called dynamic-mode decomposition (DMD) (Schmid, 2010). It is also
possible to obtain ROMs using DMD (Schmid, 2010), which is also based on
concepts from linear algebra and assumes that the system’s state can be
advanced in time via a linear operator $\mathbf{A}$. While the POD modes are
orthogonal in space, the DMD ones are orthogonal in time, and each mode is
associated with a particular frequency and a growth rate. Therefore, DMD may
help to identify temporal patterns in the data more clearly than POD. In
contrast, POD may lead to a more compact low-order representation of the
original system due to its optimality property. See an illustrative example in
Fig. 3.9.

Dynamic-mode decomposition (DMD) (Schmid, 2010) is a technique used to
approximate the normal modes and eigenvalues of a linear system. Additionally,
these modes can be associated with a damped or driven sinusoidal behaviour in
time. DMD is useful for identifying a system’s frequency and decay/growth
rate. Let us define a dynamical process formulated as
$\frac{d{\mathbf{x}}}{dt}=f({\mathbf{x}},t,\mu),$ where ${\mathbf{x}}$ defines
a measurement, $t$ is a time, $\mu$ is a parametric dependence, and $f$
indicates an unspecified system but from which we obtain many data. Therefore,
the complex dynamical system $f$ can be approximated as follows
$\frac{d{\mathbf{x}}}{dt}\approx{\bf A}{\mathbf{x}}$, where
$x\in\mathbb{R}^{n}$, $n\gg 1$ and ${\bf A}$ defines a linear dynamical
system. Then its general solution is the ‘exponential solution’ defined as
${\mathbf{x}}={\mathbf{v}}e^{\lambda t}$, where ${\mathbf{v}}$ and $\lambda$
are eigenvectors and eigenvalues of the linear system ${\bf A}$. The problem
of finding the eigenvectors ${\mathbf{v}}$ and the eigenvalues $\lambda$ is a
eigenvalue problem defined as $\lambda{\mathbf{v}}={\bf A}{\mathbf{v}}.$ Yet,
we are interested in obtaining ${\bf A}$, not its eigendecomposition. This is
what the so-called ‘exact DMD’ does. DMD uses observations/measurements
$x_{j}={\mathbf{x}}(t_{j})$, defined at a time point $j$ to construct two
matrices: the first concatenating the data from the first snapshot to
$(m-1)$-th snapshot, and the second with the shifted-by-1-time-step samples,
${\bf X}$ and ${\bf Y}$, respectively. The goal is thus building a linear
dynamical system $A$ fitted with $\frac{d{\mathbf{x}}}{dt}={\bf
A}{\mathbf{x}}$, and thus learn the linear dynamical system ${\bf A}$ that
takes the data ${\mathbf{x}}$ from current state $(j-1)$ to future state
$(j)$, that is ${\bf Y}={\bf A}{\bf X}$. The linear dynamical system $A$ can
be extracted using a pseudo-inverse ${\bf X}^{\dagger}$ of ${\bf X}$, that is
${\bf A}={\bf Y}{\bf X}^{\dagger}.$ Intuitively, the linear dynamical system
${\bf A}$ performs a least-square fitting from the current state ${\bf X}$ to
the future state ${\bf Y}$.

Over the last decades, different numerical methods have been introduced:
Ulam’s method (Ulam, 1960), extended dynamic-mode decomposition (EDMD)
(Williams et al., 2015a, b; Klus et al., 2016), and the variational approach
of conformation dynamics (VAC) (Noé and Nüske, 2013; Nüske et al., 2014). The
advantage of purely data-driven methods is that they can be applied to
simulation and observational data. Hence, information about the underlying
system itself is not required. An overview and comparison of such methods can
be found in (Klus et al., 2018). Applications and variants of these methods
are also described in (Rowley et al., 2009b; Tu et al., 2014; McGibbon and
Pande, 2015), while kernel-based reformulations of the methods above have been
proposed before in (Williams et al., 2015b; Schwantes and Pande, 2015). Note
that the framework of higher-order dynamic-mode decomposition (HODMD) (Le
Clainche and Vega, 2017) enables relaxing the linear assumption by including
several temporal snapshots to build the operator by exploiting Takens’ delay-
embedding theorem (Takens, 1981). The HODMD approach requires additional
hyper-parameter tuning, but it has led to very insightful results, for
instance, in the context of complex turbulent flows, where this method has
enabled identifying the coherent structures responsible for the concentration
of pollutants in cities (Lazpita et al., 2022). Another relevant application
of HODMD includes cardiovascular flows (Groun et al., 2022).

(a) Space | (b) Time | (c) Data | (d) PCA | (e) DMD
---|---|---|---|---
$x_{1}(s,\cdot)$ | $x_{1}(\cdot,t)$ | $x(s,t=10)$ | Mode 1 | Mode 1
|  |  |  |
$x_{2}(s,\cdot)$ | $x_{2}(\cdot,t)$ |  | Mode 2 | Mode 2
|  |  |  |
Figure 3.9: Comparison between DMD and PCA with synthetic spatio-temporal
data. The signal under analysis $x(s,t)$ (c) is the sum of two generative
signals (a,b): $x_{1}(s)$ is a Gaussian that decays exponentially, and
$x_{2}(s,t)$ is a square that oscillates at a lower frequency. The projections
onto the two top components of PCA (d) and DMD (e) show that DMD extracts
cleaner spatial coherence patterns from the data.

####  3.2.3 Dynamic modes in neural-network latent spaces

Despite the interesting properties of POD and DMD, their inherent linearity
typically leads to the requirement of very large numbers of modes to
reconstruct most of the variance of the original signal, for example, in
three-dimensional turbulent flows (Baars and Tinney, 2014). Neural networks,
especially autoencoders (AEs), have been proposed to obtain a reduced-order
nonlinear representation of the original data. AEs exploit non-linear
activation functions to produce significantly more compact representations in
the latent space than those with, e.g. POD (Hinton and Salakhutdinov, 2006).
Figure 3.10[a-b] shows the use of AEs for learning (dynamic) feature
representations.

AEs have been used in fluid mechanics to obtain compact modal decompositions
of the flow around a two-dimensional cylinder (Murata et al., 2020) and in
more complex turbulent flows, e.g. the flow in a simplified urban environment
(Eivazi et al., 2022). Interestingly, when restricting neural networks to
linear activation functions, one recovers the POD modes, as shown by Milano
and Koumoutsakos (2002) with a multilayer perceptron (MLP) in turbulent
channel flow. Shallow NNs have been used for flow reconstruction, in this case
from sparse measurements, as illustrated by Erichson et al. (2020) for several
flow cases. A more general illustration of the potential of AEs based on
convolutional neural networks (CNNs) was presented by Lee and Carlberg (2020),
and an application to spectral submanifolds was developed by Cenedese et al.
(2022). The reader is referred to Refs. (Benner et al., 2015; Carlberg et al.,
2017) for a survey of classical methods applicable to linear subspaces.

Despite the superior compression performance of AEs compared with POD, the
former does not have two very interesting properties of the latter, namely the
optimality and orthogonality of the resulting modes. These are important
properties due to their connection with interpretable and parsimonious ROMs.
Regarding optimality, Fukami et al. (2020) proposed an interesting approach
based on hierarchical autoencoders (HAEs). They first trained a CNN-based AE
fixing the dimension of the latent space to just one, obtaining one latent
vector. Then, they trained another CNN-AE with a latent dimension of two and
fixed the first latent variable to the one obtained in the previous NN, thus
obtaining a second latent vector. Through this recursive strategy, they
obtained a sequence of latent vectors exhibiting progressively less
contribution to the reconstruction of the original signal, allowing them to
establish a ranking in the resulting modes. This approach was tested in the
flow around a two-dimensional cylinder, although it is important to note that
the resulting modes were not orthogonal. This was addressed by Eivazi et al.
(2022), who used $\beta$-variational autoencoders ($\beta$-VAEs), which enable
introducing stochasticity in the latent space to impose orthogonality in the
resulting AE modes, a phenomenon that was explained in Rolinek et al. (2019)
among a connection of $\beta$-VAEs to PCA. Also in the case of the
$\beta$-VAEs, the modes can also be ranked in terms of their contribution to
the reconstruction.

Figure 3.10: Approximating Koopman operator with explicit nonlinear mappings
using autoencoders. (a) An autoencoder neural network learns a mapping
$\bm{\Phi}$ compressed latent representation ${\bf y}$ from input data ${\bf
x}$ by minimising the reconstruction error $L$. (b) One can incorporate the
Koopman linear operator ${\bf K}$ operating in the latent representation ${\bf
y}_{t}$, which can then be used for prediction ${\bf x}_{t+1}$ from the
transformed ${\bf y}_{t+1}$ (Brunton et al., 2017). (c) Those representations
are not necessarily physically consistent. This can be addressed by enforcing
equation dictionaries using SINDy in the loss for the simultaneous discovery
of coordinates and parsimonious dynamics (Champion et al., 2019). The loss now
accounts for the reconstruction of the input data, as in a regular
autoencoder, and the temporal dynamics (gradients $\nabla_{x},\nabla_{y}$) of
${\bf x}$ and ${\bf y}$ projected onto SINDy bases.

####  3.2.4 Equation discovery in latent representations

Perhaps the biggest challenge in data-driven model discovery is balancing
model efficiency with descriptive capabilities. Parsimonious models with the
fewest terms required to capture essential interactions promote
interpretability and generalizability. However, obtaining parsimonious models
is linked to the coordinate system in which the dynamics are measured. The
previous methods based on dimensionality reduction, e.g. ROM, DMD or AE,
extract expressive components without simultaneously discovering coordinates.
In (Champion et al., 2019), AEs were trained for data reconstruction and to
recover a parsimonious dynamical system model through sparse regression using
SINDy (see Sec. 3.1.2). See Fig. 3.10[c]. The joint goal of discovering models
and coordinates is critical for understanding many modern systems. Using SINDy
as an explicit equation discovery regulariser in the latent space balances
simplifying coordinate transformations and nonlinear dynamics to identify
coordinate transformations where only a few nonlinear terms are present.

####  3.2.5 Discovering fundamental variables

Despite advances in equation discovery (either through implicit or explicit
representations), the main core problem is identifying state variables. The
discovery typically refers to the identification of the governing equations,
not the identification of the physical forces or variables. The vast majority
of data-driven models of discovery rely, however, on pre-existing knowledge of
the state variables, e.g., the position and velocity of a rigid body object.
This relies on deep domain knowledge and strong assumptions. In addition, such
assumptions cannot work properly for new physical systems or when those state
variables cannot be measured. The work (Chen et al., 2022) proposed a
principle for determining the number and identity of state variables in a
system from high-dimensional data and demonstrated high effectiveness using
video recordings of physical systems. The algorithm discovered the intrinsic
dimension of the observed dynamics and could identify candidate sets of state
variables without prior knowledge of the underlying physics. Alternatively,
other studies sought to identify the fundamental state variables via manifold
learning in ambient RKHS (termed Diffusion maps, e.g., (Kemeth et al., 2022;
Thiem et al., 2020), see also Section 4.1.2).

In short, the field of variable discovery is filled with many opportunities in
the physical, biological and chemical sciences (Pukrittayakamee et al., 2009;
Kemeth et al., 2022; Schütt et al., 2017), as well as many challenges (Chen et
al., 2022). Finally, and interestingly, we want to emphasise that variable
discovery is intimately related to revealing latent confounders in the field
of causal inference (Stegle et al., 2010; Monti et al., 2020; Diaz et al.,
2023).

###  3.3 Perspectives

Let us indulge ourselves with a brief overview of the main challenges (both
conceptual and technical) and the opportunities for future research in the
field of equation discovery for the physical sciences.

####  3.3.1 Challenges

The field of equation discovery from data is very prolific and is situated at
the intersection of many communities: statistics, machine learning,
computational fluid dynamics, Bayesian inference, dynamical systems and
control theory, functional analysis and causal inference. The field has
occupied scientists for centuries at all levels. The quest for optimal and
automated solutions has traditionally considered moving in one or several
subspaces in the sparsity-extrapolation-generalisation space, i.e., models
should be simple, generalisable/robust, and capable of extrapolating outside
the sample space (Fig. 3.11). These are very ambitious goals, implying both
theoretical and practical challenges.

Figure 3.11: On the quest for the optimal model in the sparsity-extrapolation-
generalisation space.

##### Theoretical challenges.

From a more theoretical perspective, both the identifiability of the system’s
equations (Fajardo-Fontiveros et al., 2023; Antonelli et al., 2022a) and the
role (or preference) for sparsity (simplicity) have been questioned (Fuentes
et al., 2021; Guimerà et al., 2020). In addition, there is a long-standing
debate on the evaluation of the obtained solution, where many criteria can be
adopted. Is it only about invariances and robustness in space and time? Is
Nature always simple and compositional, such that compactness and sparsity
rule in natural systems? The issue of model’s (i.e. hypotheses)
intercomparison and evaluation also speaks to the more elusive question of how
to reconcile solutions offered by different equation discovery methods.

Another important theoretical challenge is related to the fact that, very
often, one (1) assumes that all involved state variables are given/observed,
which resembles the sufficiency assumption in causal inference, and (2)
selects a subset of representative states, assumes a particular basis to
express the solution, or can operate on a manifold subspace (Champion et al.,
2019; Brunton and Kutz, 2022; Daniels and Nemenman, 2015; Waltz and Buchanan,
2009). Both are strong assumptions that challenge the process of discovering
equations and raise many questions. How do we choose the right variables to
include in the equation discovery method? How much does the solution change
when a variable is omitted or added? Here, identifiability issues and latent
confounders play a substantial role. And foremost, what if we cannot measure
the underlying state variables? Can we identify them automatically? Several
methods have arguably proposed to discover the latent variables (Chen et al.,
2022), and other efforts exploit the link between RKHS techniques and a
plausible master equation underlying the spatiotemporal evolution of the data
probability density function (Coifman and Lafon, 2006) to automatically learn
a set of fundamental coordinates of the system, cf. Sec. 3.2.5. The inferred
subspace has shown effectiveness in identifying the dynamics of (partial)
differential equations (Kemeth et al., 2022; Thiem et al., 2020), see examples
in Sec. 4.1.2. However, like other explicit methods discussed in this section,
it requires prior guidelines on the differential equation to model and
incorporates a range of heuristics during its processing pipeline (Kemeth et
al., 2022).

##### Practical challenges.

Important practical challenges are related to the model development and data
characteristics: (1) high dimensionality, (2) nonlinear relationships and (3)
risk of overfitting. In many cases, the number of variables and parameters can
be very large, making it difficult to find the most relevant features and
relationships. This is the scenario where non-identifiability arises. Another
challenge is that real-world systems often exhibit highly nonlinear
relationships between input and output variables. This makes finding a good
representation of the underlying dynamics difficult, and the search space for
equations can become very large. While several nonlinear methods that capture
complex dynamics exist (kernels, neural networks), performance evaluation and
hyperparameter tuning are still important challenges. Finally, symbolic
regression approaches can also suffer from overfitting, where the model fits
the training data well but performs poorly on new data, which speaks to the
trade-off between the model’s accuracy and complexity.

When working with spatiotemporal data, it is relatively unclear how to
incorporate information about time-lagged relations and interventional data.
Both the explicit and the implicit approximations show particular challenges,
though: On the one hand, symbolic regression techniques (such as SINDy) face
significant problems in defining the basis functions, working in high-
dimensional problems, and the impact of (even a limited) amount of noise.
Other methods, like those based on AI Feynman, even if they incorporate
sensible criteria to guide the equation discovery (like compositionally,
reversibility or physical unit consistency), almost predict completely
different equations when changing constants in the true equation (Suseela et
al., 2022). On the other hand, implicit methods (like DMD or Koopman
operators) do not provide an explicit equation but a latent feature
representation to explain system dynamics. These methods struggle with
nonstationarities, nonlinearities, and gaps noise, which remain unresolved
problems in the literature (Wu et al., 2021). Note that similar challenges to
those in causal discovery remain here, cf. Sec. 2.2. Besides, DMD methods fail
to generalise outside the training data and violate basic physical laws. To
alleviate this, integrating domain knowledge (such as symmetries, invariances
and conservation laws) in DMD has been recently introduced as an effective,
robust approach (Baddoo et al., 2023).

####  3.3.2 Opportunities

While symbolic regression presents challenges, it offers exciting research
opportunities for the physical sciences. Three main opportunities can be
identified: (1) model interpretability, (2) model compression and evaluation,
and (3) model selection. Equation discovery is a step forward in the system’s
understanding. The field leverages fully interpretable models, and unlike
causal models, equation discovery (symbolic regression) models are directly
applicable predictive models. The discovered equations from data can provide
insights into the underlying dynamics of a system, thus helping researchers
better understand how different factors interact and contribute to a
particular outcome. Even with implicit latent representations,
interpretability can be accomplished with interventional analysis. Another
interesting opportunity is model compression, as symbolic regression models
can also be used to compress large datasets into simple equations that capture
the system’s essential features; this can make it easier to analyse and
visualise the data and make it more computationally efficient to work with.
Another practical opportunity of symbolic regression models is that they offer
a reduced set of possible solutions typically ranked in amenable Pareto
fronts, which of course, trigger difficulties in choosing the right model but
also fruitful scientific discussions about the plausibility of identified
relations. This can save researchers time and effort and help identify
unexpected patterns and relationships in the data.

##### Model interpretability and intervention analysis.

The discovered explicit models are interpretable in nature. However, when
complexity cannot be traded for accuracy or whenever an implicit feature
representation is learned, intervention (or sensitivity) analysis offers
opportunities for interpretability. For example, one can (1) ignore or
simplify the problem by performing small perturbations away from real-world
dynamics, which might help identify the proper relationship between variables;
(2) intervene on exogenous variables (e.g. wind or solar irradiation, mixing
coefficients, initial conditions in climate sciences, or targeted, direct
brain stimulation in neurosciences) which is equivalent to collecting more
data; (3) create a library of trajectories under different conditions and
select those which match the desired intervention; and (4) intervene in the
learned latent space and decode the intervention back to input space, thus
allowing us to generate interventions that follow the system’s natural
trajectories. Further analysis methods, e.g. for studying interventions of the
learned ODE, might be fruitful, as they can access long-term dynamical
properties, such as how distortion in the eigenvalues affects the system’s
stability as the phase space changes.

##### Model compression and evaluation.

Scientists frequently use metrics to evaluate new ideas or distinguish between
competing hypotheses. As we have seen before, a governing equation should be
simple but not simpler, accurate for prediction, robust under distortions and
changes, and invariant in space and time. Equation discovery offers a direct
way to learn plausible models and an indirect way to contrast and evaluate
derived models. For this, one typically assesses (1) the predictive accuracy
when answering how well the (simplest) hypothesis explains the data; (2) the
model’s invariance under distributional shift to account for the causal
mechanisms; and (3) the robustness under interventions to study how the
proposed process (or descriptive equation –representation) is consistent with
interventions on the model dynamics, such as deactivation of components or
targeted modification of exogenous variables. The latter is the rarest form of
validation due to its high computational cost and difficulty in experimental
design.

##### Model selection.

Enforcing sparsity in model selection can lead to unrealistically too simple
models. That is why methods that can provide solutions along the Pareto line
(Fig. 3.2) are needed to capture complex relationships and offer subsets of
plausible model solutions. Alternative regularisation schemes will likely be
important alongside profound estimates of uncertainty and extrapolation
indicators.

##  Chapter 4 Case studies in the physical sciences

This section gives concrete examples of applying different data-driven causal
and equation discovery in important fields of the physical sciences:
neuroscience, Earth and climate sciences, and fluid and mechanical dynamics,
cf. Table 4.1.

Table 4.1: Case studies presented and the main methods used in this section. | Neuroscience | Earth & climate | Fluid dynamics
---|---|---|---
Causal discovery |  Causal connectivity (DCM, GC, TE, SCM) |  Carbon-water interactions Climate model comparison (CCM, PCMCI) | –
Equation discovery |  Learning trajectories (kFDA, GP, Variational Bayes RNN, Diffusion Maps) |  Ocean Mesoscale closures (RVM, DMD, SINDy) |  Turbulence understanding Vortex shedding (SINDy, Genetic Programming)

###  4.1 Neuroscientific applications of physics-based machine learning

####  4.1.1 Overview of parsimonious models for neural population dynamics

Neuroscientific modelling falls within the remit of the field known as
computational or theoretical neuroscience, which studies the transmission of
information in the nervous system at multiple spatiotemporal scales (ranging
from neuronal to whole-brain levels) in relation to perception, cognition, and
behaviour (e.g., (Gerstner et al., 2014; Koch, 2004)). Thus, an ongoing
challenge in computational neuroscience is to link biophysically detailed
models operating at microscopic levels with meso/macroscopic theories of
cortical processing (Rabinovich and Varona, 2018). This enterprise is often
addressed by deducting low-dimensional systems of partial differential
equations or maps describing _coarse-grained_ variables derived from
collective neural responses. Such different neurobiologically plausible
simplifications are commonly termed ensemble, population, neural-mass, or
simply _firing-rate_ models (see, for instance, (Byrne et al., 2021; Brunel,
2000)).

These synthesis efforts are intimately connected with empirically discovering
a reduced dynamical system generating the observed neuronal activity. However,
neurocomputational modelling traditionally focused on analytical, deductive
approaches mapping realistic cortical networks to _tissue-level_ descriptions,
as opposed to data-driven model discovery, reviewed in Sec. 3. Thus, ensemble
models are typically principles-based, often hinged on assumptions about
dynamical interactions arising within homogeneous pools of neurons (e.g.,
(Amari, 1977; Byrne et al., 2020; Tabas et al., 2019; Wilson and Cowan,
2021)).

Early neural ensemble models stemmed from applying statistical mechanical
principles to the interaction of homogeneous pools of (excitatory and
inhibitory) populations (Wilson and Cowan, 1972; Potthast, 2013). Later,
physics formalisms like the Fokker-Planck approach for describing the
spatiotemporal evolution of the probability distribution of neuronal activity
enabled theoretical neuroscientists to take a more holistic approach to
identify mean-field approximations of networks of spiking neurons (e.g.,
(Gerstner et al., 2014; Mattia and Del Giudice, 2002; Brunel, 2000)). These
and other nonlinear dynamical systems tools (Rabinovich et al., 2008) provided
closed-from, exact solutions for the collective behaviour of neural
populations (Montbrió et al., 2015; Mattia et al., 2019; Byrne et al., 2020)
capable of an extensive dynamical repertoire, although strongly dependent on
universal theoretical assumptions, given their deductive nature (see Sec. 1).
Alternatively, a Laplacian assumption on this probability distribution
resulted in neural mass descriptions (Marreiros et al., 2010), recently
proposed as building blocks for whole-brain models with translational
applications (Schirner et al., 2022).

Overall, these chiefly deductive approaches rendered compact models of
differential equations based on _a priori_ assumptions about neural and
synaptic variables, fostering the interpretability of high-complex neuronal
networks. By contrast, inferential approaches in neuroscience have been
typically utilised to empirically identify neural dynamics underlying
cognition and behaviour, as discussed next.

####  4.1.2 Empirical reconstruction of neuronal trajectories

There is an increasing focus on applying classic and deep machine learning
approaches to reconstruct attracting and transient dynamics of cerebral cortex
responses (Duncker and Sahani, 2021). A neural trajectory $T\;(n\times d)$ is
often defined as the sequence of $n$ neural response vectors $\bm{x}(t)$
embedded in a $d$-dimensional state-space (the _ambient_ space), spanned by
neural ensemble activity or proxies thereof (e.g., firing-rates,
electromagnetic potentials), their lags and nonlinear transformations (Galgali
et al., 2023; Duncker and Sahani, 2021; Balaguer-Ballester et al., 2011).

Traditionally, standard dimensionality reduction techniques (e.g., PCA, Multi-
dimensional Scaling, Discriminant analysis etc., see Sec. 3.2) were directly
applied for the visualisation of high-dimensional neural trajectories,
showcasing coarse-grained aspects of firing-rate dynamics concerning, for
example, cognitive decisions or motor functions (Cunningham and Yu, 2014;
Hyman et al., 2012). More recently, Gaussian processes provided a flexible
approach to derive a low-dimensional manifold representing the dynamical
systems generating the observed activity. They just require a reasonable
hypothesis on temporal correlations between observations (the prior covariance
function (Rasmussen and Williams, 2006)). For instance, Gaussian process-based
factor analysis (GPFA) (Yu et al., 2009), and other latent-variable methods
(Aoi et al., 2020), provide such low-dimensional subspace while simultaneously
approximating the probability of spiking -without the compelling need for
probability density estimation (Yu et al., 2009; Gokcen et al., 2022). These
and related approaches can identify latent neural trajectory manifolds in
prefrontal and motor cortices underlying decision-making (Rutten et al., 2020;
Aoi et al., 2020; Duncker and Sahani, 2018). Specifically, recent GPFA
variants were able to discern between competing models for the contribution of
upstream areas to recurrent dynamics supporting decision-making in the monkey
prefrontal cortex (Galgali et al., 2023).

In general, new developments in deep auto-encoders such as Latent Factor
Analysis via Dynamical Systems (LFADS) enabled the successful reconstruction
of single-trial spiking activity (Pandarinath et al., 2018). Moreover, related
approaches such as preferential subspace identification (PSI) incorporate
behavioural labels to inform the dimensionality reduction algorithm;
effectively discerning _behaviourally relevant_ dynamics from the general
neuro-dynamical landscape Sani et al. (2021). Remarkably, recent self-
supervised methods, also guided by behavioural observations, leverage
contrastive learning and nonlinear ICA to produce consistent sub-spaces across
experimental sessions and subjects Schneider et al. (2023). In addition, like
previous approaches, they can operate in a supervised fashion, fostering
decoding with respect to competing methods in a range of calcium and
electrophysiological recordings in different species Schneider et al. (2023).

Covariance (kernel) function methods can also recreate salient facets of
cortical dynamics like attracting sets. To this end, they leverage delay-
embedding techniques in RKHS spanned by neuronal correlations and their
temporal structure (Balaguer-Ballester et al., 2011, 2014; Lapish et al.,
2015; Balaguer-Ballester et al., 2020) for identifying compact manifolds
mapping animal’s choice with attracting sets of ensemble trajectories (Lapish
et al., 2015; Balaguer-Ballester et al., 2020). Figure 4.1 presents an
illustrative example of these RKHS techniques for recreating neuronal
trajectories underlying the effect of dopamine at the circuit level. This
approach facilitated the evaluation of mechanistic theories of dopamine
modulation during decision-making, which was challenging given the limitations
of direct experimental manipulations (Lapish et al., 2015). The figure shows
the flow field of trajectories derived from the activity of neuronal
constellations in the rodent anterior cingulate cortex. Interestingly, the
dynamic landscape depicted during working memory tasks in an optimal RHKS can
be approximately described as transients connecting multiple attracting sets
mapping spatial choices (Fig. 4.1a). This robust multi-stable scenario is
completely disrupted by high doses of amphetamine (a well-known trigger of
dopamine release, Figure 4.1b), while it is enhanced by low doses (see (Lapish
et al., 2015)), in line with long-standing theoretical predictions of
biophysical models (Durstewitz and Seamans, 2008).

Figure 4.1: Example of the reconstruction of neural trajectories in an RKHS
derived from rodent anterior cingulate cortex (ACC) multi-array recordings,
taken from Lapish et al. (2015). a. The flow field stems from projecting an
ACC ensemble firing rates onto the three main coordinates of a discriminant
subspace (DC1-DC3, computed here by kernel-fisher discriminant analysis
orthogonalised for a faithful representation, details in, e.g., (Balaguer-
Ballester et al., 2011, 2020, 2014; Lapish et al., 2015)). The DC1-DC3 is
embedded into an optimal ambient high-dimensional RKHS (top schematics),
spanned by neuronal ensemble firing rate and its higher-order correlations (up
to 3rd order in this example) operating on a delay-coordinate map. The colour
code (bottom) corresponds to rat choices in this experiment (the schematic of
the radial-arm maze used in this experiment is taken from (Balaguer-Ballester
et al., 2011, 2014)), which occupy distinct regions of the subspace. The flow
field indicates faster shifts (large vector lengths) during the transition
points, while it slows downs nearer the centroids of the clusters, suggesting
an attracting-like dynamic landscape. The inset quantifies this uneven
distribution of flow field speeds as a function of the distance to the
centroid, further supporting this observation on flow convergence. b. This
ordered phase-space structure cannot be achieved in the original feature space
or with delay-coordinate maps. c. It is also destroyed when the animal
receives a high dose of amphetamine, even in an optimal ambient space. Figure
adapted from (Lapish et al., 2015) and from (Balaguer-Ballester et al., 2011)
with publishers’ permission (the Society for Neuroscience and the Public
Library of Science).

Further facets of brain dynamics, such as chaotic attractors, have been
recently addressed with recurrent, piecewise-linear architectures, amenable to
optimisation via back-propagation variants or Bayesian variational inference
(Durstewitz, 2017; Brenner et al., 2022). In these approaches, tractability is
promoted by leveraging units’ linearisation for approximating trajectory
inference in a parsimonious, transparent fashion. Empowered by these
characteristics, such linearised recurrent networks could infer pathological
whole-brain dynamics from functional magnetic resonance imaging (fMRI)
recordings (Koppe et al., 2019). Moreover, recent developments of these
methods embody biophysically-inspired computations such as dendritic
processing, fostering their reconstruction capabilities of nonlinear dynamical
systems (Brenner et al., 2022).