# Informed Equation Learning

Matthias Werner corresponding author<EMAIL_ADDRESS>Learning Group at ETAS GmbH, Bosch GroupMax Planck Institute for Intelligent
Systems, Tübingen, GermanyUniversity of Tübingen Andrej
Junginger11footnotemark: 1 Philipp Hennig33footnotemark: 322footnotemark: 2
Georg Martius22footnotemark: 2

###### Abstract

Distilling data into compact and interpretable analytic equations is one of
the goals of science. Instead, contemporary supervised machine learning
methods mostly produce unstructured and dense maps from input to output.
Particularly in deep learning, this property is owed to the generic nature of
simple standard link functions. To learn equations rather than maps, standard
non-linearities can be replaced with structured building blocks of atomic
functions. However, without strong priors on sparsity and structure,
representational complexity and numerical conditioning limit this direct
approach. To scale to realistic settings in science and engineering, we
propose an informed equation learning system. It provides a way to incorporate
expert knowledge about what are permitted or prohibited equation components,
as well as a domain-dependent structured sparsity prior. Our system then
utilizes a robust method to learn equations with atomic functions exhibiting
singularities, as e.g. logarithm and division. We demonstrate several
artificial and real-world experiments from the engineering domain, in which
our system learns interpretable models of high predictive power.

##  1 Introduction

Mathematical models formed of equations are key elements in natural sciences
to describe phenomena and their underlying principles. Physics, in particular,
uses equations to build a coherent description of natural laws. For example,
Newton’s and Kepler’s laws or Schrödinger’s equation describe and encode the
phenomena of motion, celestial dynamics (Feynman et al., 1965), or quantum
mechanics (Ballentine, 1970), respectively. In the engineering domain,
mathematical equations are used, e.g., in model-predictive control García et
al. (1989) or as components describing complex system. The designed equations
then correspond to a hypothesis for the inner system’s behavior and represent
its relations and properties. Implementing such control strategies into
embedded control units for industrial-scale applications raises various
additional constraints like low computational effort, minimal memory usage,
low effort of calibration, and generalization to similar systems.

Learning mathematical models in an automated fashion is referred to as
_equation learning_ or _symbolic regression_. The desired mathematical
expressions are usually compositions of _atomic units_ consisting of a set of
functions $\\{f_{i}\\}_{i\leq F}$ and a set of operations that connect the
single functions. The former can be, e.g., functions like
$\\{\sin,\,\mathrm{sqrt},\,\log,\,\dots\\}$ and they can be connected by
operations such as $\\{+,-,*,/,\circ,\dots\\}$. Out of these atomic units,
equation learning finally tries to infer an analytic expression to describe
the relations between the quantities of interest. In this work, we consider a
Gaussian regression setting with a parameterized analytical function
$\psi_{\theta}$, to map a $D$-dimensional input ${\bm{x}}$ to a
$D^{\prime}$-dimensional output ${\bm{y}}$

${\bm{y}}_{i}=\psi_{\theta}({\bm{x}}_{i})+\bm{{\epsilon}}\;,\qquad\bm{{\epsilon}}\sim{\mathcal{N}}(0,\gamma^{2})\;.$
(1)

Here, $\theta$ denotes the parameters, the dataset
$({\bm{X}},{\bm{Y}})\in\mathcal{D}$ is assumed to be sampled iid., and it
contains $N$ data points each of which has Gaussian noise $\bm{\epsilon}$ with
variance $\gamma^{2}$. The structure of $\psi_{\theta}$ as well as its
parameters shall then be inferred from the data. Many contemporary supervised
machine learning algorithms optimize for prediction performance alone. Thus,
they typically produce unstructured and dense equations, which do not reveal
insight into the underlying principles of the data-generating system.
Following the principle of Occam’s razor (MacKay, 1992), an equation learning
system must trade-off complexity against expressivity to overcome this issue.
This typically results in a regularization of the objective function. A
generic symbolic regression method using machine learning, named
$\mathrm{EQL}$ (Martius and Lampert, 2016), replaces standard link functions
in a multi-layer feed-forward neural network with atomic units
($\sin,\,\cos,\,\mathrm{identity},\,*$) and a Lasso regularization
(Tibshirani, 1996). In principle, this approach can be extended to all
continuously differentiable functions with unbounded domains. However, for
common atomic units, like $\mathrm{division}$ or $\mathrm{logarithm}$ with a
half-bounded domain $\mathbb{D}$ or with a singularity, continuous
optimization of an equation-learning neural network can break down for two
reasons: First, the cascading transformation of intermediate results in a deep
architecture can project values outside of the constrained input domains.
Applying domain-limiting mapping functions, such as $\mathrm{softplus}$, is
not ideal since it compromises the final symbolic equation. Second,
singularities in the atomic functions or their derivative can lead to large
gradients and make optimization highly unstable. A first approach
incorporating divisions in the final layer in $\mathrm{EQL}$ has been
presented in Sahoo et al. (2018). Their algorithm uses a hard-coded
curriculum, which depends on the number of epochs itself. This makes it hard
to apply to datasets with different sample sizes. In contrast, our algorithm
overcomes those issues and introduces a learnable parameter instead. It is a
robust and simple algorithm for atomic units with singularities. We show that
this unlocks new classes of functions, as e.g. division, logarithm, etc. in
all layers of the $\mathrm{iEQL}$, which is not covered by previous
publications.

In general, a large expressivity of the equation learner is desired. The
$\mathrm{EQL}$ architecture is, in principle, well suited for scaling to large
problems. However, the set of possible solutions also grows drastically. Thus,
in general, one cannot expect to obtain a symbolic expression that is easily
interpretable or that captures the true correspondences. Therefore, we propose
to incorporate expert knowledge into the system such that suitable solutions
are preferred. Our approach is based on excluding certain combinations of
functions or on providing a user-dependent weighting scheme. It was developed
in close cooperation with applied engineers. Technically, this is achieved by
designing a generalized sparsity regularization, that allows one to choose
specific complexities for each atomic unit, and that can flexibly prohibit
certain, user-defined function compositions.

In this paper we advance the state-of-the-art in symbolic regression by

  1. (a)

proposing an informed equation learner ($\mathrm{iEQL}$) allowing to
incorporating expert knowledge,

  2. (b)

a robust training method for atomic units with singularities (e.g. logarithm
and division), and

  3. (c)

showing that interpretable expressions can be discovered from real-world
datasets from the engineering domain.

The paper is structured as follows: First, we discuss related work in section
2. In section 3, we introduce the architecture of the $\mathrm{iEQL}$, outline
the instance selection criteria and a domain specific complexity measure with
its sparsity inducing regularization. Moreover, we present the robust training
method for atomic units with singularities. In section 4, we present
applications of our method to artificial datasets as well as two real-world
applications in the engineering domain, and we conclude in section 5.

##  2 Related Work

Equation learning studies the relation between input and output data. As a
thread within the wider area of explainability, it is of increasing importance
to machine learning, which mostly produces black box models. A major challenge
in equation learning is to efficiently search the space of possible
expressions, which increases exponentially with the number of atomic units
required to describe the data relation. Current machine learning methods for
equation learning can be separated in bottom up and top down search
strategies.

Bottom up methods improve and learn based on individually sampled mathematical
expressions until the final expression is simple and accurate. For this
purpose, machine learning offers efficient search strategies (Zaremba et al.,
2014; Kusner et al., 2017; Li et al., 2019; Lample and Charton, 2020). A
physics-inspired method for equation learning that uses neural networks to
identify symmetries and separability within the data was introduced by Udrescu
and Tegmark (2020).Common methods to overcome the exponential growth of the
search space are genetic programming and evolutionary algorithms (Langdon et
al., 2008; Searson et al., 2010; McRee, 2010). Algorithms in this vein have
been used to automate the discovery of natural laws (Dubčáková, 2011; Schmidt
and Lipson, 2009). However, none of these methods use gradient information
about the learned expressions during training. This was addressed in a
reinforcement learning formulation by Petersen (2019). They transformed the
discrete optimization to a continuous gradient based optimization scheme via
risk-seeking policy gradients, which is a modification of the Reinforce
(Williams, 1992) policy gradient.

Top down search strategies are in contrast to these bottom up search
strategies. They start with a highly general expression and try to omit
irrelevant parts during training until a simple but accurate expression is
found. This idea is also applied by the $\mathrm{EQL}$ (Martius and Lampert,
2016), which defines a neural network architecture as general expression. It
represents all possible symbolic expressions within its architecture. During
training irrelevant parts get pruned until the $\mathrm{EQL}$ itself
represents a simple but accurate equation. A similar architecture has been
studied by Long et al. (2019) on differential equations, and the integration
of an $\mathrm{EQL}$-like architecture within other deep learning frameworks
has been studied by Kim et al. (2020). Also in physics it has been used to
obtain analytical expressions of classical free energy functionals in Lin et
al. (2020). Generally, these top down methods have the advantage that they
permit the use of efficient continuous gradient based optimization. In
contrast to (Martius and Lampert, 2016; Sahoo et al., 2018), which used
$L_{1}$ regularization, we use $L_{0}$ regularization (Louizos et al., 2018),
which does not shrink the weights during training and allows to naturally
align a domain specific complexity measure with the regularization. A first
approach to use divisions in $\mathrm{EQL}$ has been presented by Sahoo et al.
(2018).Their architecture only considers division units in the final layer. We
propose a robust training method for atomic units with singularities in
general that does not require a predefined curriculum and allows to using
division units in all layers.

##  3 Method

###  3.1 Architecture of the $\mathrm{iEQL}$

$x_{1}$$x_{2}$$x_{3}$$x_{4}$$y_{1}$$y_{2}$$B$$B$$B$$B$$B$$\cos$$\exp$$\log$$*$$\nicefrac{{\;}}{{\;}}$$\mathbf{copy}$$\cos$$\exp$$\log$$*$$\nicefrac{{\;}}{{\;}}$$\mathbf{copy}$$\cos$$\exp$$\log$$*$$\nicefrac{{\;}}{{\;}}$$\mathbf{copy}$$\cos$$\exp$$\log$$*$$\nicefrac{{\;}}{{\;}}$$\mathbf{copy}$lin
Figure 1: Exemplary setup of the $\mathrm{iEQL}$ with 4 hidden layers, three
types of unary atomic units $\\{\cos,\,\exp{\;},\,\log\\}$ and two types of
binary atomic units $\\{*,\,/\\}$ per layer. The $\mathrm{copy}$ units
contains all units from all previous layers. The red dashed arrows mark
forbidden connections (see equation 4) which are removed, also from the
$\mathrm{copy}$ units. After successful training, the $\mathrm{iEQL}$
represents the desired equation. The green connections visualize an exemplary
model for the equation $y_{2}=0$ and $y_{1}=(x_{2}\cdot x_{1})/(\cos
x_{1}+1.5)$.

Our goal is to identify a compact equation describing data accurately. A
standard neural network as a universal function approximator is representing a
large uninterpretable equation. Replacing units of a neural network with
various building blocks of equations and forcing a sparse connection in the
network can yield simple and interpretable mathematical equations. The
proposed architecture, see Figure 1, of a multi-layered feed-forward equation
learning network is similar to previous work by Martius and Lampert (2016);
Sahoo et al. (2018). It consists of $L$ layers, whose $L-1$ hidden layers
correspond to the maximum number of function compositions. Each hidden layer
applies a non-linear transformation, which consists of $u$ unary atomic units
$f_{i\leq u}:\mathbb{R}\to\mathbb{R}$ and $v$ binary atomic units $g_{j\leq
v}:\mathbb{R}\times\mathbb{R}\to\mathbb{R}$. The latter can be, e.g., division
$a/b$ or multiplication $a\cdot b$. We apply an affine transformation on the
output of all preceding layers, including the input data
(${\bm{y}}^{0}:={\bm{x}}$), instead of just the previous one (similarly used
in dense-nets (Huang et al., 2016))

$\displaystyle{\bm{z}}^{l}$
$\displaystyle={\bm{W}}^{l}\tilde{\bm{y}}^{l-1}+{\bm{b}}^{l}\,,\,\text{with:
}\quad\tilde{\bm{y}}^{l-1}=({\bm{y}}^{0},\dots,{\bm{y}}^{l-1})$ (2)
$\displaystyle y^{l}$
$\displaystyle=\big{(}f_{1}(z_{1}^{l}),..,f_{u}(z_{u}^{l}),g_{1}(z_{u+1},z_{u+2}),\dots\big{)}.$
(3)

Those skip connections allow identified sub-solutions (terms in the equation)
to be available with no additional cost to all subsequent layers. Figure 1
visualizes a possible layout of the architecture, in which all previous’
layers output is marked by the $\mathbf{copy}$ units. Naturally, the skip
connections favor less nested expressions and allow to ignoring unnecessary
layers. Moreover, the search for appropriate hyperparameters is simplified as
the initial network can be chosen larger and more complex. Depending on the
application, certain combinations of atomic units like

$\displaystyle\cos(\cos(\cdot))$ $\displaystyle\cos(\exp(\cdot))$
$\displaystyle\exp(\exp(\cdot))$ $\displaystyle\log(\log(\cdot)$ (4)

with respective arguments can be undesirable (i.e. they should only be chosen
rarely), or they might not make sense at all. _Expert_ and _domain knowledge_
can provide such information on possible equation structures. Therefore, all
forbidden connections are removed from the architecture as shown in figure 1.
That way, an expert finally decides which combinations to use and which to
exclude. Such information could not be used in previous work by Martius and
Lampert (2016); Sahoo et al. (2018), since the respective architecture relied
on $\mathrm{linear}$ units, which accumulate output of all previous layer in
an affine transformation. That made it impossible to remove certain
combinations. However, with the present framework based on $\mathrm{copy}$
units, it is easily possible to prohibit certain combinations of atomic units
in the $\mathrm{iEQL}$.

###  3.2 Model Selection Criteria

We apply model selection criteria V${}^{\text{int}}$, V${}^{\text{int}}$-S and
V${}^{\text{int\&ex}}$ proposed by Sahoo et al. (2018) to select a good
equation amongst all found equations. They are based on a normalized
validation error $\tilde{\nu}_{\text{in}}$, normalized complexity
$\tilde{s}^{2}$ and a normalized extrapolation validation error
$\tilde{\nu}_{\text{ex}}$, which requires some points (here 40) from the
extrapolation domain

$\displaystyle\text{V${}^{\text{int}}$}=\operatorname*{arg\,min}[\tilde{\nu}_{\text{int}}^{2}]$
(5)
$\displaystyle\text{V${}^{\text{int}}$-S}=\operatorname*{arg\,min}[\tilde{\nu}_{\text{int}}^{2}+\tilde{s}^{2}]$
(6)
$\displaystyle\text{V${}^{\text{int\&ex}}$}=\operatorname*{arg\,min}[\tilde{\nu}_{\text{int}}^{2}+\tilde{\nu}_{\text{ex}}^{2}]\,.$
(7)

###  3.3 Complexity Measure

The “complexity” of an equation — not the computational complexity per se, but
whether an expert would consider it “basic” or “involved”, depends both on the
number of terms and a domain specific complexity cost of each atomic
operation. To capture this aspect, we use a weighted sum of the number of
parameters as a measure of complexity. This is a commonly accepted way to
introduce domain knowledge in symbolic regression, e.g. in evolutionary search
(Dubčáková, 2011). An expert defines domain-specific complexity factors
$c_{u}$ for each atomic unit type $u$. The $c_{u}$ can depend on soft
criteria, like domain or personal experience, but also quite concretely on
computational cost. A specific choice of the complexity factors $c_{u}$ for
our experiments is shown in table 1.

Table 1: Domain specific complexity factors $c_{u}$ for different atomic
units. No preference is denoted as plain. For the real world applications
(power loss of an electric machine (section 4.3) and torque model of a
combustion engine (section 4.4)). We use values suggested by domain experts
denoted as motor.

$\nicefrac{{+}}{{-}}$ $*$ $/$ $x^{2}$ $\log$ $\sqrt{\phantom{x}}$ $\exp$
$\cos$ plain 1 1 1 1 1 1 1 1 motor 1 2 5 2 5 3 5 10

###  3.4 Sparsity Inducing Method

Sparsity inducing regularization is necessary to extract interpretable and
compact equations from the $\mathrm{iEQL}$. We are able to directly optimize
for the complexity measure, described in section 3.3, as regularization
$\mathcal{L}_{\text{C}}$ and the mean squared error as data loss
$\mathcal{L}_{\text{D}}$:

$\displaystyle\mathcal{L}_{\text{D}}=\frac{1}{N}\sum\limits_{i=1}^{N}\lVert{\bm{y}}_{i}-\mathrm{iEQL}({\bm{x}}_{i},{\bm{W}})\rVert_{2}^{2}\,,\quad\mathcal{L}_{\text{C}}=\sum\limits_{u=1}^{U}c_{u}\lVert
W_{u}\rVert_{0}\,.$ (8)

$W_{u}$ represents all weights, which correspond to an atomic unit type.
Opposed to the $L_{1}$ regularization scheme used by Martius and Lampert
(2016); Sahoo et al. (2018), we use a differentiable version of
$L_{0}$-regularization (Louizos et al., 2018). A final retraining of the
pruned network is not necessary anymore, as it does not induce shrinkage on
the actual values of the weights. Minor modifications are required to apply
the original implementation of Louizos et al. (2018) to weights instead of
nodes. Each weight is multiplied by a non-negative stochastic Bernoulli
distributed gate $g_{j}\sim\mathrm{Ber}(\pi_{j})$. The dropout rate of each
weight is thus given by $1-\pi_{j}$. The expectation of a weight being
relevant is given by $\pi_{j}$. Those gates learn collectively which weights
are relevant. The expected objective is then

$\displaystyle\mathcal{L}=\mathbb{E}_{q({\bm{g}}\mid\pi)}\left[\frac{1}{N}\sum\limits_{i=1}^{N}\lVert{\bm{y}}_{i}-\mathrm{iEQL}({\bm{x}}_{i},{\bm{W}})\rVert^{2}_{2}\right]+\lambda\sum\limits_{u=1}^{U}c_{u}\sum\limits_{j=1}^{\lvert
W_{u}\rvert}\pi_{j}^{u}\,.$ (9)

$\lvert W_{u}\rvert$ is the total number of weights that correspond to an
atomic unit type ($u$). Further details on the variational optimization scheme
are given in appendix A.3.

$\displaystyle-1$$\displaystyle 0$$\displaystyle 1$$\displaystyle
2$$\displaystyle 3$$\displaystyle 4$$\displaystyle 5$$\displaystyle
0$$\displaystyle 2$$\displaystyle 4$$\displaystyle 6$$\displaystyle
8$$\displaystyle
z$$f(z)=1/z$$f$$\alpha\\!=\\!0.10$$\alpha\\!=\\!0.36$$\alpha\\!=\\!0.56$$\alpha\\!=\\!1.00$

$\displaystyle-1$$\displaystyle 0$$\displaystyle 1$$\displaystyle
2$$\displaystyle 3$$\displaystyle 4$$\displaystyle
5$$\displaystyle-4$$\displaystyle-2$$\displaystyle 0$$\displaystyle
2$$\displaystyle z$$f(z)=\log{z}$

Figure 2: Relaxation of the division (left) and logarithm (right) for
different values of $\alpha$.

###  3.5 Atomic Units with Singularities

Our algorithm is designed for atomic units with singularities111we refer to
them as singular units $f:\;\mathbb{D}\to\mathbb{R}$ on a half-bounded domain
$\mathbb{D}\equiv(a,\infty)$ exhibiting a singularity at $a$ (e.g.
$\log:(0,\infty)\to\mathbb{R}$ with a singularity at $a=0$). It also applies
to $\mathrm{division}$ units $\nicefrac{{c}}{{d}}$ under the assumption that
real systems do not diverge in the domain of application. Then, it is
sufficient to only consider the positive branch of the hyperbola
$\nicefrac{{1}}{{d}}$ and let $c$ choose the sign. Cascading transformation of
intermediate results in the deep architecture of the $\mathrm{iEQL}$ can
project values outside of the constrained input domains of singular units.
Therefore, the domain $\mathbb{D}$ of a singular unit is continuously extended
to $\mathbb{R}$ to avoid forbidden inputs during training (see equation 11),
similar to Sahoo et al. (2018). During training, an additional penalty
$\mathcal{L}_{\text{su}}$ is necessary to constrain the solution space of the
$\mathrm{iEQL}$ to networks that respect the domain $\mathbb{D}$ of all its
singular units $f_{u}(z_{u})$

$\mathcal{L}_{\text{su}}=\frac{1}{N}\sum\limits_{i=1}^{N}\sum\limits_{u\in\mathcal{SU}}\mathcal{L}_{\text{su}}^{\text{u}}(z_{i}^{u})\,,\qquad\mathcal{L}_{\text{su}}^{\text{u}}(z_{i}^{u})=\begin{cases}0,&\text{for
}z_{i}^{u}>a,\\\ |a-z_{i}^{u}|,&\text{for }z_{i}^{u}\leq a\,.\end{cases}$ (10)

The set of all singular units in the $\mathrm{iEQL}$ is given by
$\mathcal{SU}$. Still, training can be corrupted by unconstrained values and
exploding gradients in the neighborhood of the singularities. Therefore, we
propose a learnable relaxation of the singular units that does not compromise
the final symbolic equation. It is shown in figure 2 for the logarithm and
division. The input $z$ is shifted $\hat{z}=z+\alpha$ by a learnable positive
relaxation-parameter $\alpha=\log(1+e^{\hat{\alpha}})>0$

$\hat{f}(\hat{z})=\begin{cases}f(z+\alpha),&\text{for }z>a,\\\ 0,&\text{for
}z\leq a\end{cases}\;.$ (11)

Since the input $z=Wy+b$ is an affine transformation of the previous layer
with bias $b$ the relaxation-parameter $\alpha$ can be safely added to the
bias and does not affect the structure of the final equation. This method
assures that, during and after training, the maximum absolute value in the
neighborhood of the singularity is smaller than $|f(a+\alpha)|$ and the
maximum value of the derivative is smaller than $|\partial f(a+\alpha)|$ on
the training dataset. The final objective function is thus a sum of data loss
$\mathcal{L}_{\text{D}}$, regularization $\mathcal{L}_{\text{C}}$ and domain
penalty $\mathcal{L}_{\text{su}}$ with domain penalty strength $\delta$

$\displaystyle\mathcal{L}=\mathcal{L}_{\text{D}}+\delta\,\mathcal{L}_{\text{su}}+\lambda\,\mathcal{L}_{\text{C}}\,.$
(12)

In order to guarantee that the domain constraints on singular units also hold
on the extrapolation area and that the order of magnitude does not change
drastically on extrapolation data, intrinsic penalty epochs are necessary as
proposed in Sahoo et al. (2018). The bound-penalty
$\mathcal{L}_{\text{bound}}$ is given by

$\displaystyle\mathcal{L}_{\text{bound}}=\frac{1}{N}\sum\limits_{i=1}^{N}\operatorname*{arg\,max}(\lvert\mathrm{iEQL}({\bm{x}}_{i},{\bm{W}})\rvert-B,0)$
(13)

For an intrinsic penalty epoch, we randomly sample $N$ input data points from
the test area without labels. Then the $\mathrm{iEQL}$ is trained on the loss
function

$\displaystyle\mathcal{L}_{\text{penalty}}=\delta\,\mathcal{L}_{\text{su}}+\delta\,\mathcal{L}_{\text{bound}}\,.$
(14)

##  4 Experiments

We demonstrate the application of the $\mathrm{iEQL}$ to four different use
cases. The first one is to learn complex equations on simulated data of
several equations. The second use case is a simulated ambiguity dataset, with
which we demonstrate how to incorporate expert knowledge to influence the
outcome of the equation learner. Additionally, we study the relative frequency
of selected atomic units in found equations with and without expert knowledge
on a simulated dataset and a real world dataset. Use cases three and four both
address real-world applications in industry determining expressions for the
power loss of an electric machine and a torque model of a combustion engine.
All experiments are performed using an $\mathrm{iEQL}$ with four hidden
layers. Each hidden layer has
$\\{\cos,\,\exp{\;},\,\log,\,\sqrt{\phantom{x}},\,x^{2},\,*,\,/\\}$ as atomic
units, and each atomic unit is applied four times in each layer. Combinations
from equation 4 are prohibited. We compare the $\mathrm{iEQL}$ to five
different algorithms:

  * •

EQL÷, a state-of-the-art method from Sahoo et al. (2018) with atomic unit
types $\\{\sin,\cos,*,\text{identity}\\}$ in each hidden layer and division in
the final layer,

  * •

a multi-layer perceptron (MLP) with $\tanh$ activation functions and five
hidden layers with 50 neurons each,

  * •

a genetic algorithm (PySR,Cranmer (2020)) with two different configurations
$\mathrm{GA1}$ and $\mathrm{GA2}$,

  * •

a Gaussian Process (GP) for the real world datasets, calculated with ascmo
(Hoffmann et al., 2015), which is a standard tool from the engineering domain,

  * •

the mean predictor (MP) on the train set.

Further details on training and parameter settings are outlined in appendix
A.1. All experiments are executed five times, and we report median, minimum
and maximum (in sub and superscript) of root mean squared error on test
datasets.

###  4.1 Learning Complex Equations

We evaluate the $\mathrm{iEQL}$ on two different sets of ground truth
equations shown in table 2: (S0–S6) with $\\{\pm,\sin,\cos,*,/\\}$ as
operations and (A1–A4) also including $\\{\log,\exp\\}$. A two dimensional
contour plot of all datasets is shown in appendix figure 6. Most functions are
gathered from different papers on symbolic regression. A star (∗) indicates
that input and output values were scaled. We closely follow the evaluation
strategy proposed in Sahoo et al. (2018) with their model selection criteria
V${}^{\text{int\&ex}}$ with 40 extrapolation points (see section 3.2) on the
found equations of $\mathrm{iEQL}$, EQL÷ and $\mathrm{GA}$. Training datasets
consist of $10^{4}$ randomly sampled data points in the train domain. Outputs
are corrupted with standard normal noise of standard deviation $\gamma=0.01$.
For validation $10\%$ of the training dataset is used. Test datasets consist
of $10^{3}$ randomly sampled data points from the test domain (see table 2).
Note that for the equations S0–S6 and A1 with their four dimensional input,
the extrapolation area of the test dataset is 15 times larger than its
interpolation area. The included extrapolation domain is a good indicator for
whether or not the ground truth has been found.

Table 2: Complex equations on which the $\mathrm{iEQL}$ is evaluated. Most
functions are gathered from different papers on symbolic regression. A star
(∗) indicates that input and output were scaled.

Data- Ground truth motivated Domain set equation from train test S0
$y=(1-x_{2}^{2})/(\sin(2\pi\,x_{1})+1.5)$ $[-1,1]^{4}$ $[-2,2]^{4}$ S1
$y=\left[\sin(\pi x_{1})+\sin(2\pi
x_{2}+\nicefrac{{\pi}}{{8}})+x_{2}-x_{3}x_{4}\right]/3$ Sahoo et al. (2018)
$[-1,1]^{4}$ $[-2,2]^{4}$ S2 $y=\left[\sin(\pi
x_{1})\\!+\\!x_{2}\\!\cos(2\\!\pi
x_{1}+\nicefrac{{\pi}}{{4}})\\!+\\!x_{3}\\!-\\!x_{4}^{2}\right]/3$ Sahoo et
al. (2018) $[-1,1]^{4}$ $[-2,2]^{4}$ S3 $y=\left[(1+x_{2})\sin(\pi
x_{1})+x_{2}x_{3}x_{4}\right]/3$ Sahoo et al. (2018) $[-1,1]^{4}$ $[-2,2]^{4}$
S4∗
$y=(3.0375\,x_{1}x_{2}+5.5\sin{(\nicefrac{{9}}{{4}}\,(x_{1}-\nicefrac{{2}}{{3}})(x_{2}-\nicefrac{{2}}{{3}}))})/5$
Jin et al. (2019) $[-1,1]^{4}$ $[-2,2]^{4}$ S5∗
$y=\frac{(5x_{1})^{4}}{(5x_{1})^{4}+1}+\frac{(5x_{2})^{4}}{(5x_{2})^{4}+1}$
Trujillo et al. (2016) $[-1,1]^{4}$ $[-2,2]^{4}$ S6∗
$y=((1-x_{1})^{2}+(1-x_{3})^{2}+100(x_{2}-x_{1}^{2})^{2}$ 4D-Rosenbrock
$[-1,1]^{4}$ $[-2,2]^{4}$ $\quad\quad+100(x_{4}-x_{3}^{2})^{2})/1500$ A1∗
$y=(1.5\mathrm{e}^{1.5x_{1}}+5\cos{(3\,x_{2})})/10$ Jin et al. (2019)
$[-1,1]^{4}$ $[-2,2]^{4}$ A2∗ $y=\log(2x_{2}+1)-\log(4x_{1}^{2}+1)$ Trujillo
et al. (2016) $[0,1]^{2}$ $[0,2]^{2}$ A3∗
$y=\frac{\exp(-(4x_{1}-0.4)^{2})}{1.2+(4x_{2}-1.9)^{2}}$ Trujillo et al.
(2016) $[0,1]^{2}$ $[0,2]^{2}$ A4∗
$y=\frac{x_{3}+1}{3\pi^{2}}\frac{(5x_{1}+1)^{3}}{\exp\left(\frac{(x_{3}+1)(5x_{1}+1)}{(x_{4}+1)(0.5x_{2}+1)}\right)-1}$
Udrescu and Tegmark (2020) $[0,1]^{4}$ $[0,2]^{4}$

Table 3: Learning complex equations S0-S6: median, minimum and maximum (in sub
and superscript) of root mean squared error (RMSE) on test datasets. Baselines
are the EQL÷ from Sahoo et al. (2018), a genetic algorithm (GA) from Cranmer
(2020), a multi-layer perceptron (MLP) and the mean predictor (MP) on the
train dataset.

S0 S1 S2 S3 S4 S5 S6 MP 1.57 0.66 0.75 0.64 1.17 0.44 0.82 MLP 1.11
$\,{}^{\scriptscriptstyle 1.17}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
1.04$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.45 $\,{}^{\scriptscriptstyle 0.54}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.35$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.41 $\,{}^{\scriptscriptstyle 0.42}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.39$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.36 $\,{}^{\scriptscriptstyle 0.37}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.34$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.65 $\,{}^{\scriptscriptstyle 0.70}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.57$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.07 $\,{}^{\scriptscriptstyle 0.15}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.02$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.65 $\,{}^{\scriptscriptstyle 0.66}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.65$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$\mathrm{GA1}$ 1.07 $\,{}^{\scriptscriptstyle 1.58}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.91$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.77 $\,{}^{\scriptscriptstyle 0.77}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.77$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.33 $\,{}^{\scriptscriptstyle 0.35}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.47}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>1.20 $\,{}^{\scriptscriptstyle 1.23}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
1.07$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.29 $\,{}^{\scriptscriptstyle 0.29}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.29$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.82 $\,{}^{\scriptscriptstyle 0.82}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.82$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$\mathrm{GA2}$ 0.39 $\,{}^{\scriptscriptstyle 1.47}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.66 $\,{}^{\scriptscriptstyle 0.77}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.37 $\,{}^{\scriptscriptstyle 0.77}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.60 $\,{}^{\scriptscriptstyle 1.78}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.58$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.01}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.08 $\,{}^{\scriptscriptstyle 0.42}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>4.26 $\,{}^{\scriptscriptstyle\text{inf}}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.74$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>EQL÷ 0.01 $\,{}^{\scriptscriptstyle 0.12}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.17}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.05}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.01}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 1.38}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.01}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.02}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$\mathrm{iEQL}$ 0.01 $\,{}^{\scriptscriptstyle 0.01}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.01}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.01}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.01}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.01}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.01}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.01}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>
All three network-based models (MLP, EQL÷, $\mathrm{iEQL}$) can fit the
training data perfectly. Their prediction error on training and validation
data is at noise level for all datasets S0–S6, and A1–A4. Table 3 shows the
results on the test datasets for S0–S6. Unsurprisingly, the MLP is not able to
capture the data-generating function (ground truth) by any means. Therefore,
it is not able to perform well on the test domain, which includes large
extrapolation domains, except for equation A3, which is mainly zero in the
extrapolation domain. In contrast, $\mathrm{GA}$, EQL÷, and $\mathrm{iEQL}$
can learn meaningful equations, which can extrapolate on the test domain. If
the correct equation has been identified by the V${}^{\text{int\&ex}}$
criteria, the test RMSE is expected to be at noise level. The genetic
algorithm ($\mathrm{GA1}$) captures the data generating function on the
datasets S2 and S3 at least in one experiment and $\mathrm{GA2}$ on the
datasets S0, S1, S2, S4 and S5. On datasets S0–S6 the $\mathrm{iEQL}$ captures
the data generating function reliably for all five runs as shown in table 3.
Despite its larger expressivity, with
$\\{\cos,\,\exp{\;},\,\log,\,\sqrt{\phantom{x}},\,x^{2},\,*,\,/\\}$ atomic
units, it even outperforms the EQL÷ architecture, which uses just the
necessary atomic units $\\{\cos,\,\sin,\,*,\,/\\}$. Equations A1–A4 are more
difficult to learn. $\mathrm{iEQL}$ outperforms all baselines as shown in
table 4.

Table 4: Learning complex equations (A1–A4): median, minimum and maximum (in
sub and superscript) of root mean squared error (RMSE) on test datasets.
Baselines are the EQL÷ from Sahoo et al. (2018), a geneatic algorithm (GA)
from Cranmer (2020), a multi-layer perceptron (MLP) and the mean predictor
(MP) on the train dataset.

A1 A2 A3 A4 MP 0.88 1.09 0.19 0.66 MLP 0.67 $\,{}^{\scriptscriptstyle
0.70}_{{\hbox to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.64$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.24 $\,{}^{\scriptscriptstyle 0.25}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.22$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.02}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.39 $\,{}^{\scriptscriptstyle 0.40}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.37$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$\mathrm{GA1}$ 0.82 $\,{}^{\scriptscriptstyle 1.11}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.28$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.35 $\,{}^{\scriptscriptstyle 0.62}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.26$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.14 $\,{}^{\scriptscriptstyle 0.19}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.14$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.60 $\,{}^{\scriptscriptstyle 0.66}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.56$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$\mathrm{GA2}$ 0.28 $\,{}^{\scriptscriptstyle 0.45}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.04$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.10 $\,{}^{\scriptscriptstyle 0.51}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.06$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.10 $\,{}^{\scriptscriptstyle 0.77}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.04$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.35 $\,{}^{\scriptscriptstyle 0.60}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.28$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>EQL÷ 0.09 $\,{}^{\scriptscriptstyle 0.11}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.09$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.09 $\,{}^{\scriptscriptstyle 0.19}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.06$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.02 $\,{}^{\scriptscriptstyle 0.33}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.36 $\,{}^{\scriptscriptstyle 1.89}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.22$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$\mathrm{iEQL}$ 0.02 $\,{}^{\scriptscriptstyle 0.06}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.03 $\,{}^{\scriptscriptstyle 0.05}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.02$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.01 $\,{}^{\scriptscriptstyle 0.02}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.01$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>0.07 $\,{}^{\scriptscriptstyle 0.13}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.07$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>
###  4.2 Expert Knowledge

$\displaystyle 0.10$$\displaystyle 0.11$$\displaystyle 0.12$$\displaystyle
0.14$$\displaystyle 0.17$$\displaystyle 0.20$$\displaystyle
0.25$$\displaystyle 0.33$$\displaystyle 0.50$$\displaystyle
1.00$$\displaystyle 2.00$$\displaystyle 3.00$$\displaystyle
4.00$$\displaystyle 5.00$$\displaystyle 6.00$$\displaystyle
7.00$$\displaystyle 8.00$$\displaystyle 9.00$$\displaystyle
10.00$$\displaystyle 0$$\displaystyle 20$$\displaystyle 40$$\displaystyle
60$$\displaystyle 80$$\displaystyle 100$ratio $\displaystyle
r=c_{\cos}/c_{x^{2}}$[%]$\displaystyle x^{2}$$\displaystyle\cos$ (a) relative
frequency of atomic units

$\displaystyle 0.1$$\displaystyle 0.3$$\displaystyle 1.0$$\displaystyle
3.0$$\displaystyle 10$$0$$1$$2$ratio $\displaystyle
r=c_{\cos}/c_{x^{2}}$$\displaystyle\\#$ units$\cos$$x^{2}$ (b)
V${}^{\text{int}}$-S with complexity factors

$\displaystyle 0.1$$\displaystyle 0.3$$\displaystyle 1.0$$\displaystyle
3.0$$\displaystyle 10$$0$$1$$2$ratio $\displaystyle
r=c_{\cos}/c_{x^{2}}$$\displaystyle\\#$ units$\cos$$x^{2}$ (c)
V${}^{\text{int}}$-S without complexity factors

Figure 3: The effect of expert knowledge is quantified by means of an
ambiguous dataset with its ground truth given by equation 15. For each
complexity cost ratio $r=\nicefrac{{c_{\cos}}}{{c_{x^{2}}}}$ the relative
frequency of atomic units is calculated for each found equation. Panel (a)
shows the average of all those relative frequencies for $\cos$ and $x^{2}$
units. The number of $\cos$ and $x^{2}$ units appearing in the selected
equation depending on the complexity cost ratio is shown in panel (b) selected
with V${}^{\text{int}}$-S with complexity factors and in panel (c) with
V${}^{\text{int}}$-S without complexity factors.

Expert knowledge is especially important if there are several possible
equations for an ambiguous dataset. Incorporating this knowledge into the
complexity measure as well as the regularization can decide whether the right
equation is found. The following experiment studies how different sets of
complexity factors affect the relative frequency of atomic units in the set of
found equations as well as the importance of complexity factors for the
selection criteria V${}^{\text{int}}$-S, which is applied to the set of found
equations. The selection criteria V${}^{\text{int\&ex}}$ is not suited for
this experiment, since it ignores the model complexity. In order to construct
an ambiguous dataset we choose the ground truth

$\displaystyle y$ $\displaystyle=8\cos(0.5\,x)-4$ (15)

which, close to $x=0$, can be modeled equivalently by a $\cos$ function or a
polynomial of degree 2. The output is corrupted with standard normal noise
${\mathcal{N}}(0,0.01^{2})$. The training dataset consists of $10^{4}$
randomly sampled data points on the train domain $[-1,1]$. It resembles a
parabola. We study 9 configurations that prefer a periodic structure
$c_{\cos}\in\\{\nicefrac{{1}}{{2}},\nicefrac{{1}}{{3}},\dots,\nicefrac{{1}}{{10}}\\}$
and 9 configurations that prefer a polynomial structure
$c_{x^{2}}\in\\{\nicefrac{{1}}{{2}},\nicefrac{{1}}{{3}},\dots,\nicefrac{{1}}{{10}}\\}$
as well as the plain configuration. For each complexity cost ratio
$r=\nicefrac{{c_{\cos}}}{{c_{x^{2}}}}$ we calculate the relative frequency of
atomic units for each of the 78 found equations (number of different
regularization strength, see appendix A.1). Figure 3(a) shows the average of
all those relative frequencies for $\cos$ and $x^{2}$ units. As desired the
relative frequencies follow the prior. Ratios with $r<1$ prefer $\cos$ units
and ratios with $r>1$ prefer $x^{2}$ units. Without any further information
($r=1$) a polynomial solution seems very plausible. It turns out that
$\mathrm{iEQL}$ does not converge to all kinds of simple ambiguous solutions
with the same probability. Often, certain solutions are more likely to
converge to. In this experiment, for instance, $\mathrm{iEQL}$ predominantly
uses $x^{2}$ units without a biased preference ($r=1$). Figure 3(b) shows the
number of $\cos$ and $x^{2}$ units appearing in the selected equation
depending on the complexity cost ratio. Ratios with $r<1$ prefer a periodic
equation and ratios with $r>1$ prefer a polynomial equation. As desired, the
prior is followed in case of such a strong ambiguity. Without complexity
factors a polynomial solution is mostly selected as shown in figure 3(c).
Hence, complexity factors are crucial to select the correct equation with
respect to the prior. Since in this case both types of found equations have
the same complexity we found that the Adam optimizer converges better for
$x^{2}$ atomic units. This leads to the preference of polynomial solutions.

Figure 4(a) shows the impact of motor specific complexity factors (shown in
table 1) on the relative frequencies of atomic units evaluated on all found
equations (5 experiemnts, each with 78 different regularization strength, see
appendix A.1) on the real world dataset ‘torque of an internal combustion
engine‘, which is introduced in section 4.4. The motor specific complexity
factors clearly reduce the relative frequency of $\\{\exp,\cos,/\\}$ atomic
units and enhances polynomial structures with $\\{x^{2},*\\}$ units. The
$\sqrt{\phantom{x}}$ unit has been barely chosen at all. This coincides with
the expected behavior. Figure 4(b) shows the impact of penalizing $\cos$ units
with $(c_{\cos}=5)$ on a more sophisticated synthetic dataset (S5), as
expected the relative frequency of the cosine units is reduced.

logsqrtexpcosdiv$\displaystyle x^{2}$mul$\displaystyle 0$$\displaystyle
10$$\displaystyle 20$$\displaystyle 30$[%]plainmotor53510522$c_{u}$ (a)
combustion engine

logsqrtexpcosdiv$\displaystyle x^{2}$mul$\displaystyle 0$$\displaystyle
10$$\displaystyle 20$$\displaystyle 30$$\displaystyle
40$[%]plain$c_{\cos}=5$1115111$c_{u}$ (b) S5 dataset

Figure 4: Relative frequency of atomic units with and without domain specific
factors over all found equations (5 experiemnts, each with 78 different
regularization strength, see appendix A.1) of the combustion engine
experiments in (a) and all found equations of the S5 dataset experiments in
(b). Domain specific complexity factors $c_{u}$ are shown in the upper x-axis.

###  4.3 Power Loss of an Electric Machine

Table 5:  Data properties for the power loss of an electric machine. The train domain covers just $80\%$ the range of operation. Except for rotor temperature $T_{\text{rot}}$ which has just three different operation points.  Quantity | Test Domain | Train Domain | Description | Type
---|---|---|---|---
$I_{\text{D}}[\text{A}]$ | $[1,525]$ | $[76,451]$ | direct current | input
$I_{\text{Q}}[\text{ A}]$ | $[1,525]$ | $[76,451]$ | quadratic current | input
$T_{\text{rot}}[^{\circ}\text{C}]$ | $[-20,150]$ | $[-20,150]$ | rotor temperature | input
$n_{\text{mtr}}[\text{rpm}]$ | $[1,16000]$ | $[2001,14001]$ | motor speed | input
$P_{\text{mod}}[\text{W}]$ | $[0,4558]$ | $[95,3364]$ | power loss | output

In this section, we address with the power loss of an electric machine a
relevant problem in the industrial domain (Buchner et al., 2020). Increasing
requirements on control of electric machines in automotive power trains lead
to the necessity of more and more sophisticated and complex models to describe
the system. Those models have to be able to capture non-linear effects at
different components in the system, as in the dataset for the power loss of an
electric machine. Further, those models have to be minimal in computational
effort and memory demand due to embedded hardware and latency constraints in
the range of micro-seconds.

In this industrial dataset, the power loss of a winding head is measured
depending on the direct current $I_{\text{D}}$, quadratic current
$I_{\text{Q}}$, rotor temperature $T_{\text{rot}}$, and motor speed
$n_{\text{mtr}}$. The data was measured at stationary operation points. It
contains $24684$ data points measured at equidistant variations of the
quantities listed in table 5 within their range of operation. The test dataset
contains unseen data from the entire domain and training dataset contains data
only from $80\%$ of its range of operation. Further details on data
preparation are given in appendix A.2.

We use the model selection criteria V${}^{\text{int}}$-S on the found
equations of $\mathrm{iEQL}$, EQL÷ and $\mathrm{GA}$.

Table 6: Results on real-world datasets. Reported are median, minimum and
maximum (in sub and superscript) of root mean squared error (RMSE) on the real
world test datasets. Domain knowledge is used for
$\mathrm{iEQL_{\text{motor}}}$ as given in table 1. We use the model selection
criteria V${}^{\text{int}}$-S and state the number of active parameters. For
the combustion engine dataset we also present the best validation models
selected with V${}^{\text{int}}$.

electric machine [W] combustion engine [Nm] V${}^{\text{int}}$-S #param.
V${}^{\text{int}}$-S #param. V${}^{\text{int}}$ #param. MP 1042.70 60.17 60.17
GP 0.92 1.79 1.79 $\mathrm{GA1}$ $188.40$ ${\ {}^{\scriptscriptstyle
347.93}_{{\hbox to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
155.84$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$9$ ${\ {}^{\scriptscriptstyle 9}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
3$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$16.79$ ${\ {}^{\scriptscriptstyle 22.37}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
13.35$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$7$ ${\ {}^{\scriptscriptstyle 11}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
5$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$5.91$ ${\ {}^{\scriptscriptstyle 7.66}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
5.00$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$33$ ${\ {}^{\scriptscriptstyle 40}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
23$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$\mathrm{GA2}$ $230.23$ ${\ {}^{\scriptscriptstyle 230.23}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
230.23$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$8$ ${\ {}^{\scriptscriptstyle 8}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
8$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$22.37$ ${\ {}^{\scriptscriptstyle 22.37}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
16.79$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$5$ ${\ {}^{\scriptscriptstyle 7}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
5$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$5.33$ ${\ {}^{\scriptscriptstyle 6.29}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
4.86$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$35$ ${\ {}^{\scriptscriptstyle 39}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
34$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>EQL÷ $0.03$ ${\ {}^{\scriptscriptstyle 0.04}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.03$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$10$ ${\ {}^{\scriptscriptstyle 14}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
7$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$1.75$ ${\ {}^{\scriptscriptstyle 1.96}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
1.55$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$95$ ${\ {}^{\scriptscriptstyle 129}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
70$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$1.46$ ${\ {}^{\scriptscriptstyle 1.61}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
1.35$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$382$ ${\ {}^{\scriptscriptstyle 427}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
225$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$\mathrm{iEQL}$ $0.03$ ${\ {}^{\scriptscriptstyle 0.03}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.02$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$6$ ${\ {}^{\scriptscriptstyle 6}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
6$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$2.48$ ${\ {}^{\scriptscriptstyle 3.04}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
2.18$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$65$ ${\ {}^{\scriptscriptstyle 79}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
42$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$1.44$ ${\ {}^{\scriptscriptstyle 1.71}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
1.40$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$448$ ${\ {}^{\scriptscriptstyle 470}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
395$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$\mathrm{iEQL_{\text{motor}}}$ $0.02$ ${\ {}^{\scriptscriptstyle 0.02}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.02$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$6$ ${\ {}^{\scriptscriptstyle 6}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
6$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$3.17$ ${\ {}^{\scriptscriptstyle 6.16}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
2.90$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$32$ ${\ {}^{\scriptscriptstyle 48}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
15$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$1.60$ ${\ {}^{\scriptscriptstyle 1.79}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
1.39$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$339$ ${\ {}^{\scriptscriptstyle 386}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
170$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>
Results are shown in table 6. $\mathrm{iEQL}$ and EQL÷ outperform the genetic
algorithm (GA) and the Gaussian Process (GP) even on the training data and
provide very accurate predictions on the test set. The GP does not capture the
underlying relationship as can be seen at the large test RMSE. A closer look
at the results of the genetic algorithm revealed that the internal selection
criteria of $\mathrm{GA}$ led to competitive results: $\mathrm{GA1}$ with
$0.47\,\,^{\scriptscriptstyle 159.11}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.02$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>and $19\,\,^{\scriptscriptstyle 23}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
9$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>and $\mathrm{GA2}$ with $0.02\,\,^{\scriptscriptstyle 0.16}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
0.02$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>and $17\,\,^{\scriptscriptstyle 19}_{{\hbox
to0.0pt{\leavevmode\hbox{\set@color$\scriptscriptstyle
15$}\hss}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>$\mathrm{iEQL}$ needs even fewer parameters than EQL÷. This is due to the use
of copy units, which allow to reusing features from previous layers directly.
We emphasize that, out of the large set of possible functions, the
$\mathrm{iEQL}$ reliably extracted a simple quadratic equation that suitably
describes the dataset. To highlight the simplicity of the $\mathrm{iEQL}$’s
result, we state in the following the structure of the selected equation,
which is the same for all 5 experiments

$\displaystyle y=$ $\displaystyle
w_{1}I_{Q}+w_{2}I_{D}+w_{5}(w_{3}I_{D}+b_{1})^{2}+w_{6}(w_{4}I_{Q}+b_{2})^{2}+b_{3}\,.$
(16)

Weights are indicated with $w$ and bias with $b$. All selected equations can
be simplified to the same equation
$y=0.4I_{D}^{2}-1.12I_{D}+0.41I_{Q}^{2}+1.13I_{Q}-0.31$. Numbers were rounded
to three figures and input and output dimensions were anonymized. In addition,
we emphasize that the algorithm learned that only two of the input variables,
namely $[I_{D},\,I_{Q}]$, are relevant to the output value, while the others
have negligible influence.

###  4.4 Torque Model of an Internal Combustion Engine

Table 7: Data properties for the torque model of an internal combustion
engine.

Quantity Unit Min. Max. Description Type $\phi_{\text{ex}}$ ∘Crank $-20$ $20$
exhaust camshaft input $\phi_{\text{in}}$ ∘Crank $-4$ $36$ intake camshaft
input $r_{\text{l}}$ $\%$ $13$ $86$ relative load input $\phi_{\text{ign}}$
∘Crank $-27$ $61$ ignition angle input $n_{\text{eng}}$ rpm $597$ $6000$
engine speed input $M_{\text{eng}}$ Nm $-38$ $261$ engine torque output

A second task from the industrial domain is the modeling of the torque of a
combustion engine in dependence of the control parameters. Ongoing
improvements in power train technologies require more and more precise models
for control strategies. This is commonly correlated with a complex structure
of the models as well as the large number of calibration parameters. For
utilization of those models in embedded systems it is essential to have fast
development cycles, i.e. an automated adaptation to a new system, and low
production costs. The unit costs grow with the computational needs on an
embedded controller, thus computationally efficient at evaluation time is
important. In practice, this leads often to a trade-off between accuracy of
the model and its complexity.

In this dataset the engine torque is measured depending on exhaust camshaft
$\phi_{\text{ex}}$, intake camshaft $\phi_{\text{in}}$, relative load
$r_{\text{l}}$, ignition angle $\phi_{\text{ign}}$ and engine speed
$n_{\text{eng}}$. The data was measured at stationary operation points. It
contains $1775$ data points measured at variations of the quantities listed in
table 7 within their range of operation. The dataset is split into $80\%$
training and $20\%$ testing data and $10\%$ of the train data set is used for
validation.

$\displaystyle 10^{1}$$\displaystyle 10^{2}$$\displaystyle
10^{3}$$\displaystyle 2$$\displaystyle 4$$\displaystyle 10$complexityRMSE
[Nm]validations.eqeq-1 (a) pareto plot

$\displaystyle-10$$\displaystyle-5$$\displaystyle 0$$\displaystyle
5$$\displaystyle 10$residuals [Nm]s.eqeq-1 (b) histogram of residuals

Figure 5: Equation learning with $\mathrm{iEQL_{\text{motor}}}$ on a real
world dataset for torque of a combustion engine. Here, we show the run with
median performance with the V${}^{\text{int}}$-S criteria. Panel (a) shows the
pareto plot of the $\mathrm{iEQL}$ and (b) shows the histogram of residuals of
the equation (s.eq) selected with the V${}^{\text{int}}$-S criteria and a more
complex equation (eq-1).

Instance selection can be done with the V${}^{\text{int}}$-S criteria or
solely based on validation performance as long as the equation is
computationally efficient at evaluation time. We present both solutions in
table 6 alongside with the number of active parameters. The GP is supposed to
perform best on this dataset, but the EQL-networks, selected with
V${}^{\text{int}}$, perform similar or even better. With about 170-470
parameters they are fast to compute and small enough to fit on embedded
systems. The EQL÷ finds the best test RMSE results. This might indicate that
for this dataset it is sufficient to deal with
$\\{\sin,\cos,*,\text{identity}\\}$ units. Despite its larger expressive
power, $\mathrm{iEQL}$ can compete with the EQL÷ architecture. $\mathrm{GA}$
is not able to learn equations that compete on the performance level, but it
can identify simple equations with about 33 parameters. Yet,
$\mathrm{iEQL_{\text{motor}}}$ also finds equations of similar complexity, but
with a twice as good median test RMSE. Table 6 shows a noticeable relation
between the performance of the models and their complexity for
$\mathrm{iEQL}$, $\mathrm{iEQL_{\text{motor}}}$ and EQL÷. The more complex
the model the betters its test RMSE. This effect occurs also in figure 5(a),
which shows the pareto plot for the $\mathrm{iEQL_{\text{motor}}}$ experiment
(median performance and V${}^{\text{int}}$-S criteria). The histograms of
residuals of the selected equation alongside with the best performing equation
show no significant distribution shifts, see figure 5(b).

We analyze the equations of median performance and the equations with the
smallest number of parameters selected with the V${}^{\text{int}}$-S criteria
for $\mathrm{iEQL}$ and $\mathrm{iEQL_{\text{motor}}}$ in appendix A.5. Their
test RMSE is not as good as a GP, but their structure is interpretable. The
motor specific complexity factors lead to equations of mainly polynomial
structure with higher degree than without motor specific complexity factors.
The latter combine polynomials with nonlinear units like $\cos,\exp,\log$ or
division units.Using more nonlinear units reduces the degree of the used
polynomials. This is in agreement with the analysis of relative frequencies of
atomic units shown in figure 4(a), which clearly shows that the use of motor
specific complexity factors leads to a reduction of the use of
$\\{\log,\cos,\exp\\}$ and division units and an increase of $x^{2}$ and
multiplication units.

##  5 Conclusion

We introduced the informed equation learner $\mathrm{iEQL}$. It aims to learn
compact and interpretable models in form of concise mathematical equations
from data. Our method differs from previous work in two key aspects: First, it
is able to handle functional units with singularities in all hidden layers.
Secondly, it can incorporate expert knowledge about the underlying system. To
do so, prior knowledge is encoded in a complexity cost for individual atomic
units. We evaluated our method by finding compact functional expressions for
four different use cases: The $\mathrm{iEQL}$ performs well on established
reference datasets. Customized priors can be applied to prefer certain
functional terms in an ambiguous dataset supporting several possible
explanations. Finally, we evaluated the algorithm on two real industrial
problems. In both settings, the underlying relations of the system were
unknown a priori, and the $\mathrm{iEQL}$ was able to extract simple,
explainable models describing the respective output variable. Depending on the
needs it also provides a high-performance solution with low computational
cost.

Funding:

Philipp Hennig and Georg Martius are members of the Machine Learning Cluster
of Excellence, funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy – EXC number 2064/1 –
Project number 390727645. We acknowledge the support from the German Federal
Ministry of Education and Research (BMBF) through the Tübingen AI Center (FKZ:
01IS18039B).

Authors’ contributions:

All authors contributed to the conception and design of the study. Code
preparation, data collection and analysis were performed by MW and supervised
by GM and PH. AJ has supported the work as domain expert. The manuscript was
drafted by MW and edited and reviewed by all authors.

## References

  * Akaike [1998] H. Akaike.  Information theory and an extension of the maximum likelihood principle.  In _Selected papers of hirotugu akaike_ , pages 199–213. Springer, 1998.
  * Ballentine [1970] L. E. Ballentine.  The Statistical Interpretation of Quantum Mechanics.  _Reviews of Modern Physics_ , 42(4):358–381, Oct. 1970.
  * Buchner et al. [2020] J. S. Buchner, S. Boblest, P. Engel, A. Junginger, and H. Ulmer.  An artificial-intelligence-based method to automatically create interpretable models from data targeting embedded control applications.  _IFAC_ , 2020.  21th IFAC World Congress.
  * Cranmer [2020] M. Cranmer.  Pysr: Fast & parallelized symbolic regression in python/julia, Sept. 2020\.  URL http://doi.org/10.5281/zenodo.4041459.
  * Dubčáková [2011] R. Dubčáková.  Eureqa: software review.  _Genetic Programming and Evolvable Machines_ , 12(2):173–178, Jun 2011.  ISSN 1573-7632.
  * Feynman et al. [1965] R. P. Feynman, R. B. Leighton, and M. Sands.  The feynman lectures on physics; vol. i.  _American Journal of Physics_ , 33(9):750–752, 1965.
  * García et al. [1989] C. E. García, D. M. Prett, and M. Morari.  Model predictive control: Theory and practice—a survey.  _Automatica_ , 25(3):335 – 348, 1989.  ISSN 0005-1098.
  * He et al. [2016] K. He, X. Zhang, S. Ren, and J. Sun.  Identity mappings in deep residual networks.  In _European Conference on Computer Vision (ECCV 2016)_ , pages 630–645. Springer International Publishing, 2016.
  * Hoffmann et al. [2015] S. Hoffmann, M. Schrott, T. Huber, and T. Kruse.  Modellbasierte Methoden zur Applikation moderner Verbrennungsmotoren.  _MTZ_ , 76(4):46–51, 2015.
  * Huang et al. [2016] G. Huang, Z. Liu, L. van der Maaten, and K. Q. Weinberger.  Densely connected convolutional networks.  _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 2261–2269, 2016.
  * Jin et al. [2019] Y. Jin, W. Fu, J. Kang, J. Guo, and J. Guo.  Bayesian symbolic regression.  _arXiv preprint arXiv:1910.08892_ , 2019.
  * Kim et al. [2020] S. Kim, P. Y. Lu, S. Mukherjee, M. Gilbert, L. Jing, V. Čeperić, and M. Soljačić.  Integration of neural network-based symbolic regression in deep learning for scientific discovery.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2020\.
  * Kingma and Ba [2015] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization. international conference on learning representations (2015), 2015.
  * Kusner et al. [2017] M. J. Kusner, B. Paige, and J. M. Hernández-Lobato.  Grammar variational autoencoder.  In _Proceedings of the 34th International Conference on Machine Learning-Volume 70_ , pages 1945–1954. JMLR. org, 2017.
  * Lample and Charton [2020] G. Lample and F. Charton.  Deep learning for symbolic mathematics.  In _International Conference on Learning Representations_ , 2020.
  * Langdon et al. [2008] W. B. Langdon, R. Poli, N. F. McPhee, and J. R. Koza.  Genetic programming: An introduction and tutorial, with a survey of techniques and applications.  In _Computational intelligence: A compendium_ , pages 927–1028. Springer, 2008.
  * Li et al. [2019] L. Li, M. Fan, R. Singh, and P. Riley.  Neural-guided symbolic regression with semantic prior.  _arXiv preprint arXiv:1901.07714_ , 2019.
  * Lin et al. [2020] S.-C. Lin, G. Martius, and M. Oettel.  Analytical classical density functionals from an equation learning network.  _The Journal of Chemical Physics_ , 152(2):021102, 2020.
  * Long et al. [2019] Z. Long, Y. Lu, and B. Dong.  PDE-Net 2.0: Learning PDEs from data with a numeric-symbolic hybrid deep network.  _Journal of Comp. Physics_ , 399:108925, 2019.  ISSN 0021-9991.
  * Louizos et al. [2018] C. Louizos, M. Welling, and D. P. Kingma.  Learning sparse neural networks through $l_{0}$ regularization.  In _International Conference on Learning Representations_ , 2018.
  * MacKay [1992] D. J. MacKay.  Bayesian interpolation.  _Neural computation_ , 4(3):415–447, 1992.
  * Maddison et al. [2016] C. J. Maddison, A. Mnih, and Y. W. Teh.  The concrete distribution: A continuous relaxation of discrete random variables, 2016.
  * Martius and Lampert [2016] G. Martius and C. H. Lampert.  Extrapolation and learning equations, 2016.  arXiv https://arxiv.org/abs/1610.02995.
  * McRee [2010] R. K. McRee.  Symbolic regression using nearest neighbor indexing.  In _Proceedings of the 12th Annual Conference Companion on Genetic and Evolutionary Computation_ , GECCO ’10, pages 1983–1990, New York, NY, USA, 2010. ACM.  ISBN 978-1-4503-0073-5.
  * Mitchell and Beauchamp [1988] T. J. Mitchell and J. J. Beauchamp.  Bayesian variable selection in linear regression.  _Journal of the American Statistical Association_ , 83(404):1023–1032, 1988.
  * Petersen [2019] B. K. Petersen.  Deep symbolic regression: Recovering mathematical expressions from data via policy gradients.  _arXiv preprint arXiv:1912.04871_ , 2019.
  * Sahoo et al. [2018] S. S. Sahoo, C. H. Lampert, and G. Martius.  Learning equations for extrapolation and control.  In _International Conference on Machine Learning, ICML 2018_ , volume 80, pages 4442–4450. PMLR, 2018.
  * Schmidt and Lipson [2009] M. Schmidt and H. Lipson.  Distilling free-form natural laws from experimental data.  _Science_ , 324(5923):81–85, 2009.  ISSN 0036-8075.
  * Schwarz et al. [1978] G. Schwarz et al.  Estimating the dimension of a model.  _The annals of statistics_ , 6(2):461–464, 1978\.
  * Searson et al. [2010] D. P. Searson, D. E. Leahy, and M. J. Willis.  GPTIPS: an open source genetic programming toolbox for multigene symbolic regression.  In _International Multiconference of Engineers and Computer scientists_ , volume 1, pages 77–80. IMECS Hong Kong, 2010.
  * Tibshirani [1996] R. Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):267–288, 1996.
  * Trujillo et al. [2016] L. Trujillo, L. Muñoz, E. Galván-López, and S. Silva.  neat genetic programming: Controlling bloat naturally.  _Information Sciences_ , 333:21–43, 2016.  ISSN 0020-0255.
  * Udrescu and Tegmark [2020] S.-M. Udrescu and M. Tegmark.  Ai feynman: A physics-inspired method for symbolic regression.  _Science Advances_ , 6(16):eaay2631, 2020.
  * Williams [1992] R. J. Williams.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  _Machine learning_ , 8(3-4):229–256, 1992.
  * Zaremba et al. [2014] W. Zaremba, K. Kurach, and R. Fergus.  Learning to discover efficient mathematical identities.  In _Advances in Neural Information Processing Systems_ , page 1278, 2014.

##  Appendix A Appendix

(a) S0 (b) S1 (c) S2 (d) S3 (e) S4 (f) S5 (g) S6 (h) A1 (i) A2 (j) A3 (k) A4

Figure 6: Two dimensional slices through the input space of simulated
equations S0–S6 and A1–A4 from table 2. The black rectangle indicates the
train domain.

###  A.1 $\mathrm{iEQL}$ Training

All experiments are performed using an $\mathrm{iEQL}$ with four hidden
layers. Each hidden layer has
$\\{\cos,\,\exp{\;},\,\log,\,\sqrt{\phantom{x}},\,x^{2},\,*,\,/\\}$ as atomic
units, and each atomic unit is applied four times in each layer. We prohibit
the following combinations $\cos(\cos)$, $\cos(\exp)$, $\exp(\exp)$,
$\log(\log)$. The $\mathrm{iEQL}$ has thus $6405$ learnable weights. Training
is executed in two phases. In phase 1 we train for $T_{1}=2000$ epochs without
regularization. This phase overcomes bad initialization and assures that the
model is close to a minimum when pruning starts. In phase 2 we train for
$T_{2}=10000$ epochs with regularization strength $\lambda$. Just for the
combustion engine dataset, which has about 5 to 6 times fewer data points, we
increased the number of epochs by a factor of 8. We skipped train phase 1 for
the ambiguous dataset since the $\mathrm{iEQL}$ is already sufficiently close
to a minimum. After each epoch, an intrinsic penalty epoch is calculated with
$100$ randomly sampled data points from the test domain (without labels). The
maximum desired output value is set to $B=10$ with $\eta=1$. We use the Adam
optimizer [Kingma and Ba, 2015] with learning rate $0.001$, moving average
$\beta_{1}=0.4$ and ${\epsilon}=10^{-8}$ for numerical stability. The initial
dropout rate of the Bernoulli gates is set to $0.5$ and the domain/bound
penalty strength is set to $\delta=1$. Further information on dropout rate is
given in appendix 3.4. Since the optimal regularization scale parameter
$\lambda$ is not known in advance, we train several models with different
regularization strengths $\lambda=10^{k}$ where $k$ is in the range from
$-5.0$ to $0.0$ with 78 equally spaced steps. This results in 78 equations
with different complexities and root mean squared errors (RMSE).

#### Baselines:

We calculate the mean predictor (MP) on the train set and a multi-layer
perceptron (MLP) with $\tanh$ activation functions and five hidden layers with
50 neurons each. It is trained with batch size $100$ for $5000$ epochs and the
Adam optimizer with a learning rate of 0.001 and $\beta_{1}=0.9$. A grid
search on the learning rate revealed that it is pretty robust in the range
$[0.001,0.0001]$. We select the best validation model to avoid overfitting.

We compare to EQL÷, a state-of-the-art method from Sahoo et al. [2018] with
atomic unit types $\\{\sin,\cos,*,\text{identity}\\}$ and division in the
final layer. It is sufficient to use the hyperparameters proposed in Sahoo et
al. [2018] since the datasets on which we compare have similar properties to
the ones in Sahoo et al. [2018]. So we applied, Adam optimizer with learning
rate $0.001$ and $\epsilon=10^{-4}$, mini-batch size $20$, domain penalty
$\eta=10$ with $B=10$, and 10 atomic units per type in each layer. The number
of total epochs is given by $T=(L-1)\cdot 10,000$, where $L$ is the number of
hidden layers. Just for the combustion engine dataset, which has about 5 to 6
times fewer data points, we had to increase the number of total epochs to
$T=(L-1)\cdot 80,000$. We perform model selection amongst the following
parameters: regularization strengths $\lambda=10^{k}$ where $k$ is in the
range from $-6.0$ to $-3.5$ with 26 equally spaced steps and
$L\in\\{2,3,4\\}$. This results in 78 equations with different complexities
and root mean squared errors (RMSE).

Another baseline is (PySR, Cranmer [2020]) a genetic algorithm for symbolic
regression with hyperparameters shown in table 8 for two different
configurations ($\mathrm{GA1}$, $\mathrm{GA2}$). Both settings turned out to
perform well on different datasets. We had to restrict the size of the
datasets to 1000 data points in order to avoid exploding memory size. We do
not compare to Eureqa, the current state-of-the-art tool for symbolic
regression [Dubčáková, 2011], since it has become proprietary and was merged
into an online service.

Table 8: Tuned hyperparameters for the genetic algorithm [Cranmer, 2020].  | $\mathrm{GA1}$ | $\mathrm{GA2}$
---|---|---
niterations | 40 | 10
npop | 1000 | 1000
populations | 60 | 30
binary operators | $\\{\pm,*,/\\}$ | $\\{\pm,*,/\\}$
unary operators | $\\{\cos,\exp,\log\\}$ | $\\{\cos,x^{2},\exp,\log\\}$
maxsize | 40 | 40
parsimony | 0 | 0
warmupMaxsizeBy | 0.5 | 0
useFrequency | False | True
annealing | False | False
optimizer algorithm | BFGS | BFGS
optimizer iterations | 100 | 100
procs | 30 | 30

For the real world datasets a Gaussian Process (GP) is calculated with ASCMO
[Hoffmann et al., 2015], a standard tool from the engineering domain.

###  A.2 Power Loss of an Electric Machine

Here, we present details on the preparation of training and test dataset.
First, $20\%$ of the whole dataset is used for testing and excluded from
training..Then, the train domain is restricted to $80\%$ of its range of
operation ($[a,b]\rightarrow[0.8a,0.8b]$ for each dimension). Rotor
temperature $T_{\text{rot}}$ is not restricted, since it consists just of
three different operation points. Further, $10\%$ of the train dataset is used
for validation. Restricting the train domain assures that the test dataset
contains samples from an extrapolation domain as well as samples from the
train domain.

###  A.3 Sparsity Inducing Method

In this section we outline the $L_{0}$ regularization by Louizos et al. [2018]
in more detail, which is used in our sparsity inducing method. For simplicity
the complexity factors $c_{u}$ are neglected, since it is straightforward to
integrate them. The objective function we want to optimize is given by

$\displaystyle\mathcal{L}=\frac{1}{N}\sum\limits_{i=1}^{N}\lVert{\bm{y}}_{i}-\mathrm{iEQL}({\bm{x}}_{i},{\bm{W}})\rVert_{2}^{2}+\lambda\lVert
W\rVert_{0}\,.$ (17)

The $L_{0}$-norm counts the number of non-zero weights, which total number is
denoted by $\lvert W\rvert$. Specific choices of the scale parameter $\lambda$
refer to well-known model selection criteria as the Bayesian Information
Criterion (BIC, Schwarz et al. [1978]) and the Akaike Information Criterion
(AIC, Akaike [1998]). Each weight of the $\mathrm{iEQL}$ is multiplied by a
non-negative stochastic Bernoulli distributed gate

$\displaystyle w_{i}=\tilde{w}_{i}\cdot g_{j}$ $\displaystyle
q(g_{j}\mid\pi_{j})=\mathrm{Ber}(\pi_{j})\,.$ (18)

The dropout rate of each weight is thus given by $1-\pi_{j}$. The expectation
of a weight being relevant is given by $\pi_{j}$. Those gates learn
collectively which weights are relevant. The expected objective is then

$\displaystyle\mathcal{L}=\mathbb{E}_{q({\bm{g}}\mid\pi)}\left[\frac{1}{N}\sum\limits_{i=1}^{N}\lVert{\bm{y}}_{i}-\mathrm{iEQL}({\bm{x}}_{i},{\bm{W}})\rVert^{2}_{2}\right]+\lambda\sum\limits_{j\leq\lvert
W\rvert}\pi_{j}\,.$ (19)

This objective can be described as a special case of a variational bound over
weights with spike and slab [Mitchell and Beauchamp, 1988] priors and
approximate posteriors. We refer the interested reader to appendix A of
Louizos et al. [2018] for further details.

####  A.3.1 Reparameterization Trick

In order to apply gradient based optimization the objective is smoothed.
Therefore, the gates $z$ are hard-sigmoid rectifications $b(\cdot)$ of a
continuous random variable $s$ with variables $\phi$

$\displaystyle s\sim q(s\mid\phi)\,,$ $\displaystyle z=b(s)\,,$
$\displaystyle\text{with: }b(\cdot)=\min(1,\max(0,\cdot)).$ (20)

The expected loss is then

$\displaystyle\mathcal{L}\\!$
$\displaystyle=\\!\mathbb{E}_{q(s\mid\phi)}\\!\left[\\!\frac{1}{N}\\!\sum\limits_{i=1}^{N}\lVert{\bm{y}}_{i}\\!-\\!\mathrm{iEQL}\\!({\bm{x}}_{i};{\bm{W}}\\!\odot\\!b(s)\rVert^{2}_{2}\\!\right]\\!+\\!\lambda\\!\sum\limits_{j\leq\lvert
W\rvert}\\!(1\\!-\\!Q(s_{j}\\!\leq\\!0\\!\mid\\!\phi_{j}))$ (21)

with the cumulative distribution (CDF) $Q(\cdot)$ of $s$. This penalizes the
probability of a gate being non-zero. By choosing a suitable continuous
distribution $q(s)$ we can apply the reparameterization trick with parameter
free noise distribution $\epsilon$ and a deterministic and differentiable
transformation $f(\cdot)$

$\displaystyle\mathcal{L}\\!$
$\displaystyle=\\!\mathbb{E}_{p(\epsilon)}\\!\left[\\!\frac{1}{N}\\!\sum\limits_{i=1}^{N}\lVert{\bm{y}}_{i}\\!-\\!\mathrm{iEQL}\\!({\bm{x}}_{i};{\bm{W}}\\!\odot\\!b(f(\phi,\epsilon))\rVert^{2}_{2}\\!\right]\\!+\\!\lambda\\!\sum\limits_{j\leq\lvert
W\rvert}\\!(1\\!-\\!Q(s_{j}\\!\leq\\!0\\!\mid\\!\phi_{j}))$ (22)

The final reparameterization is given by the hard concrete distribution, which
is outlined in the following.

####  A.3.2 Hard Concrete Distribution

$s$ is a binary concrete random variable (Maddison et al. [2016], He et al.
[2016]) by choice, distributed over $(0,1)$ with probability density
$q(s\mid\phi)$ and cumulative density $Q_{\beta}(s\mid\phi)$. The distribution
has two parameters. Its location is denoted by $\log\alpha$ and the degree of
approximation is controlled by $\beta$. A sampling method for the stretched
version for the interval $(\gamma,\zeta)$ is given by

$\displaystyle u$ $\displaystyle\sim U(0,1)$ (23) $\displaystyle s$
$\displaystyle=\mathrm{Sigmoid}\big{(}(\log
u-\log(1-u)+\log\alpha)/\beta\big{)}$ (24) $\displaystyle\bar{s}$
$\displaystyle=s(\zeta-\gamma)+\gamma$ (25) $\displaystyle z$
$\displaystyle=b(\bar{s})$ (26)

with more details

$\displaystyle q_{s}(s\mid\phi)$ $\displaystyle=\frac{\beta\alpha
s^{-\beta-1}(1-s)^{-\beta-1}}{(\alpha s^{-\beta}+(1-s)^{-\beta/2})^{2}}$ (27)
$\displaystyle Q_{s}(s\mid\phi)$ $\displaystyle=\mathrm{Sigmoid}\left((\log
s-\log(1-s))\beta-\log\alpha\right)$ (28) $\displaystyle
q_{\bar{s}}(\bar{s}\mid\phi)$
$\displaystyle=\frac{1}{\lvert\zeta-\gamma\rvert}q_{s}\left(\frac{\bar{s}-\gamma}{\zeta-\gamma}\mid\phi\right)$
(29) $\displaystyle Q_{\bar{s}}(\bar{s}\mid\phi)$
$\displaystyle=Q_{s}\left(\frac{\bar{s}-\gamma}{\zeta-\gamma}\mid\phi\right)$
(30)

and finally the hard sigmoid $z$ we obtain the distribution

$\displaystyle q(z\mid\phi)=$ $\displaystyle
Q_{\bar{s}}(0\mid\phi)\delta(z)+\big{(}1-Q_{\bar{s}}(1\mid\phi)\big{)}\delta(z-1)$
$\displaystyle+\big{(}Q_{\bar{s}}(1\mid\phi)-Q_{\bar{s}}(0\mid\phi)\big{)}\,q_{\bar{s}}(z\mid\bar{s}\in(0,1),\phi)$
(31)

Thus the complexity loss is given by

$\displaystyle\mathcal{L}_{\text{C}}$
$\displaystyle=\mathrm{Sigmoid}\left((\log\alpha-\beta\log\frac{-\gamma}{\zeta}\right)$
(32)

and for final estimation of the parameters under a hard concrete gate:

$\displaystyle\hat{z}=\min\left(1,\max\left(0,\mathrm{Sigmoid}\left(\log\alpha\right)(\zeta-\gamma)+\gamma\right)\right)$
(33)

In our implementation we use the same hyperparameters as Louizos et al. [2018]

$\displaystyle\zeta=1.1,$ $\displaystyle\gamma=-0.1,$
$\displaystyle\beta=2/3\,.$ (34)

$\displaystyle 0.0$$\displaystyle 0.2$$\displaystyle 0.4$$\displaystyle
0.6$$\displaystyle 0.8$$\displaystyle 1.0$$\displaystyle 0$$\displaystyle
1$$\displaystyle 2$$\displaystyle 3$$\displaystyle z$$\displaystyle
p(z)$concretehard concrete

$\displaystyle-10$$\displaystyle-5$$\displaystyle 0$$\displaystyle
5$$\displaystyle 10$$\displaystyle 0.0$$\displaystyle 0.25$$\displaystyle
0.5$$\displaystyle 0.75$$\displaystyle
1.0$$\displaystyle\log\alpha$$\displaystyle z$concretehard
concrete$\mathbb{E}_{p(\epsilon)}\\!\left[\text{concrete}\right]$$\mathbb{E}_{p(\epsilon)}\\!\left[\text{hard
concrete}\right]$

Figure 7: The probability density function for the concrete gate (equation 28)
and hard concrete gate (equation 30) are shown in (a) and their expectation
(equation 24 and equation 26 respectively) for $10,000$ samples are shown in
(b). We added the concrete and hard concreate without noise as well. All plots
are produced with the hyperparameters mentioned in equation 34.

###  A.4 Model Selection Criteria

Both criteria V${}^{\text{int\&ex}}$ and V${}^{\text{int}}$-S perform
similarly good on equations (S0–S6), but the V${}^{\text{int\&ex}}$ methods
performs better on the more sophisticated equations (A1–A4). However, it
requires additional extrapolation data points, which might not be available in
real world datasets. Therefore, we apply the V${}^{\text{int}}$-S method on
the two real world datasets.

To illustrate the importance of the equation selection process we evaluate the
selected equation alongside with a simpler and a more complex solution on the
S0 dataset, see figure 8. If the selected equation is too complex its
validation loss is at noise level, but it is not able to extrapolate, as seen
e.g. in the left panel of figure 8. If the selected equation is too simple
validation and extrapolation loss are not at noise level.

$\displaystyle-2$$\displaystyle-1$$\displaystyle 0$$\displaystyle
1$$\displaystyle 2$$\displaystyle-2$$\displaystyle-1$$\displaystyle
0$$\displaystyle 1$$\displaystyle 2$$\displaystyle
x_{1}=x_{2}=x_{3}=x_{4}=x$$\displaystyle y$gteq-1s.eqeq-2

$\displaystyle 10$$\displaystyle 30$$\displaystyle 100$$\displaystyle
10^{-2}$$\displaystyle 10^{-1}$complexityRMSEvalidation

Figure 8: The panel in the left column displays a slice through the input
space of equation S0. The blue shading indicates the train domain.The selected
equation (dotted line) was determined with the V${}^{\text{int\&ex}}$
selection criteria. Additionally, a too simple (dashdotted line) and too
complex equation (dashed line) are plotted for comparison. The right column
shows the corresponding pareto plot. Each point corresponds to an independent
run of the $\mathrm{iEQL}$ with a different complexity regularization. The
three equations shown on the left are marked with corresponding vertical
lines.

###  A.5 Different Combustion engine Equations

We present the equations of median performance and the equations with the
smallest number of parameters selected with the V${}^{\text{int}}$-S criteria
for $\mathrm{iEQL}$ and $\mathrm{iEQL_{\text{motor}}}$. An overview of all
four analyzed equations on the combustion engine dataset is shown in table 9.

Table 9: Overview of analysed $\mathrm{iEQL}$ expressions on the combustion engine dataset. The RMSE on the test dataset with the number of active parameters is shown.  $\mathrm{iEQL}$ | RMSE [Nm] | #parameters | equation
---|---|---|---
motor simple | 6.16 | 15 | 35
motor median | 3.17 | 32 | 36
plain simple | 3.04 | 42 | 37
plain median | 2.48 | 79 | 38

Numbers were rounded to three figures and input and output dimensions were
anonymized. The motor specific complexity factors lead to equations of mainly
polynomial structure with higher degree than without motor specific complexity
factors. The latter combine polynomials with simple nonlinear units like
$\cos,\exp,\log$ or division units. Using more nonlinear units reduces the
degree of the used polynomials.

#####  $\mathrm{iEQL_{\text{motor}}}$ simple:

  * •

compact notation of a polynomial of degree 8

$y=2.48x_{2}+0.66\left(0.3x_{2}+0.79\right)\left(-0.8x_{1}+1.86x_{3}-0.08\right)-0.52\left(-2.13x_{3}+0.57\left(0.05-1.19x_{3}\right)^{2}-0.62\left(-0.89x_{1}+0.85\left(0.46x_{1}-0.32\right)\left(1.28x_{3}+0.34x_{5}+0.71\right)+0.42\right)^{2}+0.56\right)^{2}+0.94$
(35)

#####  $\mathrm{iEQL_{\text{motor}}}$ median:

  * •

compact notation of a polynomial of degree 7

  * •

$\cos$ units in final layer with previous polynomial terms as arguments.

$y=-0.3x_{1}+2.32x_{2}+0.27x_{3}+0.6c_{2}+0.3c_{3}+0.35\cos{\left(1.12x_{3}+1.42c_{2}+0.59c_{3}-0.23\right)}+0.51$
(36)

substitutions:

$\displaystyle c_{1}=$
$\displaystyle\left(0.05-0.69x_{1}\right)\left(1.0x_{1}+0.64\right)$
$\displaystyle c_{2}=$
$\displaystyle\left(-0.95x_{2}-0.69\right)\left(-1.06x_{3}+0.96\left(-0.21x_{1}+0.66x_{3}-1.01\right)^{2}-1.2\right)$
$\displaystyle c_{3}=$
$\displaystyle\left(-0.74x_{3}-0.76\left(0.6x_{1}+0.73c_{1}+0.47\left(0.76x_{2}+0.19x_{4}-0.19x_{5}-0.6\right)^{2}-0.9\right)^{2}+1.36\right)$
$\displaystyle\left(2.0x_{3}+0.88\left(-0.73x_{1}-1.17c_{1}+0.26\right)\left(0.6x_{4}-0.72x_{5}-0.46\right)-0.31\right)$

#####  $\mathrm{iEQL}$ simple:

  * •

polynomial of degree 4 with simple $\cos$ and $\exp$ units

  * •

$\cos$ units with previous polynomial terms as arguments.

$y=1.49x_{2}+0.53\left(0.94x_{2}+0.8\right)\left(1.89x_{3}-0.29\cos{\left(4.79x_{3}+2.83\right)}+0.46\right)+0.81e^{-0.27x_{1}+0.18x_{2}}+0.94e^{-0.32x_{1}+0.25x_{2}+0.55x_{3}}-0.05c_{1}+0.22\cos{\left(1.45x_{3}-1.02c_{2}+0.07\right)}-0.29\cos{\left(0.66x_{2}-1.13c_{2}-0.74c_{3}+2.73\right)}-0.23\cos{\left(-1.21x_{1}+1.99x_{3}-2.26\left(-0.94x_{2}-0.52\right)\left(0.97x_{3}+0.02\right)-0.97c_{3}+3.71\right)}-1.71$
(37)

substitutions:

$\displaystyle c_{1}=$ $\displaystyle\cos{\left(4.0x_{1}+2.38\right)}$
$\displaystyle c_{2}=$
$\displaystyle\left(-0.49x_{1}+1.1x_{3}+0.76\left(1.02x_{1}-0.02\right)^{2}-0.85\right)^{2}$
$\displaystyle c_{3}=$
$\displaystyle\left(-0.41x_{2}-0.53x_{3}-0.29x_{4}+0.34x_{5}+0.59\left(0.07-0.97x_{1}\right)\left(0.28-1.02x_{2}\right)+0.23c_{1}+0.6\right)^{2}$

#####  $\mathrm{iEQL}$ median:

  * •

polynomial of degree 4 with simple $\cos,\log$ and division units

  * •

division units occur also in intermediate layers

  * •

$\cos$ units with previous polynomial terms as arguments

$y=2.0x_{2}+0.6x_{3}+0.06x_{4}-0.05x_{5}-0.49\left(0.83x_{2}+0.58\right)\left(-0.81x_{3}-0.42\right)+\frac{0.39\left(0.99x_{2}+0.79\right)}{6.75\left(1.77c_{1}-0.75\right)^{2}+1.96}-0.24\left(-0.77x_{2}-0.36c_{1}-0.42\right)\left(1.11x_{3}+0.32c_{2}-0.89\left(-1.1x_{3}+0.73c_{1}+0.41\right)^{2}-0.47\cos{\left(5.98x_{3}+3.27\right)}+0.09\right)-0.29\left(-0.78x_{2}+0.35\left(1.77c_{1}-0.75\right)^{2}+0.35\cos{\left(1.72x_{2}-2.6\right)}-0.67\right)^{2}-0.28c_{3}+0.18\log{\left(1.08-0.88x_{1}\right)}-0.2\cos{\left(-1.21x_{2}+1.28c_{3}+4.61\right)}-0.15\cos{\left(1.95x_{2}+1.58x_{3}+1.0c_{4}+3.51\right)}+0.18\cos{\left(0.77x_{1}-1.71x_{2}-0.76c_{4}+0.71c_{3}+1.51\right)}+0.74-\frac{0.39\left(0.85x_{1}-0.58\right)}{3.29-2.85x_{2}}$
(38)

substitutions:

$\displaystyle c_{1}=$
$\displaystyle\cos{\left(2.07x_{1}+0.21x_{2}-0.02\right)}$ $\displaystyle
c_{2}=$
$\displaystyle\left(0.33x_{1}+1.09c_{1}+0.51\right)\left(-0.74x_{4}+0.61x_{5}+1.02c_{1}-0.62\right)$
$\displaystyle c_{3}=$
$\displaystyle\left(0.29x_{1}-1.65x_{3}-0.09x_{4}+0.09x_{5}-0.27\left(-0.69x_{1}-0.14\right)^{2}+0.35\left(0.68x_{1}-1.27x_{3}+0.3\right)^{2}+0.4\right)^{2}$
$\displaystyle c_{4}=$
$\displaystyle\left(-0.43x_{4}+0.36x_{5}+0.8+\frac{0.86\left(-1.22x_{2}-0.94x_{3}-0.66\right)}{2.11x_{1}+2.35}\right)$
$\displaystyle\left(0.6x_{2}+1.28x_{3}+0.89x_{4}-1.18x_{5}+0.77c_{2}-0.05\right)$