1 Introduction↩︎

SR is a machine learning paradigm that distills symbolic expressions from data. Other machine learning methods assume a fixed functional form and seek to identify its parameters. In contrast, SR searches for both the functional form and its associated parameters, aiming for more concise representations. Symbolic expressions are human interpretable and have long been used to describe and study phenomena in various disciplines. As such, SR is a uniquely powerful tool for scientific discovery, system identification, and explainable artificial intelligence.

Koza [1] introduced SR itself and proposed an evolutionary-inspired method, GP. Putting GP in simple terms, a population of randomly generated expressions is randomly combined and changed. Inferior expressions are discarded in favor of better ones, as determined by the optimization objectives that balance fit quality and expression complexity. Beyond [1], GP methods have been continually improved in the last \(30\) plus years and are still among the best approaches to SR [2], [3]. There are numerous other approaches for SR including, but not limited to, exhaustive search [4], [5], reinforcement-learning-based methods [6]–[8], and amortized approaches [9], [10]. Although some are well-poised to do so, these have not yet eclipsed GP-based methods in all respects.

SR has successfully been applied to scientific discovery in many domains, , astrophysics [11]–[13], material science [14]–[16], and thermodynamics [17]–[21]. For some problems, it is beneficial to describe similar systems using the same functional form and allow it to differ only in the parameter values. For example, [22] developed a single, shared expression linking the residual entropy to the viscosity of \(124\) fluids. Some of the parameters are fluid-agnostic and are shared among all fluids, while other parameters are fluid-specific (non-shared) and change to accommodate the different behavior. Intuitively, finding a functional form capable of describing multiple similar phenomena increases the chances of it being meaningful and of generalizing beyond the provided category values. The chances of the functional form that is valid for \(124\) fluids being valid for a \(125^{\mathrm{th}}\) fluid are high. Notably, this viscosity correlation was developed without the use of SR. This approach essentially introduces one categorical variable, with the various fluids serving as its different category values.¹

Thus, categorical variables allow different continuous datasets to be combined, either when describing the same phenomenon measured under different conditions or when describing different but similar phenomena. In some cases, categorical variables may hide unwanted complexity or compress not-yet-understood interrelations of missing independent variables capable of explaining the distinct categorical values, as is also discussed in [23]. For example, the ideal gas law is \(p = \rho \frac{R}{M}T\), where \(R\) is the ideal gas constant, \(M\) is the substance-specific molar mass, and \(p\), \(\rho\), \(T\) are the variables pressure, density, and temperature. If the substance-specific molar masses were unknown, we could develop a shared expression defining each substance as a category value. For the viscosity correlation mentioned, it remains unclear which latent continuous variables account for the variance between different fluids.

Sharing a functional form across category values regularizes the functional form [24]. This may prevent “symbolic overfitting”, improve generalization performance, and increase robustness to noise. Some of the prior approaches use a fully distinct set of parameters for each category value, while others incorporate category-agnostic or “shared” parameters. Sharing some of the parameters between the category values reduces the number of parameters and potentiates all mentioned benefits, as discussed in [25]. Furthermore, this approach also leads to transfer learning to some degree, where knowledge is reused across phenomena [26], [27].

In recent years, several works have explored incorporating a single categorical variable into SR in such a (or a similar) manner with some differences in their approaches [23]–[25], [27], [28]; see also the implementations  [29] and  [30]. However, no prior work considered two or more categorical variables, or proposed a method capable of taking the effects of multiple categorical variables and their interactions into account separately.

Problems with two or more categorical variables can be approached using category-value combinations as a single, condensed categorical variable. However, this ignores a crucial part of the problem and fails to fully exploit the potential for parameter sharing and reduction. To the best of our knowledge, all prior works but one have explored approaches capable of incorporating only one categorical variable. Only Kronberger  [23] could, in principle, handle multiple categorical variables separately. However, their approach cannot account for the interactions. In the present work, we demonstrate a novel approach that uses three types of parameter sharing, “fully-shared” (global), “partially-shared” (category-value-specific), and “non-shared” (category-value-combination-specific).

The proposed approach is introduced in the following section. In 3, we discuss existing methods and compare them to our approach. Finally, the approach and its benefits are demonstrated in 4, before a conclusion is drawn in 5.

2 Method↩︎

In this section, for simplicity’s sake, the proposed approach is shown for a two-category example. However, the approach is applicable to arbitrarily many categorical variables.

2.1 Concept↩︎

To explain the concept, we introduce an example problem. The problem has two categorical variables \(U\) (uppercase) and \(L\) (lowercase). Category \(U\) has four category values (A, B, C, D), while Category \(L\) has three values (a, b, c). Thus, in the datasets, there are four times three possible category-value combinations.

The core of our approach is the use of three distinct levels of parameter sharing: fully-shared, partially-shared, and non-shared, which depend on neither, either, or both category values, respectively. The highest level of sharing is represented by fully-shared parameters, which are the same across all category-value combinations. Non-shared parameters are unique to each category-value combination. The intermediate sharing level is represented by partially-shared parameters. These parameters are shared across one of the two categories, but depend on the values of the other category. This approach allows us to distinguish category-agnostic effects, effects of either category, and also category-value-combination interactions. 1 illustrates the sharing concepts for the example at hand.

As an intuitive example, if \(U\) were patient groups and \(L\) were medications, a parameter partially shared across the patient types may account for the effects of a particular medicine on all patient types. A parameter partially shared across medications may account for a particular patient-type’s reaction to medication in general. Finally, non-shared parameters would correct for binary interactions of particular medications and particular patient types, , allergies. The approach is reminiscent of mixed-effect models or group-contribution methods.

2.2 Evaluation and Parameter Identification↩︎

Our implementation shares some similarity with Kronberger  [23], with some extra considerations regarding the new types of parameter sharing. Shared parameters \(p_{\mathrm{shared}}\) consist of a single real number minimized across the entire dataset. Partially-shared parameters depending on category \(U\), denoted as \(p_{\mathrm{partial},U}\), contain one real value for each unique level \(u \in U\); the same logic applies to other categories. Non-shared parameters \(p_{\mathrm{nonshared}}\) are optimized for each category-value combination \((u,l)\) individually, containing as many real numbers as there are unique category-value combinations \((u,l)\).

For each category pair \((u,l)\) and each data point \(i\) belonging to that pair, the observed value \(y_{u,l,i}\) is predicted with \[\label{eq:predict} \hat{y}_{u,l,i} = f\paren[big](){ x_{u,l,i};\, p_{\mathrm{shared}},\, p_{\mathrm{partial},U,u},\, p_{\mathrm{partial},L,l},\, p_{\mathrm{nonshared},(u,l)} } .\tag{1}\] The least-squares fitting objective is formulated as a nested loop over the categorical variables and the data points \(\cI_{u,l}\) within each category-value combination: \[\label{eq:minimize} \text{Minimize} \quad \sum_{u \in U} \; \sum_{l \in L} \; \sum_{i \in \cI_{u,l}} \paren[big](){ y_{u,l,i} - \hat{y}_{u,l,i} }^{2} \quad \text{\wrt\;} P ,\tag{2}\] where \(P = \paren[big](){p_{\mathrm{shared}},\, p_{\mathrm{partial},U},\, p_{\mathrm{partial},L},\, p_{\mathrm{nonshared}}}\).

To illustrate, we introduce the shared expression \[\begin{gather} \label{eq:example1} \hat{y}_{u,l,i} = p_{\mathrm{shared}} \cdot v_1 + p_{\mathrm{partial},U,u} \cdot v_1^2 \@ifclassloaded{acmart}{\\}{} + p_{\mathrm{partial},L,l} \cdot v_1^3 + p_{\mathrm{nonshared},(u,l)} \cdot v_1^4, \end{gather}\tag{3}\] where \(v_1\) is a continuous, independent variable, \(p_{\mathrm{shared}} = 100\) is a fully-shared parameter, \(p_{\mathrm{partial},U} = \set{\texttt{A} \mapsto 10,\, \texttt{B} \mapsto 20,\, \texttt{C} \mapsto 30,\, \texttt{D} \mapsto 40}\) is a partially-shared parameter depending on the categorical variable \(U\), \(p_{\mathrm{partial},L} = \set{\texttt{a} \mapsto 1,\, \texttt{b} \mapsto 2,\, \texttt{c} \mapsto 3}\) is a partially-shared parameter depending on \(L\), and \(p_{\mathrm{nonshared}} = \set{\texttt{Aa} \mapsto 0.01,\, \texttt{Ab} \mapsto 0.02,\, \texttt{Ac} \mapsto 0.03,\, \texttt{Ba} \mapsto 0.04,\, \texttt{Bb} \mapsto 0.05,\, \texttt{Bc} \mapsto 0.06,\, \texttt{Ca} \mapsto 0.07,\, \texttt{Cb} \mapsto 0.08,\, \texttt{Cc} \mapsto 0.09,\, \texttt{Da} \mapsto 0.1,\, \texttt{Db} \mapsto 0.11,\, \texttt{Dc} \mapsto 0.12}\) is a non-shared parameter different for each category-value combination \((u,l)\).

1 illustrates which values each parameter assumes depending on the category-value combination using 3 for a simple dataset. The sparsity pattern of the Jacobian (consisting of the partial derivatives of the predictions \(\hat{y}_{u,l,i}\) the parameters \(P\)) required for identifying the parameters of this example is shown in 2. This sparsity pattern can be exploited during differentiation to avoid a time and memory complexity of \(O(nk)\), where \(n\) is the number of data points and \(k\) denotes the number of “individual” parameters² (\(1 + 4 + 3 + 12 = 20\) in this case), to achieve a time and memory complexity of \(O(nm)\), where \(m\) is the number of parameters (\(1 + 1 + 1 + 1 = 4\) in this case).

2.3 Incorporation into Symbolic Regression↩︎

Regardless of the SR algorithm this approach is implemented in, the terminal set has to be expanded to accommodate the three types of parameters. Partially-shared parameters are subdivided based on the categorical variables they reference, as their behavior and the number of individual parameters they contain depend on the specific category. This may affect several parts of the respective algorithm or implementation. Here, we briefly describe a selection of the changes required for a GP method.

First, the new parameter types can be included during generation of new expressions. However, they may also be omitted, leaving their introduction to subsequent mutations. We set the probability of including any but fully-shared parameters into new expression to zero.

Next, a new string representation is required for the different parameter types. For partially-shared parameters, we use C followed by the category index they refer to (, C1 for parameters depending on categorical variable \(1\)). For non-shared parameters depending on the combination of all categories, we use CI (I for interaction).

Prior research [23], [31] has shown that fewer parameter identification iterations are required if parameter values are cached between generations. Thus, in our approach, the individual values within a new parameter node are preserved separately, ensuring that genetic operators treat each scalar as a distinct unit.

Mutations which may concern terminals should be adapted for the new terminal set. We extended standard “point mutations” to allow nodes to transition between existing and new parameter types. Also, we modified parameter perturbations to independently adjust each scalar within non-shared or partially-shared parameters.

Finally, because lower levels of sharing increase the number of individual parameters, this should be either directly reflected in the complexity measure or penalized using an additional objective. We define the number of individual parameters as an additional objective in a multiobjective GP approach (see 4).

Many SR implementations allow reusing parts of the expression within itself, often referred to as directed acyclic graphs. This is particularly useful for shared expressions, as it allows a single shared parameter to be referenced in multiple locations, potentially reducing the total count of individual parameters.

2.4 More Categorical Variables and Limitations↩︎

The concept is explained and demonstrated for two categories, but is generalizable to arbitrarily many categories with increasingly more levels of intermediate sharing. For two categorical variables, the terminal set is extended by four new types of parameters. This increases the search space and does thus impact the search efficiency. However, the number of additional parameter types increases exponentially with the number of categories. Thus, including all emerging, intermediate sharing types does not scale well with the number of categorical variables.

3 Related Work↩︎

The method developed by Kronberger  [23] is closest to what we propose in this paper. Focusing on finding common behavior across similar systems, they introduce one categorical variable using category-value-specific parameters in the manner already discussed. Both fully-shared and category-value-specific parameters are allowed. Multiple categorical variables are mentioned, but are deemed unfeasible, as interactions could not be accounted for and datasets would become too sparse. They also compare this varying-parameter-sets-approach to one-hot-encoding, which is another method for including categorical variables. One-hot-encoding is deemed less interpretable because the categorical influences are distributed across the resulting expressions.

In [23], [24], [27], [28] and the implementations [29], [30], various approaches are introduced each of which incorporates only one categorical variable. While [24], [27], [30] allow for category-value-specific parameters only, [23], [28], [29] also include parameters shared across all category values.

Russeil  [25] compares [24], [28]–[30] using a common benchmark with real-world problems containing one categorical variable. They discuss benefits and drawbacks, and offer recommendations for future research in this promising field. The option of using more than one categorical variable is briefly mentioned in [25], but they do not discuss adapting their methods to accommodate for them.

4 Application↩︎

To demonstrate the proposed approach, first, we further expand on the previously introduced example 3 with the two categories \(U\) and \(L\). The presented approach may be uniquely useful for cases, where there are very few data for some of the category combinations, while there are an abundance of data for the others. We design an example to gauge the parameter sharing capabilities, and thus, the data sharing and transfer-learning limits in terms of data reduction. In essence, we seek to understand to what degree data distributed across the category combinations may be used between the categories to correctly learn the parameters for category combinations with very limited data.

3 , along with its parameters, remains the same. For \(v_1\) of each of the twelve category-value combinations, eight data points are randomly and uniformly sampled in the range of \(\interval[][]{-20}{20}\), resulting in a total of \(96\) data points. In each step of a “procession”, a data point is randomly transferred from the training set to the test set. All individual parameters \(p\) are then randomly perturbed with \[\label{eq:perturb} p_{\mathrm{perturbed}} = p + 0.1 \cdot p \cdot r \quad \text{with } r \sim \cN(0,1)\tag{4}\] before they are re-identified and the prediction is compared to the test set. For one procession, this is repeated until the minimum data requirements are no longer satisfied. Excerpts of four of a total of \(100\) conducted processions are shown in 2.

As can be seen in 2, the test data are predicted well, despite very few data points in many of the category-value combinations. For many category-value combinations, only one data point is required, as long as there are sufficient data in other category-value combinations. For example, in ID 3:26, seven out of twelve category-value combinations require only a single data point.

In this example, there are \(20\) individual parameters, which act like \(48\) parameters. Thus, in the best case, as can be seen in procession \(2\), only a minimum of \(20\) data points are required to identify the parameters. However, as can be seen in processions \(1\), \(3\), and \(4\), the minimum data requirements are no longer fulfilled despite the availability of more than \(20\) points. This stems from the fact that data points only contribute to the identification of parameters they specifically reference.

For each non-shared parameter, one point is required in each category-value combination. For a partially-shared parameter, at least one additional point is required for each value of the category upon which it depends. And finally, one additional point in any category-value combination is required to identify a shared parameter.

Whether the data distribution across the category-value combinations is sufficient according to those conditions is indicated in the “req.” column of 2. As evident in the table, as soon as the minimum data requirements are no longer satisfied, the test data can no longer be predicted well, which is reflected in all of the \(100\) conducted processions. Of course, sufficient data for the parameter identification alone does not guarantee a successful refitting, as it also depends on the distribution of the \(v_1\) values in the training dataset. To conclude this example, although there are four parameters in the expression, the parameters can be identified correctly for category combinations containing only one point. This highlights the strengths and limitations of our approach’s transfer learning capabilities.

The proposed approach incorporates multiple categorical variables in a novel manner. Although we expect a wide applicability of the proposed approach to many scientific and engineering domains, publicly available datasets containing multiple categorical variables remain scarce. The methods introduced by prior studies are limited to one categorical variable, and most of their examples do only contain one. One exception is the astrophysics example introduced by [24], [25], which we utilize for the present study.

This dataset contains two categorical variables, which the authors flattened to one by considering their combination as a single categorical variable. It describes the change of radiation flux of supernovae (three datasets) over time using two types of photometric filters (two bands), which results in three times two (equals six) unique category-value combinations.

In [24], the authors show three resulting expressions for this problem with three, four, and five parameters and \(R^2\)-values of \(0.99\), \(0.987\), and \(0.992\), respectively. However, as they only contain the present-studies-equivalent of non-shared parameters, each of those contribute six individual parameters, resulting in \(18\), \(24\), and \(30\) individual parameters. The expressions and the relevant information are shown in 3.

We implement the proposed approach on the basis of the open-source SR package [32], [33]. leverages NSGA-II to provide a multi-objective approach to GP, featuring a highly modular architecture for adding custom objectives and configurations. To evaluate different approaches to handling categorical variables, we implement five variants inside the same algorithmic framework. Four of the five variants (conventional, separate, flattened, predictive) represent established methodologies, while one (novel) is novel.

The first variant, conventional, utilizes standard SR on the entire dataset, effectively ignoring categorical distinctions. Within our framework, this corresponds to generating expressions using only fully-shared parameters. At the opposite extreme, the separate variant treats each category-value combination independently. The two categories are collapsed into one by considering the combinations of category values as a single categorical variable. While the expressions share a functional form, their parameters are entirely separate and thus non-shared. This variant aligns with methods proposed by [24], [27], [30]. The flattened variant also uses only a single flattened categorical variable, but allows both fully-shared and non-shared parameters, consistent with the approaches in [28], [29]. The predictive variant exclusively utilizes full and partial parameter sharing, omitting non-shared parameters. This variant is especially useful for zero-shot generalization; by excluding parameters tied to specific category interactions, the resulting expressions can be applied to out-of-distribution category combinations entirely absent from the training dataset. The approach [23] is aligned with this approach. However, for one categorical variable, their approach also allows realizing the flattened variant. Lastly, the novel variant incorporates all three parameter-sharing types discussed in the present work. The variants and their sharing modalities are summarized in 4.

For all experiments, we adopt the dataset and preprocessing of [24], and choose the function set as \[\set{ \texttt{+} ,\, \texttt{-} ,\, \texttt{*} ,\, \texttt{/} ,\, \texttt{\^{}} ,\, \texttt{exp} ,\, \texttt{log} , \texttt{square} , \texttt{sqrt} } .\] The selection objectives are \(1 - R^2\), the complexity (number of operators and operands), and the number of individual parameters. The maximum allowed complexity is limited to \(12\) to allow direct comparisons to [24]. Each of the five variants is executed five times, with each run lasting \(15\) minutes. The \(1 - R^2\) and number of individual parameters of the expressions found using each variant are shown in 3.

As can be seen in 3, the conventional variant only produces solutions with a small number of individual parameters. Among the tested variants, its expressions are Pareto-optimal only for two or fewer parameters. The separate variant generates expressions with parameter counts in multiples of six; this is because it relies exclusively on non-shared parameters, each of which contributes six individual parameters to the total. Compared to other variants, separate only reaches the Pareto frontier at zero, \(18\), and \(24\) individual parameters. At the zero-parameter baseline, all studied approaches perform identically. For \(18\) and \(24\) parameters, approaching upper extreme of number of parameters of this study, the performance of separate and flattened converges. However, because flattened can utilize shared parameters, it possesses greater flexibility in allocating its parameter budget. It identifies superior solutions at six and twelve individual parameters compared to separate, and provides more granular trade-offs at counts that are not multiples of six.

A similar improvement is observed when comparing the separate and novel variants. The novel approach leverages greater degrees of freedom by optimizing both the allocation of individual parameters and the choice of sharing modality. Consequently, its solutions lie almost exclusively on the Pareto frontier. Across a wide range of parameter counts, it produces expressions on par with the best-performing variants, while often uniquely identifying the best trade-offs. Finally, the predictive variant produces expressions limited to ten or fewer parameters. However, it successfully reaches the Pareto frontier at three, six, seven, and nine individual parameters.

Lastly, an excerpt of the expressions found using the variants novel and predictive are shown in 5. Although the novel variant has the flexibility to use any sharing type, it also identifies Pareto-optimal expressions that omit non-shared parameters, particularly when the search is incentivized by the parameter-count objective. These solutions are highly desirable as they provide the best trade-off between accuracy and the ability to generalize to unseen category combinations.

All the shown expressions have fewer individual parameters than the ones reported by [24], while they retain a similar \(R^2\)-value. For example, the expression in the third row of 5 uses one parameter depending on the type of band, and two non-shared parameters, which results in only \(14\) individual parameters. Although it uses four fewer individual parameters, its \(R^2\)-value is the same as that of the three-parameter expression from [24] (\(R^2 = 0.99\)), shown in the first row of 3. Comparing the best expressions of the present study and [24], the proposed approach requires half of the parameters (\(15\) versus \(30\)) to achieve a similar \(R^2\) value (\(R^2 = 0.993\) versus \(0.992\)).

An additional benefit of the proposed approach is that more information is can be extracted about the problem and the categories. By inspecting the resulting expressions, domain experts can distinguish between effects that are category-agnostic, category-specific, or unique to particular category-value interactions.

5 Conclusion↩︎

This paper presents a novel framework for incorporating multiple categorical variables into SR. By introducing fully-shared, partially-shared, and non-shared parameters, our approach significantly extends the flexibility of SR in the presence of multiple categorical variables. These sharing modalities allow the model to simultaneously capture category-agnostic trends, isolate category-specific effects, and account for unique category-value interactions. This enables the identification of additional problem structure, reduces the number of individual parameters, and lowers data requirements. Through a synthetic benchmark, we characterize the limits of data sharing and transfer learning within this framework. The results show that the proposed sharing mechanisms dramatically reduce per-combination data requirements, provided that sufficient global data exist to identify the shared parameters. This framework can be implemented in any SR architecture that allows for the expansion of the terminal set with new parameter types and the necessary adaptations to identify those parameters.

We implement the framework in a multiobjective GP-based SR algorithm and demonstrate it on a real-world astrophysics dataset used in a related study [24], [25]. Within the same framework, we implement five algorithmic variants that incorporate categorical variables in different ways. Four variants align with prior methods, while the remaining one is based on our proposed method. Comparative analysis shows that our proposed variant consistently outperforms established methods by identifying a more diverse Pareto frontier of solutions with superior trade-offs between fit quality and parameter count. The proposed approach will be merged into the publicly available branch of in the future.

References↩︎