How much can a picture tell us? The image above shows an initial disturbance slowly breaking up into many small waves, which themselves merge and split as time progresses. Such complex behaviour is certainly beautiful to behold, but as scientists we wish take it a step further. Can we understand the underlying system and capture it in a (quantitative) model? This seems a daunting task considering figure 1.1, yet this process is at the heart of science. The classical approach to construct such models is known as First-Principle modelling. It starts, as the name implies, from basic princi-ples and relies on simple abstractions and laws such as conservation of mass and energy to derive a model. FP modelling is an iterative process, tweaking and improving the model until it describes the observed data sufficiently, making it very time-consuming. Additionally, it requires expert knowledge; without any knowledge about the underlying system, it is next to impossible to come up with or improve the model. Its biggest drawback however is its ability to only model fairly simple systems. Figure 1.1 is a one-dimensional solution of the so-called Kuramoto-Shivashinsky equation,

𝜕𝑢 𝜕𝑡 = 𝜕 2 𝑢 𝜕𝑥 2 - 𝜕 4 𝑢 𝜕𝑥 4 -𝑢 𝜕𝑢 𝜕𝑥 .
(1.1) While this is a fairly complex equation due to its fourth-order derivatives and non-linear terms, it does consist of only four terms. Faced with complex systems such as the airflow over a wing, the behaviour of crowd or the movement of impurities in a plasma, first-principle modelling becomes impossible; it would simply take too much time to derive a (effective) model describing all the interactions and processes. In these cases one resorts to data-driven modelling. In this paradigm we seek only to mathematically describe the data we observe. Instead of concerning itself with why the data looks like it does, it restricts itself to accurately modelling the data without knowledge of the underlying system. This does not make data-driven modelling more trivial than first-principle modelling: constructing a model which accurately describes the noisy, often high-dimensional data is no easy feat1 .

While data-driven modelling thus allows us to work with very complex systems, it lacks the key feature which makes first-principle modelling so powerful and popular: interpretability. This powerful property has various connotations across various fields, and we consider it here that a given dataset can be described by a fairly compact equation. Interpretable models in this form are implicitly tied to trust (can we prove mathematically a certain state can or cannot be reached?), abstraction (we can decompose complex behaviour as the interaction between subsystems), and exploitation (if we would suppress this interaction, we could use it), and thus forms a key property of science. The ideal approach would thus marry the interpretability of first-principle modelling to the ability of data-driven modelling to handle complex data. Simply put, we wish for an algorithm which, given figure 1.1, would be capable of autonomously discovering eq. (1.1)2 .

Model discovery

The process of using an algorithm to autonomously discover a systems governing equation is known as model discovery 3 . Before I introduce the typical approaches to model discovery, I want to discuss what these models typically look like. In the context of model discovery we typically recognize three categories. Given some spatio-temporal data 𝑢(𝑥, 𝑡),

• explicit models give an explicit solution for the data we're measuring, 𝑢(𝑥, 𝑡) = 𝑓(𝑥, 𝑡). Having an expression for the observable makes these among the most powerful, but in unfortunately these models can usually only be derived for a very restricted subset of systems. • A much more common approach is to derive a model in terms of its temporal evolution, i.e. in the form of a differential equation 𝜕𝑢/𝜕𝑡 = 𝑓(…). It is usually much simpler and more natural to specify a model in this way, as and such they admit a much larger class of problems than explicit models. We can recognize two distinct subclasses:

-With ordinary differential equations4 the propagation function 𝑓 is a function of the state of the system, i.e. 𝜕𝑢/𝜕𝑡 = 𝑓(1, 𝑢, 𝑢 2 , …), -With partial differential equations the propagation function also contains higher-order partial derivatives with respect to other (spatial) features, i.e. 𝜕𝑢 𝜕𝑡 = 𝑓(1, 𝑢, 𝑢 𝑥 , …). The distinction between ODEs and PDEs might seem minor, but, as I will discuss later on, the spatial derivatives make discovering PDEs a much more challenging problem.

These categories are ofcourse an oversimplification, but they correspond well to the different approaches and main problems encountered in model discovery. No matter the type of model, a first, naive approach to model discovery would be to simply fit all possible models and compare them. However, there are simply too many models to construct and test -such an approach would take prohibitively long, even on small datasets with powerful computers. Luckily, this brute-force approach is not necessary: all we need is a method to smartly search the space of all models.

amounts. While he initially thought this to be wrong, he inadvertently (and unwillingly) opened the door to quantum physics!

Symbolic regression

Symbolic regression covers a large class of approaches to model discovery, which all, to a varying degree, are based on searching the space of all possible models. Before discussing various approaches to this search, we must first understand how to computationally construct a symbolic equation. Equations can be considered so-called expression trees. For example, the expression tree for the equation (𝑎 + 𝑏) ⋅ 𝑐 + 7 is given by figure 1.2. The expression tree approach allows to computationally work with symbolic equations -for example, to represent the equation (𝑎 + 𝑏) ⋅ 𝑐 + 7 ⋅ 𝑑, simply add '*' and 'd' nodes to the '7' node. In other words, expression trees also easily allow to change a given equation. This idea forms the core of genetic programming (GP), one of the most popular approaches to symbolic regression. Starting with a random population of these trees, GP randomly mutates the trees, conserving only the well-performing trees and repeat this until the convergence, analogous to how a population of genes evolves [START_REF] Bongard | Automated reverse engineering of nonlinear dynamical systems[END_REF][START_REF] Schmidt | Distilling Free-Form Natural Laws from Experimental Data[END_REF][START_REF] Maslyaev | Data-driven PDE discovery with evolutionary approach[END_REF]. Alternatively, Petersen [2020] train a recurrent neural network to generate the trees corresponding to the unknown equations. A completely different approach to symbolic regression is taken by the team behind AI-Feynman [Udrescu andTegmark, 2019, Udrescu et al., 2020]. They consider specifically equations originating from physics, and note that these equations often possess symmetry, separability or other various simplifying properties. By recursively applying these simplifications, the problem is split into a number of much easier to solve sub-problems, discovering a variety of seminal equations from physics. A final approach worth mentioning is that of the Bayesian model scientist [START_REF] Guimerà | A Bayesian machine scientist to aid in the solution of challenging scientific problems[END_REF]. While the entire model space might be (infinitely) large, in practice many models and equations look alike or share much of their expression. Having mined all the equations on Wikipedia, the Bayesian model scientist generates various hypotheses for the form of the equation by combining fitting with this knowledge of all previous equations. Symbolic regression is an active area of research, with the field seemingly focusing on discovering complex, explicit equations from synthetic, low-noise datasets [START_REF] Silviu | AI Feynman: a Physics-Inspired Method for Symbolic Regression[END_REF]. However, most real, physical systems are governed by (coupled) differential equations, and the observations are noisy and sparse. Symbolic regression struggles in these settings, and as I will discuss later, the difficulty lies not in finding the equation, but in accurately calculating the features (i.e. the derivatives). An alternative approach based on sparse regression is much more for these systems and settings.

Sparse regression

Model discovery can be greatly simplified if we consider only differential equations -this doesn't reduce the generality of our approach, as many physical systems are expressed as DEs. Considering again equation (1.1), note that it essentially is a linear combination of features, just like many other PDEs. This implies that many models can be written as5 , 𝜕𝑢 𝜕𝑡 =𝑤 0 𝑋 0 +…+𝑤 𝑀 𝑋 𝑀 (1.2)

=𝑤𝑋 𝑇 (1.3)
where 𝑤 𝑖 is the coefficient corresponding to feature 𝑋 𝑖 and 𝑋 the full feature matrix. Even within this restricted model space (i.e. the model is a linear combination of features), a combinatorial approach remains prohibitively expensive. However, by constructing a large set of possible features and combining these into a single matrix Θ, model discovery can be approached as a regression problem,

̃𝑤 = min 𝑤 ∥ 𝜕𝑢 𝜕𝑡 -𝑤Θ∥ 2 .
(1.4)

where ̃𝑤 determines the model: ̃𝑤 𝑖 =0if feature 𝑖 is not a part of the equation. More precisely, this has transformed model discovery into a variable selection problem; finding the correct model involves selecting the correct (i.e. non-zero) components of 𝑤. The number of candidate terms 𝑁 is typically much larger than the number of terms making up the equation 𝑀 , implying that most components of ̃𝑤 should be zero. To promote such solutions, a sparsity promoting penalty is added -hence this approach is known as sparse regression. It was pioneered by [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF] and has since become the de-facto method of performing model discovery.

Contributions

The sparse regression approach is elegant, flexible and widely applicable, but struggles on noisy and sparse data. When working with differential equations, the features making up the equation (𝜕𝑢/𝜕𝑡 for both ODEs and PDEs and higher order spatial derivatives for PDEs) are not directly observed. Rather, the observed data 𝑢(𝑥, 𝑡) must first be numerically differentiated to calculate the features. Numerical differentiation fundamentally is build on approximations, and when faced with sparse and noisy data it produces highly inaccurate results6 . This in turn fundamentally limits model discovery to low-noise, densely sampled datasets; if the features are corrupted, no approach will be able to select the right ones. As higher-order derivatives are more inaccurate, this issue specifically affects model discovery of PDEs, which often contains several higher-order spatial derivatives. If the calculation of the features is the limiting factor, this implies that model discovery should be approached as a distinct two-step process: calculating the features (step 1) is just as important, if not more, as selecting the right variables (step 2).

Motivated by this line of reasoning, we realized that if our goal was to apply model discovery on experimental data, improving the accuracy of the features would have much more impact than constructing more robust sparse regression algorithms. A common approach to improve the accuracy of the features is to approximate the data using some data-driven model, for example using polynomials or a fourier series, and perform sparse regression based on the features calculated from that approximation. My work is build on the idea that by integrating these two steps, model discovery can be made much more robust: the approximation makes it easier to discover the governing equation, while this equation in turn can be used to improve the approximation. In other words, I argue that data-driven and first-principle modelling are not opposites, but that they can mutually improve each other7 .

ORGANIZATION OF THIS THESIS

Specifically, my work focuses on 1) constructing accurate representations of (noisy) data using neural networks and 2) how to incorporate the knowledge of the governing equation back into these networks. Our overarching goal was to improve model discovery of partial differential equations on noisy and sparse datasets, opening the way for use on real, experimental datasets.

In Both et al. [2019] we show that such neural networks significantly improve the robustness of model discovery compared to unconstrained neural networks or classical approaches. Improving on this, in Both et al. [2021b] we present a modular framework, showing how these constrained networks can utilize any sparse regression algorithm and how to iteratively apply the sparse regression step, boosting performance further. In our latest work [Both and Kusters, 2021], we build upon the first two articles to construct a fully differentiable model discovery algorithm by combining multitask learning and Sparse Bayesian Learning. Besides these methodological papers, we show in Both et al. [2021a] that neural network based model discovery needs much less samples to find the underlying equation when they are randomly sampled, as opposed to on a grid. This marks a fundamental difference with classical approaches, which struggle with off-grid sampled data. In Both and Kusters [2019] we introduce Conditional Normalizing Flows. CNFs allow to learn a fully flexible, time-dependent probability distribution, for example to construct the density of single particle data. In Both and Kusters [2021] we build on this to discover a population model directly from single particle data. This series of articles was released under the name DeepMoD -Deep Learning based Model Discovery. Most of this work is collected in an open-source software package called DeePyMoD. Taken together, our work strongly establishes the argument for neural network based surrogates, constrained by the underlying physics, for model discovery on real, experimental data. 8

Organization of this thesis

I've written my thesis using Distill for R Markdown, a format based on the excellent online journal Distill. Distills typesetting and typography is strongly based on the Tufte style and is optimized for digital viewing: I've included some interactive plots, and hover over footnotes to see their content 9 or over citations to see the full one -often with a link to an Arxiv version [Both and Kusters, 2021]. All this means that while you can read it on paper, you will miss out on these features, so I highly recommend using this ing providing better features, while first-principle modelling improves data-driven modelling by infusing it with physical biases.

8 I've worked on several other topics, but I have chosen to focus on my main projects. 9 No more scrolling! website 10 . This thesis itself is also a small experiment. Instead of writing a long static document which is likely to only be read by the committee, can I write something a little more modern which will be useful to others, perhaps those not familiar with the field? Instead of writing a comprehensive review examining all the papers in the field, can I write something which gives a broader overview, discusses key papers in-depth and introduce the key challenges? The website in front of you is the answer I came up with.

• I have split the background of my work into two pages: Regression shortly discusses the regression approach to model discovery and presents the key papers of the field. As I argue in this thesis, surrogates should be a key component of model discovery algorithms. In Surrogates I discuss various classical approaches, but will spend the bulk of the chapter on Physics Informed Neural Networks (PINNs). • DeepMoD is the umbrella under which I've collected all of my work done during my PhD on model discovery from January 2019 until July 2021. I start by explaining the main idea behind our approach, and then summarize our papers and list their respective contributions. • In Discussion I reflect on our work and put it in context. Finally, I consider the open questions and challenges for model discovery.

Each chapter article can be read on its own, but I suggest following the order in which I've introduced them above. I hope this thesis is able to convey all that I've learned and done the last few years -happy reading! Chapter 2

Regression

Figure 2.1: One-dimensional solution of the Kuramoto-Shivashinsky equation.

Consider again the Kuramoto-Shivashinsky equation,

𝜕𝑢 𝜕𝑡 = 𝜕 2 𝑢 𝜕𝑥 2 - 𝜕 4 𝑢 𝜕𝑥 4 -𝑢 𝜕𝑢 𝜕𝑥 . (2.1)
It is a fourth order, non-linear partial differential equation giving rise to the beautifully rich behaviour of figure 2.1. This complexity -a fairly unusual 11 fourth-order derivative and a non-linear term -also poses a challenge for model discovery. Is it really possible to discover an equation this complex autonomously from data? In this chapter I introduce how sparse regression makes this possible, using the Kuramoto-Shivashinksy equation as an example1 . Before reviewing regression and various sparsity-promoting regularizations, I shortly consider why the sparse regression technique can be used. This chapter ends with an overview of the key extensions which have been developed.

Structure of differential equations

While the Kuramoto-Shivashinksy equation appears complex, its structure is not particularly unique. In fact, it shares many features with other models.

Consider for example the Cahn-Hilliard equation, which describes phase separation of a fluid with concentration 𝑐, diffusion coefficient 𝐷 and a length scale √ 𝛾:

𝜕𝑐 𝜕𝑡 =𝐷∇ 2 (𝑐 3 -𝑐-𝛾∇ 2 𝑐) , (2.2)
or the Korteweg-de Vries equation, which describes the speed 𝑢 of waves in shallow water,

𝜕𝑢 𝜕𝑡 =- 𝜕 3 𝑢 𝜕𝑥 3 + 6𝑢 𝜕𝑢 𝜕𝑥 .
(2.3)

All three equations describe completely different systems, but share a similar structure: each consists of 1) a linear combination of 2) relatively simple features2 . Indeed, considering the list of non-linear partial differential equations on Wikipedia, only very few equations have features composed of more than two terms, such as 𝑢⋅𝑢 𝑥 ⋅𝑢 𝑥𝑥𝑥 . Property 2 implies that most models can be constructed from a small, finite set of elements consisting of (2.8)

This least-squares problem is trivial to solve6 , but will, unfortunately, not yield the correct result in most cases7 . Since the number of features 𝑀 is likely to be much larger than the number of terms making up the unknown equation, solving eq. (2.8) will overfit the model: very little if any components of ξ will be zero. Luckily, equation (2.8) can be adapted to yield results ξ win which many components are exactly 0, a property known as sparsity. This is known as sparse regression and this approach to model discovery was pioneered by [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF] for ODEs and extended to PDEs in [START_REF] Samuel | Data-driven discovery of partial differential equations[END_REF] . Various ways to promote the sparsity of solutions to eq. (2.8) exist, and I discuss those in the next section.

Regularized regression

To promote sparsity of the resulting vector, we augment the objective function eq. (2.8) with a regularization function 𝑅(𝜉) on the unknown coefficients 𝜉8 ,

ξ=min 𝜉 1 𝑁 ‖𝑢 𝑡 -Θ𝜉‖ 2 2 + 𝜆𝑅(𝜉) (2.9)
where 𝜆 sets the strength of the regularization. The choice of 𝑅(𝜉) affects which behaviour is penalized. Three common choices are,

• 𝑅(𝜉) = ‖𝜉‖ 0 , also known as ℓ 0 regularization, penalizes the number of non-zero components, promoting sparsity in the final solution. • 𝑅(𝜉) = ‖𝜉‖ 1 , known as ℓ 1 regularization, penalizes the magnitude of 𝜉 linearly, also promoting sparsity in the final solution. • 𝑅(𝜉) = ‖𝜉‖ 2 2 , also called ℓ 2 regularization, penalizes the magnitude of 𝜉 quadratically, preventing large magnitude coefficiens in the final solution.

These penalties can also be combined: for example, the elastic net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] combines the benefits of both ℓ 1 and ℓ 2 regularization by defining

𝑅(𝜉) = 𝛼 ‖𝜉‖ 1 + (1 -𝛼) ‖𝜉‖ 2 2 .
(2.10)

Since these regularizations are central to sparse regression and model discovery, I want to discuss each of these three regularizations in-depth.9 2.3.1 ℓ 0 ℓ 0 regularization penalizes the number of non-zero components through the penalty function, (2.11) This makes it an obvious choice for promoting sparsity: if the amount of terms is penalized, solutions with fewer terms are favoured, and additionally, it does not affect the magnitude of the components, making it an unbiased estimator. Unfortunately, by penalizing the number of terms (also known as the support) the objective function becomes non-differentiable and nonconvex. Simply put, this regularization still effectively involves solving a combinatorial problem. [START_REF] Maddu | Stability selection enables robust learning of partial differential equations from limited noisy data[END_REF] show that an approximation known as Iterative Hard Thresholding can be an effective method in model discovery. A different avenue of research tackles these issues head-on by taking a probabilistic approach. For example, [START_REF] Louizos | Learning Sparse Neural Networks through $L_0$ Regularization[END_REF] consider the penalty as the sum of Bernouilli random variables, allowing gradientbased optimization. Nonetheless, these approaches are still computationally expensive.

‖𝜉‖ 0 =∑ 𝑗 𝜉 * 𝑗 ,𝜉 * 𝑗 ={ 1 if 𝜉 𝑗 ≠0 0 if 𝜉 𝑗 =0

ℓ 1 -Lasso

Given these issues, a common approach is to relax the problem and apply ℓ 1 regularization,

‖𝜉‖ 1 =∑ 𝑗 |𝜉 𝑗 |.
(2.12)

Compared to an ℓ 0 regularized problem, this too yields sparse solutions (I will discuss why it does later on), but is convex and can be efficiently solved using proximal method such as FISTA [START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF], despite its non-differentiability at 𝑥=0 . Usually this approach is called Lasso (Least Absolute Shrinkage and Selection operator), as it the sparsest norm still yielding a convex problem [START_REF] Tibshirani | Regression shrinkage and selection via the lasso: a retrospective: Regression Shrinkage and Selection via the Lasso[END_REF].

The Lasso is a popular approach when sparsity is required, but it does not come without its drawbacks. It yields a biased estimate of the coefficients, but a more significant issue is that it is not guaranteed to be consistent10 . An estimator is consistent if the estimate approaches the true value with increasing sample size:

ξ→𝜉 true as 𝑁→∞ (Consistency)
While this might seem like a trivial point, its implications are not: an inconsistent estimator will not be able to recover the truth, even with an infinite number of samples. The Lasso is consistent only if the Irrepresentable condition (IRC) is satisifed [Zhao and Yu]. In practice the IRC is often violated, implying the Lasso is unable to recover the ground truth. The violation is caused by correlation between variables, and the lasso can be made consistent by scaling the components in the penalty. For example, the adaptive Lasso [START_REF] Zou | The Adaptive Lasso and Its Oracle Properties[END_REF] scales the penalty as

∑ 𝑗 𝑤 𝑗 |𝜉 𝑗 | (Adaptive Lasso)
where 𝑤 𝑗 = 1/ |𝜉 * 𝑗 | 𝛾 , with 𝜉 * 𝑗 the result of a different consistent estimator and 𝛾 a scaling factor.

ℓ 2 -Ridge

The third popular regularization is an ℓ 2 penalty,

∑ 𝑗 |𝜉 𝑗 | 2 , (Ridge)
also known as ridge regression. Ridge regression can be analytically solved, thus having essentially the same computational cost as least squares11 , and can handle highly correlated variables. Unfortunately, ridge regression favours results where none of the coefficients are large, but not exactly zero either -in other words, it does not yield a sparse coefficient vector. To see why Lasso yields sparse results while ridge regression does not, consider figure 2.4. The ovals show the level sets of the loss due the fitting term (i.e. the first term of eq. (2.9)) and the loss due the regularization is represented by the circle (ridge) and square (lasso). The solution of the regularized problem is the intersection of these two losses. For the ℓ 1 regularization, the intersection happens at the corners, leading to sparse coefficients, whereas for the ℓ 2 regularization it happens in the middle, yielding small but non-zero coefficients. Model discovery often relies on heuristics and some approaches based on ridge regression show good performance, despite it not mathematically yielding sparse coefficients. For example, SINDy [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF], one of the most popular model discovery approaches, uses ridge regression.

Heuristics

Besides these regularizations, various other heuristic are often applied to improve the sparsity of the obtained coefficients. I discuss two commonly used ones: thresholding and sequential regression.

Thresholding

It is quite common that we recover a sparse vector which is approximately right: the required terms stand out, but several other small but non-zero components exist, despite regularization. Performance can be strongly improved by thresholding this sparse vector. Since all components have different dimensions, the features are first normalised by their ℓ 2 norm12 ,

Θ * 𝑗 = Θ 𝑗 ∥Θ 𝑗 ∥ 2 , (2.13)
allowing us to define normalized coefficients as,

𝜉 * 𝑗 =𝜉 𝑗 ⋅ ∥Θ 𝑗 ∥ 2 . (2.14)
Given some threshold value 𝜂, the threshold operation is then defined as

𝜉 𝑗 ={ 𝜉 𝑗 if 𝜉 * 𝑗 ≥𝜂 0 if 𝜉 * 𝑗 <𝜂 (2.15)
The value of 𝜂 is usually preset (i.e. a hyperparameter), but it can also be learned [START_REF] Samuel | Data-driven discovery of partial differential equations[END_REF].

Sequential regression

A second trick is to sequentially refine the estimate by iteratively performing sparse regression on only the active features selected by the previous estimator,

𝜉 𝑖+1 = min 𝜉 𝑖 1 𝑁 ∥𝑢 𝑡 -Θ 𝑖 𝜉 𝑖 ∥ 2 2 + 𝜆𝑅(𝜉 𝑖 ) (2.16) Θ 𝑖+1 =Θ 𝑖 [𝜉 𝑖+1 ≠ 0]
(2.17) until 𝜉 converges. While this does bias the results as coefficients can only be removed, not added, in practice it works well.

Figure 2.5: Illustration of sequential regression.

Extensions

Sparse regression forms a very flexible and general approach to model discovery and has been adapted, extended and improved in various ways since its introduction. In the next section I discuss the key papers which have either solved important open questions, or make significant progress towards it.

Stability selection

In the previous section we considered the effect of different regularizations on sparsity. An additional important factor is the strength of this regularization, denoted by the factor 𝜆 in eq. (2.9). Figure 2.6 shows the solution ξ as a function of the regularization strength 𝜆. The resulting sparse vector is strongly dependent on 𝜆; the lower the amount of regularization, the more features are selected. The default method to choose hyper-parameters such as 𝜆 is cross-validation (CV). CV consists of splitting the data into 𝑘 folds 13 , training the model on 𝑘-1 folds and testing model performance on the remaining fold. This process is repeated for all folds and averaging over the results yields a data efficient but computationally expensive approach to hyperparameter tuning. CV optimizes for predictive performance -how well does the model predict unseen data? -, but model discovery is interested in finding the underlying structure of the model. Optimizing for predictive performance might not be optimal; a correct model certainly generalizes well, but optimizing for it on noisy data will likely bias the estimator to include additional, spurious, terms. An alternative metric which optimizes for variable selection is stability selection [START_REF] Meinshausen | Stability Selection[END_REF].

The key idea is that while a single fit probably will not be able to discriminate between required and unnecessary terms we expect that among multiple fits on sub-sampled datasets the required terms will consistently be non-zero. Contrarily, the unnecessary terms are essentially modelling the error and noise, and will likely be different for each subsample. More specifically, stability selection bootstraps a dataset of size 𝑁 into 𝑀 subsets 𝐼 𝑀 of 𝑁 /2 samples, and defines the selection probability Π 𝜆 𝑗 as the fraction of subsamples in which a term 𝑗 is non-zero

Π 𝜆 𝑗 = 1 𝑀 ∑ 𝑀 1(𝑗 ∈ 𝑆 𝜆 (𝐼 𝑀 )) (2.18)
where 𝑆 𝜆 𝑗 is the recovered support, i.e. whether a term is non-zero. The support can be determined using any sparse regression method: [START_REF] Meinshausen | Stability Selection[END_REF] use Lasso, but [START_REF] Maddu | Stability selection enables robust learning of partial differential equations from limited noisy data[END_REF] use a relaxed ℓ 0 penalty. Calculating the selection probability for a range of regularization strengths Λ yields the stability path (figure 2.7). As the active components are determined from a range of 𝜆, stability selection does not require making a point estimate for the regularization strength as with CV. Stability selection has been shown to outperform cross validation in model discovery, especially at higher noise levels [START_REF] Maddu | Stability selection enables robust learning of partial differential equations from limited noisy data[END_REF][Maddu et al., , 2020]].

Multi-experiment datasets

A dataset rarely consists of a single experiment -they are comprised of several experiments, each with different parameters, initial and boundary conditions. Despite these variations, they all share the same underlying dynamics. To apply model discovery on these datasets we need a mechanism to exchange information about the underlying equation among the experiments. Specifically, we need to introduce and apply the constraint that all experiments share the same support (i.e. which terms are non-zero).

Consider 𝑘 experiments, all taken from the same system, {Θ ∈ ℛ 𝑁×𝑀 ,𝑢 𝑡 ∈ ℛ 𝑁 } 𝑘 , with an unknown sparse coefficient vector 𝜉∈ℛ 𝑀×𝑘 . Applying sparse regression on each experiment separately is likely violate this constraint; while it will yield a sparse solution for each experiment, the support will not be the same (see figure 2.8). Instead of penalizing single coefficients, we must penalize groups -an idea known as group lasso14 [START_REF] Huang | The Benefit of Group Sparsity[END_REF]. Group lasso first calculates a norm over all members in each group (typically an ℓ 2 norm), and then applies an ℓ 1 norm over these group-norms (see figure 2.8),

𝑅(𝜉) = 𝑘 ∑ 𝑗=1 | √ √ √ ⎷ 𝑀 ∑ 𝑖=1 |𝜉 𝑖𝑗 | 2 | (Group sparsity)
This approach drives groups of coefficients to zero, rather than single coefficients, and thus satisfies the constraint (see figure 2.8). Group sparsity has been successfully used for model discovery on multiple experiments [START_REF] De Silva | Discovery of Physics from Data: Universal Laws and Discrepancy Models[END_REF], or with space-or time-dependent coefficients (see next section), but can also be exploited to impose symmetry [Maddu et al., 2020].

Figure 2.8: Lasso applies the penalty vertically, while group lasso applies the penalty horizontally.

Time/space dependent coefficients

So far we have implicitly assumed that all the governing equations have fixed coefficients. This places a very strong limit on the applicability of model discovery; usually coefficients are fields depending on space or time (or even both!), so that 𝜉 = 𝜉(𝑥, 𝑡). As a first attempt at solving this, [START_REF] Rudy | Data-driven identification of parametric partial differential equations[END_REF] learn the spatial or temporal dependence of 𝜉 by considering each corresponding slice as a separate experiment with different coefficients, and applying the group sparsity approach we introduced in the previous section.

A Bayesian version of this approach was studied by [START_REF] Chen | Robust data-driven discovery of partial differential equations with time-dependent coefficients[END_REF].

It is important to note that in both these works 𝜉(𝑥, 𝑡) is not parametrized; it is not the functional form which is inferred, but simply its values at the locations of the samples.15

Information-theoretic approaches

The sparse regression approach tacitly assumes that the sparsest equation is also the correct one. While not an unreasonable assumption, it does require some extra nuance, as there often is not a single correct model. Many (effective) models are constructed as approximations of a certain order depending on the accuracy required. A first order approximation yields the sparsest equation, while the second order model could describe the model better. Another example is that a system can be locally modelled by a simpler, sparser equation. In both cases neither model is incorrect -they simply make a different trade-off. Information-theoretic approaches give a principled way to balance sparsity of the solution with accuracy.

Consider two models 𝐴 and 𝐵, each with likelihood ℒ 𝐴,𝐵 and 𝑘 𝐴,𝐵 terms. Model B has a higher likelihood, ℒ 𝐵 >ℒ 𝐴 , indicating a better fit to the data, but also consists of more terms, 𝑘 𝐵 >𝑘 𝐴 . The Bayesian Information Criterion (BIC) and the closely related Akaike Information Criterion (AIC) are two metrics to decide which balance these two objectives,

BIC =𝑘ln 𝑛-2ln ℒ (2.20)
and the Akaike Information Criterion, AIC = 2𝑘 -2 ln ℒ.

(2.21)

If model 𝐵 has a lower AIC or BIC than model 𝐴 it is a better, even if 𝑘 𝐵 >𝑘 𝐴 . In other words, these criteria give a principled way to decide the trade-off between how well the model fits fits the data and the complexity of the model. This approach has been used by [START_REF] Mangan | Model selection for dynamical systems via sparse regression and information criteria[END_REF] to decide between various discovered models, and a similar approach based on minimum description length was used by Udrescu and Tegmark [2019].

As we discussed, it is likely for multiple models to be correct simultaneously, all with similar BIC or AIC. Each of these models is essentially the model best describing the data at a certain accuracy, and all of the models together describe a curve in model space known as the Pareto Frontier16 . All models on this curve are correct, but of varying complexity, and we merely need to pick the model corresponding to the accuracy to which we wish to describe the data [Udrescu et al., 2020].

Bayesian approaches

The requirement that the resulting coefficient vector 𝜉 be sparse is essentially a form of prior knowledge. A principled way to include such prior knowledge is given by Bayes' theorem, which states for some data 𝑋 and parameters

𝑤 posterior ⏞ 𝑝(𝑤 | 𝑋) = likelihood ⏞ 𝑝(𝑋 | 𝑤) prior ⏞ 𝑝(𝑤) 𝑝(𝑋)
.

(2.22)

The likelihood describes the probability of the data 𝑋 for a given parameter 𝑤, while the prior describes the probability distribution of the parameter 𝑤. The posterior then is the distribution of 𝑤 given the data and the prior. Contrary to maximum likelihood approaches which yield point estimates, Bayesian approaches yield distributions, allowing to quantify uncertainty. Additionally, the inclusion of the prior effectively regularizes the problem, making Bayesian regression typically much more robust than MLE approaches.

The sparsity of 𝜉 can be encoded in the prior in various ways. A common approach known as Sparse Bayesian Learning [START_REF] Tipping | tipping01a[END_REF] takes a zerocentred Gaussian as prior,

𝑝(𝜉 𝑗 |𝛼 𝑗 ) = 𝒩(𝑤 𝑖 | 0, 𝛼 -1 𝑗 ).
(2.23)

The SBL assumes all terms are zero, but the confidence of that decision (𝛼 𝑗 ) is inferred from data. As 𝛼 𝑗 →∞ , the prior becomes infinitely sharp at zero, and the component is pruned from the model. The Gaussian prior is not particularly sparse, but [START_REF] Yuan | Machine Discovery of Partial Differential Equations from Spatiotemporal Data[END_REF] nonetheless show good performance with model discovery. Sparsity can be more strongly promoted by using a Laplacian prior to yield a Bayesian Lasso [Helgøy and Li, 2020], or a so-called spike-and-slab-prior [START_REF] Nayek | On spike-and-slab priors for Bayesian equation discovery of nonlinear dynamical systems via sparse linear regression[END_REF].17

Chapter 3

Differentiation and surrogates

Surrogates 1 are a widely used approach to approximate a dataset or model by a different, data-driven model with certain desirable properties. For example, a computationally cheap surrogate can be used to replace a computationally expensive simulation, or an easily-differentiable surrogate can be used to approximate a dataset we wish to differentiate. I discuss various surrogates used in the model discovery community and make the case for neural networks as surrogates. A neural network-based surrogate can be strongly improved by constraining it to the underlying physics using an approach known as Physics Informed Neural Networks (PINNs), and I end the chapter with a short introduction to Normalizing Flows (NFs). Before all of this though, I want to discuss why surrogates are needed in model discovery.

The need for surrogates

Using the sparse regression approach to discover PDEs requires calculating various higher-order derivatives. For example, a diffusive term is given by a second order derivative, ∇ 2 𝑢, while the Kuramoto-Shivashinksy equation requires a fourth order derivative. Consider the definition of the derivative,

𝑑𝑢 𝑑𝑥 = lim ℎ→0 𝑢(𝑥 + ℎ) -𝑢(𝑥) ℎ . (3.1)
1 also known as digital twins [START_REF] Rasheed | Digital Twin: Values, Challenges and Enablers From a Modeling Perspective[END_REF] 25 Since data is always sampled at a finite spacing, in practice a derivative must always be approximated,

𝑑𝑢 𝑑𝑥 ≈ 𝑢(𝑥 + ℎ) -𝑢(𝑥) ℎ . (3.2)
This equation forms the basis of numerical differentiation. In what follows we assume the data has been sampled on a grid with spacing Δ𝑥, allowing the use of a more natural notation in terms of grid index 𝑖. Equation (3.2) takes into account only the 'following' sample 𝑢 𝑖+1 at position +ℎ and is hence known as a forward-difference scheme. A more accurate estimate can be obtained by taking into account both the previous sample 𝑢 𝑖-1 and following sample 𝑢 𝑖+1 , yielding the popular second-order central difference scheme2 ,

𝑑𝑢 𝑖 𝑑𝑥 ≈ 𝑢 𝑖+1 -𝑢 𝑖-1 2Δ𝑥 . (3.3)
Applying eq. (3.3) recursively yields higher order derivatives,

𝑑 2 𝑢 𝑖 𝑑𝑥 2 ≈ 𝑢 𝑖-1 -2𝑢 𝑖 +𝑢 𝑖+1 Δ𝑥 2 (3.4) 𝑑 3 𝑢 𝑖 𝑑𝑥 3 ≈ -𝑢 𝑖-2 + 2𝑢 𝑖-1 -2𝑢 𝑖+1 +𝑢 𝑖+2 2Δ𝑥 3 (3.5)
Numerical differentiation through finite differences has various attractive properties: it's easy to implement and computationally cheap. However, these advantages come with large drawbacks:

• The approximation is only valid for small Δ𝑥. If the samples are spaced far apart, the resulting derivative will be inaccurate. Higher order schemes can alleviate this only slightly. • Higher order derivatives are calculated by successively applying the approximation (see (3.5)). Consequently, small errors in lower orders get propagated and amplified to higher orders, making estimates of these higher-order derivatives highly inaccurate. • At the edges of the dataset (e.g. 𝑢 0 ), the derivatives cannot be calculated. Doing so would require extrapolating, which is highly inaccurate, so typically the data at the edge is simply discarded. Since higher order derivatives utilize large stencils, this can amount to a significant part of the data when working with small datasets.

However, the dominant issue is numerical differentiation's inability to deal with noisy data. Even low noise levels of i.i.d. white noise (<10%) can corrupt the features to such an amount that model discovery becomes impossible. We illustrate these issues in figure 3.1, where we show with a toy function 𝑓(𝑥) = 𝑐𝑜𝑠(𝑥 2 ) the effect of data sparsity and 10% applied noise on the first and third derivative. Note that while the effect on the first derivative is still relatively limited, the third derivative is strongly corrupted. For the sparse data we can still qualitatively observe the behaviour, albeit with much lower peaks, but the noise has corrupted the data almost completely. Denoising the to-be-differentiated data and regularizing the differentiation process increases accuracy [START_REF] Van Breugel | Numerical Differentiation of Noisy Data: A Unifying Multi-Objective Optimization Framework[END_REF], but remains inadequate for very noisy data. Luckily, surrogates offer a way out. By fitting an easily-differentiable model such as a polynomial to the data, the problem of numerical differentiation is circumvented; the derivative of a polynomial can be trivially calculated symbolically. This does not mean that calculating the features is easy; the difficulty of determining the derivatives is now displaced to the task of accurately fitting the surrogate model the noisy data. In the next section I discuss the surrogate models which are often used in model discovery.

Surrogates: local versus global

Model discovery requires surrogates to be easily differentiable, but also to be flexible. The surrogate model must be able to accurately model the data: given some data 𝑢(𝑥), the surrogate 𝑝(𝑥) ought to have enough degrees of freedom to approximate 𝑢 reasonably, i.e. 𝑝(𝑥) ≈ 𝑢(𝑥). A classical approach is to represent 𝑝(𝑥) as a series expansion3 ,

𝑝(𝑥) = ∑ 𝑛 𝑎 𝑛 ℎ 𝑛 (𝑥) (3.6)
where 𝑎 𝑛 are picked such that 𝑝(𝑥) ≈ 𝑢(𝑥). The derivatives are now easily (and accurately!) calculated -it's simply the sum over the derivatives of the basis functions:

𝑑𝑢 𝑑𝑥 ≈ 𝑑𝑝(𝑥) 𝑑𝑥 =∑ 𝑖 𝑎 𝑖 𝑑ℎ 𝑖 (𝑥) 𝑑𝑥 (3.7)
The flexibility of the surrogate strongly depends on the choice of basis functions. A choice well-known to physicists is to use the Fourier basis,

ℎ 𝑛 (𝑥) = 𝑒 2𝜋𝑖𝑛𝑥 . (3.8)
Using Fourier series has several attractive properties: they're computationally efficient, well established, and are particularly useful for calculating derivatives. The fourier transform of a 𝑝th order derivative is4 ,

ℱ( 𝑑 𝑝 𝑓 𝑑𝑥 𝑝 ) = (𝑖𝑘) 𝑝 f(𝑘), (3.9)
which allows to efficiently calculate all derivatives5 . Additionally, the Fourier representation allows natural denoising of the data by applying a low-pass filter. This approach has been successfully applied to model discovery [START_REF] Schaeffer | Learning partial differential equations via data discovery and sparse optimization[END_REF], but is limited by its issues with sharp transitions (Gibbs ringing) and its requirement for periodic boundary condition.

A different choice of basis would be polynomials, but these do no have enough expressivity to model most data. However, data can locally be approximated polynomially, and spline interpolation exploits this idea by locally fitting a polynomial in a sliding window, and ensuring continuity at the edges6 . Various types of splines exist, depending on 1) the type of interpolating polynomial used, 2) the order of the interpolating polynomials and 3) the continuity conditions at the edges. They are computationally cheap, have well known properties, and can be used to smooth data through with a smoothness parameter. These properties make it a popular approach to calculate derivatives, but when data is extremely noisy, sparse or highdimensional results are inaccurate

Neural networks as surrogates

The most basic multilayer perceptron (MLP) neural networks cconsist of a basic matrix multiplication with an applied non-linearity,

𝑧 = 𝑓(𝑥𝑊 𝑇 + 𝐵) = ℎ(𝑥) (3.10)
where 𝑥∈ℛ 𝑁×𝑀 is the input, 𝑊∈ℛ 𝐾×𝑀 the kernel or weight matrix, 𝑏∈ℛ 𝐾 the bias and 𝑓 a non-linearity such as tanh. These layers are composed to increase expressive power, yielding a deep neural network,

𝑧=ℎ 𝑛 ∘ℎ 𝑛-1 …∘ℎ 0 (𝑥) = 𝑔 𝜃 (𝑥) (3.11)
in the rest of this thesis, we refer to a neural network of arbitrary depth and width as 𝑔 𝜃 (𝑥), where 𝜃 are the networks weights and biases.

Neural networks are excellent function interpolators, making them very suitable for surrogates: they scale well to higher dimensions, are computationally efficient and are very flexible7 . Given a dataset 𝑢, 𝑥 𝑖 , we typically train neural networks by minimizing the mean squared error,

ℒ= 𝑁 ∑ 𝑖=1 |𝑢 𝑖 -𝑔 𝜃 (𝑥 𝑖 )| 2 (3.12)
A neural networks flexibility is also its weakness; it is very sensitive to overfitting. To counter this, the dataset is split into a test and train set and train until the MSE on the test set starts increasing. While this works well in practice for many problems, the approach so far is completely datadriven; the only information used to train the network are the samples in the dataset. We often have more information available about a dataset -for example, its underlying (differential) equation8 . How can we include such knowledge in a neural network?

Physics Informed Neural Networks

One of the most popular ways to include physical knowledge in neural networks is using so-called Physics Informed Neural Networks (PINNs) [Raissi et al., 2017a,b]. Consider a (partial) differential equation

𝜕𝑢 𝜕𝑡 = 𝑓(…) (3.13)
and a neural network ̂𝑢 = 𝑔 𝜃 (𝑥, 𝑡). PINNs train a neural network by minimizing

ℒ(𝜃) = 1 𝑁 𝑁 ∑ 𝑖 | 𝜕û 𝑖 𝜕𝑡 -𝑓( û 𝑖 , …), | 2 . (3.14)
As ℒ→0, the output of the network ̂𝑢 satisfies the given differential equation9 , and the network learns the solution of a differential equation, without explicitly solving it.

PINNs can also be used to perform parameter inference. Given a dataset {𝑢, 𝑡, 𝑥} 𝑖 and a differential equation with unknown coefficients 𝑤, PINNs minimize,

ℒ(𝜃, 𝑤) = 1 𝑁 𝑁 ∑ 𝑖 |𝑢 𝑖 -û 𝑖 | 2 + 𝜆 𝑁 𝑁 ∑ 𝑖 | 𝜕û 𝑖 𝜕𝑡 -𝑓(𝑤, û 𝑖 , …), | 2 .
(3.15)

Note that we minimize both the neural network parameters 𝜃 and the unknown parameters 𝑤. Here the first term ensures the network learns the data, while the second term ensures the output satisfies the given equation.

PINNs have become a very popular method for both solving differential equations and for parameter inference, for various reasons:

• PINNs are very straightforward to implement, requiring no specialized architecture or advanced numerical approaches. The parametrization of the data or solution as a neural network makes PINNs a meshless method. Not having to create meshes is a second simplifying property. • Automatic differentiation can be used to calculate the derivatives, yielding machine-precision derivatives. • When used for parameter inference, PINNs essentially act as physically-consistent denoisers, making them particularly robust and useful when working with noisy and sparse data.

Notwithstanding their popularity and ease of use, PINNs are not without flaws -specifically, they can unpredictably yield a low accuracy solution, or even fail to converge at all [START_REF] Vincent Sitzmann | Implicit Neural Representations with Periodic Activation Functions[END_REF], Wang et al., 2020b,a].

Extensions & alternatives

PINNS have been extended and modified along several lines of research. To boost performance on complex datasets, it has proven fruitful to decompose the domain and apply a PINN inside each subdomain, while matching the boundaries [START_REF] Hu | When Do Extended Physics-Informed Neural Networks (XPINNs) Improve Generalization?[END_REF][START_REF] Shukla | Parallel Physics-Informed Neural Networks via Domain Decomposition[END_REF]. A different line of research instead focuses on appropriately setting the strength of the regularization term, adapting ideas from multitask learning 𝜆 [Maddu et al., 2021, McClenny and[START_REF] Mcclenny | Self-Adaptive Physics-Informed Neural Networks using a Soft Attention Mechanism[END_REF]. A final avenue worth mentioning is the inclusion of uncertainty, either through Bayesian neural nets [START_REF] Yang | B-PINNs: Bayesian Physics-Informed Neural Networks for Forward and Inverse PDE Problems with Noisy Data[END_REF] or using dropout [Zhang et al.].

Various other approaches have been proposed to include physical knowledge in neural networks. Champion et al. [2019b] use a variation on a PINN with auto-encoders to discover an underlying equation. Other works highlight the connection between numerical differentiation and convolutions and include the equation through this [START_REF] Long | PDE-Net: Learning PDEs from Data[END_REF].

A very promising set of approaches are differentiable ODE and PDE solvers [START_REF] Ricky | Neural Ordinary Differential Equations[END_REF][START_REF] Rackauckas | Universal Differential Equations for Scientific Machine Learning[END_REF], and are also often referred to as differentiable physics [START_REF] Ramsundar | Differentiable Physics: A Position Piece[END_REF]. Here, instead of modelling the solution 𝑢 by a neural network, we model the time-derivative 𝑢 𝑡 and apply a solver. The upside of this approach is that the solution can be guaranteed to satisfy the physics, contrary to the 'soft-constraint' or regularization approach of PINNs, but it does reintroduce all the nuances involved in solving PDEs -discretisation, stability issues and possible stiffness. Calculating the gradient involves solving the adjoint problem, which itself can also be unstable.

Kernel-based methods have also been used to incorporate physical knowledge -specifically Gaussian Processes [START_REF] Särkkä | The Use of Gaussian Processes in System Identification[END_REF][START_REF] Atkinson | Bayesian Hidden Physics Models: Uncertainty Quantification for Discovery of Nonlinear Partial Differential Operators from Data[END_REF]. GPs have build-in uncertainty estimation and can be made to respect a given equation by using it as a kernel function. Unfortunately, this requires linearisation of the equation as GPs cannot accommodate non-linear kernels.10

A final approach worth mentioning explores including more abstract principles. [START_REF] Cranmer | Lagrangian Neural Networks[END_REF] present an architecture which parametrizes the Lagrangian of a system, while [START_REF] Greydanus | Hamiltonian Neural Networks[END_REF] present a Hamiltonianpreserving architecture.

Normalizing Flows

We discussed already that neural networks are universal function approximators and can thus theoretically approximate any function, including probability distributions. However, for a function 𝑝(𝑥) to be a proper distribution we require Simply put, this results in a PINN which discovers its constraining equation during training1 . This approach has several benefits:

∫ +∞ -∞ 𝑝(𝑥) = 1, ( 3
• The neural network is a very flexible surrogate, and we can use automatic differentiation to calculate machine-precision derivatives, circumventing the numerical differentiating problem.

• By constraining the network to solutions of the constraining equation, it acts as a physically-consistent denoiser and interpolator. In other words, the data is denoised and interpolated with knowledge of some underlying equation -in contrast to for example splines, which are purely data-driven. • The integration of the two steps of model discovery -calculating features and sparse regression -into a single, differentiable model starts a 'virtuous cycle'. Accurate features improve the model discovery task, which in turn more strongly constrains the network, leading to a better estimate of the data, finally yielding more accurate features.

Equation (4.2) only describes the core idea, and in a series of three articles we explore and refine this idea under the name DeepMoD -Deep Learning for Model Discovery. Additionally, we worked on two tangential projects: one compares splines with DeepMoD and studies the effect of sampling strategies, while the second introduces Conditional Normalizing Flows (CNFs) to infer time-dependent probabilities and densities. In this chapter, I shortly discuss the contributions and methodology. followed by the results of each paper.

The result sections are adapted from the corresponding paper.

DeepMoD: Deep Learning for Model discovery in noisy data

Gert-Jan Both, Subham Choudhury, Pierre Sens, Remy Kusters -PDF This paper lays out and introduces the ideas I've sketched above, forming the basis for the rest of our work. We propose to train a neural network with the loss function

ℒ(𝜃, 𝜉) = ∑ 𝑖 ‖𝑢 𝑖 -û 𝑖 ‖ 2 2 ⏟⏟⏟⏟⏟⏟⏟ ℒ fit +∑ 𝑖 ‖( û 𝑡 ) 𝑖 -Θ 𝑖 𝜉‖ 2 2 ⏟⏟⏟⏟⏟⏟⏟ ℒ reg +𝜆∑ 𝑗 |𝜉 𝑗 | ⏟ ℒ sparse (4.3)
where ̂𝑢 𝑖 =𝑔 𝜃 (𝑥 𝑖 ,𝑡 𝑖 ) is the output of the neural network, and the features ̂𝑢 𝑡 and Θ are calculated from the output of the neural network. The first term in equation eq. ( 4.3) ensures the network learns the data, the second term constrains it to the given equation and the third term promotes the sparsity of 𝜉 by applying an ℓ 1 penalty. Training the network involves both optimizing the neural network parameters 𝜃 and the sparse coefficient vector 𝜉. After convergence, the found 𝜉 is normalized and thresholded to yield the final equation (see figure 4.1 for a schematic overview).

Results

We test the performance of this approach, which we call DeepMoD, on a set of case studies: the Burgers' equation with and without shock, the 2D advection-diffusion equation and the Keller-Segel model for chemotaxis. These examples show the ability of DeepMoD to handle (1) non-linear equations, (2) solutions containing a shock wave, (3) coupled PDEs and finally (4) higher dimensional and experimental data.

Non-linear PDEs

The Burgers equation occurs in various areas of gas dynamics, fluid mechanics and applied mathematics and is evoked as a prime example to benchmark model discovery [START_REF] Samuel | Data-driven discovery of partial differential equations[END_REF][START_REF] Long | PDE-Net: Learning PDEs from Data[END_REF] and coefficient inference algorithms [Raissi et al., 2017a,b], as it contains a non-linear term as well as second order spatial derivative,

𝜕 𝑡 𝑢 = -𝑢𝑢 𝑥 +𝜈𝑢 𝑥𝑥 . (4.4)
Here 𝜈 is the viscosity of the fluid and 𝑢 its velocity field. We use the dataset produced by Rudy et al. [START_REF] Samuel | Data-driven discovery of partial differential equations[END_REF], where 𝜈 = 0.1. The numerical simulations for the synthetic data were performed on a dense grid for numerical accuracy. DeepMoD requires significantly less datapoints than this grid and we hence construct a smaller dataset for DeepMoD by randomly sampling the results through space and time. From now on, we will refer to randomly sampling from this dense grid simply as sampling. Also note that this shows that our method does not require the data to be regularly spaced or stationary in time. For the data in Fig. 4.2 we add 10% white noise and sampled 2000 points for DeepMoD to be trained on.

We train the neural network using an Adam optimizer and plot the different contributions of the cost function as a function of the training epoch in Fig. 4.2a and we show the value of each component of 𝜉 as a function of the training epoch in Fig. 4.2b. Note that for this example, after approximately 2000 epochs, the MSE is converged, while at the same time we observe the components of 𝜉 only start to converge after this point. We can thus identify three 'regimes': in the initial regime (0 -2000 epochs), the MSE is trained. Since the output of the neural network is far from the real solution, so is Θ, and the regression task cannot converge (See first 2000 epochs in Fig. 4.2b). After the MSE has converged, ̂𝑢 is sufficiently accurate to perform the regression task and 𝜉 starts to converge. After this second regime (2000 -6000 epochs), all components of the cost converged (>6000 epochs) and we can determine the solution. From this, we obtain a reconstructed solution (see Fig. 4.2b) and at the same time recover the underlying PDE, with coefficients as little as 1% error in the obtained coefficients. We show the impact of including the regression term in the cost function in Fig. 4.2c, where the obtained solution of DeepMoD is compared with a neural network, solely trained on the MSE, to reconstruct the data. Including the regression term regularizes the network and prevents overfitting, despite the many terms in the library. We conclude that it is the inclusion of the regression in the neural network which makes DeepMoD robust to noisy data and prevents overfitting.

Next, we characterize the robustness of DeepMoD in Fig. 4.2d, where we run DeepMod for five times (differently sampled data set) for a range of sample sizes and noise levels. The color in Fig. 4.2d shows how many of the five runs return the correct equation and the value in the grid displays the mean error over all correct runs. Observe that at vanishing noise levels, we recover the Burgers equation with as little as 100 data-points, while for 5000 data points we recover the PDE with noise levels of up to 75%. Between the domain where we recover the correct equation for all five runs and the domain where we do not recover a single correct equation, we observe an intermediate domain where only a fraction of the runs return the correct equation, indicating the importance of sampling.

To benchmark DeepMoD, we can directly compare the performance of our algorithm with respect to two state-of-the-art methods, (i): PDE-Find by Rudy et al. [START_REF] Samuel | Data-driven discovery of partial differential equations[END_REF] and (ii) PDE-Stride by Maddu et al. [START_REF] Maddu | Stability selection enables robust learning of partial differential equations from limited noisy data[END_REF]. Considering an identical Burgers' data set with 10 5 data points, approach (i) recovers the correct equation for up to 1% Gaussian noise [START_REF] Samuel | Data-driven discovery of partial differential equations[END_REF] while method (ii) discovers the correct equation up to 5% noise [START_REF] Maddu | Stability selection enables robust learning of partial differential equations from limited noisy data[END_REF]. Compared to the results in Fig. 4.2d we note that even for two order of magnitude fewer samples points, 10 3 w.r.t. 10 5 , DeepMoD recovers the correct equation up to noise levels > 50% Gaussian noise. Deep-MoD allows up to two orders of magnitude higher noise-levels and smaller sample sizes with respect to state-of-the-art model discovery algorithms.

Shock wave solutions

If the viscosity is too low, the Burgers' equation develops a discontinuity called a shock (See Fig. 4.3). Shocks are numerically hard to handle due to divergences in the numerical derivatives. Since DeepMoD uses automatic differentiation we circumvent this issue. We adapt the data from Raissi et al. [START_REF] Raissi | Physics Informed Deep Learning (Part II): Data-driven Discovery of Nonlinear Partial Differential Equations[END_REF] which has 𝜈 = 0.01/𝜋, sampling 2000 points and adding 10% white noise (See Fig. 4.3). We recover ground truth solution of the Burgers' equation as well as the corresponding PDE, 𝜕 𝑡 𝑢 = -0.99𝑢𝑢 𝑥 + 0.0035𝑢 𝑥𝑥 , (4.5)

with a relative error of 5% on the coefficients. In Fig. 4.3 we show the inferred solution for 𝑡 = 0.8.

Coupled differential equations

Next, we apply DeepMod to a set of coupled PDE's in the form of the Keller-Segel (KS) equations, a classical model for chemotaxis [START_REF] Painter | Mathematical models for chemotaxis and their applications in self-organisation phenomena[END_REF]. Chemotactic attraction is one of the leading mechanisms that accounts for the morphogenesis and self-organization of biological systems. The KS model describes the evolution of the density of cells 𝑢 and the secreted chemical 𝑤,

𝜕 𝑡 𝑢 = ∇ ⋅ (𝐷 𝑢 ∇𝑢 -𝜒𝑢∇𝑤) 𝜕 𝑡 𝑤=𝐷 𝑤 Δ𝑤 -𝑘𝑤 + ℎ𝑢. (4.6)
Here the first equation represents the drift-diffusion equation with a diffusion coefficient of the cells, 𝐷 𝑢 and a chemotactic sensitivity 𝜒, which is a measure for the strength of their sensitivity to the gradient of the secreted chemical 𝑤. The second equation represents the reaction diffusion equation of the secreted chemical 𝑤, produced by the cells at a rate ℎ and degraded with a rate 𝑘. For a 1D system, we sample 10000 points of 𝑢 and 𝑤 for parameter values of 𝐷 𝑢 = 0.5, 𝐷 𝑣 = 0.5, 𝜒 = 10.0, 𝑘 = 0.05 and ℎ = 0.1 and add 5% white noise. We choose a library consisting of all spatial derivatives (including cross terms) as well as first order polynomial terms, totalling 36 terms. For these conditions we recover the correct set of PDEs, 𝜕 𝑡 𝑢 = 0.50𝑢 𝑥𝑥 -9.99𝑢𝑤 𝑥𝑥 -10.02𝑢 𝑥 𝑤 𝑥 𝜕 𝑡 𝑤 = 0.48𝑤 𝑥𝑥 -0.049𝑤 + 0.098𝑢, (4.7)

as well as the reconstructed fields for 𝑢 and 𝑤 (See Fig. 4.4). Note that even the coupled term, 𝑢 𝑥 𝑤 𝑥 , which becomes vanishingly small over most of the domain, is correctly identified by the algorithm, even in the presence of considerable noise levels.

Experimental data

To showcase the robustness of DeepMoD on high-dimensional and experimental input data, we consider a 2D advection diffusion process described by,

𝜕 𝑡 𝑢 = -∇ ⋅ (-𝐷∇𝑢 + ⃗ 𝑣 ⋅ 𝑢) , (4.8)
where ⃗𝑣 is the velocity vector describing the advection and 𝐷 is the diffusion coefficient. In the SI of the paper we apply DeepMod on a simulated dataset of eq. 4.8, with as initial condition, a 2D Gaussian with 𝐷 = 0.5 and ⃗ 𝑣 = (0.25, 0.5). For as little as 5000 randomly sampled points we recover the correct form of the PDE as well as the vector ⃗𝑣 for noise levels up to ≈ 25%. In the absence of noise the correct equation is recovered with as little as 200 sample points through space and time. This number is surprisingly small considering this is an 2D equation.

Finally, we apply DeepMoD on a time-series of images from an electrophoresis experiment, tracking the advection-diffusion of a charged purple loading dye under the influence of a spatially uniform electric field. In Fig. 4.5a we show time-lapse images of the experimental setup where we measure the time-evolution of two initial localised purple dyes. Fig. 4.5b shows the resultant 2D density field for three separate time-frames (in arbitrary units), corresponding to the red square in Fig. 4.5a by substracting the reference image (no dye present). The dye displays a diffusive and advective motion with constant velocity 𝑣, which is related to the strength of the applied electric field.

We apply DeepMoD on 5000 sampled data-points sampled through space and time and consistently recover the advection term 𝑢 𝑦 as well as the two diffusive components (𝑢 𝑥𝑥 and 𝑢 𝑦𝑦 ). In Fig. 4.5c we show the scaled coefficients as function of the number of epochs. After thresholding the scaled coefficients (|𝜉| < 0.1), we obtain for the unscaled coefficients, the resultant advection diffusion equation, 0.31𝑢 𝑦 + 0.011𝑢 𝑥𝑥 + 0.009𝑢 𝑦𝑦 = 0.

(4.9)

Analysing the second diffusing dye (dashed box in Fig. 4.5) result in nearly identical values for the drift velocity, 𝑣 ≈ 0.3, and the diffusion coefficients, 𝐷 ≈ 0.01 indicating the robustness of the obtained value of 𝐷 and 𝑣. In contrast to the artificial data presented in previous paragraphs, some higherorder non-linear terms, in particular 𝑢𝑢 𝑦𝑦 and 𝑢𝑢 𝑥𝑥 remain small, yet nonzero. This suggests that an automatic threshold strategy may not guarantee the desired sparse solution. Fixing a threshold of the scaled coefficients (|𝜉| < 0.1 in this particular case) or other thresholding strategies such as coefficient cluster detection would be better suited for this task.

Sparsely constrained neural networks for model discovery of PDEs

Gert-Jan Both, Gijs Vermarien, Remy Kusters -PDF

SPARSELY CONSTRAINED NEURAL NETWORKS FOR MODEL DISCOVERY OF PDES43

This paper was conceived and written mostly during the first confinement in Paris -plenty of time to think and wonder!

The original DeepMoD outperformed other model discovery approaches, but left several questions unanswered. Specifically, we wondered:

• How can the optimal strength of ℓ 1 penalty, 𝜆, be determined? Retraining the network would be prohibitively expensive, and while performance was not extremely sensitive to the choice of 𝜆, it is a hyperparameter which much be correctly chosen. • Lasso has several undesirable properties and perhaps suboptimal performance. Can other model discovery algorithms such as SINDy be integrated in this framework? • We never remove terms from the equation; rather, they approach or become zero due to the ℓ 1 penalty. Can performance be improved further removing the zero components and iteratively refining the estimate, similar to the iterative regression technique?

Next to answering these questions, this paper also represented a deepening of our understanding of the DeepMoD framework. We realized the coefficient vector 𝜉 essentially performs two tasks: 1) selecting which terms are active, i.e. the 'model discovery' task and 2) constraining the neural network. The main innovation in this paper is the introduction of a mask 𝑔 defining the support of 𝜉, allowing us to separate these two tasks:

ℒ(𝜃, 𝜉) = ∑ 𝑖 ‖𝑢 𝑖 -û 𝑖 ‖ 2 2 +∑ 𝑖 ‖𝑢 𝑡 -Θ(𝑔 ⋅ 𝜉)‖ 2 2 . (4.10)
This mask 𝑔 is calculated non-differentiably, i.e. it effectively acts as an external 'forcing variable'. In ML applications, non-differentiability is typically considered a short-fall (see next section), but here it forms a bridge between neural-network based approaches and classical approaches. It is through the mask 𝑔 that any sparse regression method can be utilized to perform variable selection on a neural-network based library. For example, lasso with cross validation can be used, but we also show an example with SINDy. We thus show that deep-learning based and classical methods can be synthesized, benefitting from both approaches. Together with this innovation we present a modular version of our framework, reprinted in figure 4.6.

In this framework, (I) A function approximator constructs a surrogate model of the data, (II) from which a library of possible terms and the time derivative is constructed using automatic differentiation. (III) A sparsity estimator constructs the sparsity mask 𝑔 to select the active terms in the library using some sparse regression algorithm and (IV) a constraint constrains the function approximator to solutions allowed by the active terms obtained from the sparsity estimator.

We released this paper together with a pytorch-based framework, DeePy-MoD. Overall the paper and the software package made our approach significantly more flexible: each module can be easily changed, independent of another.

Results

Constraint

The sparse coefficient vector 𝜉 in eq. (4.10) is typically found by optimising it concurrently with the neural network parameters 𝜃. Considering a network with parameter configuration 𝜃 * , the problem of finding 𝜉 can be rewritten as min 𝜉 |𝑢 𝑡 (𝜃 * ) -Θ(𝜃 * )𝜉| 2 . This can be analytically solved by least squares under mild assumptions; we calculate 𝜉 by solving this problem every iteration, rather than optimizing it using gradient descent. In figure 4.7 we compare the two constraining strategies on a Burgers data-set, by training for 5000 epochs without updating the sparsity mask. Panel A) shows that the least-squares approach reaches a consistently lower loss. More strikingly, we show in panel B) that the mean absolute error in the coefficients is three orders of magnitude lower. We explain the difference as a consequence of the random initialisation of 𝜉: the network is initially constrained by incorrect coefficients, prolonging convergence. The random initialisation also causes the larger spread in results compared to the least squares method. The least squares method does not suffer from sensitivity to the initialisation and consistently converges.

Sparsity estimator

Implementing the sparsity estimator separately from the neural network allows us to use any sparsity promoting algorithm. Here we show that a classical method for PDE model discovery, PDE-find [START_REF] Samuel | Data-driven discovery of partial differential equations[END_REF], can be used together with neural networks to perform model discovery in highly sparse and noisy data-sets. We compare it with the thresholded Lasso approach in figure 4.8 on a Burgers data-set with varying amounts of noise. The PDE-find estimator discovers the correct equation in the majority of cases, even with up to 60% -80% noise, whereas the thresholded lasso mostly fails at 40%. We emphasise that the modular approach we propose here allows to combine classical and deep learning-based techniques. More advanced sparsity estimators such as SR3 [Champion et al., 2019c] can easily be included in this framework.

Function approximator

We show in figure 4.9 that a tanh-based NN fails to converge on a data-set of the Kuramoto-Shivashinksy (KS) equation (panel A and B). Consequently, the coefficient vectors are incorrect (Panel D). As our framework is agnostic to the underlying function approximator, we instead use a SIREN [START_REF] Vincent Sitzmann | Implicit Neural Representations with Periodic Activation Functions[END_REF], which is able to learn very sharp features in the underlying dynamics. In panel B we show that a SIREN is able to learn the complex dynamics of the KS equation and in panel C that it discovers the correct equation 2 . This example shows that the choice of function approximator can be a decisive factor in the success of neural network based model dis-Figure 4.8: Fraction of correct discovered Burgers equations (averaged over 10 runs) as function of the noise level for the thresholded lasso and PDE-find sparsity estimator. covery. Using our framework we can also explore using RNNs, Neural ODEs [START_REF] Rackauckas | Universal Differential Equations for Scientific Machine Learning[END_REF] or Graph Neural Networks [START_REF] Seo | Differentiable Physics-informed Graph Networks[END_REF].

Fully differentiable model discovery

Gert-Jan Both, Remy Kusters -PDF

This is my favourite paper. I think it's an elegant approach solving many of the issues and questions of our previous papers, and as a bonus it also works in practice!

The modular approach in the second paper improved flexibility and performance, but left several things to be desired:

• It required the tuning of a fairly large set of hyperparameters: when to first update the mask, when to update it after, and which sparse regression approach to use (which itself introduces various hyperparameters). • The non-differentiability of the mask allows to easily incorporate various sparse regression techniques, but also presents a weakness. If it is set incorrectly (for any reason), the training fails as the constraint is incompatible with the data. Additionally, due to this nondifferentiability, model selection is effectively a side-effect of the training, and we hypothesize a fully differentiable method would increase performance and stability.

In this paper we create a fully-differentiable model discovery algorithm, which, considering eq. ( 4.10), requires making the mask 𝑔 differentiable. Differentiable masking is challenging due to the binary nature of the problem, and instead we relax the application of the mask to a regularisation problem. Specifically, we propose to use Sparse Bayesian Learning [START_REF] Tipping | tipping01a[END_REF] to select the active features and act as constraint.

PINNs constrained by Sparse Bayesian Learning

Sparse Bayesian Learning (SBL) [START_REF] Tipping | tipping01a[END_REF] using an iterative method (see [START_REF] Tipping | tipping01a[END_REF]). The marginal likelihood also offers insight how the SBL provides differentiable masking. Considering only the first two terms of eq. ( 4.16), 𝛽 ‖𝑢 𝑡 -Θ𝜇‖ 2 +𝜇 𝑇 𝐴𝜇 (4.17)

we note that the SBL essentially applies a coefficient-specific ℓ 2 penalty to the posterior mean 𝜇. If 𝐴 𝑗 →∞ , the corresponding coefficient 𝜇 𝑗 →0 , pruning the variable from the model. Effectively, the SBL replaces the discrete mask 𝑚 𝑗 ∈ {0, 1} by a continuous regularisation 𝐴 𝑗 ∈ (0, ∞], and we thus refer to our approach as continuous relaxation.

Model

To integrate SBL as a constraint in PINNs, we place a Gaussian likelihood on the output of the neural network,

û ∶ 𝑝(𝑢; û , 𝜏 ) = 𝑁 ∏ 𝑖=1 𝒩(𝑢 𝑖 ;û 𝑖 ,𝜏 -1 ), (4.18)
and define a Gamma hyper prior on 𝜏 , 𝑝(𝜏 ) = Γ(𝜏 ; 𝑒, 𝑓), yielding the loss function,

3 Neglecting the hyper-prior, this loss function can also written more compactly as Our approach does not rely on any specific property of the SBL, and thus generalises to other Bayesian regression approaches.

ℒ SBL (𝐴, 𝛽) = log |𝐶| + 𝜕 𝑡 𝑢 𝑇 𝐶 -1 𝜕 𝑡 𝑢, 𝐶 = 𝛽 -1 𝐼 + Θ𝐴 -1 Θ 𝑇 , ( 4

Training

The loss function for the SBL-constrained PINN contains three variables which can be exactly minimised, and denote these as  , τ and β . With these values, we introduce LSBL-PINN (𝜃) ≡ ℒ SBL-PINN (𝜃,  , τ , β ) and note that we can further simplify this expression as the gradient of the loss with respect to these variables is zero. For example, ∇ 𝜃 ℒ(  ) = ∇ 𝐴 ℒ⋅∇ 𝜃 𝐴| 𝐴=  =0, as ∇ 𝐴 ℒ| 𝐴=  =0 . Thus, keeping only terms directly depending on the neural network parameters 𝜃 yields,

LSBL-PINN (𝜃) = ̂𝜏 2 ‖𝑢 -û ‖ 2 + β 2 ‖𝑢 𝑡 -Θ𝜇‖ 2 +𝜇 𝑇 Â 𝜇 -log|Σ| = 𝑁τ 2 [ℒ fit (𝜃) + β̂𝜏 ℒ reg (𝜃, 𝜇)] ⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟ =ℒ PINN (𝜃,𝜇) +𝜇 𝑇 Â 𝜇 -log|Σ| (4.21)
where in the second line we have rewritten the loss function in terms of a classical PINN with relative regularisation strength 𝜆= β / τ and coefficients 𝜉=𝜇 . Contrary to a PINN however, the regularisation strength is inferred from the data, and the coefficients 𝜇 are inherently sparse.

An additional consequence of ∇ 𝜃 ℒ( Â , β , τ ) = 0 is that our method does not require backpropagating through the solver. While such an operation could be efficiently performed using implicit differentiation [START_REF] Shaojie Bai | Deep Equilibrium Models[END_REF], our method requires solving an iterative problem only in the forward pass. During the backwards pass the values obtained during the forward pass can be considered constant.

Considering eq. ( 4.21), we note the resemblance to multitask learning using uncertainty, introduced by [START_REF] Cipolla | Multi-task Learning Using Uncertainty to Weigh Losses for Scene Geometry and Semantics[END_REF]. Given a set of objectives, the authors propose placing a Gaussian likelihood on each objective so that each task gets weighed by its uncertainty. The similarity implies that we are essentially reinterpreting PINNs as Bayesian or hierarchical multi-task models.

Physics Informed Normalizing Flows

Having redefined the PINN loss function in terms of likelihoods (i.e. eq. (4.21)) allows to introduce a PINN-type constraint to any architecture with a probabilistic loss function. In this section we introduce an approach with normalizing flows, called Physics Informed Normalizing Flows (PINFs). As most physical equations involve time, we first shortly discuss how to construct a time-dependent normalizing flow. We show in the experiments section how PINFs can be used to directly infer a density model from single particle observations.

Conditional Normalizing Flows

Normalizing flows construct arbitrary probability distributions by applying a series of 𝐾 invertible transformations 𝑓 to a known probability distribution 𝜋(𝑧), Most physical processes yield time-dependent densities, meaning that the spatial axis is a proper probability distribution with ∫ 𝑝(𝑥, 𝑡)𝑑𝑥 = 1. Contrarily, this is not valid along the temporal axis, as ∫ 𝑝(𝑥, 𝑡)𝑑𝑡 = 𝑓(𝑥). To construct PINFs, we first require a Conditional Normalizing Flow capable of modelling such time-dependent densities. Instead of following the method of Both and Kusters [2019], which modifies the Jacobian, we employ time-dependent hyper-network. This hyper-network ℎ outputs the flow parameters 𝜃, and is only dependent on time, i.e. 𝜃 = ℎ(𝑡), thus defining a time-dependent normalizing flow as 𝑧=𝑔 ℎ(𝑡) (𝑥).

PINFs

Conditional normalizing flows yield a continuous spatio-temporal density, and the loss function of a PINF is defined as simply adding the SBL-loss to that of the normalizing flow, yielding

LPINF (𝜃) = ℒ NF (𝜃) + 𝑁 β 2 ℒ reg (𝜃, 𝜇) + 𝜇 𝑇 Â 𝜇 -log|Σ|. (4.23)

Results

We now show several experiments illustrating our approach. We start this section by discussing choice of hyperprior, followed by a benchmark on several datasets and finally a proof-of-concept with physics-informed normalizing flows.

Choosing prior

The loss function for the SBL constrained approach contains several hyperparameters, all defining the (hyper-) priors on respectively 𝐴, 𝛽 and 𝜏 .W e set uninformed priors on 𝐴 and 𝛽, 𝑎 = 𝑏 = 𝑒 = 𝑓 = 1𝑒 -6 , but those on 𝛽, the precision of the constraint, must be chosen more carefully. Figure 4.10 illustrates the learning dynamics on a dataset of the Korteweg-de Vries equation4 when the 𝛽 hyperprior is uninformed, i.e. 𝑐 = 𝑑 = 1𝑒 -6 . Observe that the model fails to learn the data, while almost immediately optimising the constraint. We explain this behaviour as a consequence of our assumption that the likelihoods factorise, which implies the two tasks of learning the data and applying the constraint are independent. Since the constraint contains much more terms than required, it can fit a model with high precision to any output the neural network produces. The two tasks then are not independent but conditional: a high precision on the constraint is warranted only if the data is reasonably approximated by the neural network. To escape the local minimum observed in figure 4.10, we couple the two tasks by making the hyper-prior on 𝛽 dependent on the performance of the fitting task.

Figure 4.10: Regression loss as a function of fitting loss during training, comparing an uninformed prior with a dynamic prior.

Our starting point is the update equation for 𝛽 (see [START_REF] Tipping | tipping01a[END_REF] for details), β= 𝑁-𝑀+∑ 𝑖 𝛼 𝑖 Σ 𝑖𝑖 + 2𝑐 𝑁ℒ reg + 2𝑑 (4.24)

We typically observe good convergence of normal PINNs with 𝜆=1 , and following this implies β≈ τ , and similarly ℒ reg →0as the model converges.

Assuming 𝑁≫𝑀+∑ 𝑖 𝛼 𝑖 Σ 𝑖𝑖 , we have

̂𝜏 ≈ 𝑁 + 2𝑐 2𝑑 , (4.25)
which can be satisfied with 𝑐 = 𝑁 /2, 𝑑 = 𝑛/ τ . Figure 4.10 shows that with this dynamic prior the SBL constrained PINN does not get trapped in a local minimum and learns the underlying data. We hope to exploit multitask learning techniques to optimize this choice in future work.

Experiments

We present three experiments to benchmark our approach. We first study the learning dynamics in-depth on a solution of the Korteweg-de Vries equation, followed by a robustness study of the Burgers equation, and finally show the ability to discover the chaotic Kuramoto-Shivashinsky equation from highly noisy data. Reproducibility details can be found in the appendix of the paper.

Korteweg-de Vries

The Korteweg-de Vries equation describes waves in shallow water and is given by 𝑢 𝑡 =𝑢 𝑥𝑥𝑥 -𝑢𝑢 𝑥 . Figure 4.11a shows the dataset: 2000 samples with 20% noise from a two-soliton solution. We compare our approach with I) Sparse Bayesian Learning with features calculated with numerical differentiation, II) a model discovery algorithm with PINNs, but non-differentiable variable selection called DeepMoD [Both and Kusters, 2021] and III) PDE-find [START_REF] Samuel | Data-driven discovery of partial differential equations[END_REF], a popular model discovery method for PDEs based on SINDy [START_REF] Steven | Discovering governing equations from data by sparse identification of nonlinear dynamical systems[END_REF]. The first two benchmarks also act as an ablation study: method I uses the same regression algorithm but does not use a neural network to interpolate, while method II uses a neural network to interpolate but does not implement differentiable variable selection.

In figure 4.11b and c we show that the differentiable approach recovers the correct equation after approximately 3000 epochs. Contrarily, DeepMoD recovers the wrong equation. Performing the inference 10 times with different seeds shows that the fully-differentiable approach manages to recover the Kortweg-de Vries equation nine times, while DeepMoD recovers the correct equation only twice -worse, it recovers the same wrong equation the other 8 times. Neither PDE-find nor SBL with numerical differentiation is able to discover the Korteweg-de Vries equation from this dataset, even at 0% noise due to the data sparsity.

Burgers

We now explore how robust the SBL-constrained PINN is with respect to noise on a dataset of the Burgers equation, 𝑢 𝑡 =𝜈 𝑢 𝑥𝑥 -𝑢𝑢 𝑥 (figure 4.12a). We add noise varying from 1% to 100% and compare the equation discovered by benchmark method II (DeepMoD, panel b) and our approach (panel c) -the bold orange and blue lines denote 𝑢 𝑥𝑥 and 𝑢𝑢 𝑥 respectively, and the black dashed line their true value. Observe that DeepMoD discovers small additional terms for >50% noise, which become significant when noise >80%. Contrarily, our fully differentiable approach discovers the same equation with nearly the same coefficients across the entire range of noise, with only very small additional terms (𝒪(10 -4 )). Neither PDE-find nor SBL with numerical differentiation is able to find the correct equation on this dataset at 10% noise or higher.

Kuramoto-Shivashinsky

The Kuramoto-Shivashinksy equation describes flame propagation and is given by 𝑢 𝑡 = -𝑢𝑢 𝑥 -𝑢 𝑥𝑥 -𝑢 𝑥𝑥𝑥𝑥 . The fourth order derivative makes it challenging to learn with numerical differentiation-based methods, while its periodic and chaotic nature makes it challenging to learn with neural network based methods [Both et al., 2021b]. We show here that using the SBL-constrained approach we discover the KS-equation from only a small slice of the chaotic data (256 in space, 25 time steps), with 20% additive noise. We use a tanh-activated Figure 4.12: Exploration of robustness of SBL-constrained model for model discovery for the Burgers equation (panel a). We show the discovered equation over a range of noise for DeepMoD (panel b) and the approach presented in this paper (panel c). The bold orange and blue lines denotes 𝑢 𝑥 𝑥 and 𝑢𝑢 𝑥 , and black dashed line their true value. network with 5 layers of 60 neurons each, and the library consists of derivatives up to 5th order and polynomials up to fourth order for a total of thirty terms. Additionally, we precondition the network by training without the constraint for 10k epochs.

Training this dataset to convergence takes significantly longer than previous examples, as the network struggles with the data's periodicity (panel b). After roughly 70k epochs, a clear separation between active and inactive terms is visible in panel c, but it takes another 30k epochs before all inactive terms are completely pruned from the model. Panels d and e show the corresponding posterior and the maximum likelihood estimate of the coefficients using the whole library. Remarkably, the MLE estimate recovers the correct coefficients for the active terms, while the inactive terms are all nearly zero. In other words, the accuracy of the approximation is so high, that least squares identifies the correct equation. Figure 4.13: Recovering the Kuramoto-Shivashinsky equation. We show the chaotic data and a cross section in panels a and b. The periodicity makes this a challenging dataset to learn, requiring 200k iterations to fully converge before it can be recovered (panel c). Panels d and e show that the posterior and MLE of the coefficients yield nearly the same coefficients, indicating that the network was able construct an extremely accurate approximation of the data.

Model discovery with normalizing flows

Consider a set of particles whose movement is described by a micro-scale stochastic process. In the limit of many of such particles, such processes can often be described with a deterministic macro-scale density model, determining the evolution of the density of the particles over time. For example, a biased random walk can be mapped to an advection-diffusion equation. The macro-scale density models are typically more insightful than the corresponding microscopic model, but many (biological) experiments yield single-particle data, rather than densities. Discovering the underlying equation thus requires first reconstructing the density profile from the particles' locations. Classical approaches such as binning or kernel density estimation are either non-differentiable, non-continuous or computationally expensive. Normalizing Flows (NFs) have emerged in recent years as a flexible and powerful method of constructing probability distribution, which is similar to density estimation up to a multiplicative factor. In this section we use physics informed normalizing flows to learn a PDE describing the evolution of the density directly from unlabelled single particle data.

Figure 4.14: Using a tempo-spatial Normalizing Flow constrained by Sparse Bayesian Learning to discover the advection-diffusion equation directly from single particle data. Panel a shows the true density profile, and in panels b, c and d we show the density inferred by binning (blue bars), inferred by NF (red) and the ground truth (black, dashed) at 𝑡 = 0.1, 2.5, 4.5. Note that although the estimate of the density is very good, we see in panel e that we recover two additional terms (bold blue line: 𝑢 𝑥 , bold orange line 𝑢 𝑥𝑥 .

Since the conditional normalizing flow is used to construct the density, a precision denoting the noise level does not exist, and instead we set as prior for 𝛽 (𝑎 = 𝑁 , 𝑏 = 𝑁 ⋅ 10 -5 ). We consider a flow consisting of ten planar transforms [Rezende and Mohamed, 2015] and a hyper-network of two layers with thirty neurons each. The dataset consists of 200 walkers on a biased random walk for 50 steps, corresponding to an advection-diffusion model, with an initial condition consisting of two Gaussians, leading to the density profile shown in figure 4.14a. The two smallest terms in panel e correspond to the advection (bold green line) and diffusion (bold red line) term, but not all terms are pruned. Panels b, c and compare the inferred density (red line) to the true density (dashed black line) and the result obtained by binning. In all three panels the constrained NF is able to infer a fairly accurate density from only 200 walkers. We hypothesise that the extra terms are mainly due to the small deviations, and that properly tuning the prior parameters and using a more expressive transformation would prune the remaining terms completely. Nonetheless, this shows that NF flows can be integrated in this fully differentiable model discovery framework.

Temporal Normalizing Flows

Gert-Jan Both, Remy Kusters. PDF This is another of my favourite papers. Although the execution leaves much to be desired, I believe the core idea is quite simple and elegant and can be of wide use.

In the biological sciences, single particle tracking (SPT) has become the method of choice to investigate many cellular processes. The obtained trajectories can be either analyzed as random walks, or from the perspective of walker densities. This density model is deterministic and typically more insightful than the stochastic random walk. However, in experiments the position of single particles, and from these positions a time-dependent density must be reconstructed. Classical methods such as binning or kernel density estimation are computationally expensive, non-differentiable 5 , but above all unable to explicit model a time-dependent density. Rather, they construct snapshots of the density at given times. In this paper we introduce conditional normalizing flows (CNFs)6 , a flexible and differentiable method which models the entire continuous spatio-temporal density.

We construct a Conditional Normalizing Flow by constraining the Jacobian of the transformation. Consider a two dimensional normalizing flow in spatial coordinate 𝑥 and time 𝑡:

𝑝(𝑥, 𝑡) = 𝜋(𝑧,

Results

We now demonstrate tNFs on three datasets:

• A multi-scale toy problem to show tNFs can accommodate different length scales in a single distribution; • A dataset of Brownian motion to show how tNFs enhance density estimation for sparse datasets; • A dataset of chemotactic walkers to show that tNFs can correctly estimate a multi-modal, non-Gaussian density.

Multi-scale density estimation

A key problem in density estimation is inferring an accurate distribution when vastly different length scales are present within a single dataset. Classical approaches such as binning and KDE require a single characteristic length scale, prohibiting an accurate estimate of a multi-scale distribution. We now show that normalizing flows, and by extension tNFs, are capable of accurately inferring such a distribution.

We build an artificial distribution consisting of three normal distributions with standard deviations 𝜎 of 0.01, 0.1 and 1.0 (thus spanning three orders of magnitude) and respective weights 0.013, 0.13 and 0.85. Figure 4.15 shows the inferred distribution from 5000 samples for the NF and the KDE with Scott's rule determining the lengthscale. Observe that, as expected, the KDE is unable to accommodate the different scales and that due the different weighting of each peak, the widest is dominating the lengthscale estimation. Contrarily, the NF provide an accurate density estimate for all lengthscales present in the problem, independent of their weights. We stress that tNF do not require any prior information of the different lengthscales present in the data-set.

Brownian motion

Brownian motion is the most basic and ubiquitous random walk and thus an ideal test case to assess the performance of tNFs, comparing them to time independent NFs and classical binning. We generate a single trajectory for a Brownian random walker by the recursive relation, ⃗𝑥 𝑛+1 = 𝒩( ⃗ 𝑥 𝑛 , √ 2𝐷Δ𝑡). Here 𝑛 is the step number with 𝑥 0 the initial position, 𝐷 the diffusive coefficient and Δ𝑡 the time step. In the limit of an infinite number of walkers, the walker density 𝑐 is described by the diffusion equation, 𝜕 𝑡 𝑐=𝐷∇ 2 𝑐.

Our dataset consists of 𝑀 = 500 walkers with 𝐷 = 2.0, with snapshots being taken every Δ𝑡 = 0.1 for 𝑁 = 100 frames. The initial positions were sampled from a Gaussian centered at 𝑥 = 1.5 with width 𝜎 = 0.5; in this case, the diffusion equation can be solved exactly and the solution behaves as a spreading Gaussian in time. We show the estimated density at 𝑡 = 0.75 and 𝑡 = 4.25 in figure 4.16 (a) and (b) for the tNF, the time independent NF and binning. The tNF provides a significantly better density estimate than the time independent NF, illustrated by the difference in ℓ 2 error; 2.6 ⋅ 10 -5 for the tNF and 7.4 ⋅ 10 -5 for the NF, averaged over frame 15 and 85.

Normalizing flows are based on neural networks and hence prone to overfitting. We analyze the effect of overfitting in Appendix I of the paper and show that NFs overfit more strongly than tNFs and perform worse in terms of the ℓ 2 error. We mainly attribute this improvement to the temporal correlations in the dataset, which suppresses the natural frame-to-frame variations in the density estimate. Nonetheless, tNFs are not immune to overfitting and we speculate performance could be enhanced by applying techniques such as early stopping. For the diffusion equation the true mapping can be trivially derived. We compare it to the learned mapping in figure 4.16c. It shows perfect agreement at 𝑡 = 4.25, but deviates from the true curve for 𝑥 < -5 at 𝑡 = 0.75. As can be seen in figure 4.16a, no samples were present in this domain, explaining the deviance. Nonetheless, it implies that the network does not generalize well outside the sampling domain. We speculate that techniques such as batch normalization or a different architecture for the network (a recurrent network, for example) might further improve performance.

Chemotaxis

The Brownian motion presented in the previous section was a linear problem with a uni-modal, Gaussian solution. We now apply tNFs to so-called chemotactic walkers, a non-linear problem with a multi-modal solution. Bacteria and other micro-organisms sense gradients of chemicals throughout their environment and use this to guide their motion towards a food source. This effect is known as chemotaxis and is typically modelled by a random walker with a superimposed drift; ⃗𝑥 𝑛+1 = 𝒩( ⃗ 𝑥 𝑛 + 𝜒∇𝑝( ⃗ 𝑥 𝑛 )𝑑𝑡, √ 2𝐷𝑑𝑡), where 𝑝 is the chemical density and 𝜒 is the chemotactic sensitivity, which controls the interaction between the chemical and the bacteria. In the infinite walker limit, the walker and chemical density are given by the Keller-Segel model: 𝜕 𝑡 𝑐 = ∇ ⋅ (𝐷 𝑐 ∇𝑐 -𝜒𝑐∇𝑝) and 𝜕 𝑡 𝑝=𝐷 𝑝 ∇ 2 𝑝-𝐾𝑝. Here 𝐷 𝑐 and 𝐷 𝑝 are the diffusion coefficients of the bacteria and the chemical respectively and a decay set by 𝐾 has been added to the chemical density.

Our dataset consisted of 𝑀 = 500 walkers with 𝐷 = 0.5 and we sampled the initial position from a Gaussian centred at 𝑥 = -2.5. The food source was modelled by a Gaussian with diffusion coefficient 𝐷 = 0.25, centred at 𝑥 = 2.5; the walkers will thus drift towards food source over time. Figure 4.17 shows a comparison of the time independent NF, tNF and the binning method. In figure 4.17(a) and (b) we find that the tNF leads to a significantly more accurate density estimation, illustrated by the difference in ℓ 2 error (1.45⋅10 -4 for the NF versus 1.8⋅10 -5 for the tNF, averaged over 𝑡 = 6.0 and 19.0). The tNF captures the multi-modal distribution at 𝑡 = 6.0 excellently, without overfitting, contrarily to the time independent NF. The mapping, as shown in figure 4.17c, is non-linear, in contrast to the mapping obtained for the Brownian motion.

Model discovery in the sparse sampling regime

Gert-Jan Both, Georges Tod, Remy Kusters. PDF This paper started with a curious observation: between our first and second paper we noticed a curious regression in performance. We eventually traced it back to a change in sampling strategy: when the data was randomly sampled instead of on a grid, DeepMoD was able to recover the underlying equation with much less data.

DeepMoD can handle much sparser data than other classical interpolation techniques. In this paper we compare our method in-depth to spline interpolation, how it fails as data becomes sparser and link it to the data's characteristic lengthscale.

Results

Synthetic data -Burgers equation

We consider a synthetic data set of the Burgers equation 𝑢 𝑡 =𝜈𝑢 𝑥𝑥 -𝑢𝑢 𝑥 , with a delta peak initial condition 𝑢(𝑥, 𝑡 = 0) = 𝐴𝛿(𝑥) and domain 𝑡∈ [0.1, 1.1], 𝑥 ∈ [-3, 4]. This problem can be solved analytically to yield a solution dependent on a dimensionless coordinate 𝑧 = 𝑥/ √ 4𝜈𝑡. We recognize the denominator as a time-dependent length scale: a Burgers data set sampled with spacing Δ𝑥 thus has a time-dependent dimensionless spacing Δ𝑧(𝑡). We are interested in the smallest characteristic length scale, which for this data set is 𝑙 𝑐 = √4𝜈𝑡 0 , where 𝑡 0 = 0.1 is the initial time of the data set. We set 𝐴=1and 𝜈 = 0.25, giving 𝑙 𝑐 = √ 0.1 ≈ 0.3.

Splines do not scale effectively beyond a single dimension, making it hard to fairly compare across both the spatial and temporal dimensions. We thus study the effect of spacing only along the spatial axis and minimize the effect of discretization along the temporal axis by densely sampling 100 frames, i.e. Δ𝑡 = 0.01. Along the spatial axis we vary the number of samples between 4 and 40, equivalent to 0.5 < 4.5.1.1.1 Off-grid sampling Whereas higher-order splines are constrained to interpolating along a single dimension, DeepMoD uses a neural network to interpolate along both the spatial and temporal axis. This releases us from the constraint of on-grid sampling, and we exploit this by constructing an alternative sampling method. We observe that for a given number of samples 𝑛, DeepMoD is able to interpolate much more accurately if these samples are randomly drawn from the sampling domain.

We show in figure 4.18b and c (Red) that the relative error 𝜖 in the sparse regime, can be as much as two orders of magnitude lower compared to the grid-sampled results at the same number of samples. We hypothesize that this is due to the spatio-temporal interpolation of the network. By interpolating along both axes, each sample effectively covers its surrounding area, and by randomly sampling we cover more of the spatial sampling domain. Contrarily, sampling on a grid leaves large areas uncovered; we are effectively sampling at a much lower resolution than when using random sampling.

To test whether or not this improvement is intrinsically linked to the randomness of sampling, we also construct an alternative sampling method called shifted-grid sampling. Given a sampling grid with sensor distance Δ𝑥, shifted-grid sampling translates this grid a distance Δ𝑎 every frame, leading to an effective sample distance of Δ𝑎 ≪ Δ𝑥. This strategy, similarly as random sampling varies the sensor position over time, but does so in a deterministic and grid-based way. We show this sampling strategy in figure 4.18a, while panels b and c confirm our hypothesis; shifted grid sampling (Blue) performs similarly to random sampling. Shifted-grid sampling relies on a densely sampled temporal axis 'compensating' for the sparsely sampled spatial axes. This makes off-grid sampling beneficial when either the time or space axis, but not both, can be sampled with a high resolution.

In the experimental section we show that if both axes are sparsely sampled, we do not see a strong improvement over grid sampling.

Experimental data -2D Advection-Diffusion

In an electrophoresis experiment, a charged dye is pipetted in a gel over which a spatially uniform electric field is applied (see Figure 4.19a)). The dye passively diffuses in the gel and is advected by the applied electric field, giving rise to an advection-diffusion equation with advection in one direction: 𝑢 𝑡 = 𝐷(𝑢 𝑥𝑥 +𝑢 𝑦𝑦 )+𝑣𝑢 𝑦 . Both et al. [2019] showed that DeepMoD could discover the correct underlying equation from the full data-set (size 120 x 150 pixels and 25 frames). Here, we study how much we can sub-sample this data and still discover the advection-diffusion equation.

In figure 4.19 c) and d) we perform grid based as well as a random subsampling of the data. The neural network-based method discovers the advection-diffusion equation on as few as 6 x 8 spatial sampling points with 13 time-points, or with 20 x 25 pixels on only 3 time-points. The minimum number of required samples is similar for grid and random sampling, confirming that when both axes are poorly sampled, there is no benefit to sample randomly.

The smallest characteristic length scale in the problem is the width of the dye at 𝑡=𝑡 0 , which we estimate as 𝑙 𝑐 ≈ 10 pixels. For the data presented in figure 4.19c) and 4.19d), at a resolution below 10 × 13 sensors classical approaches would inherently fail, even if no noise was present in the data set. This is indeed what we observe: using a finite difference-based library we were unable to recover the advection-diffusion equation, even after denoising with SVD. However, we do observe that below this value we occasionally discover the correct PDE, despite the experimental noise present in the data. The use of a neural network and random sampling lead to non-deterministic behaviour: the neural network training dynamics depend on its initialization and two randomly sampled datasets of similar size might not lead to similar results. In practice this leads to a 'soft' decision boundary, where a fraction of a set of runs with different initialization and datasets fail.

Noiseless synthetic data set

To further confirm our results from the previous section, we simulate the experiment by solving the 2D advection-diffusion with a Gaussian initial condition and experimentally determined parameters (𝐷 = 0.025 and and 𝑣 = (0, 0.2). Figure 4.19e) shows the selected terms and their magnitude as functions of the applied noise levels for a highly subsampled data-set (grid sampling, spatial resolution of 11x11 and temporal resolution 14). The correct AD equation is recovered up to noise levels of 100% (See figure 4.19e), confirming the noise robustness of the NN-based model discovery. In panel f) we compare the deep-learning based model discovery using grid and random sampling with classical methods for various noise levels and sensor spacing with a fixed temporal resolution of 81 frames (Data for the FD was pre-processed with a SVD filter, see SI for further details). We confirm that, similarly to the Burgers example of the previous section, the correct underlying PDE is discovered even below the smallest characteristic length-scale in the problem (indicated by a red dashed line in figure 4.19f). This figure confirms our three main conclusions: 1) In the noiseless limit, classical approaches are only slightly less performing than NN-based model discovery for grid sampling. 2) Increasing the noise level dramatically impacts classical model discovery while barely impacting NN-based model discovery and 3) random sampling over space considerably improves performance, allowing model discovery with roughly 4-8 times fewer sample points for this particular data-set (depending on the noise level).

Experimental data -Cable equation

Applying a constant voltage to a RC-circuit with longitudinal resistance (see figure 4.20 a) result in time-dependent voltage increase throughout the circuit due to the charging of the capacitors. This rise is modeled by the cable equation, which is essentially a reaction-diffusion equation $ u_t = u_{xx}/ (R_l C) + u/(R_m C)$ with 𝐶 the capacitance, 𝑅 𝑙 the longitudinal resistance and 𝑅 𝑚 the parallel resistance of the circuit. The discrete nature of the experiment automatically gives Δ𝑥 = 𝑂(𝑙 𝑐 ). We consider an extreme case where we only have seven sensors throughout the circuit (i.e. spatial axis), but take 2500 samples along the time axis. Figure 4.20b shows the measured voltage at these seven elements. Initially, all the capacitors are uncharged and we observe a sharp voltage increase at the first element. As the capacitors charge, this further propagates through the circuit, charging the capacitors and resulting in the curves shown in the figure. We apply both a classical approach with the library generated with splines and DeepMoD to a varying amount of elements. Chapter 5

Conclusion

Despite the massive progress in the last few years, model discover of (partial) differential equations is still a field in its infancy. At the start of this thesis in early 2019, both the sparse regression approach for PDEs and PINNs were novel; to our knowledge, no one had attempted to perform model discovery on noisy experimental data. Considering the experiments on synthetic data that would have been a futile attempt -they all showed model discovery required densely sampled datasets with very little noise (<10%, but more often <2% noise). With this we set out to construct methods able to handle noisy and sparse datasets originating from experiments. The line of work we have developed with DeepMoD makes large strides towards this goal. In each paper we show that our methods are able to handle >50% noise and work with an order of magnitude less data than classical methods on increasingly challenging (synthetic) datasets such as the Kuramoto-Shivashinsky equation. DeepMoD also recovered the underlying equation when applied to data generated by simple experiments, validating our approach.

A second, more implicit goal was to create accessible methods. Complex methods are less likely to be adapted than simple ones, and if the goal is to make model discovery a valuable addition to scientists' toolboxes, complexity should be strongly limited. Ideally, scientists without a background in numerical methods should be able to understand and apply our methods. Deep learning is an ideal vehicle for this, as many scientists these days have at least a basic understanding. Automatic differentiation, the key but complex mechanism at the core of DL, is often abstracted away to frameworks such as Pytorch, making it possible to construct powerful approaches without in-depth numerical knowledge. This leads to the peculiar situation where a neural network-based approach is regarded as more accesible than splines. Indeed, our work shows that an integrated combination of simple features (a basic MLP, a simple Lasso and a constraint) can easily and strongly outperform much more complex traditional approaches -a 69 testament to the power of differentiable programming.

Our work also carries strong implications for future work and for model discovery as a field. First and foremost, it shows that the limiting factor is the accuracy of the features, implying that accurately modelling and denoising data is just as important as the sparse regression. Perhaps paradoxically then, significant progress can be made by focussing on data-driven modelling. In all cases, neural networks (especially PINNs) should feature prominently; the benefits of automatic differentiation and excellent interand extrapolation only become more pronounced in high-dimensional data.

That is not to say that non-DL-based approaches should be neglected. A main thread in our work has been to show how 'classical' methods such as sparse and Bayesian regression can be integrated in DL approaches, and most progress can be made by synthesizing these approaches. The Bayesian regression and model selection literature is especially rich, and combining DL-based modelling with these methods should prove fruitful.

Taken together, our work strongly establishes the argument for physicsconstrained, neural network-based surrogates for model discovery of PDEs on experimental data.

Challenges and questions unanswered

Model discovery is a young and exciting field, and as with all young fields, many limitations, challenges and questions remain. As such I list these in no particular order below -I've mentioned some of them before, others might have occured to you while reading this thesis and a few of them are of a more philosophical nature.

• Perhaps the single biggest barrier to apply model discovery on novel, experimental data is the lack of methods able to handle data with spatially and temporally varying coefficient fields. Initial work [Rudy et al., 2018, Chen and[START_REF] Chen | Robust data-driven discovery of partial differential equations with time-dependent coefficients[END_REF] has focussed on leveraging group sparsity, but so far this approach works only for a single varying dimension -sufficient for temporal dependence, but not for spatial dependence (2D at least). Scaling this approach to spatio-temporal dependence yields a compressed sensing problem: every sample becomes a separate regression problem (𝑛=1 ) with many features (𝑝≫𝑛 ). An alternative approach would be to extend on DeepMoD, and use a neural network to model the field. This implicitly encodes a smoothness bias into the fields and initial results have been promising. However, to perform actual model discovery would require some criterion for deciding when a field can be considered inactive.

• Contrarily, perhaps one of the biggest opportunities for model discovery would be the synthesis of data from multiple sources. Rarely does a single experiment tell the whole story: we need multiple experiments to capture all of a systems dynamics. Initial work using group sparsity [START_REF] De Silva | Discovery of Physics from Data: Universal Laws and Discrepancy Models[END_REF], Tod et al., 2021] only scratches the surface of what is possible. How can data from multiple experiments with different length and time scales be integrated? Initial work on combining timescales [Champion et al., 2019a] on ODEs is promising, but has not been expanded to PDEs, nor to (multiple) lengthscales. Even more interesting would be to combine data from different experimental setups and instruments, also known as multimodal fusion.

• Related to the data synthesis challenge is that of sampling, active learning and experimental design. Sampling has been studied in-depth from the perspective of signal reconstruction [START_REF] Brunton | Optimal Sensor Placement and Enhanced Sparsity for Classification[END_REF], but not from the viewpoint of model discovery -a subtly different problem.

Is there an optimal sampling strategy for model discovery, and if so, what would it be? Having more knowledge about this could form the basis for an active learning approach, where the model itself suggests where to sample to optimize model discovery. The most ambitious form of this would be experiment design, with an agent suggesting novel experiments designed to aid model discovery.

• All of the neural networks used in our work were simple MLPs -no normalization, batching, or attention mechanisms. How much do the implicit biases of the function approximator impact the discovered equation? Various works show that incorporating underlying symmetries and invariances in the network improves performance [START_REF] Cranmer | Lagrangian Neural Networks[END_REF][START_REF] Greydanus | Hamiltonian Neural Networks[END_REF][START_REF] Finzi | A Practical Method for Constructing Equivariant Multilayer Perceptrons for Arbitrary Matrix Groups[END_REF], but in the context of model discovery these would have to be discovered from data.

• Model discovery is usually applied directly on the observable data; we observe some quantity 𝑢 and also wish to find a model for 𝑢. In some cases the observable is not fully observed -the model is described by an 𝑛-dimensional state, and we only observe some of those dimensions -, or not at all what we wish to model -we observe single particles but want to model a density-. This significantly complicates model discovery, as it requires first constructing or inferring these 'latent' dynamics. This is likely to yield biased estimates and, in our experience, model discovery deals poorly with bias in data. Making model discovery more robust at dealing with bias is a sorely understudied subject.

• A reflection on in what kind of settings model discovery can be useful is needed. Most papers (including our own) use model discovery on systems where the dynamics are already known, and relatively simple.

In what kind of data could model discovery be truly useful? Can the complex, high dimensional data associated with the modern world by described by relatively simple, physics-inspired PDEs? Is model discovery only useful to discover effective models, or would it be able to discover new, more fundamental dynamics?

  Figure 1.1: One-dimensional solution of the Kuramoto-Shivashinsky equation.

Figure 1 . 2 :

 12 Figure 1.2: Example of binary expression tree. Figure taken from https: //en.wikipedia.org/wiki/Binary_expression_tree.

Figure 2 . 3 :

 23 Figure 2.3: Penalty functions of ℓ 0 ,ℓ 1 and ℓ 2 regularization.

Figure 2 . 4 :

 24 Figure 2.4: Comparison of Lasso and Ridge losses. Figure from https:// www.astroml.org/book_figures/chapter8/fig_lasso_ridge.html

Figure 2

 2 Figure 2.6: Lasso path. Figure from https://scikit-learn.org/0.18/auto_ examples/linear_model/plot_lasso_lars.html

Figure 2

 2 Figure 2.7: Stability paths. Figure from https://thuijskens.github.io/2018/ 07/25/stability-selection Using the stability path, the final support is determined by selecting all components with a selection probability above a given threshold 𝜋 thr ,

Figure 3 . 1 :

 31 Figure 3.1: Numerical differentiation of noisy and sparse data. Blue are the samples, black dashed line the ground truth.

Figure 4 . 1 :

 41 Figure 4.1: Schematic representation of the first DeepMoD iteration.

  Figure 4.2: a) Cost functions and b) coefficient values as function of the number of epochs for the Burgers' dataset consisting of 2000 points and 10% white noise. Initially, the neural network optimizes the MSE and only after the MSE is converged the coefficient vector is optimized by the network. c) Velocity field 𝑢 for 𝑡=5obtained after training with (no overfitting) and without (overfitting) the regression regularization ℒ 𝑅𝑒𝑔 . d) The values in the grid indicate the accuracy of the algorithm tested on the Burgers' equation, defined as the mean relative error over the coefficients, as function of the sample size of the data set and level of noise. The coloring represents the fraction of correct runs (Red indicates that in none of the five iterations the correct PDE is discovered).

Figure 4 . 3 :

 43 Figure 4.3: Ground truth, Noisy and Inferred data at 𝑡 = 0.8 for the Burgers' equation with a shock wave (10% noise and 2000 sample points).

Figure 4 . 4 :

 44 Figure 4.4: Ground truth, noisy and reconstructed solutions for the density of cells, 𝑢 (top row) and the density of secreted chemicals 𝑤 (bottom row) in the Keller Segel model for 5% white noise and 10000 samples.

Figure 4

 4 Figure 4.5: a) Time-series images of the electrophoresis essay. The two red boxes indicate the analysed region. b) Region indicated in the solid red box of (a) showing the density of the dye at three different time-frames (in pixels). c) Scaled coefficients values of all the candidate for a single run. The pink region indicates the terms with scaled coefficient |𝜉| < 0.1.

Funct

  Figure 4.6: Schematic representation of the modular DeepMoD approach.

Figure 4

 4 Figure 4.7: A) Loss and B) mean absolute error of the coefficients obtained with the gradient descent and the least squares constraint as a function of the number of epochs. Results have been averaged over twenty runs and shaded area denotes the standard deviation.

Figure 4

 4 Figure 4.9: A) Solution of the KS equation. Lower panel shows the cross section at the last time point: 𝑡 = 44. B) MSE as function of the number of epochs for both the tanh-based and SIREN NN. C) Coefficients as function of number of epochs for the SIREN and D) the tanh-based NN. The bold curves in panel C and D are the terms in the KS equation components; green: 𝑢𝑢 𝑥 :, blue: 𝑢 𝑥𝑥 and orange: 𝑢 𝑥𝑥𝑥𝑥 . Only SIREN is able to discover the correct equation.

  𝑧=𝑓 𝐾 ∘…∘𝑓 0 (𝑥) ≡ 𝑔 𝜃 (𝑥) log 𝑝(𝑥) = log 𝜋(𝑧)

Figure 4 .

 4 Figure 4.11: Comparison of a differentiable SBL-constrained model and an non-differentiable OLS-constrained model on a Korteweg-de Vries dataset (panel a) with a library consisting of 4th order derivatives and 3rd order polynomials, for a total of 20 candidate features. In panel b and c we respectively plot the inferred prior  and the posterior coefficients 𝜇. In panel d we show the non-differentiable DeePyMod approach. In panels b and c we see that the correct equation (bold blue line: 𝑢 𝑥𝑥 , bold orange line: 𝑢𝑢 𝑥 ) is discovered early on, while the non-differentiable model (panel d) selects the wrong terms.

Figure 4 .

 4 Figure 4.15: Comparison of density estimation with different scales. 𝑁= 5000 samples were taken from a density consisting of three Gaussians with widths 0.01, 0.1, and 1, at respective locations 0.1, 2 and 10, and weights 0.013, 0.13 and 0.85. The KDE used a Gaussian kernel with the kernel width set by Scott's rule.

Figure 4 .

 4 Figure 4.16: Results of inferring a Brownian walker density from 𝑀 = 500 walkers with 𝐷=2and where a snapshot was taken every Δ𝑡 = 0.1 for 𝑁 = 100 frames. The initial position was sampled from a Gaussian placed at 𝑥 = 1.5 with width 𝜎 = 0.5. Panel a and b compare the inferred density by binning, time independent normalizing flow (NF) and time dependent normalizing flow (tNF) at two different times. Panel c compares the learned mapping with the true mapping.

Figure 4 .

 4 Figure 4.17: Results of inferring a chemotactic walker density from 𝑀 = 500 walkers with 𝐷 = 0.5 and chemotactic sensitivity 𝜒 = 10, where a snapshot was taken every Δ𝑡 = 0.1 for 𝑁 = 100 frames. The initial position was sampled from a Gaussian placed at 𝑥 = -2.5 with width √ 0.5, while the food source is a Gaussian located at 𝑥 = 2.5 with width √ 0.5, diffusing with 𝐷 = 0.25 and decaying at a rate of 0.05. Panel a and b compare the inferred density by binning, time independent normalizing flow (NF) and time dependent normalizing flow (tNF) at two different times. Panel c compares the learned mapping with the true mapping.

Figure 4 .

 4 Figure 4.18: a) The three sampling strategies considered. b) and c) Error in the function library as function of the distance between the sensors Δ𝑥, normalized with 𝑙 𝑐 = √4𝜈𝑡 0 , for b) noise-less data and c) 20% of additive noise. The yellow symbols correspond to a spline interpolation and the green, blue and red correspond to the NN-based model discovery with various sampling strategies. The circles indicate that model discovery was successful while the triangles indicate that the incorrect model was discovered. The horizontal dashed line indicates the smallest characteristic length-scale in the problem: Δ𝑥/𝑙 𝑐 =1.

Figure 4 .

 4 Figure 4.19: a) Experimental setup of the gel electrophoresis. b) Three time frames of the density with a spatial resolution of 20x25. c) and d) Success diagram for the experimental data indicating correct model discovery (Yellow indicates the correct AD equation 𝑢 𝑡 = 𝐷(𝑢 𝑥𝑥 +𝑢 𝑦𝑦 )+𝑣𝑢 𝑦 is discovered) as function of the spatial and temporal resolution for c) grid sampling and d) random sampling. e) Obtained mask and coefficients (𝐷 = 0.025 and and 𝑣 = (0, 0.2)) for the artificial data-set as function of the noise level (11x11 spatial resolution). Hereby, yellow indicates the terms selected by the algorithm and the red dashed box the terms that are expected in the PDE. f) Success diagrams for various levels of additive noise, comparing the result of DeepMoD with a grid and random sampling strategy and the classical LassoCV algorithm on a Finite Difference (FD)-based library (after SVD filtering of the different frames).

Figure 4

 4 Figure 4.20: a) Schematic overview of the electronic setup to generate the cable equation. b) The voltage drop 𝑢 as function of time for various positions along the circuit for a circuit with 7 elements. The mask obtained for c) NN-based and d) cross validated Lasso with spline based library model discovery (Yellow indicates the term was recovered by the algorithm). The red boxes indicate the two expected terms in the equation.

  𝑁 is the time derivative and 𝜉∈ℛ 𝑀 the coefficient vector describing the equation. If a component 𝜉 𝑗 is zero, the corresponding feature Θ 𝑗 is not a part of the equation, and vice versa. For example, if Θ= [1𝑢𝑢 𝑥 𝑢𝑢 𝑥 ],

	Θ= ⎡ ⎢ ⎣ corresponds to the equation || | 1𝑢𝑢 𝑥 …𝑢 2 𝑢 𝑥𝑥 | | || | | | 𝜕𝑢 𝜕𝑡 = Θ𝜉 𝜕𝑢 𝜕𝑡 = 1 2 𝑢-2 ⎤ ⎥ ⎦ 𝜕 2 𝑢 𝜕𝑥 2 . With eq. (2.5), instead of comparing all models which can be build with the (2.4) (2.5) 𝑇 (2.6) (2.7) candidate features, now we only need to find a single vector 𝜉 by minimizing 5 , 𝜕𝑡 ∈ℛ 𝜉=[ 0.0 0.5 0.0 -2.0] where 𝜕𝑢 ξ=min 𝜉 2 ‖𝑢 𝑡 -Θ𝜉‖ 2

3

, 1. Polynomials of the observable 𝑢 up to order 𝑝: 1, 𝑢, 𝑢 2 ,…. 2. Derivatives up to order 𝑞: 𝑢 𝑥 ,𝑢 𝑥𝑥 ,𝑢 𝑥𝑥𝑥 ,… 3. Combinations of these two features: 𝑢⋅𝑢 𝑥 ,𝑢 2 ⋅𝑢 𝑥𝑥𝑥 ,… For most known equations, 𝑝≤3and 𝑞≤4 , typically yielding 20-30 candidate features. Even at this size, a combinatorial search is computationally too expensive, as it involves solving the PDE. However, by exploiting property 1, i.e. that the unknown equation is a linear combination of the candidate features, model discovery can be expressed as a regression problem 4 . Figure 2.2: Solution of the Cahn-Hilliard equation. Image taken from https: //www.hindawi.com/journals/mpe/2016/9532608/

2.2 Regression for model discovery

Expressing model discovery as a regression problem starts by combining all 𝑀 candidate features into a single library matrix Θ, i.e. each column of Θ∈ℛ 𝑁×𝑀 contains a single feature. Assuming the equation is a linear combination of these features, the unknown equation can be written as,

  .16) which cannot be calculated when using a normal neural network to model 𝑝(𝑥); the integral is intractable as it runs from +∞ to -∞. While the integral can be approximated, this property is never guaranteed and will always be an approximation. Other approaches to density estimation rely for example on a mixture of Gaussians, which guarantee the normalization but are neither flexible enough to model any arbitrary density, but also too computationally expensive. Normalizing Flows[START_REF] Jimenez | Variational Inference with Normalizing Flows[END_REF] have become the prime method for modelling densities, being flexible, computationally cheap, and guaranteed to yield a proper probability distribution.Various types of flows 𝑓 have been proposed, from the relatively simple planar flows of[START_REF] Jimenez | Variational Inference with Normalizing Flows[END_REF] to more complicated ones such as the Sylvester normalizing flows[START_REF] Van Den | Sylvester Normalizing Flows for Variational Inference[END_REF] or monotonic neural networks[START_REF] Wehenkel | Unconstrained Monotonic Neural Networks[END_REF]. Recent work has focused on flows on different geometries and satisfying certain equivariances.In the previous two chapters I have introduced two, largely separate ideas: using sparse regression to perform model discovery of differential equations, and using surrogates, specifically PINNs, to accurately interpolate data, consistent with a given equation. The central premise of my work is that these two ideas can be combined, yielding a robust model discovery algorithm. Integrating sparse regression in PINNs is easily achieved. Considering again the loss of a PINN with known equation 𝑢 𝑡 =𝑋𝑤,

	Chapter 4
	DeepMoD
	ℒ=∑ 𝑖	‖𝑢 𝑖 -𝑔 𝜃 (𝑥 𝑖 ,𝑡 𝑖 )‖ 2 2 +∑ 𝑖	‖(𝑢 𝑡 ) 𝑖 -𝑋 𝑖 𝑤‖	2 2 .	(4.1)
	we simply replace this equation by the candidate library Θ and the sparse vector 𝜉,
	ℒ=∑ 𝑖	𝑥=𝑓 𝐾 ∘𝑓 𝐾-1 ∘…𝑓 0 (𝑧) ‖𝑢 𝑖 -𝑔 𝜃 (𝑥 𝑖 ,𝑡 𝑖 )‖ 2 2 +∑ 𝑖 ‖(𝑢 𝑡 ) 𝑖 -Θ 𝑖 𝜉‖ 2 2 .	(4.2)	(3.18)
	To train a normalizing flow, we typically minimize the negative log likelihood over the samples,
		ℒ=-	𝑁 ∑	log 𝑝(𝑥).		(3.19)
			𝑖=1

NFs are based on the change-of-variable theorem for probability distributions,

𝑝(𝑥) = 𝜋(𝑧) | 𝑑𝑧 𝑑𝑥 | , 𝑥 = 𝑓(𝑧).

(3.17)

Essentially, this allows us to express the unknown distribution 𝑝(𝑥) to be expressed in terms of a known distribution 𝜋(𝑧) and a transformation 𝑥 = 𝑓(𝑧).

NFs can thus learn any distribution by learning the transformation 𝑓. The only requirement is for 𝑓 to be a bijective function -in other words, it needs to one-to-one and invertible. The expressive power of NFs depends on the transformation 𝑓, but constructing arbitrary invertible functions is extremely challenging 11 . Normalizing flows work around this by constructing a series of simple transformations 12 ,

  is a Bayesian approach to regression yielding sparse results. SBL defines a hierarchical model, starting with a Gaussian likelihood with noise precision 𝛽≡𝜎 -2 , and a zero-mean Gaussian with precision 𝛼 𝑗 on each component 𝜉 𝑗 as prior, 𝜕 𝑡 ̂𝑢 ∈ ℛ 𝑁 , Θ∈ℛ 𝑁×𝑀 , 𝜉∈ℛ 𝑀 , and we have defined 𝐴≡diag(𝛼). The posterior distribution of 𝜉 is a Gaussian with mean 𝜇 and covariance Σ, Many of the terms in 𝐴 will go to infinity when optimised, and correspondingly the prior for term 𝑗 becomes a delta peak. We are thus certain that that specific term is inactive and can be pruned from the model. This makes SBL a very suitable choice for model discovery, as it gives a rigorous criterion for deciding whether a term is active or not. Additionally it defines hyper-priors over 𝛼 and 𝛽,The inference of 𝐴 and 𝛽 cannot be performed exactly, and SBL uses type-II maximum likelihood to find the most likely values of  and β by minimising the negative log marginal likelihood 3

			𝑝(𝛼) =	𝑀 ∏ 𝑗=1	Γ(𝛼 𝑗 ; 𝑎, 𝑏)	(4.14)
			𝑝(𝛽) = Γ(𝛽; 𝑐, 𝑑)
	ℒ SBL (𝐴, 𝛽) =	1 2	[𝛽 ‖𝑢 𝑡 -Θ𝜇‖
			𝑝(𝜕 𝑡 û ; Θ, 𝜉, 𝛽) =	𝑁 ∏ 𝑖=1	𝒩(𝜕 𝑡 ̂𝑢 𝑖 ;Θ 𝑖 𝜉, 𝛽 -1 ),	(4.11)
					𝑝(𝜉; 𝐴) =	𝑀 ∏ 𝑗=1	𝒩(𝜉 𝑗 ; 0, 𝛼 -1 𝑗 ),	(4.12)
	with Σ = (𝛽Θ 𝑇 Θ + 𝐴) -1 𝜇 = 𝛽ΣΘ 𝑇 𝜕 𝑡 û .	(4.13)

2 +𝜇 𝑇 𝐴𝜇 -log|Σ| -log|𝐴| -𝑁 log 𝛽]-𝑀 ∑ 𝑗=1 (𝑎 log 𝛼 𝑗 -𝑏𝛼 𝑗 )-𝑐 log 𝛽+𝑑𝛽,

(4.16)

  𝑧 and 𝜏 are thus some sort of latent coordinates. Since only the spatial axis represents a valid probability distribution, i.e. ∫ 𝑝(𝑥, 𝑡)𝑑𝑥 = 1, but ∫ 𝑝(𝑥, 𝑡)𝑑𝑡 = 𝑓(𝑥), only 𝑥 can be validly transformed to 𝑧. We enforce this by simply requiring 𝑡=𝜏 , i.e. the latent time 𝜏 is the real time 𝑡.

	transformation then reduces to				The
	𝑝(𝑥, 𝑡) = 𝜋(𝑧, 𝑡) |	𝑑𝑧(𝑥, 𝑡) 𝑑𝑥	| , 𝑧 = 𝑓(𝑥, 𝑡)	(4.27)
	i.e. the bijective transform depends on 𝑥 and 𝑡, but the coordinate 𝑡 is not transformed.
	𝜏 ) |	𝑑𝑧 𝑑𝑥	𝑑𝜏 𝑑𝑡	-	𝑑𝜏 𝑑𝑥	𝑑𝑧 𝑑𝑡	| , [𝑧, 𝜏 ] = 𝑓(𝑥, 𝑡).	(4.26)

  𝑖 sums the derivatives and 𝑗 runs over the frames. The derivatives are normalised every frame by the 𝑙 2 norm of the ground truth to ensure 𝜖 is independent of the magnitude of 𝑢. 𝜖 does not take into account the nature of noise (e.g. if it is correlated and non-gaussian), nor if the correct equation is discovered. However, taken together with a success metric (i.e if the right equation was discovered), it does serve as a useful proxy to the quality of the interpolation.Figure4.18b) shows 𝜖 as a function of the relative spacing Δ𝑥/𝑙 𝑐 and whether the correct equation was discovered. The error when using splines (yellow) increases with Δ𝑥 and, as expected, we are unable to discover the correct equation for Δ𝑥 > 0.8𝑙 𝑐 (dots indicate the correct equation is discovered and triangles indicates it failed to do so). Considering the NN-based DeepMoD method, sampled on a grid (green), we observe that it is able to accurately interpolate and discover the correct equation up to Δ𝑥 ≈ 1.2𝑙 𝑐 . The reason for this is that NN-based interpolation constructs a surrogate of the data, informed by both the spatial and the temporal dynamics of the data set, while classical interpolation is intrinsically limited to a single time frame.In figure4.18c) we consider the same graph with 20% white noise on the data. Despite smoothing, the spline is unable to construct an accurate library and fails to discover the correct equation in every case. DeepMoD stands in stark contrast, discovering the correct equation with comparable relative error as in the 0% noise case.

	a)		b)
	Grid sampling	x		Spline grid NN grid NN random NN shifted	Failed Succes	No noise
	Shifted sampling	x	c)	Spline grid NN grid	Failed Succes	20% noise
				NN random
	Δ𝑥 𝑙 𝑐 𝜖 as the sum of the relative errors for all the derivatives, normalized over <5 . We study the relative error every frame, 𝜖= 3 ∑ 𝑖=1 ⟨ ∥𝜕 𝑖 𝑥 𝑢 𝑗 -𝜕 𝑖 𝑥 ̂𝑢 𝑗 ∥ 2 ∥𝜕 𝑖 𝑥 𝑢 𝑗 ∥ 2 ⟩ 𝑗 (4.28) time x where Random sampling NN shifted

A nice example of the time it takes for theories to develop is the anomalous perihelion procession of Mercury. It was already known in 1859 that Mercury's orbit could not be completely described by Newtonian mechanics, but it remained unexplained until Einstein's postulation of general relativity in 1915.

2 A great example of the extrapolation power of models is Planck's law. To explain the blackbody emission spectrum, Planck derived a model which exchanged energy in discrete

This term is mostly used by natural scientists working on physical datasets. Several other terms are in use, depending on the field and the exact goal. In engineering, it is often referred to as system identification, while in machine learning it is often known as causal inference.

In this setting also often referred to as dynamical systems.

I'm neglecting a huge group of models here: those with spatially or temporally varying coefficients! I'll return to this issue in the next chapter.

I'll illustrate later how bad numerical differentiation really is.

More generally put, data-driven modelling improves first-principle modelling by yield-

For more information on the KS equation, see Wikipedia or Encyclopedia of Math.

Although they all share a diffusive term -we can't seem to escape entropy.

An interesting question is if this is a fundamental fact or a consequence of our approach to physics.

I'm only considering non-coupled, one-dimensional cases here, but the argumentation extends to higher dimensions. The library size will be much bigger, and the combinatorial search will be even more expensive, but most features will still be a combination of two components.

The way I like to think about this, is that with this approach we're essentially fitting all models at the same time.

It has an analytical solution but for numerical stability one often uses QR decomposition or SVD. Addtionally, least squares is a so-called BLUE -which means that the estimate ξ is unbiased and has the lowest variance as long as the errors are uncorrelated.

7 Technically this norm is only appropriate for Gaussian, white noise (i.e. zero-mean, uncorrelated). Extensions to other noise models are called generalized linear models, but in practice only white noise is

used.8 The 1/𝑁 factor is not strict required, but has the benefit of making the regularization strength 𝜆 sample size-independent.

More generally, 𝜉 can be penalized by any 𝑝-norm, but in practice only 𝑝 = 1, 2 are used.

In literature this sometimes also referred to as the oracle property.

Simply add √ 𝜆 along the diagonal of the gram matrix, or augment the data and solve it like a standard least squares problem

We can also normalize the target vector, i.e. the time derivative, to make everything dimensionless.

Technically this is group sparsity, but since almost always lasso is used, these terms are used interchangeably.

A interesting approach not based on sparse regression is taken by[START_REF] Long | PDE-Net 2.0: Learning PDEs from Data with A Numeric-Symbolic Hybrid Deep Network[END_REF] . They note that since derivatives can be approximated as a convolution, learnable filters can be used to discover the underlying equation, even with spatially varying sources.

Pareto optimality is a very wide applicable concept, as it essentially describes how to make the trade-off between multiple conflicting objectives.

The field of Bayesian approaches to model selection is so rich it deserves a book on its own. I restrict myself here mainly to some approaches used in the model discovery of PDEs.

Actually, we can exactly calculate the derivative, but the trade-off is that we need to approximate the data using some other model -a surrogate!

Note that the surrogate model doesn't need to be interpretable -we're merely trying to represent and interpolate the data in such a a way that we can easily take derivatives.

Transforming the entire PDE to Fourier space can indeed get rid off all derivatives, but unfortunately non-linear terms correspond to convolutions in fourier space, making everything much more complicated.

This is also known as the spectral method to calculating derivatives, and is often used when solving PDEs[START_REF] Schaeffer | Learning partial differential equations via data discovery and sparse optimization[END_REF]

I'm seeing this as basically Taylor expanding around every sample, and making sure everything is continuous.

In fact, neural networks are known to be universal function approximators, meaning that a network with a single hidden layer of infinite width can approximate any continuous function.

but not the solution!

Similar terms can be added for initial and boundary conditions, but we omit these here for clarity.

In my opinion kernel methods such as GPs have been vastly underused. Kernel methods typically perform well in small datasets and can be made non-local through the use of neural networks, so-called Deep GPs.

Invertible neural networks are an active field -ICML has a recurring workshop on them.

If every single layer is invertible, then so will the whole flow.

In this interpretation model discovery almost comes about as a side-effect; we simply train a PINN with many terms and it simply decides itself which to use.

In bold; 𝑢𝑢 𝑥 : green, 𝑢 𝑥𝑥 : blue and 𝑢 𝑥𝑥𝑥𝑥 : orange

We choose to plot the losses of the original PINN loss ℒ data and ℒ reg because these are more easily interpreted than the likelihood-based losses we have introduced.

meaning that the inferred density is a given and cannot be changed.

When we wrote this paper, we only considered dependence on time -hence the name temporal normalizing flow. Calling them CNFs generalizes our approach and brings it in line with literature.