11institutetext: Gianforte School of Computing, Montana State University, USA
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Univariate Skeleton Prediction in Multivariate Systems Using
Transformers††thanks: Paper accepted at ECML PKDD 2024

Giorgio Morales  John Sheppard

###### Abstract

Symbolic regression (SR) methods attempt to learn mathematical expressions
that approximate the behavior of an observed system. However, when dealing
with multivariate systems, they often fail to identify the functional form
that explains the relationship between each variable and the system’s
response. To begin to address this, we propose an explainable neural SR method
that generates univariate symbolic skeletons that aim to explain how each
variable influences the system’s response. By analyzing multiple sets of data
generated artificially, where one input variable varies while others are
fixed, relationships are modeled separately for each input variable. The
response of such artificial data sets is estimated using a regression neural
network (NN). Finally, the multiple sets of input–response pairs are processed
by a pre-trained Multi-Set Transformer that solves a problem we termed Multi-
Set Skeleton Prediction and outputs a univariate symbolic skeleton. Thus, such
skeletons represent explanations of the function approximated by the
regression NN. Experimental results demonstrate that this method learns
skeleton expressions matching the underlying functions and outperforms two GP-
based and two neural SR methods.

###### Keywords:

Symbolic regression Transformer networks Symbolic skeletons Multivariate
regression Explainable artificial intelligence.

##  1 Introduction

Symbolic regression (SR) aims to identify mathematical equations or symbolic
expressions that capture the underlying relationships and dynamics of the
studied phenomena [19]. One of the main advantages of the expressions learned
by an SR model is that they are interpretable by humans and allow for the
identification of cause-effect relationships between the inputs and outputs of
a system [11]. These techniques not only capture the behavior of empirical
data into analytical equations but also reduce the computational complexity
during the inference phase, and have more powerful extrapolation ability than
black-box models [23].

A limitation of existing SR approaches lies in their primary focus on
minimizing prediction errors rather than distilling the underlying equations
that govern the system dynamics [4]. Consequently, the generated equations may
exhibit high complexity, effectively approximating the observed data but
failing to correspond to the underlying equations [19]. This limitation poses
challenges when performing out-of-sample inference, where generalization may
be ineffective.

SR methods are often based on population-based algorithms, especially genetic
programming (GP) [25]. Nevertheless, a notable drawback of GP-based SR methods
is that they suffer from slow computation. Their computational inefficiency is
related to the inherent complexity of the search space and the expensive
iterative calls to numerical optimization routines after each generation [25].
In addition, these methods do not consider past experiences, as they require
learning each problem from scratch. As such, the obtained models do not
benefit from additional data or insights from different equations, hindering
their capacity for improvement and limiting their generalization capabilities
[4, 18].

In response to these limitations, recent advancements have seen the emergence
of neural SR methods as a promising alternative. These methods may utilize
general-purpose pre-trained models that generate symbolic expressions through
a straightforward forward pass and, possibly, a single call to a numerical
optimization routine [4, 18, 31]. Thus, neural SR methods offer a substantial
time speedup compared to GP-based approaches. However, there remains a gap in
prediction accuracy between the two paradigms [18].

Although recent neural SR methods may achieve low prediction errors when
dealing with multivariate systems, they struggle to represent the dependency
of the system’s response and each independent variable correctly. As a first
step in addressing this, our proposed method aims to generate univariate
symbolic skeletons that estimate functional relationships for each variable
with respect to the response. A symbolic skeleton expression is an abstract
representation of a mathematical expression that captures its structural form
without setting specific numerical values. In other words, we propose to
generate univariate skeletons that constitute mathematical “explanations”
describing the interaction between each independent variable and the system’s
response.

Given a multivariate regression problem that can be expressed in terms of a
mathematical equation, our method identifies univariate symbolic skeleton
expressions for each variable. The skeleton prediction process begins by
training a regression model, such as a neural network (NN), that approximates
the system’s function. This regression model is used to estimate the response
of multiple randomly generated sets of points where only the variable of
interest is allowed to vary while the remaining variables are held fixed. We
then define a novel Multi-Set Transformer model, which is pre-trained on a
large dataset of synthetic symbolic expressions, that receives as inputs the
sets generated by the NN and is used to find the univariate skeletons.

####  1.0.1 Contributions

We hypothesize that our method will generate univariate skeletons that are
more similar to those corresponding to the underlying equations in comparison
to other SR methods. The generated skeletons represent explanations of how
each variable is related to the system’s response. Thus, from an
explainability standpoint, producing more faithful univariate skeletons means
that we can offer better explanations about the system’s behavior. In
addition, the generated skeletons may be taken as building blocks to be used
to estimate the overall function of the system. Our main contributions are
summarized as:

  * •

We introduce a method that learns univariate symbolic skeletons that explain
the functional form between each independent variable and the system’s
response using regression NNs and pre-trained transformer models.

  * •

We introduce a new SR problem called Multi-Set symbolic skeleton prediction
(MSSP). It receives multiple sets of input–response pairs, where all sets
correspond to the same functional form but use different equation constants,
and outputs a common skeleton expression.

  * •

We present a novel transformer model called “Multi-Set Transformer” to solve
the MSSP problem. The model is pre-trained on a large dataset of synthetic
symbolic expressions.

##  2 Related Work

SR is commonly tackled using GP-based methods. Such methods evolve a
population of tree-like individuals using operations like selection,
crossover, and mutation to improve their fitness over multiple generations.
Each individual represents a symbolic expression that maps the inputs and the
output, and its fitness function determines how well it fits the data set
being modeled. Variations of this approach attempt to design improved
operators and fitness functions to reduce the complexity of the search [13,
25, 28].

Two significant challenges of GP for SR are code growth (i.e., bloat) and the
huge search space [1]. Code growth leads to computationally expensive
evolution of large programs and hinders their generalization ability [8]. The
huge search space is attributed to the variability in program lengths
permitted during the evolutionary process. This allows for the generation of
multiple solution trees representing mathematically equivalent functions [9].
However, GP tends to generate a greater number of large solution trees
compared to smaller ones [21]. The high complexity of large solutions entails
poor generalization performance. In most GP-based methods, the program fitness
is determined by its overall output and not by the intermediate outputs of its
subexpressions. Thus, a program’s subexpressions are optimized indirectly.
Arnaldo et al. [2] pointed out that the indirect optimization approach may
result in suboptimal fitness optimization because it does not focus on
evolving and learning suitable building blocks.

Furthermore, Martius and Lampert [23], Sahoo et al. [27], and Werner et al.
[33] proposed the use of NN architectures whose units represent conventional
mathematical operators. The objective is to reduce the prediction error and
gradually prune irrelevant parts of the network until a simple equation can be
extracted from the network. However, similar to the case of GP-based methods,
the learned expressions tend to be large and, as a consequence, difficult to
explain. Another limitation is that each problem is learned from scratch.

As a recent alternative, a few methods based on large language models have
been proposed. Biggio et al. [4] introduced the use of pre-trained transformer
models for SR. A dataset of multivariate equations is generated to pre-train a
model based on a Set Transformer [22], and this transformer acts as a general-
purpose model to predict the symbolic skeleton from a corresponding set of
input-output pairs. The skeleton’s constants are then fit using a non-convex
optimizer, such as the BFGS algorithm [12].

Kamienny et al. [18] pointed out that the loss function minimized by the BFGS
algorithm can be highly non-convex and the correct constants of the skeletons
are not guaranteed to be found. They avoided performing skeleton prediction as
an intermediary step and proposed an end-to-end (E2E) transformer model that
estimates the full mathematical expression directly. The learned constants can
then be refined using a non-convex optimizer. Experiments showed that E2E
performs better than previous NN-based methods. Since previous methods have
shown scalability issues when dealing with equations having many variables,
Chu et al.[6] proposed a method that decomposes multi-variable SR into a
sequence of single-variable SR problems, combined in a bottom-up manner. The
process involves learning a data generator using NNs from observed data,
generating variable-specific samples with controlled input variables, applying
single-variable SR, and iteratively adding variables until completion.

##  3 Multivariate Skeleton Prediction

Consider a system whose response $y\in\mathbb{R}$ is sensitive to variables
$\textbf{x}=\\{x_{1},\dots,x_{t}\\}$ ($\textbf{x}\in\mathbb{R}^{t}$). The
underlying function that maps the feature value space and the response value
space is denoted as $f:\mathcal{X}^{t}\rightarrow\mathcal{Y}$, such that
$y=f(\mathbf{x})=f(x_{1},\dots,x_{t})$. In addition, let $\kappa(\cdot)$
represent a skeleton function that replaces the numerical constants of a given
symbolic expression by placeholders $c_{i}$; e.g.,
$\kappa(3x^{2}+e^{2x}-4)=c_{1}\,x^{2}+e^{c_{2}\,x}+c_{3}$. We assume $f$ can
be expressed as a mathematical expression with unary (e.g., sin, cos, and log)
and binary operators (e.g. $+$, $-$, $*$, and $/$). We aim to obtain the
univariate skeletons $\hat{\mathbf{e}}(x_{1}),\dots,\hat{\mathbf{e}}(x_{t})$
that describe the functional form between each variable and the system’s
response. In the following, we lay out the steps proposed in this method.

###  3.1 Neural Network Training

The underlying function $f$ can be approximated based on observed data using
any regression model. Let
$\textbf{X}=\\{\textbf{x}_{1},\dots,\textbf{x}_{N_{R}}\\}$ be a data set with
$N_{R}$ samples, where each sample is denoted as
$\textbf{x}_{j}=\\{x_{j,1},\dots,x_{j,t}\\}$, and
$\textbf{y}=\\{y_{1},\dots,y_{N_{R}}\\}$ is the set of corresponding target
observations. A NN regression model, whose compute function is denoted as
$\hat{f}(\cdot;\boldsymbol{\theta}_{NN})$ ($\boldsymbol{\theta}_{NN}$
represents the weights of the network), is constructed to capture the
association between X and y. Thus, a target estimate for a given input
$\textbf{x}_{j}$ is computed as
$\hat{y}_{j}=\hat{f}(\textbf{x}_{j},\boldsymbol{\theta}_{NN})$ or, simply,
$\hat{y}_{j}=\hat{f}(\textbf{x}_{j})$. The parameters
$\boldsymbol{\theta}_{NN}$ of the network $\hat{y}$ are obtained by minimizing
the mean squared error of the predictions; that is,
$\boldsymbol{\theta}_{NN}^{*}=\;{\text{argmin}_{\boldsymbol{\theta}_{NN}}}\
\frac{1}{{N_{R}}}\sum_{j=1}^{N_{R}}(\hat{y}_{j}-y_{j})^{2}$. Note that we
selected an NN to generate the function $\hat{f}$ due to its ease of training
and high accuracy; however, other regression methods could be applied.

###  3.2 Multi-Set Symbolic Skeleton Prediction

Figure 3.1: $x_{1}$ vs. $y$ curves when $x_{2}=4.45$, $0.2$, and $1.13$.

We tackle the SR problem by decomposing it into single-variable sub-problems.
To do so, we deviate from the Symbolic Skeleton Prediction (SSP) problem that
has been explored previously by existing research [4, 26, 31]. We explain the
rationale behind this decision through the use of an example: Consider the
function $y={x_{1}}/{\log(x_{1}^{2}+x_{2})}$. We analyze the relationship
between $x_{1}$ and $y$. Similar to the work presented by Chu et al. [6], the
variables that are not being analyzed are held fixed. Fig. 3.1 depicts that
the function’s behavior differs when using three different $x_{2}$ values. As
such, conventional SSP solvers may produce expressions with different
functional forms depending on the fixed values of the variables not currently
under analysis.

What is more, the fixed values of the remaining variables could project the
function into a space where the functional form is not easily identifiable,
which could be worsened due to the limited range of values that the analyzed
variable can take. Therefore, the SSP problem would benefit from the injection
of additional context data. When analyzing the variable $x_{v}$, we could use
multiple sets of input–response pairs, each of which is constructed using
different fixed values for the variables $\mathbf{x}\setminus\\{x_{v}\\}$. The
key idea is to process the information from these sets simultaneously to
produce a skeleton that is common to all input sets. We refer to this problem
as Multi-Set Symbolic Skeleton Prediction (MSSP).

More formally, assume we are given a data set of $N_{R}$ input–response pairs
$(\mathbf{X},\mathbf{y})$ where $\mathbf{X}\in\mathbb{R}^{N_{R}\times t}$ and
$\mathbf{y}\in\mathbb{R}^{N_{R}\times 1}$. Suppose we analyze the relationship
between the $v$-th variable of the system, $x_{v}$ (i.e., $v\in[1,\dots,t]$),
and the response variable $y$. We construct a collection of $N_{S}$ sets,
denoted as $\mathbf{D}=\\{\mathbf{D}^{(1)},\dots,\mathbf{D}^{(N_{S})}\\}$.
Each set $\mathbf{D}^{(s)}$ comprises $n$ input–response pairs such that
$\mathbf{D}^{(s)}=(\mathbf{X}_{v}^{(s)},f(\mathbf{X}^{(s)}))=(\mathbf{X}_{v}^{(s)},\mathbf{y}^{(s)})$,
where $\mathbf{X}^{(s)}\in\mathbb{R}^{n\times t}$,
$\mathbf{y}^{(s)}\in\mathbb{R}^{n}$, and $\mathbf{X}_{v}^{(s)}$ denotes the
$v$-th column of $\mathbf{X}^{(s)}$ (i.e., the data corresponding to the
$x_{v}$ variable). $\mathbf{X}^{(s)}$ is constructed in such a way that each
of the variables in $\mathbf{x}\setminus\\{x_{v}\\}$ are assigned random
values and held fixed for all samples.

$\mathbf{X}^{(s)}$ can be constructed by selecting $n$ samples from
$\mathbf{X}$ that meet the criterion that the variable $x_{v}$ is allowed to
vary while the variables in $\mathbf{x}\setminus\\{x_{v}\\}$ are fixed. In
this case, the set $\mathbf{y}^{(s)}$ is known as it is the corresponding set
of observed response values of $\mathbf{X}^{(s)}$. Conversely, if the set
$\mathbf{X}$ is not large enough, $\mathbf{X}^{(s)}$ can be generated
artificially using the regression model, and the corresponding set of response
values would be obtained as $\mathbf{y}^{(s)}=f(\mathbf{X}^{(s)})$. Since $f$
is not known a priori, the black-box model $\hat{f}$ that has been trained to
approximate $f$ (i.e., $\hat{f}(\mathbf{x})\approx f(\mathbf{x})$) is used to
estimate the response of $\mathbf{X}^{(s)}$ (cf. Section 3.1).

Since the variables in $\mathbf{x}\setminus\\{x_{v}\\}$ have been fixed to
constant values to construct each set $\mathbf{D}^{(s)}$, the underlying
function that explains the relationship between $\mathbf{X}_{v}^{(s)}$ and
$\mathbf{y}^{(s)}$ could be expressed solely in terms of variable $x_{v}$. As
such, the underlying function of the $s$-th set is denoted by
$f^{(s)}(x_{v})$. It is important to note that functions
$f^{(1)}(x_{v}),\dots,f^{(N_{S})}(x_{v})$ have been derived from the same
function $f(\mathbf{x})$ and only differ in their coefficient values due to
the selection of different values for the variables in
$\mathbf{x}\setminus\\{x_{v}\\}$ for each constructed set. Thus, if we apply
the skeleton function $\kappa(\cdot)$ to functions
$f^{(1)}(x_{v}),\dots,f^{(N_{S})}(x_{v})$, they should produce the same target
symbolic skeleton $\mathbf{e}(x_{v})$, in which the constant values have been
replaced by placeholders; i.e.,
$\mathbf{e}(x_{v})=\kappa\left(f^{(1)}(x_{v})\right)=\dots=\kappa\left(f^{(N_{S})}(x_{v})\right)$.

Then, the collection $\mathbf{D}$ is fed as an input to the MSSP problem. The
objective is to generate an skeleton $\hat{\mathbf{e}}(x_{v})$ that
characterizes the functional form of all input sets, and approximates the
target skeleton $\mathbf{e}(x_{v})$; i.e.,
$\hat{\mathbf{e}}(x_{v})\approx\mathbf{e}(x_{v})$. For the sake of generality,
we define the MSSP problem as follows:

###### Definition 1

The input of the Multi-Set Symbolic Skeleton Prediction (MSSP) problem
consists of a collection of $N_{S}$ sets, denoted as
$\mathbf{D}=\\{\mathbf{D}^{(1)},\dots,\mathbf{D}^{(N_{S})}\\}$. Each set
$\mathbf{D}^{(s)}$ comprises $n$ input–response pairs such that
$\mathbf{D}^{(s)}=(\mathbf{X}_{v}^{(s)},\mathbf{y}^{(s)})$, where
$\mathbf{X}_{v}^{(s)}\in\mathbb{R}^{n}$ and
$\mathbf{y}^{(s)}\in\mathbb{R}^{n}$. The underlying function of the $s$-th set
is denoted by $f^{(s)}$; i.e.,
$\mathbf{y}^{(s)}=f^{(s)}(\mathbf{X}_{v}^{(s)})$. The underlying functions of
all input sets are assumed to share a common unknown symbolic skeleton,
denoted as $\mathbf{e}$. Thus, the objective of the MSSP problem is to
generate a symbolic skeleton $\hat{\mathbf{e}}\approx\mathbf{e}$ that
characterizes the functional form of all input sets $\mathbf{D}$.

####  3.2.1 Multi-Set Transformer

Our method for solving the MSSP problem draws inspiration from the Set
Transformer [22], which is an attention-based neural network that is derived
from the transformer model [32]. The Set Transformer is designed for set-input
problems and, as such, it is capable of processing input sets of varying sizes
and exhibits permutation invariance. We propose a Multi-Set Transformer that
presents modifications to address the limitations of the Set Transformer and
adapt it to the specific requirements of the proposed research. Details about
the differences between the Set Transformer and our Multi-Set Transformer are
given in the Supplementary Materials111Please visit https://github.com/NISL-
MSU/MultiSetSR for the supplementary file.

Let the $s$-th input set be denoted as
$\mathbf{D}^{(s)}=(\mathbf{X}_{v}^{(s)},\mathbf{y}^{(s)})=\\{(x_{v,i}^{(s)},y_{i}^{(s)})\\}_{i=1}^{n}$.
The first step involves arranging $\mathbf{D}^{(s)}$ in a manner analogous to
the input structure of the Set Transformer, which consists of a matrix where
each row represents a $d_{in}$-dimensional element of the input set.
$\mathbf{S}^{(s)}\in\mathbb{R}^{n\times d_{in}}$ denotes the $s$-th input of
our proposed Multi-Set Transformer such that its $i$-th row,
$\mathbf{s}_{i}^{(s)}$, consists of the concatenation of the input value
$x_{v,i}^{(s)}$ and its corresponding output $y_{i}^{(s)}$; i.e.,
$\mathbf{s}_{i}^{(s)}=[x_{v,i}^{(s)},y_{i}^{(s)}]$. Hence, $\mathbf{S}^{(s)}$
is defined as a matrix with $d_{in}=2$ columns.

Our Multi-Set Transformer comprises two primary components: an encoder and a
decoder. The encoder maps the information of all input sets into a unique
latent representation $\mathbf{Z}$. To do so, an encoder stack $\phi$
transforms each input set $\mathbf{S}^{(s)}$ into a latent representation
$\mathbf{z}^{(s)}\in\mathbb{R}^{d}$ (where $d$ is context vector length or the
“embedding size”) individually. Our encoder, denoted as $\Phi$, comprises the
use of the encoder stack $\phi$ to generate $N_{S}$ individual encodings
$\mathbf{z}^{(1)},\dots,\mathbf{z}^{(N_{S})}$, which are then aggregated into
a unique latent representation:

$\mathbf{Z}=\Phi(\mathbf{S}^{(1)},\dots,\mathbf{S}^{(N_{S})},\boldsymbol{\theta}_{e})=\rho(\phi(\mathbf{S}^{(1)},\boldsymbol{\theta}_{e}),\dots,\phi(\mathbf{S}^{(N_{S})},\boldsymbol{\theta}_{e})),$

where $\rho(\cdot)$ is a pooling function, $\boldsymbol{\theta}_{e}$
represents the trainable weights of the encoder stack, and
$\mathbf{z}^{(s)}=\phi\left(\mathbf{S}^{(N_{S})},\boldsymbol{\theta}_{e}\right)$.
We define $\phi$ as a stack of $\ell$ induced set attention blocks (ISABs)
[22] to encode high-order interactions among the elements of an input set in a
permutation-invariant way. ISAB is a multi-head attention layer with low
complexity ($\mathcal{O}(mn)$) that acts as a universal approximator of
permutation invariant functions. Furthermore, we include a pooling by multi-
head attention (PMA) layer [22] in $\phi$ to aggregate the features extracted
by the ISAB blocks, whose dimensionality is $n\times d$, into a single
$d$-dimensional latent vector. Finally, the function $\rho(\cdot)$ that is
used to aggregate the latent representations $\mathbf{z}^{(s)}$ is implemented
using an additional PMA layer.

The decoder $\psi$ generates sequences conditioned on the representation
$\mathbf{Z}$ generated by $\Phi$. This objective is aligned with that of the
transformer decoder [32] and, thus, the same architecture is used for our
model. $\psi$ consists of a stack of $M$ identical blocks, each of which is
composed of three main layers: a multi-head self-attention layer, an
encoder–decoder attention layer, and a position-wise feedforward network. Let
$\hat{\mathbf{e}}=\\{\hat{e}_{1},\dots,\hat{e}_{N_{out}}\\}$ denote the output
sequence produced by the model, which represents the skeleton as a sequence of
indexed tokens in prefix notation. Each token in this sequence is transformed
into a numerical index according to a pre-defined vocabulary that contains all
unique symbols. During inference, each element $\hat{e}_{i}$ is generated
auto-regressively; that is, $\psi$ produces a probability distribution over
the elements of the vocabulary:
$\sigma\left(\psi\left(\mathbf{Z},\boldsymbol{\theta}_{d}|\hat{e}_{1},\dots,\hat{e}_{i-1}\right)\right)=P\left(\hat{e}_{i}|\hat{e}_{1},\dots,\hat{e}_{i-1},\mathbf{Z}\right)$,
where $\boldsymbol{\theta}_{d}$ represents the weights of $\psi$. This
distribution is obtained by applying a softmax function $\sigma(\cdot)$ to the
decoder’s output. $\hat{e}_{i}$ is thus selected as the token with the highest
probability.

####  3.2.2 Multi-Set Transformer Training

We train our model on a large dataset of artificially generated MSSP problems.
Let $\mathbf{D}^{b}=\left\\{\mathbf{D_{1}},\dots,\mathbf{D}_{B}\right\\}$
denote a training batch with $B$ samples, each of which represents a
collection of $N_{S}$ input sets; i.e.,
$\mathbf{D}_{j}=\\{\mathbf{D}_{j}^{(1)},\dots,\mathbf{D}_{j}^{(N_{S})}\\}$
($j\in[1,\dots,B]$). In addition,
$\mathbf{E}^{b}=\\{\mathbf{e}_{1},\dots,\mathbf{e}_{B}\\}$ is the
corresponding set of target skeletons, each of which represents a sequence of
variable length; i.e.,
$\mathbf{e}_{j}=\left\\{e_{j,1}\dots,e_{j,N_{j}}\right\\}$ and
$N_{j}=|\mathbf{e}_{j}|$. The function computed by the model is denoted as
$g(\cdot)$, and $\boldsymbol{\Theta}$ denotes its weights. Note that
$\boldsymbol{\Theta}=[\boldsymbol{\theta}_{e},\boldsymbol{\theta}_{d}]$
contains the weights of the encoder and the decoder stacks. Given an input set
collection $\mathbf{D}_{j}$, $g(\mathbf{D}_{j},\boldsymbol{\Theta})$ computes
the estimated skeleton $\hat{\mathbf{e}}_{j}$ with length $N_{out,\,j}$.
$g(\cdot)$ is trained to generate estimated skeletons so that
$\hat{\mathbf{e}}_{j}\approx\mathbf{e}_{j}$.

During training, the model is provided with past elements of the target
skeleton sequence as inputs for generating subsequent tokens. Hence, the
normalized probability distribution produced by the decoder $\psi$ for the
$i$-th element of the predicted skeleton sequence of the $j$-th sample would
be expressed as
$P\left(\hat{e}_{j,i}|e_{j,1},\dots,e_{j,i-1},\mathbf{Z}_{j}\right)$. In cases
where $\hat{\mathbf{e}}_{j}$ and $\mathbf{e}_{j}$ differ in length, we use
padding and masking to ensure that the loss is only calculated for valid
tokens. Thus, our optimization objective is defined as the cross-entropy loss
between the padded target and predicted skeleton sequences, which is
calculated using their corresponding probability distributions over the set of
possible tokens:

$\mathcal{L}=-\frac{1}{B}\sum_{j=1}^{B}\sum_{i=1}^{T}\omega_{j,i}P(e_{j,i})\log
P(\hat{e}_{j,i}|e_{j,1},\dots,e_{j,i-1},\mathbf{Z}_{j}),$

where $w_{j,i}$ is 0 for padding positions and 1 otherwise, and $T_{j}$ is the
length of the $j$-th sequence. Hence, the optimization problem is expressed
as: $\boldsymbol{\Theta}=\text{argmin}_{\boldsymbol{\Theta}}\ \mathcal{L}.$

We generate a training dataset of expressions stored in prefix notation. The
generation process is explained in Section B of the Supplementary Materials1.
Our method differs from the generation method used by Biggio et al. [4] and
Lample et al. [20]. Their approach generates several expressions that contain
binary operators exclusively and, thus, are limited to simple expressions such
as $\hat{e}(x)=c_{1}x$. Conversely, we use a method that builds the expression
tree recursively in preorder, which allows us to enforce certain conditions
and constraints effectively. That is, we forbid certain combinations of
operators and set a maximum limit on the nesting depth of unary operators
within each other.

The training routine of the Multi-Set Transformer $g$ takes as inputs the set
of pre-generated expressions, denoted as $\mathbf{Q}$. The skeleton
$\mathbf{e}_{j}$ and its associated data collection $\mathbf{D}_{j}$,
corresponding to the $j$-th expression in $\mathbf{Q}$, are generated using a
function called $\texttt{generateSets}(\mathbf{Q}[j],N_{S},n)$ (Algorithm 1).
The estimated skeletons $\hat{\mathbf{e}}_{j}$ are obtained processing the
input $\mathbf{D}_{j}$ and target $\mathbf{e}_{j}$ through the network $g$
using a teacher forcing strategy.
$\mathcal{L}(\mathbf{E}^{B},\hat{\mathbf{E}}^{B})$ represents the loss
function while update($g,L$) encompasses the conventional backpropagation and
stochastic gradient descent processes used to update the weights of model $g$.

Algorithm 1 describes the function generateSets, which takes as input a pre-
generated expression $\mathbf{ex}$. Function
$\texttt{getConstants}(\mathbf{ex})$ returns $\mathbf{c}$, the list of
constant placeholders in $\mathbf{ex}$, and $n_{c}=|\mathbf{c}|$. Function
$\texttt{selectConstants}(\mathbf{ex},\mathbf{c},n_{f})$ retrieves an
expression with $n_{f}$ constant placeholders selected randomly ($2\leq
n_{f}\leq n_{c}$). Then, we generate the underlying functions $f^{(s)}$
corresponding to each of the $N_{S}$ sets using the function
$\texttt{sampleConstants}(\mathbf{ex})$, which samples the values of each
constant in $\mathbf{ex}$ from a uniform distribution $\mathcal{U}(-10,10)$
independently. The $n$ input points $\mathbf{X}_{v}^{(s)}$ are then sampled
using the function $\texttt{sampleSupport}(n)$ from a uniform distribution
$\mathcal{U}(-x^{\text{limit}},x^{\text{limit}})$, where
$x^{\text{limit}}\sim\mathcal{U}(1,10)$ (i.e., the domain range is randomly
sampled each time). The function
$\texttt{avoidNaNs}(\mathbf{X}_{v}^{(s)},f^{(s)})$ may modify the coefficients
in $f^{(s)}$ or sample additional $\mathbf{X}_{v}^{(s)}$ values to avoid
undefined results (see Section B of the Supplementary Materials1).

1:function generateSets($\mathbf{ex},N_{S},n$)

2: $\mathbf{c},n_{c}\leftarrow\texttt{getConstants}(\mathbf{ex})$

3:
$\mathbf{ex}\leftarrow\texttt{selectConstants}(\mathbf{ex},\mathbf{c},n_{f}=\texttt{randInt}(2,n_{c}))$

4:
$\mathbf{D},s,\mathbf{e}\leftarrow[\,],1,\kappa(\texttt{sampleConstants}(\mathbf{ex}))$

5: while $s\leq N_{S}$ do

6: $\mathbf{X}_{v}^{(s)}\leftarrow\texttt{sampleSupport}(n)$

7:
$\mathbf{X}_{v}^{(s)},f^{(s)},\\_\\_\leftarrow\texttt{avoidNaNs}(\mathbf{X}_{v}^{(s)},\texttt{sampleConstants}(\mathbf{ex}))$

8: if $s>1$ and $\kappa(f^{(s)})\neq\mathbf{e}$ then

9: continue $\triangleright$ Verify that all sets correspond to the same
skeleton

10: $\mathbf{y}^{(s)}\leftarrow f^{(s)}(\mathbf{X}_{v}^{(s)})$

11: $\mathbf{D}.\texttt{append}((\mathbf{X}_{v}^{(s)},\mathbf{y}^{(s)}))$

12: $s\leftarrow s+1$

13: return $\mathbf{D},\mathbf{e}$

Algorithm 1 Multi-Set Transformer Data Generation

####  3.2.3 Univariate Symbolic Skeleton Prediction

Suppose we are currently analyzing the $v$-th variable, $x_{v}$. As an MSSP,
we generate $N_{S}$ artificial sets of points
$\\{\tilde{\mathbf{X}}^{(1)},\dots,\tilde{\mathbf{X}}^{(N_{S})}\\}$ where the
variable $x_{v}$ is allowed to vary while the variables
$\mathbf{x}\setminus\\{x_{v}\\}$ are fixed to random values. Specifically, the
$s$-th artificial set is denoted as
$\tilde{\mathbf{X}}^{(s)}=\\{\tilde{\mathbf{x}}^{(s)}_{1},\dots,\tilde{\mathbf{x}}^{(s)}_{n}\\}$.
The value of the $v$-th dimension of the $j$-th sample is obtained by sampling
from the distribution $\mathcal{U}(x_{v}^{\min},x_{v}^{\max})$ whose lower and
upper bounds, $x_{v}^{\min}$ and $x_{v}^{\max}$, respectively, are calculated
from the observed data. The values assigned to the remaining dimensions are
sampled independently using similar uniform distributions; however, the same
value is shared across all samples (i.e.,
$\tilde{\mathbf{x}}^{(s)}_{1,k}=\tilde{\mathbf{x}}^{(s)}_{2,k}=\dots=\tilde{\mathbf{x}}^{(s)}_{n,k}$,
$\forall k\in[1,\dots,t]$ and $k\neq v$).

Our approach requires processing data sets where the variables not under
analysis are held constant. However, it is not always possible to find subsets
of data that meet this condition for all variables in real-world datasets.
Even when such subsets exist, they may not be sufficiently large. Thus, we
generate sets of data $\tilde{\mathbf{X}}^{(s)}$ with the desired behavior and
estimate their response using a regression model that has been trained on
observed data. Then, our method derives univariate skeletons based on multiple
sets of input-estimated response pairs. It is important to use a prediction
model that learns a function $\hat{f}$ that is as close as possible to $f$ so
that it accurately estimates how the real system would respond to the
artificial inputs in $\tilde{\mathbf{X}}^{(s)}$. As a consequence, our
analysis can be regarded as an explainability method that generates univariate
symbolic skeletons as explanations of the function approximated by the
regression model.

The response of the samples in $\tilde{\mathbf{X}}^{(s)}$,
$\tilde{\mathbf{y}}^{(s)}$, is estimated using the network $\hat{f}$ as
$\tilde{\mathbf{y}}^{(s)}=\hat{f}(\tilde{\mathbf{X}}^{(s)})$. In addition,
$\tilde{\mathbf{D}}_{v}^{(s)}=(\tilde{\mathbf{X}}_{v}^{(s)},\tilde{\mathbf{y}}^{(s)})$
denotes the set of $n$ input–response pairs and is used to analyze the
relationship between the system’s response and the $v$-th variable. The
collection of $N_{S}$ sets
$\tilde{\mathbf{D}}_{v}=\\{\tilde{\mathbf{D}}_{v}^{(1)},\dots,\tilde{\mathbf{D}}_{v}^{(N_{S})}\\}$
is then fed into the pre-trained Multi-Set Transformer $g$ so that the
estimated skeleton obtained for variable $x_{v}$ is calculated as
$\tilde{\mathbf{e}}(x_{v})=g(\tilde{\mathbf{D}}_{v},\boldsymbol{\Theta})$.
This process is repeated for all the variables of the system to obtain their
corresponding symbolic skeleton expressions with respect to the system’s
response.

###  3.3 Performance Evaluation

Typically, SR methods are evaluated on benchmark datasets such as the Feynman
SR Benchmark (FSRB) [30]. Their performance is compared based on the mean
squared error (MSE) achieved by the learned expressions on a subset of the
available data. Note that a function that was learned by minimizing the
prediction error with respect to the system’s response does not necessarily
correspond to the underlying functional form; e.g.,
$f(x)=\cos(-\frac{x}{10})^{2}$ and $\hat{f}(x)=-0.0093x^{2}+0.9983$ produce an
MSE lower than $4\times 10^{-5}$ when $x\in[-5,5]$. As such, conventional
evaluation approaches may not align with our objectives.

Instead of minimizing prediction error, our method produces a set of
univariate skeletons aiming to describe the functional form between each
variable and the system’s response. We argue that an expression with minimum
prediction error can be produced as a consequence of identifying the correct
functional form. To the best of our knowledge, no previous work has addressed
the problem of testing how well the learned expression’s functional form
(i.e., its skeleton) matches the system’s underlying functional form. Thus, we
present a method to test the similarity between the underlying skeleton
corresponding to the variable $x_{v}$, represented as
$\textbf{e}(x_{v})=\kappa(f(\textbf{x}),x_{v})$, and the estimated skeleton
$\hat{\mathbf{e}}(x_{v})$. The function $\kappa(\cdot,x_{v})$ replaces the
numerical constants of a given symbolic expression with the placeholder
labeled with the suffix $c$. In addition, it considers the remaining variables
$\mathbf{x}\setminus x_{v}$ as constants as they are irrelevant when
describing the functional form between $x_{v}$ and the system’s response. For
example, if $f(\mathbf{x})=3x_{1}^{2}+\sqrt{x_{2}+1}/e^{2\,x_{3}}$, then
$\kappa(f(\textbf{x}),x_{1})=c_{1}\,x_{1}^{2}+c_{2}$.

We assign random numerical values to the coefficients of skeleton
$\mathbf{e}(x_{v})$ using the sampleConstants routine (Algorithm 1) in order
to obtain a function $f_{\text{target}}(x_{v})$. Let
$f_{\text{est}}(x_{v})=\texttt{setConstants}(\hat{\mathbf{e}}(x_{v}),\mathbf{c})$
denote the function obtained when replacing the $n_{c}$ constant placeholders
in $\hat{\mathbf{e}}(x_{v})$ with the numerical values in a given set
$\mathbf{c}=[c_{1},\dots,c_{n_{c}}]$. If the functional form of
$\hat{\mathbf{e}}(x_{v})$ is mathematically equivalent to that of
$\mathbf{e}(x_{v})$, then there exists a set of values $\mathbf{c}$ so that
the difference between $f_{\text{target}}(x_{v})$ and $f_{\text{est}}(x_{v})$
is zero for all values of $x_{v}$. The optimal set $\mathbf{c}^{*}$ is found
as:
$\mathbf{c}^{*}=\text{argmin}_{\mathbf{c}}\sum_{x_{v}\in\mathbf{X}_{v}^{\text{test}}}|f_{\text{target}}(x_{v})-f_{\text{est}}(x_{v})|$,
where $\mathbf{X}_{v}^{\text{test}}$ is a test set of $N_{\text{test}}$
elements whose elements are drawn from a distribution
$\mathcal{U}(2\,x_{v}^{\min},2\,x_{v}^{\max})$. Note that the domain of
$\mathbf{X}_{v}^{\text{test}}$ is larger than that used for training (i.e.,
$[x_{v}^{\min},x_{v}^{\max}]$). If the estimated skeleton matches the system’s
underlying functional form, it should do so regardless of the variable domain.

The coefficient fitting problem is solved using a simple genetic algorithm
(GA) [15]. The individuals of our GA are arrays of $n_{c}$ elements that
represent potential $\mathbf{c}$ sets. The population size is set to 500 and
the optimization algorithm is stopped when the change of the objective
function after 20 generations is less than $10^{-5}$. We utilize tournament
selection, binomial crossover, and generational replacement. We selected this
configuration as it demonstrated effective optimization results across all
experiments utilized in this work.

If $\hat{\mathbf{e}}(x_{v})$ and $\mathbf{e}(x_{v})$ are not equivalent, the
error
$r=\sum_{x_{v}\in\mathbf{X}_{v}^{\text{test}}}|f_{\text{target}}(x_{v})-\texttt{setConstants}(\hat{\mathbf{e}}(x_{v}),\mathbf{c}^{*})|$
is greater than 0. We use $r$ as a performance metric that indicates the
closeness between $\hat{\mathbf{e}}(x_{v})$ and $\mathbf{e}(x_{v})$ given the
sampled values of the constants of $\mathbf{e}(x_{v})$. Note that if
$\hat{\mathbf{e}}(x_{v})$ and $\mathbf{e}(x_{v})$ are similar, the error $r$
should be low regardless of the sampled values of the constants of
$\mathbf{e}(x_{v})$. Thus, for the sake of generality, we repeat this process
30 times; that is, we sample 30 different $f_{\text{target}}(x_{v})$ functions
and solve 30 optimization problems. Finally, we report the mean and the
standard deviation of the 30 resulting error metrics.

##  4 Experimental Results

A training dataset222The code and datasets are available at
https://github.com/NISL-MSU/MultiSetSR consisting of one million pre-generated
expressions ($|\mathbf{Q}|=10^{6}$) has been created to train the Multi-Set
Transformer. These expressions allow up to one nested operation and contain a
maximum of five unary operators. We also generated an independent validation
set consisting of $10^{5}$ expressions. For the model architecture, due to the
high computational expense associated with training a single model, we used a
one-factor-at-a-time approach to choose the following hyperparameters: the
number of ISAB encoder blocks $\ell=3$, the number of decoder blocks $M=5$, an
embedding size $d=512$, and the number of heads $h=8$. In addition, we set the
number of input sets to $N_{S}=10$ and the number of input-response pairs in
each input set to $n=3000$. Future work will examine the impact of altering
the values of $N_{S}$ and $n$ on the final performance.

###  4.1 Experiments with Synthetic Data

Table 4.1: Equations used for experiments

Eq. | Underlying equation | Reference | Domain range
---|---|---|---
E1 | $(3.0375x_{1}x_{2}+5.5\sin(9/4(x_{1}-2/3)(x_{2}-2/3)))/5$ | [16] | $[-5,5]^{2}$
E2 | $5.5+(1-x_{1}/4)^{2}+\sqrt{x_{2}+10}\sin(x_{3}/5)$ | — | $[-10,10]^{2}$
E3 | $(1.5e^{1.5x_{1}}+5\cos(3x_{2}))/10$ | [16] | $[-5,5]^{2}$
E4 | $((1-x_{1})^{2}+(1-x_{3})^{2}+100(x_{2}-x_{1}^{2})^{2}+100(x_{4}-x_{3}^{2})^{2})/10000$ | Rosenbrock-4D | $[-5,5]^{4}$
E5 | $\sin(x_{1}+x_{2}x_{3})+\exp{(1.2x_{4})}$ | — |  | $x_{1}\in[-10,10],\,x_{2}\in[-5,5],$
---
$\,x_{3}\in[-5,5],\,x_{4}\in[-3,3]$
E6 | $\tanh(x_{1}/2)+|x_{2}|\cos(x_{3}^{2}/5)$ | — | $[-10,10]^{3}$
E7 | $(1-x_{2}^{2})/(\sin(2\pi\,x_{1})+1.5)$ | [33] | $[-5,5]^{2}$
E8 | $x_{1}^{4}/(x_{1}^{4}+1)+x_{2}^{4}/(x_{2}^{4}+1)$ | [29] | $[-5,5]^{2}$
E9 | $\log(2x_{2}+1)-\log(4x_{1}^{2}+1)$ | [29] | $[0,5]^{2}$
E10 | $\sin(x_{1}\,e^{x_{2}})$ | [3] | $x_{1}\in[-2,2],x_{2}\in[-4,4]$
E11 | $x_{1}\,\log(x_{2}^{4})$ | [3] | $[-5,5]^{2}$
E12 | $1+x_{1}\,\sin(1/x_{2})$ | [3] | $[-10,10]^{2}$
E13 | $\sqrt{x_{1}}\,\log(x_{2}^{2})$ | [3] | $x_{1}\in[0,20],x_{2}\in[-5,5]$

We assessed the performance of the skeletons generated by our Multi-Set
Transformer using ten synthetic SR problems generated by equations inspired by
previous works and equations proposed in this work, as reported in Table 4.1.
Note that previous works used narrow domain ranges for all variables (e.g.,
$[-1,1]$) while we used extended ranges (e.g., $[-5,5]$ and $[-10,10]$) to
increase the difficulty of the problems. We adapted the benchmark equations
proposed by Bertschinger et al. [3] (i.e., E10–E13) to a multivariate setting.

In all cases, the training datasets consisted of 10,000 points where each
variable was sampled using a uniform distribution. For the estimated response
functions, $\hat{f}$, we trained feed-forward NNs with varying depths: three
hidden layers for problem E2; five hidden layers for problems E1, E4, E5, and
E7; and four hidden layers for the other cases. Each layer consisted of 500
nodes with ReLU activation. We used 90% of the samples for training and 10%
for validation.

We compared the skeletons produced by the Multi-Set Transformer (MST) to the
ones extracted from the expressions generated by three other methods: two GP-
based methods (PySR [7] and TaylorGP [13]) and two neural SR methods (NeSymReS
[4] and E2E [18]). For NeSymReS and E2E, we used the pre-trained models
available online. NeSymReS, E2E, and the Multi-Set Transformer were trained
using vocabularies with the same unary and binary operators: $+,\times,/$,
abs, acos, asin, atan, cos, cosh, exp, log, pow2, pow3, pow4, pow5, sin, sinh,
sqrt, tan, tanh. Thus, PySR and TaylorGP were executed using the same set of
operators. Our experimentation with the GP-based methods involved a maximum of
10,000 iterations, though convergence was consistently achieved in fewer
iterations across all cases. Population sizes of 100, 200, 500, and 1000 were
tested, with no discernible advantage observed beyond a size of 500.

The compared methods produce multivariate expressions, from which the skeleton
variable corresponding to variable $x_{v}$ is obtained using the skeleton
function $\kappa(\cdot,x_{v})$. NeSymReS could not be executed on E4 and E5
because its model was trained using expressions limited to three variables.
Table 4.2 shows the target and estimated skeletons corresponding to each
variable for problem E2. The skeletons obtained for the other problems are
presented in Section C of the Supplementary Materials1. We also evaluated
skeleton performance using the method described in Section 3.3. We set the
size of the test sets to $N_{\text{test}}=3000$. Using $N_{\text{test}}>3000$
did not vary the obtained results. Table 4.3 reports the rounded mean and the
standard deviation of the error metrics obtained after 30 repetitions of the
proposed evaluation. The bold entries indicate the method that achieved the
lowest mean error $r$ and that its difference w.r.t. the values obtained by
the other methods is statistically significant according to Tukey’s honestly
significant difference test performed at the 0.05 significance level.

Table 4.2: Comparison of skeleton prediction results for problem E2

Method | $x_{1}$ | $x_{2}$ | $x_{3}$
---|---|---|---
PySR | $c_{1}+\left|c_{2}+|c_{3}+x_{1}|\right|$ | $c_{1}$ | $c_{1}+c_{2}\,x_{3}$
TaylorGP | $c_{1}+c_{2}\,x_{1}$ | $c_{1}+c_{2}\,x_{2}$ | $c_{1}+c_{2}\,x_{3}$
NeSymReS | $c_{1}+c_{2}\,x_{1}$ | $c_{1}+\exp(\exp(c_{2}\,x_{2}))$ | $c_{1}+c_{2}\,x_{3}$
E2E | $c_{1}+c_{2}\,x_{1}+c_{3}(c_{4}+c_{5}\,x_{1})^{2}$ | $c_{1}+c_{2}\,(c_{3}+c_{4}\,x_{2})$ | $c_{1}+c_{2}\,x_{3}+c_{3}(c_{4}+c_{5}\cos(c_{6}+c_{7}\,x_{3}))$
MST | $c_{1}+c_{2}(c_{3}+c_{4}\,x_{1})^{2}$ | $c_{1}\sqrt{c_{2}\,x_{2}+c_{3}}+c_{4}$ | $c_{1}+c_{2}\sin(c_{3}\,x_{3}+c_{4})$
Target $\mathbf{e}(x)$ | $c_{1}+(c_{2}+c_{3}\,x_{1})^{2}$ | $c_{1}\sqrt{x_{2}+c_{2}}+c_{3}$ | $c_{1}+c_{2}\sin(c_{3}\,x_{3})$

Table 4.3: Skeleton evaluation performance comparison

Eq. | Var. | PySR | TaylorGP | NeSymReS | E2E | MST
---|---|---|---|---|---|---
E1 | $x_{1}$ | $1.4\pm 0.8$ | $1.4\pm 0.8$ | $0.9\pm 0.7$ | $\mathbf{0.2\pm 0.4}$ | $\mathbf{0.01\pm 0.02}$
$x_{2}$ | $1.5\pm 0.9$ | $1.5\pm 0.9$ | $1.3\pm 0.8$ | $1.5\pm 0.9$ | $\mathbf{0\pm 0}$
E2 | $x_{1}$ | $303.5\pm 167.3$ | $310.0\pm 170.1$ | $310.0\pm 170.1$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
$x_{2}$ | $5.4\pm 5.0$ | $4.2\pm 5.3$ | $4.6\pm 5.0$ | $4.2\pm 5.3$ | $\mathbf{0.02\pm 0.03}$
$x_{3}$ | $1.7\pm 1.0$ | $1.7\pm 1.0$ | $1.7\pm 1.0$ | $\mathbf{0.01\pm 0.02}$ | $\mathbf{0\pm 0}$
E3 | $x_{1}$ | $2\\!\times\\!10^{12}\pm 5\\!\times\\!10^{12}$ | $939.4\pm 1419.9$ | $1.9\pm 1.2$ | $\mathbf{0.8\pm 1.8}$ | $\mathbf{0.8\pm 1.8}$
$x_{2}$ | $1.3\pm 1.0$ | $1.3\pm 1.0$ | $0.8\pm 0.8$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
E4 | $x_{1}$ | $4576.2\pm 2695.7$ | $4581.5\pm 2697.4$ | — | $\mathbf{2.3\pm 3.6}$ | $\mathbf{1.1\pm 0.7}$
$x_{2}$ | $79.6\pm 41.3$ | $80.2\pm 40.8$ | — | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
$x_{3}$ | $3995.5\pm 2815.6$ | $4304.6\pm 2843.7$ | — | $\mathbf{2.0\pm 4.0}$ | $\mathbf{1.0\pm 0.9}$
$x_{4}$ | $74.5\pm 48.0$ | $75.5\pm 47.0$ | — | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
E5 | $x_{1}$ | $0.6\pm 0.05$ | $0.6\pm 0.05$ | — | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
$x_{2}$ | $1.5\pm 1.0$ | $1.5\pm 1.0$ | — | $1.5\pm 1.0$ | $\mathbf{0\pm 0}$
$x_{3}$ | $0.6\pm 0.05$ | $0.6\pm 0.05$ | — | $0.6\pm 0.05$ | $\mathbf{0\pm 0}$
$x_{4}$ | $2.7\pm 1.2$ | $487.4\pm 461.9$ | — | $\mathbf{0.6\pm 0.8}$ | $\mathbf{0.6\pm 0.8}$
E6 | $x_{1}$ | $0.8\pm 0.08$ | $0.8\pm 0.08$ | $0.3\pm 0.01$ | $0.04\pm 0$ | $\mathbf{0\pm 0}$
$x_{2}$ | $16.8\pm 12.2$ | $16.8\pm 12.2$ | $13.8\pm 11.1$ | $1.3\pm 0.9$ | $\mathbf{0\pm 0}$
$x_{3}$ | $2.9\pm 1.3$ | $1.9\pm 0.7$ | $1.6\pm 0.9$ | $1.6\pm 0.9$ | $\mathbf{0\pm 0}$
E7 | $x_{1}$ | $29.5\pm 1.0$ | $1.8\pm 2.0$ | $1.8\pm 2.0$ | $1.1\pm 1.3$ | $\mathbf{0\pm 0}$
$x_{2}$ | $63.6\pm 43.8$ | $1.6\pm 1.0$ | $42.4\pm 24.7$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
E8 | $x_{1}$ | $0.02\pm 0.01$ | $0.04\pm 0.01$ | $0.8\pm 1.2$ | $0.02\pm 0.02$ | $\mathbf{0\pm 0}$
$x_{2}$ | $0.02\pm 0.01$ | $0.04\pm 0.01$ | $0.9\pm 1.3$ | $0.02\pm 0.01$ | $\mathbf{0\pm 0}$
E9 | $x_{1}$ | $271.3\pm 446.8$ | $239.8\pm 428.2$ | $375.1\pm 485.5$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
$x_{2}$ | $\mathbf{0\pm 0}$ | $0.2\pm 0.09$ | $2.7\pm 1.7$ | $0.05\pm 0.01$ | $\mathbf{0\pm 0}$
E10 | $x_{1}$ | $\mathbf{0\pm 0}$ | $0.6\pm 0.2$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
$x_{2}$ | $\mathbf{0\pm 0}$ | $0.4\pm 0.06$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
E11 | $x_{1}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
$x_{2}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
E12 | $x_{1}$ | $21.8\pm 13.1$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
$x_{2}$ | $\mathbf{0\pm 0}$ | $2.4\pm 1.5$ | $2.5\pm 1.6$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$
E13 | $x_{1}$ | $\mathbf{0\pm 0}$ | $0.8\pm 0.5$ | $0.7\pm 0.6$ | $0.7\pm 0.8$ | $\mathbf{0\pm 0}$
$x_{2}$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$ | $3.8\pm 3.5$ | $\mathbf{0\pm 0}$ | $\mathbf{0\pm 0}$

###  4.2 Analysis of Nitrogen Fertilizer Response Curves

In the previous section, we performed experiments based on synthetic datasets.
Recall that the objective of this work is to demonstrate that our method can
explain the relationship between each system’s variable and the system’s
response. To do so, datasets with known underlying functions are needed.
However, in this section, we show the applicability of our approach in a real-
world application.

In precision agriculture, crop yield prediction is a critical task with a
significant societal impact [14]. Let $y(i,j)$ represent the observed yield at
a field site with coordinates $(i,j)$. Furthermore, let $\mathbf{x}(i,j)$
represent a set of multiple covariate factors, such as the nitrogen (N)
fertilizer rate applied, and other topographic and meteorological features. It
describes the state of the field at position $(i,j)$, and potentially its
neighboring areas. The underlying yield function of the field is denoted as
$f(\cdot)$ and $y(i,j)=f(\mathbf{x}(i,j))$. Thus, the early crop yield
prediction problem consists of generating predicted yield values
$\hat{y}(i,j)$ using a prediction model $\hat{f}(\cdot)$, such that
$\hat{y}(i,j)=\hat{f}(\mathbf{x}(i,j))$ and $\hat{y}(i,j)\approx y(i,j)$.

In practice, $f$ is a complex multivariate system with unknown functional
form. Nevertheless, tasks like N-rate optimization, which allows for profit
maximization and environmental impact maximization [14], do not require
estimating the full functional form of $f(\mathbf{x}(i,j))$. Instead, N-rate
optimization only analyzes the functional relationship between the N-rate
variable and the predicted yield values. This relationship is typically
represented using N fertilizer-yield response (N-response) curves, which
depict the estimated crop yield values for a specific field site in response
to various admissible fertilizer rates.

The experiments presented in this section aim to estimate the functional form
of N-response curves for a winter wheat dryland field called “Field A”,
discussed in a previous work [24]. Traditionally, N-response curves are fitted
employing a single parametric function for the entire field [5, 17]. However,
previous work [24] demonstrates that the N-response functional form varies
within the field, which can be clustered into regions known as management
zones (MZs) based on the fertilizer responsivity of the field sites (i.e., the
shape of the approximated N-response curves). For details about the dataset,
the generation of approximated N-response curves, and the clustering method,
please see [24].

Fig. 4.1 depicts Field A clustered into four MZs. The figure also shows fifty
example response curves selected randomly from each MZ. Due to their shape
similarity, we assume that all sites within an MZ share the same functional
form. Consider that an N-response curve is a set of input–response pairs
(N-rate $Nr$ vs. relative yield $ry$). Thus, the problem of estimating the
functional form of the N-response curves of an MZ can be posed as an MSSP. To
do so, we select $N_{S}=10$ random field sites within an MZ and feed their
corresponding N-response curves into the pre-trained MST model to obtain a
univariate skeleton.

Figure 4.1: MZs for Field A and corresponding N-response curves [24]. Table
4.4: Comparison of skeleton prediction results for problem E2

Method | MZ | Functional Form | $\bar{r}$ | Method | MZ | Functional Form | $\bar{r}$
---|---|---|---|---|---|---|---
Quadratic-plateau | 1 | $c_{1}+c_{2}\,(\min(Nr,c_{3})+c_{3})^{2}$ | 0.0680 | Quadratic-plateau | 3 | $c_{1}+c_{2}\,(\min(Nr,c_{3})+c_{3})^{2}$ | 0.0879
Exponential | $c_{1}(1-\exp(c_{2}+c_{3}\,Nr))+c_{4}$ | 0.2303 | Exponential | $c_{1}(1-\exp(c_{2}+c_{3}\,Nr))+c_{4}$ | 0.1965
MST | $c_{1}+c_{2}(Nr-\log(c_{3}+c_{4}\,Nr))$ | 0.0489 | MST | $c_{1}+c_{2}\,Nr+c_{3}\,cos(c_{4}+c_{5}\,Nr)$ | 0.0479
Quadratic-plateau | 2 | $c_{1}+c_{2}\,(\min(Nr,c_{3})+c_{3})^{2}$ | 0.0725 | Quadratic-plateau | 4 | $c_{1}+c_{2}\,(\min(Nr,c_{3})+c_{3})^{2}$ | 0.0615
Exponential | $c_{1}(1-\exp(c_{2}+c_{3}\,Nr))+c_{4}$ | 0.1825 | Exponential | $c_{1}(1-\exp(c_{2}+c_{3}\,Nr))+c_{4}$ | 0.2249
MST | $c_{1}+c_{2}\,\texttt{tanh}(c_{3}+c_{4}\,Nr)$ | 0.0495 | MST | $c_{1}+c_{2}\,Nr+c_{3}\,cos(c_{4}+c_{5}\,Nr)$ | 0.0237

Table 4.4 shows the skeletons derived by the MST for each MZ. The process was
repeated multiple times for each MZ, with different random sites selected as
inputs for the MSSP problem in each iteration, resulting in the generation of
equivalent skeletons. Furthermore, we evaluated the suitability of the
obtained skeletons for each MZ and compared them to two traditional N-response
models: quadratic-plateau [5] and exponential [17]. For each method and field
site, we fit the skeleton’s coefficient values to minimize the distance to the
corresponding N-response curve, following a similar procedure as the one
described in Sec. 3.3. We report the average error $\bar{r}$ obtained
considering all sites within each MZ. Results indicate that different MZs may
be modeled using different functional forms, unlike traditional approaches
that assume one for the entire field.

##  5 Discussion

Our method involves training a feedforward NN, which is then used to solve an
MSSP problem for each system’s variable using a Multi-Set Transformer. This
process produces univariate skeletons that describe the functional
relationship between each variable and the system’s response. After evaluation
of our univariate skeleton prediction method across the tested problems, we
observed that it generated skeletons that matched or were equivalent to the
target skeleton for all variables across all problems. For instance, for
problem E11, the target skeleton for variable $x_{2}$ is given by
$\mathbf{e}(x_{2})=c_{1}\,\log(x_{2}^{4})$, and our method produces the
skeleton
$\mathbf{\hat{e}}_{MST}(x_{2})=c_{1}^{\prime}+c_{2}^{\prime}\,\log(c_{3}^{\prime}\,x_{2}^{2})$.
Notice that these skeletons are equivalent if $c_{1}^{\prime}=0$,
$c_{2}^{\prime}=2\,c_{1}$, and $c_{3}^{\prime}=1$. In addition, from the
skeleton performance evaluation shown in Table 4.3, we verified that our
method consistently attained lower or comparable error metrics compared to
other SR methods. These results strongly support our hypothesis that our
method would generate univariate skeletons that are more similar to those
corresponding to the underlying equations in comparison to other SR methods.

Note that E2E produced the correct skeleton for at least one of the variables
in most of the cases. Recall that E2E generated expressions that minimize the
prediction error; thus, it did not prioritize identifying the correct
functional form of the variables that do not contribute substantially to the
overall error. In addition, in some cases, E2E generated skeletons that were
equivalent to the target skeletons but larger than those produced by MST. For
example, in problem E4, E2E generated the skeleton
$\mathbf{\hat{e}}_{E2E}(x_{1})=c_{1}+c_{2}|c_{3}+c_{4}\,x_{1}+c_{5}\,x_{1}^{2}+c_{6}\,x_{1}^{3}+c_{7}\,x_{1}^{4}|$
for variable $x_{1}$, which is equivalent to the one produced by MST,
$\mathbf{\hat{e}}_{MST}(x_{1})=c_{1}+c_{2}\,x_{1}+c_{3}\,x_{1}^{2}+c_{4}\,x_{1}^{3}+c_{5}\,x_{1}^{4}$.
Another advantage over the other neural SR methods is that MST requires 24.2
million parameters while E2E requires 93.5 million. NeSymReS requires 26.4
million parameters (i.e., 2.2 million more than MST), and Table 4.3 shows that
it failed to identify the correct functional form in most cases and is limited
to problems with up to three variables.

It is worth pointing out that in problems E5, E6, E8, E9, and E13, some
compared methods achieved low error metrics but are not comparable to the ones
achieved by our method. For example, in problem E6, E2E generated the skeleton
$\mathbf{\hat{e}}_{E2E}(x_{1})=c_{1}+c_{2}\,\text{atan}(c_{3}+c_{4}\,x_{1})$
for variable $x_{1}$, which does not coincide with the functional form of the
underlying skeleton $\mathbf{e}(x_{1})=c_{1}\\!+\\!\text{tanh}(c_{2}\,x_{1})$.
However, E2E achieved low error metrics because, during the coefficient
fitting process, the GA found appropriate values for the constant that
multiplies the argument of the atan function, stretching or compressing the
curve, making it resemble the shape of tanh and minimizing the error. Hence,
the skeleton generated by E2E produced low error metrics and is considered to
be similar to the target skeleton. Conversely, the functional form of the
skeleton generated by MST coincided with that of the target skeleton exactly
(i.e., $\mathbf{\hat{e}}_{MST}(x_{1})=c_{1}+c_{2}\,\text{tanh}(c_{3}\,x_{1})$)
and thus produced significantly lower error metrics.

We claim that our method can be regarded as a post-hoc explainability method
that generates univariate skeletons as explanations of the function
approximated by a black-box regression model. From an explainability
perspective, the generation of more accurate univariate symbolic skeletons is
crucial. By producing skeletons that align more closely with the underlying
data, we offer transparent insights into how each variable influences the
function’s response.

Finally, we applied our method to a precision agriculture problem. Given a
field that has been split into MZs, the problem consists of determining the
functional form that describes the N-response curves corresponding to all
sites within an MZ. This is modeled as an MSSP problem and the skeleton
generated for each MZ is reported in Table 4.4. Our results show that the
generated skeletons yield lower fitting errors and, thus are more suitable
when modeling the field’s N-response curves. In future work, the function
fitted for each site will serve as a surrogate yield model that may aid in
optimization tasks.

One potential limitation of our approach, as well as any neural SR method,
lies in its ability to generate skeletons whose complexity is bounded
inherently by the expressions produced during the pre-training phase of the
Multi-Set Transformer. For instance, we would not be able to identify the
skeleton $c_{1}+c_{2}\,x_{2}^{2}/\sin(c_{3}\,e^{c_{4}\,x_{2}})$ as it requires
eight operators, while our training set was limited to expressions with up to
seven operators. However, it is feasible to overcome this limitation through
transfer learning, so that the MST model can be trained on more complex tasks,
potentially enabling the recognition of such complex skeletons. Another
limitation is that we do not provide multi-variate expressions as other
methods do. This is because the objective of this work lies in discerning the
functional relationship between each variable and the system’s response
accurately. Future work will focus on using the generated univariate skeletons
as building blocks to produce multivariate expressions that approximate the
system’s entire underlying function. For instance, mathematical expressions
explaining the full behavior of real-world systems will be generated, and
their performance will be evaluated by comparing observed response values with
those predicted by the generated expressions.

##  6 Conclusions

Symbolic regression aims to find symbolic expressions that represent the
relationships within the observed data. As such, SR represents a promising
avenue for building explainable models. By seeking to uncover mathematical
equations that represent the relationships between input variables and their
response, the resulting equations offer transparency and clear insights into
model behavior.

Given a multivariate regression problem, the objective of the work reported in
this paper is to produce symbolic skeletons that describe the relationship
between each variable and the system’s response. To do so, we introduced a
problem called Multi-Set Skeleton Prediction that aims to generate a symbolic
skeleton expression that characterizes the functional form of multiple sets of
input–response pairs. To solve this problem, we proposed a novel Multi-Set
Transformer model, which was pre-trained on a large dataset of synthetic
symbolic expressions. Then, the multivariate SR problem is tackled as a
sequence of MSSP problems generated using a black-box regression model. The
generated skeletons are regarded as explanations of the black-box model’s
function. Experimental results showed that our method consistently produced
more accurate univariate skeletons in comparison to two GP-based SR methods
and two neural SR methods.

Future work will focus on merging the generated univariate skeletons into a
multivariate symbolic expression that approximates the underlying function of
the system. This task requires us to answer questions such as how to ensure
that the skeletons generated for all variables are compatible and can be
merged, or how to find the optimal order in which to merge the skeletons.

## Acknowledgments

The authors wish to thank the members of the Numerical Intelligent Systems
Laboratory (NISL) for their comments throughout the development of this work.
This research was supported by the Data Intensive Farm Management project
(USDA-NIFA-AFRI 2016-68004-24769 and USDA-NRCS NR213A7500013G021).
Computational efforts were performed on the Tempest HPC System, operated and
supported by University Information Technology Research Cyberinfrastructure at
Montana State University.

## References

  * [1] Amir Haeri, M., Ebadzadeh, M.M., Folino, G.: Statistical genetic programming for symbolic regression. Applied Soft Computing 60, 447–469 (2017)
  * [2] Arnaldo, I., Krawiec, K., O’Reilly, U.M.: Multiple regression genetic programming. In: 2014 Annual Conference on Genetic and Evolutionary Computation. p. 879â€“886 (2014)
  * [3] Bertschinger, A., Davis, Q.T., Bagrow, J., Bongard, J.: The metric is the message: Benchmarking challenges forÂ neural symbolic regression. In: Machine Learning and Knowledge Discovery in Databases. pp. 161–177 (2023)
  * [4] Biggio, L., Bendinelli, T., Neitz, A., Lucchi, A., Parascandolo, G.: Neural symbolic regression that scales. In: 38th International Conference on Machine Learning. vol. 139, pp. 936–945 (2021)
  * [5] Bullock, D.G., Bullock, D.S.: Quadratic and quadratic-plus-plateau models for predicting optimal nitrogen rate of corn: A comparison. Agronomy Journal 86(1), 191–195 (1994)
  * [6] Chu, X., Zhao, H., Xu, E., Qi, H., Chen, M., Shao, H.: Neural symbolic regression using control variables. ArXiv abs/2306.04718 (2023)
  * [7] Cranmer, M.: Interpretable machine learning for science with PySR and SymbolicRegression.jl. ArXiv abs/2305.01582 (2023)
  * [8] Dignum, S., Poli, R.: Generalisation of the limiting distribution of program sizes in tree-based genetic programming and analysis of its effects on bloat. In: 9th Annual Conference on Genetic and Evolutionary Computation. p. 1588â€“1595 (2007)
  * [9] Ebner, M.: On the search space of genetic programming and its relation to nature’s search space. In: 1999 Congress on Evolutionary Computation. vol. 2, pp. 1357–1361 Vol. 2 (1999)
  * [10] Edwards, H., Storkey, A.: Towards a neural statistician. In: International Conference on Learning Representations. Toulon, France (2017)
  * [11] Filho, R., Lacerda, A., Pappa, G.: Explaining symbolic regression predictions. In: 2020 IEEE Congress on Evolutionary Computation. pp. 1–8 (2020)
  * [12] Fletcher, R.: Practical Methods of Optimization. John Wiley & Sons, New York, NY, USA, second edn. (1987)
  * [13] He, B., Lu, Q., Yang, Q., Luo, J., Wang, Z.: Taylor genetic programming for symbolic regression. In: Genetic and Evolutionary Computation Conference. p. 946â€“954 (2022)
  * [14] Hegedus, P., Maxwell, B., Sheppard, J., Loewen, S., Duff, H., Morales, G., Peerlinck, A.: Towards a low-cost comprehensive process for on-farm precision experimentation and analysis. Agriculture 13(3) (2023)
  * [15] Holland, J.H.: Genetic algorithms. Scientific American (Jul 1992)
  * [16] Jin, Y., Fu, W., Kang, J., Guo, J., Guo, J.: Bayesian symbolic regression. ArXiv abs/1910.08892 (2020)
  * [17] Kakimoto, S., Mieno, T., Tanaka, T., Bullock, D.: Causal forest approach for site-specific input management via on-farm precision experimentation. Computers and Electronics in Agriculture 199, 107164 (2022)
  * [18] Kamienny, P.A., d'Ascoli, S., Lample, G., Charton, F.: End-to-end symbolic regression with transformers. In: Advances in Neural Information Processing Systems. vol. 35, pp. 10269–10281 (2022)
  * [19] La Cava, W., Orzechowski, P., Burlacu, B., de França, F.O., Virgolin, M., Jin, Y., Kommenda, M., Moore, J.: Contemporary symbolic regression methods and their relative performance. In: Advances in Neural Information Processing Systems (2021)
  * [20] Lample, G., Charton, F.: Deep learning for symbolic mathematics. In: International Conference on Learning Representations (2020)
  * [21] Langdon, W.B., Soule, T., Poli, R., Foster, J.A.: The Evolution of Size and Shape. In: Advances in Genetic Programming. The MIT Press (07 1999)
  * [22] Lee, J., Lee, Y., Kim, J., Kosiorek, A., Choi, S., Teh, Y.W.: Set Transformer: A framework for attention-based permutation-invariant neural networks. In: 36th International Conference on Machine Learning. vol. 97, pp. 3744–3753 (2019)
  * [23] Martius, G., Lampert, C.H.: Extrapolation and learning equations. ArXiv abs/1610.02995 (2016)
  * [24] Morales, G., Sheppard, J.W.: Counterfactual analysis of neural networks used to create fertilizer management zones. In: 2024 International Joint Conference on Neural Networks. Yokohama, Japan (Jun 2024)
  * [25] Orzechowski, P., La Cava, W., Moore, J.: Where are we now? A large benchmark study of recent symbolic regression methods. In: Genetic and Evolutionary Computation Conference. p. 1183â€“1190 (2018)
  * [26] Petersen, B.K., Landajuela, M., Mundhenk, T.N., Santiago, C.P., Kim, S.K., Kim, J.T.: Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. In: International Conference on Learning Representations (2021)
  * [27] Sahoo, S., Lampert, C., Martius, G.: Learning equations for extrapolation and control. In: 35th International Conference on Machine Learning. pp. 4442–4450 (10–15 Jul 2018)
  * [28] Schmidt, M., Lipson, H.: Distilling free-form natural laws from experimental data. Science 324(5923), 81–85 (2009)
  * [29] Trujillo, L., MuÃ±oz, L., GalvÃ¡n-LÃ³pez, E., Silva, S.: neat genetic programming: Controlling bloat naturally. Information Sciences 333, 21–43 (2016)
  * [30] Udrescu, S.M., Tegmark, M.: AI Feynman: A physics-inspired method for symbolic regression. Science Advances 6(16), eaay2631 (2020)
  * [31] Valipour, M., You, B., Panju, M., Ghodsi, A.: SymbolicGPT: A generative transformer model for symbolic regression. In: Advances in Neural Information Processing Systems (2022)
  * [32] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A., Kaiser, L., Polosukhin, I.: Attention is all you need. In: Advances in Neural Information Processing Systems. vol. 30 (2017)
  * [33] Werner, M., Junginger, A., Hennig, P., Martius, G.: Informed equation learning. ArXiv abs/2105.06331 (2021)
  * [34] Zaheer, M., Kottur, S., Ravanbakhsh, S., Poczos, B., Salakhutdinov, R.R., Smola, A.J.: Deep sets. In: Advances in Neural Information Processing Systems. vol. 30. Curran Associates, Inc. (2017)

##  Appendix 0.A Set Transformer vs. Multi-Set Transformer

Here, we provide background information about the Set Transformer, the
foundational basis for the Multi-Set Transformer. Furthermore, we enumerate
the differences between the architectures and application domains of both
models.

The Transformer [32] was designed to handle sequences with fixed orderings,
where each element’s position is crucial to understanding the data’s meaning.
In contrast, sets are collections of elements without any inherent order, and
their permutations do not alter the underlying semantics. This inherent
permutation invariance poses a significant obstacle for traditional sequence
models when processing sets [34, 10]. In order to address the limitations of
conventional approaches, Lee et al. [22] presented an attention-based neural
network module called Set Transformer that is based on the transformer model.
This method introduces modifications to the transformer architecture, enabling
it to handle sets without assuming a fixed ordering of elements.

A model designed for set-input problems must meet two essential criteria to
effectively handle sets: First, it should be capable of processing input sets
of varying sizes; second, it should exhibit permutation invariance. The latter
means the output of the function represented by the model remains the same
regardless of the order in which the elements of the input set are presented.
More formally, Zaheer et al. [34] described this type of function as
permutation equivariant:

###### Definition 0.A.1

Consider an input set $\mathbf{x}=\\{x_{1},\dots,x_{n}\\}$, where
$x_{i}\in\mathscr{X}$ (i.e., the input domain is the power set
$\mathcal{X}=2^{\mathscr{X}}$). A function
$u:2^{\mathscr{X}}\rightarrow\mathcal{Y}$ acting on sets must be permutation
invariant to the order of objects in the set, i.e. for any permutation $\pi$:

$u(\\{x_{1},\dots,x_{n}\\})=u(\\{x_{\pi(1)},\dots,x_{\pi(n)}\\}),$

where $\pi\in S_{n}$ and $S_{n}$ represents the set of all permutations of
indices $\\{1,\dots,n\\}$.

The Set Transformer consists of two main parts: an encoder $\phi$ and a
decoder $\psi$. The process starts by encoding the set elements in an order-
agnostic manner. As such, the encoding for each element should be the same
regardless of its position in the set. Then, $\psi$ aggregates the encoded
features and produces the desired output. Let us consider an input set
$\mathbf{S}=\\{\mathbf{s}_{1},\dots,\mathbf{s}_{n}\\}$, where each element is
$d_{in}$-dimensional ($\mathbf{S}\in\mathbb{R}^{n\times d_{in}}$). Therefore,
the output $T$ produced by the Set Transformer, whose computed function is
denoted as $g$, can be expressed as:

$T=g(\mathbf{S})=\psi\left(\phi\left(\left\\{\mathbf{s}_{1},\dots,\mathbf{s}_{n}\right\\}\right)\right).$
(0.A.1)

Typically, the order-agnostic property is achieved by making $\phi$ to act on
each element of a set independently (i.e.,
$\psi\left(\left\\{\phi(\mathbf{s}_{1}),\dots,\phi(\mathbf{s}_{n})\right\\}\right)$)
[34]. However, the Set Transformer uses self-attention mechanisms to encode
the entire input set simultaneously to recognize interactions among the set
instances.

For example, Lee et al. [22] used self-attention blocks (SABs), which, given
an input set $A$, perform self-attention between the elements of the set and
produce a set of equal size. This operation computes pairwise interactions
among the elements of $A$; therefore, a stack of multiple SAB operations would
encode higher-order interactions. One drawback of the use of SABs is that the
attention mechanism is performed between two identical sets with $n$ elements,
which leads to a quadratic time complexity $\mathcal{O}(n^{2})$. In order to
alleviate this issue, Zaheer et al. [34] introduced the induced set attention
block (ISAB). This could be interpreted as if the input set $A$ was projected
into a lower dimensional space and then reconstructed to produce outputs with
the desired dimensionality. Since the attention mechanism involved in the ISAB
operation is computed between a set of size $m$ and a set of size $n$, its
associated time complexity is $\mathcal{O}(mn)$. Furthermore, Lee et al. [22]
proved that the SAB and ISAB blocks are universal approximators of permutation
invariant functions.

It is worth noting that the Set Transformer was designed to produce fixed-size
outputs. This setting can be used for applications such as population
statistic estimation (e.g., retrieving a unique value that presents the median
of an input set), unique character counting (i.e., obtaining the number of
unique characters in a set of images), or $k$-amortized clustering where the
objective is to produce $k$ pairs of output parameters of a mixture of
Gaussians (i.e., mean and covariance).

Figure 0.A.1: An example of an MSSP problem using the Multi-set Transformer.

The most evident limitation of the Set Transformer, in the context of the
problem discussed in this paper, resides in its encoder structure, which is
specifically designed to process a single input set, as shown in Eq. 0.A.1.
This is important because the input of the MSSP problem is defined as a
collection of $N_{S}$ input sets ($N_{S}>1$). Hence, the encoder of our Multi-
set Transformer is designed to process multiple input sets simultaneously. The
Multi-Set Transformer architecture is depicted in Fig. 0.A.1. All input sets
shown in this example were generated using the equation
$y=\frac{1}{x\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{r}{x})^{2}}$ and a different
$r$ value was used for each set $\mathbf{D}^{(s)}$. Its architecture differs
from the one used in the SSP method proposed by Biggio et al. [4], whose
encoder only processes single input sets.

Note that the Set Transformer’s output is $k$-dimensional; that is, its size
is fixed depending on the problem. Conversely, the objective of the MSSP
problem is to generate a symbolic skeleton string whose length is not known a
priori and depends on each input collection. Therefore, the decoder of the
Multi-set Transformer is designed as a conditional-generative structure as it
generates output sequences (i.e., the skeleton string) based on the encoded
context.

##  Appendix 0.B Data Generation

In this section, we explain how to generate the synthetic equations and
corresponding data used to train our Multi-set Transformer. In addition, we
explain the methodology we used to avoid generating functions that lead to
undefined values during the training process.

###  0.B.1 Univariate Symbolic Skeleton Dataset Generation

Table 0.B.1: Vocabulary used to pre-train the Multi-set Transformer. Token |  Meaning |  Token |  Meaning
---|---|---|---
SOS |  Start of sentence |  sin |  Sine
EOS |  End of sentence |  sinh |  Hyperbolic sine
c |  Constant placeholder |  sqrt |  Square root
x |  Variable |  tan |  Tangent
abs |  Absolute value |  tanh |  Hyperbolic tangent
acos |  Arc cosine |  -3 |  Integer number
add |  Sum |  -2 |  Integer number
asin |  Arc sine |  -1 |  Integer number
atan |  Arc tangent |  0 |  Integer number
cos |  Cosine |  1 |  Integer number
cosh |  Hyperbolic cosine |  2 |  Integer number
div |  Division |  3 |  Integer number
exp |  Exponential |  4 |  Integer number
log |  Logarithmic |  5 |  Integer number
mul |  Multiplication |  E |  Euler’s number
pow |  Power |  — |  —

Table 0.B.1 provides the vocabulary used in this work. Our method differs from
the generation method proposed by Biggio et al. [4], which, in turn, was based
on that by Lample et al. [20]. Their method consists of constructing an
expression tree with a specified number of total operators (i.e., unary and
binary operators) by selecting operators and inserting them into random
available places of the tree structure iteratively. The operators are chosen
based on a pre-defined set of probability weights. In practice, we noticed
that this approach leads to the generation of several expressions that contain
binary operators exclusively and, thus, can be simplified into simple
expressions such as $c_{1}\,x$ and $c_{1}\,x^{2}+c_{2}$.

As an alternative, we use a generation method that builds the expression tree
recursively in a preorder fashion (i.e., the root node is created before
generating its left and right subtrees recursively). At each step, the
algorithm decides whether to add a binary operator, a unary operator, or a
leaf node. Instead of assigning specific probability weights to individual
operators, we assign weights to the types of nodes to be generated, namely
unary operators, binary operators, or leaf nodes, with weights of 2, 1, and 1,
respectively. In addition, rather than specifying the exact number of
operators to include, we set a maximum limit on the total number of operators
allowed in the expression. This adds flexibility in cases where including
additional operators may incur violations of syntaxis or pre-defined
conditions.

By generating the tree in a preorder manner, we can enforce certain conditions
and constraints effectively. In particular, we can impose restrictions such as
setting a maximum limit on the nesting depth of unary operators within each
other. In addition, we can forbid certain combinations of operators. For
example, we avoid embedding the operator log within the operator exp, or vice
versa, since such composition could lead to direct simplification (e.g.,
$\texttt{log}\left(\texttt{exp}(x)\right)=x$). We can also avoid combinations
of operators that would generate extremely large values (e.g.,
$\texttt{exp}\left(\texttt{exp}(x)\right)$ and
$\texttt{sinh}\left(\texttt{sinh}(x)\right)$). Table 0.B.2 shows the forbidden
operators we considered for some specific parent operators.

Table 0.B.2: Forbidden combinations of operators Parent Node | Forbidden Operators
---|---
abs |  sqrt, pow2, pow4
exp, tan, ln |  exp, sinh, cosh, tanh, tan, ln, pow3, pow4, pow5
sinh, cosh, tanh |  exp, sinh, cosh, tanh, tan, ln, pow2, pow3, pow4, pow5
pow2, pow3, pow4, pow5 |  pow2, pow3, pow4, pow5, exp, sinh, cosh, tanh
sin, cos, tan |  sin, cos, tan
asin, acos, atan |  asin, acos, atan
Figure 0.B.1: Example of a randomly generated expression.

Each generated tree is traversed to derive a mathematical expression in prefix
notation. The expression is then transformed from prefix to infix notation,
which is simplified using the SymPy library333https://www.sympy.org/. The
simplified expression is transformed again into an expression tree. Then, each
non-numerical node of the tree is multiplied by a unique placeholder and then
added to another unique placeholder. Some exceptions apply; for instance, the
arguments of the exp, sinh, cosh, and tanh operators are not affected by
additive placeholders. This is because adding a constant to their arguments
would lead to direct simplification (e.g.,
$c_{1}e^{c_{2}\,x+c_{3}}=(c_{1}e^{c_{2}})e^{c_{3}\,x}$). Fig. 0.B.1
illustrates the random expression generation process. Note that previous
approaches also proposed to include unique placeholders in the generated
expressions; however, they included only one placeholder at a time and in some
of the nodes. In contrast, our approach allows for the generation of more
general expressions.

###  0.B.2 Avoiding Invalid Operations

In this section, we present a detailed description of the function
$\texttt{avoidNaNs}(\mathbf{x},f)$, initially introduced in Algorithm 1 of the
main paper. This function serves the purpose of modifying certain coefficients
within the function $f$ or generating supplementary support values within the
vector $\mathbf{x}$ to avoid numerical inconsistencies. We classify the
operators that may generate undefined values as follows:

  * •

Single-bounded operators: These operators have bounded numerical arguments due
to the mathematical constraints imposed on their domains. We consider the
following single-bounded operators:

    * –

Logarithm: The log operator is bounded on its left side as it cannot receive
an input lower than or equal to 0:

$\text{Domain}\left(\texttt{log}(x)\right)=\\{x\in\mathbb{R}\;|\;x>0\\}.$

    * –

Square root: The sqrt operator is also bounded on its left side to avoid
generating complex numbers:

$\text{Domain}\left(\texttt{sqrt}(x)\right)=\\{x\in\mathbb{R}\;|\;x\geq 0\\}.$

    * –

Exponential: We decided to bound the exp operator on its right side to avoid
generating extremely large values. This also applies to the sinh, cosh, and
tanh operators:

$\text{Domain}\left(\texttt{exp}(x)\right)=\\{x\in\mathbb{R}\;|\;x\leq 7\\}.$

  * •

Double-bounded operators: Unlike the single-bounded operators, the numerical
arguments of these operators are bounded on their left and right sides. We
consider the following double-bounded operators:

    * –

Arcsine: The asin operator takes a value between -1 and 1 as its input and
returns the angle whose sine is equal to that value:

$\text{Domain}\left(\texttt{asin}(x)\right)=\\{x\in\mathbb{R}\;|\;-1\leq x\leq
1\\}.$

    * –

Arccosine: Like arcsine, the acos also takes a value between -1 and 1 as its
input operator:

$\text{Domain}\left(\texttt{acos}(x)\right)=\\{x\in\mathbb{R}\;|\;-1\leq x\leq
1\\}.$

  * •

Operators with singularities: Operators like tangent or division can exhibit
singularities at specific input values. For instance, the tangent function
becomes undefined when its input equals an odd multiple of $\pi/2$, resulting
in an asymptotic behavior where the function approaches infinity.

The method proposed to address the NaN (“not a number”) values arising from
the aforementioned functions, which we refer to as “special functions”, is
shown in Algorithm 0.B.1. Here, $f.\texttt{args}$ returns a vector containing
the numerical arguments of function $f$. If $f$ is considered as an expression
tree, $f.\texttt{func}$ returns the name of the operator located at the top of
the tree. The algorithm analyzes each argument of $f$ separately (Line 3). A
given argument of $f$, denoted as $arg$, can be considered as a sub-expression
tree. Hence, the function $\texttt{containSpecialF}(arg)$ (Line 4) traverses
the sub-expression tree and returns a true value if any special function is
found in it.

If the current sub-expression tree contains a special function, we check if
the top operator of the subtree is special using function
$\texttt{isSpecialF}(arg.\texttt{func})$ (Line 5). If not, we move down to a
deeper level of the subtree using recursion (Line 6). Otherwise, the arguments
of the sub-expression $arg$ may need to be modified to avoid undefined values.
Before doing so, in Line 8, we verify if there is another special function
contained inside the current subtree, in which case we explore a deeper level
of the subtree using recursion. Note that we apply the
$\texttt{containSpecialF}(arg.\texttt{args}[0])$ function given that $arg$ is
guaranteed to represent a special function; i.e., it is a unary function with
a single argument.

The function $arg$ is evaluated on the input values $\mathbf{x}$ and the
obtained values are stored in the vector $vals$ (Line 11). The function
$\texttt{containNaNs}(vals)$ returns a true value in the event that one or
more undefined values are found within the vector $vals$. If an undefined
value is found, a modification of the inner argument, $innArg$, or the input
vector $\mathbf{x}$ is needed. To do so, we first evaluate the function
represented by $innArg$ on the input values $\mathbf{x}$ and store the
outcomes in the variable $innVals$ (Line 13). In the case that the function
$arg$ is single-bounded, the domain of function $innArg$ is modified
accordingly using $\texttt{modifySBounded}(arg,vals,innArg,innVals)$ (Line 15)
by adding a horizontal offset. Likewise, if $arg$ is double-bounded, the
domain of $innArg$ is modified accordingly using
$\texttt{modifyDBounded}(arg,innArg,innVals)$ (Line 17) by re-scaling
$innArg$. In the case that $arg$ represents a function with singularities, a
new input vector $\mathbf{x}$ with resampled values is obtained using the
function $\texttt{avoidSingularities}(\mathbf{x},arg,Xsing)$ (Line 19), where
$Xsing$ is a variable that will store all positions at which the function
produces undefined values. Fig. 0.B.2 depicts an example of how data are
generated using the underlying function
$f(x)=\frac{-3.12x}{\texttt{sin}(1.45x)}-2.2$. The figure at the bottom shows
in detail that the avoidSingularities function generates valid intervals of
values for $\mathbf{x}$ by avoiding getting too close to the singular points
(red dotted line). Furthermore, the inner argument of function $arg$ is then
replaced with the modified function $innArg$ (Line 21). Finally, the arguments
of the original function $f$ are replaced with the modified arguments stored
in the list $newArgs$ (Line 22).

1:function avoidNaNs($\mathbf{x},f,Xsing=[\,]$)

2: $args,newArgs\leftarrow f.\texttt{args},[]$

3: for each $arg\in args$ do

4: if $\texttt{containSpecialF}(arg)$ then

5: if $!\texttt{isSpecialF}(arg.\texttt{func})$ then $\triangleright$ If the
top operator is special

6: $\mathbf{x},arg,Xsing\leftarrow\texttt{avoidNaNs}(\mathbf{x},arg,Xsing)$

7: else

8: if $\texttt{containSpecialF}(arg.\texttt{args}[0])$ then $\triangleright$
If there’s a special function inside

9: $\mathbf{x},arg,Xsing\leftarrow\texttt{avoidNaNs}(\mathbf{x},arg,Xsing)$

10: $innArg\leftarrow arg.\texttt{args}[0]$

11: $vals\leftarrow arg(\mathbf{x})$

12: if $\texttt{containNaNs}(vals)$ then

13: $innVals\leftarrow innArg(\mathbf{x})$

14: if $\texttt{isSingleBounded}(arg)$ then

15: $innArg\leftarrow\texttt{modifySBounded}(arg,vals,innArg,innVals)$

16: else if $\texttt{isDoubleBounded}(arg)$ then

17: $innArg\leftarrow\texttt{modifyDBounded}(arg,innArg,innVals)$

18: else $\triangleright$ Operations with singularities

19:
$\mathbf{x},sing\leftarrow\texttt{avoidSingularities}(\texttt{length}(\mathbf{x}),arg,Xsing)$

20: $arg.\texttt{args}[0]\leftarrow innArg$ $\triangleright$ Update function

21: $newArgs.\texttt{append}(arg)$

22: $f.\texttt{args}\leftarrow newArgs$

23: return $\mathbf{x},arg,sing$

Algorithm 0.B.1 Avoiding Invalid Operations Figure 0.B.2: Generation data for
$f(x)=\frac{-3.12x}{\texttt{sin}(1.45x)}-2.2$. (Left) Generated data on the
entire domain. (Right) Detailed view of how singularities are avoided.

##  Appendix 0.C Comparison of Predicted Symbolic Skeletons

In this section, we report the symbolic skeletons predicted by two GP-based SR
methods, PySR [7] and TaylorGP [13], and three neural SR methods, NeSymRes
[4], E2E [18], and our Multi-Set Transformer.

Tables 0.C.1 and 0.C.2 report the skeletons obtained by all compared methods.
It is important to note that these experiments are designed to validate the
functional form obtained between each system’s variable and the system’s
response. That is, instead of evaluating the expressions generated by the
compared methods directly, we first extract their skeletons with respect to
each system’s variable. For example, NeSymRes produces the estimated function
$\hat{f}(x_{1},x_{2})=0.592\,x_{1}\,x_{2}+\cos(0.014\,(-0.995\,x_{1}+x_{2}-0.08)^{2})$
for problem E1. From this, the skeletons with respect to variables $x_{1}$ and
$x_{2}$ are $\hat{e}(x_{1})=c_{1}\,x_{1}+\cos(c_{2}(c_{3}+c_{4}\,x1)^{2})$ and
$\hat{e}(x_{2})=c_{1}\,x_{2}+\cos(c_{2}(c_{3}+x_{2})^{2})$, respectively.

Table 0.C.1: Comparison of skeleton prediction results (E1–E9)

Prob. |  Var. Met. | PySR | TaylorGP | NeSymReS | E2E | MST |  Target $\mathbf{e}(x)$
---|---|---|---|---|---|---|---
E1 | $x_{1}$ | $c_{1}\,x_{1}$ | $c_{1}\,x_{1}$ | $c_{1}\,x_{1}+\cos(c_{2}(c_{3}+c_{4}\,x_{1})^{2})$ | $c_{1}+c_{2}\,x_{1}+c_{3}\sin((c_{4}+c_{5}\,x_{1})^{2})$ | $c_{1}+c_{2}\,x_{1}+c_{3}\sin(c_{4}+c_{5}\,x1)$ | $c_{1}\,x_{1}+c_{2}\,\sin(c_{3}(c_{4}+x_{1}))$
$x_{2}$ | $c_{1}\,x_{2}$ | $c_{1}\,x_{2}$ | $c_{1}\,x_{2}+\cos(c_{2}(c_{3}+x_{2})^{2})$ | $c_{1}+c_{2}(c_{3}+c_{4}\,x_{2})$ | $c_{1}+c_{2}\,x_{2}+c_{3}\sin(c_{4}+c_{5}\,x_{2})$ | $c_{1}\,x_{2}+c_{2}\,\sin(c_{3}(c_{4}+x_{2}))$
E2 | $x_{1}$ | $c_{1}+\left|c_{2}+|c_{3}+x_{1}|\right|$ | $c_{1}+c_{2}\,x_{1}$ | $c_{1}+c_{2}\,x_{1}$ | $c_{1}+c_{2}\,x_{1}+c_{3}(c_{4}+c_{5}\,x_{1})^{2}$ | $c_{1}+c_{2}(c_{3}+c_{4}\,x_{1})^{2}$ | $c_{1}+(c_{2}+c_{3}\,x_{1})^{2}$
$x_{2}$ | $c_{1}$ | $c_{1}$ | $c_{1}+e^{e^{c_{2}\,x_{2}}}$ | $c_{1}+c_{2}\,(c_{3}+c_{4}\,x_{2})$ | $c_{1}\sqrt{c_{2}\,x_{2}+c_{3}}+c_{4}$ | $c_{1}\sqrt{x_{2}+c_{2}}+c_{3}$
$x_{3}$ | $c_{1}+c_{2}\,x_{3}$ | $c_{1}+c_{2}\,x_{3}$ | $c_{1}+c_{2}\,x_{3}$ | $c_{1}+c_{2}\,x_{3}+c_{3}(c_{4}+c_{5}\cos(c_{6}+c_{7}\,x_{3}))$ | $c_{1}+c_{2}\sin(c_{3}\,x_{3}+c_{4})$ | $c_{1}+c_{2}\sin(c_{3}\,x_{3})$
E3 | $x_{1}$ | $c_{1}\,e^{x_{1}}\,|\sinh(c_{2}\,x_{1})|$ | $c_{1}\,e^{c_{2}\,x_{1}}$ | $c_{1}\,e^{c_{2}\,x_{1}}$ | $c_{1}+c_{2}e^{c_{3}\,x_{1}}$ | $c_{1}+c_{2}e^{c_{3}\,x_{1}}$ | $c_{1}+c_{2}\,e^{c_{3}\,x_{1}}$
$x_{2}$ | $c_{1}$ | $c_{1}$ | $c_{1}\,\cos(c_{2}\,x_{2})$ | $c_{1}+c_{2}\,\cos(c_{3}+c_{4}\,x_{2})$ | $c_{1}+c_{2}\,\cos(c_{3}+c_{4}\,x_{2})$ | $c_{1}+c_{2}\,\cos(c_{3}\,x_{2})$
E4 | $x_{1}$ | $c_{1}+c_{2}\,x_{1}^{2}$ | $c_{1}\,x_{1}^{2}$ | — | $c_{1}+c_{2}|c_{3}+c_{4}\,x_{1}+c_{5}\,x_{1}^{2}+c_{6}\,x_{1}^{3}+c_{7}\,x_{1}^{4}|$ | $c_{1}(c_{2}\,x1+c_{3})^{4}+c_{4}$ | $c_{1}+c_{2}\,x_{1}+c_{3}\,x_{1}^{2}+c_{4}\,x_{1}^{4}$
$x_{2}$ | $c_{1}+c_{2}\,x_{2}$ | $c_{1}$ | — | $c_{1}+c_{2}|c_{3}+c_{4}\,x_{2}+c_{5}\,x_{2}^{2}|$ | $c_{1}(c_{2}\,x2+c_{3})^{2}+c_{4}$ | $c_{1}+c_{1}\,x_{2}+c_{3}\,x_{2}^{2}$
$x_{3}$ | $c_{1}+c_{2}\,x_{3}^{2}$ | $c_{1}$ | — | $c_{1}+c_{2}|c_{3}+c_{4}\,x_{3}+c_{5}\,x_{3}^{2}+c_{6}\,x_{3}^{3}+c_{7}\,x_{3}^{4}|$ | $c_{1}(c_{2}\,x3+c_{3})^{4}+c_{4}$ | $c_{1}+c_{2}\,x_{3}+c_{3}\,x_{4}^{2}+c_{4}\,x_{3}^{4}$
$x_{4}$ | $c_{1}+c_{2}\,x_{4}$ | $c_{1}$ | — | $c_{1}+c_{2}|c_{3}+c_{4}\,x_{4}+c_{5}\,x_{4}^{2}|$ | $c_{1}(c_{2}\,x4+c_{3})^{2}+c_{4}$ | $c_{1}+c_{2}\,x_{4}+c_{3}\,x_{4}^{2}$
E5 | $x_{1}$ | $c_{1}$ | $c_{1}$ | — | $c_{1}+c_{2}\cos(c_{3}+c_{4}\,x_{1})$ | $c_{1}\sin(c_{2}\,x1+c_{3})+c_{4}$ | $c_{1}+\sin(c_{2}+c_{3}\,x_{1})$
$x_{2}$ | $c_{1}$ | $c_{1}$ | — | $c_{1}+c_{2}\,x_{2}$ | $c_{1}\sin(c_{2}\,x2+c_{3})+c_{4}$ | $c_{1}+\sin(c_{2}+c_{3}\,x_{2})$
$x_{3}$ | $c_{1}$ | $c_{1}$ | — | $c_{1}+c_{2}\,x_{3}+c_{3}\,x_{3}^{2}+c_{4}\,x_{3}^{3}$ | $c_{1}\sin(c_{2}\,x3+c_{3})+c_{4}$ | $c_{1}+\sin(c_{2}+c_{3}\,x_{3})$
$x_{4}$ | $e^{c_{1}\,x_{4}}$ | $c_{1}\,e^{x_{4}}\,e^{\sin(c_{2}\,x_{4})}$ | — | $c_{1}+c_{2}\,e^{c_{3}\,x_{4}}$ | $c_{1}e^{c_{2}\,x4}+c_{3}$ | $c_{1}+e^{c_{2}\,x_{4}}$
E6 | $x_{1}$ | $c_{1}$ | $c_{1}$ | $c_{1}+c_{2}\,x_{1}$ | $c_{1}+c_{2}\,\text{atan}(c_{3}+c_{4}\,x_{1})$ | $c_{1}+c_{2}\,\text{tanh}(c_{3}\,x_{1})$ | $c_{1}+\text{tanh}(c_{2}\,x_{1})$
$x_{2}$ | $c_{1}$ | $c_{1}$ | $c_{1}+x_{2}\,\sin(c_{2}/(c_{3}+x_{2}))$ | $c_{1}+c_{2}\,x_{2}+c_{3}\,x_{2}\,\cos(c_{4}+c_{5}\,x_{2})$ | $c_{1}+c_{2}\,|x_{2}|$ | $c_{1}+c_{2}\,|x_{2}|$
$x_{3}$ | $\text{tanh}(\exp(x_{3}))$ | $c_{1}\,\sin(c_{2}\,x_{3}^{2})/(c_{3}\,\sqrt{x_{3}}+\sin(\sqrt{x_{3}}))$ | $c_{1}$ | $c_{1}$ | $c_{1}\,\cos(c_{2}\,(c_{3}+x_{3})^{2})+c_{4}$ | $c_{1}+c_{2}\,(\cos(c_{3}\,x_{3}^{2}))$
E7 | $x_{1}$ | $c_{1}+x_{1}^{2}$ | $c_{1}$ | $c_{1}+c_{2}\,x_{1}$ | $c_{1}+c_{2}\,\sin(c_{3}+c_{4}\,x_{1})^{2}+c_{5}\,\sin(c_{6}+c_{7}\,x_{1})$ | $c_{1}+c_{2}/(c_{3}+\sin(c_{4}\,x_{1}+c_{5}))$ | $c_{1}/(c_{2}+\sin(c_{3}\,x_{1}))$
$x_{2}$ | $c_{1}/\text{sinh}(\text{sinh}(\text{tanh}(e^{\text{sinh}(\sin(c_{2}\,x_{2}))})))$ | $c_{1}\,x_{2}^{2}+\sqrt{|x_{2}|}$ | $(c_{1}+x_{2}^{2})/(c_{2}+\cos(c_{3}\,(c_{4}+c_{5}\,x_{2})^{3}))$ | $c_{1}\,(c_{2}+c_{3}\,x_{2})^{2}$ | $c_{1}+c_{2}\,x_{2}^{2}$ | $c_{1}+c_{2}\,x_{2}^{2}$
E8 | $x_{1}$ | $c_{1}+\text{tanh}(c_{2}+\text{cosh}(x_{1}))$ | $c_{1}$ | $e^{c_{1}\,\cos(c_{2}/x_{1})}$ | $c_{1}+c_{2}\,e{c_{3}\,|c_{4}+c_{5}\,x_{1}|}$ | $c_{1}+c_{2}/(c_{3}\,x_{1}^{4}+c_{4})$ | $c_{1}+c_{2}\,x_{1}^{4}/(c_{3}+c_{4}\,x_{1}^{4})$
$x_{2}$ | $c_{1}+\cos(\tan(\text{tanh}(c_{2}/x_{2}^{2})))$ | $c_{1}$ | $c_{1}+\cos(1/x_{2})$ | $c_{1}+c_{2}\,e^{c_{3}\,|c_{4}+x_{2}|}$ | $c_{1}+c_{2}/(c_{3}\,x_{1}^{4}+c_{4})$ | $c_{1}+c_{2}\,x_{2}^{4}/(c_{3}+c_{4}\,x_{2}^{4})$
E9 | $x_{1}$ | $\log(c_{1}/(c_{2}+c_{3}\,x_{1}^{2}))$ | $c_{1}+\log(c_{2}/|x_{1}|)$ | $\log(c/|x_{1}|)$ | $c_{1}+c_{2}\,\log(c_{3}+c_{4}\,x_{1}+c_{5}\,x_{1}^{2})$ | $c_{1}+\log(c_{2}\,x_{1}^{2}+c_{3})$ | $c_{1}+c_{2}\,\log(c_{3}+c_{4}\,x_{1}^{2})$
$x_{2}$ | $\log(c_{1}+c_{2}\,x_{2})$ | $c_{1}\exp(c_{2}\,x_{2})$ | $\log(c\,|x_{2}|)$ | $c_{1}+c_{2}\,\log(c_{3}+c_{4}/(c_{5}+c_{6}\,x_{2}+c_{7}\,x_{2}^{2}))$ | $c_{1}+\log(c_{2}+c_{3}\,x_{2})$ | $c_{1}+\log(c_{2}+c_{3}\,x_{2})$

Table 0.C.2: Comparison of skeleton prediction results (E10–E13)

Prob. | Var. | PySR | TaylorGP | NeSymReS | E2E | MST |  Target $\mathbf{e}(x)$
---|---|---|---|---|---|---|---
E10 | $x_{1}$ | $\sin(c_{1}\,x_{1})$ | $c_{1}\,x_{1}$ | $\sin(c_{1}\,x_{1})$ | $c_{1}+c_{2}\,\sin(c_{3}+c_{4}\,x_{1})$ | $c_{1}+c_{2}\,\sin(c_{3}\,x_{1})$ | $\sin(c_{1}\,x_{1})$
$x_{2}$ | $\sin(c_{1}\,e^{x_{2}})$ | $c_{1}\,e^{c_{2}\,\sqrt{|x_{2}|}}$ | $\sin(c_{1}\,e^{x_{2}})$ | $c_{1}+c_{2}\,\sin(c_{3}+c_{4}\,e^{c_{5}\,x_{2}})$ | $c_{1}+\sin(c_{2}\,e^{c_{3}\,x_{2}})$ | $\sin(c_{1}\,e^{x_{2}})$
E11 | $x_{1}$ | $c_{1}\,x_{1}$ | $c_{1}\,x_{1}$ | $c_{1}\,x_{1}$ | $c_{1}\,x_{1}$ | $c_{1}+c_{2}\,x_{1}$ | $c_{1}\,x_{1}$
$x_{2}$ | $c_{1}\,\log(x_{2}^{4})$ | $c_{1}\,\log(|x_{2}|)$ | $c_{1}\,\log(x_{2}^{4})$ | $c_{1}+c_{2}\,\log(c_{3}+c_{4}\,|c_{5}+c_{6}\,x_{2}|)$ | $c_{1}+c_{2}\,\log(c_{3}\,x_{2}^{2})$ | $c_{1}\,\log(x_{2}^{4})$
E12 | $x_{1}$ | $c_{1}+\sin(x_{1})\,|x_{1}|$ | $c_{1}+c_{2}\,x_{1}$ | $c_{1}+c_{2}\,x_{1}$ | $c_{1}+c_{2}\,x_{1}$ | $c_{1}+c_{2}\,x_{1}$ | $c_{1}+c_{2}\,x_{1}$
$x_{2}$ | $c_{1}+c_{2}\,\sin(c_{3}/x_{2})$ | $\begin{array}[]{ccc}c_{1}\,\sin(c_{2}/x_{2})+\\\ \sqrt{x_{2}}\,\sin(c_{3}/x_{2})\end{array}$ | $\begin{array}[]{ccc}c_{1}\,\sin(1/x_{2})+\\\ x_{2}\,\sin(1/x_{2})\end{array}$ | $c_{1}\,\sin(c_{2}/(c_{3}+c_{4}\,x_{2}))$ | $c_{1}+c_{2}\,\sin(c_{3}/x_{2})$ | $\begin{array}[]{ccc}c_{1}+\\\ c_{2}\,\sin(1/x_{2})\end{array}$
E13 | $x_{1}$ | $c_{1}\,\sqrt{x_{1}}$ | $c_{1}+c_{2}\,\sqrt{e^{\sqrt{x_{1}}}}$ | $c_{1}+c_{2}\,x_{1}$ | $c_{1}+c_{2}\,\log(c_{3}+c_{4}\,x_{1})$ | $c_{1}+c_{2}\,\sqrt{c_{3}+x_{1}}$ | $c_{1}\,\sqrt{x_{1}}$
$x_{2}$ | $c_{1}\,\log(x_{2}^{2})$ | $\begin{array}[]{ccc}c_{1}+\log(|x_{2}|)+\\\ c_{2}\,\log(|x_{2}|)\end{array}$ | $c_{1}\,\log(x_{2}^{2})$ | $c_{1}+c_{2}/(c_{3}+c_{4}\,|c_{5}+c_{6}\,x_{2}|)$ | $c_{1}+c_{2}\,\log(c_{3}\,(c_{4}+x_{2})^{2})$ | $c_{1}\,\log(x_{2}^{2})$


$D$ | 0.524 | 0.529 | 0.506 | 0.242
$D\cup S_{T}$ | 0.488 | 0.458 | 0.472 | 0.141
$P$ | 0.097 | 0.033 | 0.037 | 0.089
$P\cup D$ | 0.510 | 0.543 | 0.541 | 0.140
$P\cup D\cup S_{T}$ | 0.467 | 0.363 | 0.494 | 0.120
Table 5: FMNIST $\shortrightarrow$ MNIST

###  B.10 Comparison with Traditional Augmentation Methods

We conduct the traditional data augmentation methods, including random
rotation and Gaussian noise, on the datasets and compare the classification
accuracy with our method. The random rotation is from 0 to 90 degrees. The
Gaussian noise is with kernel size (5,9) and standard deviation 0.1 and 5. All
the comparison results are in Table 6–9.

Through this experiment, we show that using random rotation to augment the
dataset does not improve accuracy much and sometimes even decreases the
accuracy in all of the eight transfer learning settings. We note that ($D$,
$P$, $P\cup D$) are settings without augmentation while ($D\cup S_{T}$, $P\cup
D\cup S_{T}$) are settings with augmentation. Comparing with the results of
our proposed gradient flow approach in Figure 3, augmentation (by random
rotation) has much lower performances than ours.

Accuracy | $F\rightarrow M$ | $K\rightarrow M$ | $M\rightarrow F$ | $M\rightarrow K$
---|---|---|---|---
$D$ | 0.535 | 0.512 | 0.483 | 0.240
$D\cup S_{T}$ | 0.542 | 0.498 | 0.489 | 0.164
$P$ | 0.097 | 0.033 | 0.037 | 0.089
$P\cup D$ | 0.510 | 0.543 | 0.541 | 0.140
$P\cup D\cup S_{T}$ | 0.551 | 0.581 | 0.454 | 0.124
Table 6: 1-shot learning augmented with random rotation Accuracy | $F\rightarrow M$ | $K\rightarrow M$ | $M\rightarrow F$ | $M\rightarrow K$
---|---|---|---|---
$D$ | 0.788 | 0.837 | 0.519 | 0.234
$D\cup S_{T}$ | 0.569 | 0.545 | 0.490 | 0.174
$P$ | 0.097 | 0.033 | 0.037 | 0.089
$P\cup D$ | 0.772 | 0.763 | 0.581 | 0.178
$P\cup D\cup S_{T}$ | 0.738 | 0.768 | 0.593 | 0.199
Table 7: 5-shot learning augmented with random rotation Accuracy | $F\rightarrow M$ | $K\rightarrow M$ | $M\rightarrow F$ | $M\rightarrow K$
---|---|---|---|---
$D$ | 0.524 | 0.529 | 0.506 | 0.242
$D\cup S_{T}$ | 0.512 | 0.494 | 0.452 | 0.167
$P$ | 0.097 | 0.033 | 0.037 | 0.089
$P\cup D$ | 0.510 | 0.543 | 0.541 | 0.140
$P\cup D\cup S_{T}$ | 0.557 | 0.590 | 0.454 | 0.127
Table 8: 1-shot learning augmented with random Gaussian noise Accuracy | $F\rightarrow M$ | $K\rightarrow M$ | $M\rightarrow F$ | $M\rightarrow K$
---|---|---|---|---
$D$ | 0.827 | 0.824 | 0.549 | 0.265
$D\cup S_{T}$ | 0.792 | 0.803 | 0.734 | 0.320
$P$ | 0.097 | 0.032 | 0.037 | 0.089
$P\cup D$ | 0.779 | 0.764 | 0.581 | 0.178
$P\cup D\cup S_{T}$ | 0.832 | 0.797 | 0.729 | 0.258
Table 9: 5-shot learning augmented with random Gaussian noise

We show that adding Gaussian noise improves the accuracy in some settings, for
example, M$\rightarrow$F and M$\rightarrow$K in 5-shot, but decreases accuracy
in many settings. Our gradient flow method improves the accuracy in all the
settings when we compare the accuracy without $S_{T}$ and with $S_{T}$.

Finally, we again note that those simple augmentations, for example (adding
Gaussian noise or using random rotation which are designed onto the target
dataset) can be used as a complement to the gradient flow method, which
transfers data from the source dataset to the target dataset.

###  B.11 Random Gaussian Noise to MNIST

In this section, we conduct an experiment that transfers from random noise
vectors to MNIST images. The noise vectors are independently drawn from the
normal distribution. Thus, the mean of the noise vectors is a zero vector and
the covariance matrix is an identity matrix. From this experiment, we want to
show that our method can transfer from any distribution to the target
distribution.

Figure 9: Our method successfully transfers from random noise vectors to MNIST
images.

###  B.12 Evaluation Metric by SSIM and Optimal Transport

We construct a new measure that connects a human perception distance and
optimal transport to measure the distance between the distribution of
transferred samples and the distribution of target samples.

The target data consists of $M$ samples $(\bar{x}_{j},\bar{y}_{j})_{j=1}^{M}$,
and the generated data consists of $N$ samples
$(x_{i}^{T},y_{i}^{T})_{i=1}^{N}$. We use the structural similarity index
measure (SSIM) wang2004image to measure image similarity between the generated
images and the target images. SSIM measures the similarity between two
grayscale images. SSIM is between 0 and 1, and $\mathrm{SSIM}=1$ means that
the two images are identical. Thus, we define a SSIM-based cost to construct a
matrix $C\in\mathbb{R}^{N\times M}$:
$C_{ij}=1-\mathrm{SSIM}(x_{i}^{T},\bar{x}_{j}),\text{for
}i=1,\ldots,N;\leavevmode\nobreak\ j=1,\ldots,M.$

We use the matrix $C$ as a ground transportation cost and solve the Earth
Mover Distance problem rubner1998metric; rubner2000earth to obtain the
Wasserstein distance between target image and generated image distributions,
which are $\frac{1}{N}\sum_{i=1}^{N}\delta_{(x_{i}^{T},y_{i}^{T})}$and
$\frac{1}{M}\sum_{j=1}^{M}\delta_{(\bar{x}_{j},\bar{y}_{j})}.$ We compute the
distance in 1-shot FashionMNIST to MNIST experiment and Figure 10 plots the
distance of each class over 10 independent runs. The plot shows our algorithm
successfully decreases the Wasserstein distance between the distribution of
generated images and the target distribution during each step.

Figure 10: This plots the OT_SSIM distance between the distribution of
transferred images and of the target images over 140 steps. The distance
decreases as we flow the samples from the source domain to the target domain
for every class.