Symbolic Regression is NP-hard
Authors: Marco Virgolin
Solon P Pissis
Venue: Accepted by TMLR ()

Abstract: Symbolic regression (SR) is the task of learning a model of data in the form of a mathematical expression. 
By their nature, SR models have the potential to be accurate and human-interpretable at the same time. 
Unfortunately, finding such models, i.e., performing SR, appears to be a computationally intensive task. 
Historically, SR has been tackled with heuristics such as greedy or genetic algorithms and, while some works have hinted at the possible hardness of SR, no proof has yet been given that SR is, in fact, NP-hard. 
This begs the question: Is there an exact polynomial-time algorithm to compute SR models? 
We provide evidence suggesting that the answer is probably negative by showing that SR is NP-hard.

Published in Transactions on Machine Learning Research (10/2022)
Symbolic Regression is NP-hard
Marco Virgolin
marco.viryokknnecui..
Centrum Wiskunde 8 Informatica, Amsterdam, the Netherlands
Solon P. Pissis
solon, solon.pissis Bowi.nd
Centrum Wiskunde & Informatico, Amsterdam, the Netherlands Vrije Universitest, Amsterdam, the Netherlands
Reviewed OpenReview: https: // openreview. net/ ptt/form?? id= id=LTaaPage2e
htdps://openreviee
Abstract
Symbolic regression (SR) is the task of learning model of data in the form of a mathematical expression. By their nature, SR models have the potential to be accurate and human- interpretable at the same time. Unfortunately, finding such models, ie., performing SR, appears to be a computationally intensive task. Historically, SR has been tackled with heuristics such as greedy or genetic algorithms and, while some works have hinted at the possible hardness of SR, proof has yet been given that SR is, in fact, NP-hard. This begs the question: Is there an exact polynomial-time algorithm to compute SR models? We provide evidence suggesting that the answer is probably negative by showing that SR is NP-hard.
1 Introduction
Symbolic regression (SR) is a sub-field of machine learning concerned with discovering a model of the given data in the form of mathematical expression (or equation) (Koza, 1994 Schmidt & Lipson 2009). For example, consider having measurements of planet masses my and m2, the distance r between them, and the respective gravitational force F. Then, an SR algorithm would ideally re-discover the well-known expression (or equivalent formulation thereof) F = G x , with G = 6.6743 x 10-11, by opportunely combining the mathematical operations (here, of multiplication and division) with the variables and constant at play.
The appeal of learning models as mathematical expressions goes beyond obtaining predictive power alone, as is commouplace in machine learning. In fact, SR models are particularly well suited for human interpretability and in-depth analysis (Ott 2013 Virgolin et al.) 2021b La Cava et al. 2021 This aspect enables a safe and responsible use of machine learning models for high-stakes socictal applications, as requested in the AI acts by the European Union and the United States (European Commission 2021 117th US Congress) 2022 Jobin et al. 2019). Moreover, it enables scientists to gain deeper knowledge about the phenomena that underlie the data. Consoquently, SR enjoys wide applicability: SR has successfully been applied to astrophysics (Lemos
ct al. 2022), chemistry (Hernandez et al. 2019)), control (Derner et al. 2020). economics (Verstyuk &
Douglas 2022), mechanical engineering (Kronberger et al.; 2018). medicine (Virgolin et. al. 2020b). space exploration (Märtens & Izzo 2022). and more (Matsubara et al. 2022).
As we will describe in Sec. 2 many different algorithms have been proposed to address SR, ranging from genetic algorithms to dlep learning ones. Existing algorithms either lack optimality guarantees heavily restrict the space of SR models to consider. In fact, there is wide belief in the community that SR is an NP-hard probleu [Lu et al. 2016 Petersen et al. 2019 Udrescu & Tegmark 2020 Li et al. 2022). However, to the best of our knowledge, this belief had yet to be solidified in the form of proof prior to the advent of this paper. Indeed, we prove that there exist instances of the SR problem for which one cannot discover the best-pos possiilos mathematical expression in polynomial time unless P-NP. Id est, SR is an NP-hard problem.
Ta 2016) that SR NP-hard hut provide reference nor proof.

---

Published in Transactions on Machine Learning Research (10/2022)
2
Background
We begin with a historical overview of how SR has been attempted from an algorithmic perspective (Sec. 2.1 and then follow with related work concerning hardness (Soc. 2.2).
2.1 SR algorithms
The introduction of SR is generally attributed to John R. Koza (c.g.. Zelinka et al. (2005) make this claim); however, the problem of finding a matheratical expression or equation that explains empirical measurements was already considered in earlier works [Gerwin] 1974 Langley [981] Falkenhainer & Michalski 1986). Such
works build mathematical expressions by iterative application of multiple heuristic tests on the data.
Koza is best known for his pioneering work on genetic programming (GP), i.e., the form of evolutionary computation where candidate solutions are
X
variable-sized and represent programs (Koza et al., 1989 Koza, 1990
1994p. Early forms of GP were proposed by Cramer [1985]; Hicklin 1986
Koza showed that GP can be used to discover SR models by encoding
mathematical expressions as computational trees (see Fig. D. In such trees,
sin
22
I3
internal nodes represent functions (e.g., +. -, x, etc.) that are drawn from pre-decided set of possibilities, and leaf nodes represent variables or constants (c.g., 21, 22, ..... -1, π, etc.). GP evolves a population of trees by init ially sampling random (rees, and then conducts the following steps: (1) stochastic replacement and recombination of their sub-trees; (2) evaluation of the fitness by exccuting the trees and asscssing their output; and (3) stochastic survival of the fittest.
Figure 1: Example of a tree that
encodes f(x) = (sin(x1) + x2) x
0s/21-
Recently, La Cava et al. (2021 proposed SRBench, a benchmarking platform for SR that includes more than 20 algorithms and more than 250 data sets. SRBench shows that several state-of- ofrrt--atttttttoe algorithms for SR are GP-based. Among these, at the time of writing, Operon by Buulaau et al. (2020) was found to perform best in terms of discovering accurate SR models; and GP- GOMEA by Virgolin et al. (2021a) was found to perform best in terms of discovering decently-accural relatively-simple SR models (i.c., shorter mathematical expressions). Other forms of GP, such as strongly- ttype GP (Montana) 1995), Haamaar-guided GP (MeKay et al., 2010), and grammatical evolution (O'Neill & Ryan 2001], are often used to tackle dimensionally-auus SR, i.c., the search of mathematical expressions with constraints to achieve meaningful combinations of units of measurement.
and
SR has been addressed with other types of algorithms than genetic ones, including, e.g., Monte-Carlo tree scarch (Cazonave] 2013 Sun et al. 2022). Morcover, several authors proposed deterministic algorithms. For example, Worm & Chiu (2013) and Kammerer et al. (2020) proposed enumeration algorithms which make SR tractable by restricting the space of possible models to consider and including dynamic programming and pruning strategies. [ozzad] (2014); Cozad & Sahinidis (2018) showed how SR can be addressed with mixed integer nonlinear programming. McConaghy (2011) proposed FFX, which generates linear combination of many functions that are linearly-indepeenentt from each other, and then fits its coefficients with the elastic net (Zou & Hastie 2005) to promote sparsity. Olivetti de França (2018) and Rivero et all (022) propose greedy algorithrs that start from small ma mathematical expressions and iteratively expand them, by replacing existing components with larger ones from a set of possibilities.
Lastly, recent ycars have seen the proposal of deep lcarning-based algorithms for SR. Peterson et al. (0220) cast the SR problem as reinforcement learning one and train a recurrent neural network to generate accurate SR. models. [ddrescu& & Tegmark (2020) leverage neural networks in order to test for symmetries and invariances in the data that are then used to prune the space of possible SR models. An end-to-end approach is taken by Kamienny et al. (2022) and Vastl et al. 2022). , who train deep neural transformers to produce SR
models directly from the data. Li et al. (202) seek SR models by proposing a convexified formulation of deep reinforcement learning.
In summary, existing SR algorithms are either heuristics, which do not guarantee optimality (e.g., genetic, greedy, or deep learning-based algorithms), they are exact algorithms that achieve optimality but only over
2

---

Published in Transactions on Machine Learning Research (10/2022)
small subset of all possible SR models, to bound the runtime (e.g., dynamic programming and mixed-integer nonlinear programming algorithms). This strongly hints to the fact that SR is NP-hard. As mentioned earlier, no proof has yet been given.
2.2 Related hardness results
SR is typically posed as an empirical risk minimization (ERM) problem or, when regularization is considered, structural risk minimization one [Vapnik 1999). There exists a multitude of theoretical results in machine learning posing the problem as an ERM one. For example, Blum & Rivest (1992) famously proved the NP-completeness of training a three node-neural network to label a given data set correctly. The loss employed in such situations commonly the 0-1 loss, i.c., the loss that returns (or prediction) equals the output that. is expected from the data (or label) and 1 otherwise. More recently, it has been shown that, under the 0-1 loss, it is NP-hard to even train linear classifier to be e-better than random (ellmman et al. 2012). In the context of coding theory, ERM with 0-1 loss has been used to prove the NP-hardness of finding a univariate polynomial of maximum degree k over a finite field (a code) that maximizes code matching (Guruswami & Vardy 2005).
the output of the model
Under different types of loss, polynomial-time solutions to ERM exist. A famous example of this is lincar regression under the squared error loss, which can be solved in polynomial time via ordinary least squares. Recently, Backurs et al. (2017) presented fine-grained complexity results for ERM with kernel-based and neural network-based approaches.
SR can be set to search in the space of polynomials and one can choose to use the 0-1 loss. Moreover, one can in principle set SR to work in finite fields rather than on real numbers to operate with polynomials for discrete codes. For example, Koza (1990) shows how to set GP to learn Boolean circuits by composition of logic gates. This means that the hardness of SR. can follow from linking back results such the one by Feldman et al. (2012) (we sketch how this can be achieved at the end of Sec. The modern connotation of SR focused regrossion (i.c., we seek model f: IRd R with d the number of features) and commonly-used losses have co-domain in IRF such as the absolute error loss or the squared error loss; rather than the 0-1 loss. We will provide a proof of NP-hardness that is general to this type of losses (Eq. 2 ). In essence, we will show that when certain basic arithmetic operations are chosen for combining (e.g,, addition), SR becomes NP-hard. This choice of basic operations allows us to reduce from a rather classical variant of the knapsack problem, namely, the unbounded subset set problem (Kellerer et al. 2004).
3 Preliminaries
We will hereon refer SR models as functions when appropriate, as this their fundamental nature. Functions take variables as arguments. One value of the variable taken as its argument. For simplicity of exposition, and not make distinction between (non-identity) functions and variables. Similarly, we will refer to functions also for constant functions, i.., functions that can only return a single numerical value, irrespective of their arguments. This said, let us recall the concept of function composition, which is central to SR.
the identity function, i..., the function that returns the will generically refer to functions
Desinition Function composition
Given two functions f A + and g IB + C, function composition, which we denote by g o f, is the operation that produces third function h A > C, such that h(x) = (xf())).
Thanks to function composition, can now define the concept of search space of an SR problem. Definition 2. Search space of SR
Let P be a set of functions. The search space of SR is the function space F. that contains all functions that can be formed by composition of the elements of P.
To better understand what Def. states, consider that can be set to contain a mix of functions that perform basic algebraic operations such as addition, subtraction, multiplication, and division; transcendental functions such sin, cos, log, exp; constant functions (or simply constants), such ca2(x) = 42 and ea (I)) æ for
3

---

Published in Transactions on Machine Learning Research (10/2022)
any 2; and identity functions that represent variables of interest for the problem at hand, such z1, 42, T3. P is typically referred to as the primiliiee sel, and its elements as primillues (Poli et al. 2008). Once P has been decided, F is determined. For example, choosing P = {+(s-),-(+-)), x (-. ), 21, 22, - 1,, +1} means that FF will contain a subset of all possible polynomials of arbitrary degree in œ1 and Z2. In particular, F is a subset because only some coefficients can be expressed, by composing constants with addition, subtraction, and multiplication.
Let us clarify point regarding constants in particular. Normally, one would include constants which are relevant to the instance of SR at hand. For example, if the unknown phenomenon for which an SR model is sought is suspected to have sinusoidal components, it may be advisable to include multiples of i in P. Moreover, P can be set to contain special elements that represent probability distributions from which constants can be sampled (see the concept of ephemeral random constant described by Koza (1994); Poli et al. (2008)). We denote one such element by R and, e-g.; R can be chosen to represent the uniform distribution between two numbers, or the normal distribution with a certain mean and variance. When an SR algorithm picks 9t from P to compose an SR model, constant is sampled from the distribution identified by S. Here (more specifically, in Corollary D we will generously assume that any constant can be sampled directly from R, and therefore that there is no need for a real-valued optimizer to be part of the SR algorithm. For example, having P = (H)))))((ss))) ((, )) æ1, 22, R} R2 21 and d22
will mean that F contains all polynomials of arbitrary degree in
We can now proceed by providing a definition of the SR problem. While this definition can be extended to other domains, we focus on handling real-valued numbers as the majority of the works takes place in this domain, and subsets thereof.
Definition Symbolic Regression (SR) problem
Given set P of fanctions, distance L: R R - R vectors x, = (.1......., xa,i) E IRd and scalars ys E R. for i = .,....n, n, the SR problem asks for finding function S such that:
1 n
f* E arg min JeF
Lyn
(1)
where F is the search space that is defined by P.
We provide some remarks concerning the proposed definition of the SR problem. Firstly, let us map the objects provided in the definition to terms familiar to a machine learning audience. The pair (x; ys) is normally what is referred to as observation, dala poinl, example, or sample, where æj, is the value of the jth feature or variable for the ith observation, and Yi the value of the label or target variable for the same observation. The set that contains the observations upon which is computed, i.e., D {(x,, y:) is called training set. It is commonly assumed that the observations in D were drawn independently and identically distributed (i.i.d.) from an unknown probability distribution. Moreover, the distance & is called loss. Losses need not be distances, but in SR they normally are. Popular choices in the literature are the absolute error loss and the squared error loss. Here, we consider losses of the form:
Cly, = ly, - f(x))),
(2)
with r E Ng: for r = 0 one gets the 0-1 loss; for 1 one gets the absolute error loss; and for r = 2 one gets the squared error loss.
The minimization of the loss function across the observations in D makes SR an ERM problem. As is generally the case for learning, one actually desires f to generalize to new (or also called unseen) observations, i.e., observations that come from the same underlying probability distribution but were not available in D. In other words, it is not sufficient that f* is a best-possible function with respect to the training set, as the loss should remain minimal also for new observations that are not available to us. To this end, a common practice is to heuristically use a soparate set of data (the validation set) to estimate the generalization to observations outside the training set. Another approach, which often used together and not alternative to adopting validation set, is to perform structural risk minimization, i.e., account. for regularization terms such
4

---

Published in Transactions on Machine Learning Research (10/2022)
À C(f), where X R. controls the regularization strength and C: R is function of the complexity of f. Typical goals of such regularization terms are improving generalization (by limiting effects akin to Runge's phenomenon (Fornberg & Zuev 2007)) or, particularly for SR, improving the interpretability of f. Implementations C' range from weighted counting of the number of primitives that constitute f (Ekart & Nemeth 2001; Hein el. al.) 2018) to machine learning models trained from human leedback to predict s's interpretabllity (Virgolin et al., 2020a 2021b). Here, for simplicity, we do not consider regularization and focus on ERM alone. Equivalently put, we consider À O: note that this choice is not limiting because this step can also be taken when constructing the proof of Theorem g We will briefly get back to how À >> 0 and C' may be used to build interesting search spaces for the hardness of SR at the end of Sec. 4
Still, considering a "pure optimization" formulation (or ERM), as given in Eq. [ can be considered to be a pre-requisite for being able to machine-leari accurate models from the data; in fact, is commonplace for literature that concerns the hardness of learning to provide results with respect to the training set (see, e.g.; [ellmaan et al.) 2012; Hu et al., 2019)). In a similar fashion, here we will consider the case of minimizing the empirical risk with respect to the training set. D and show that this alone already poses challenge for any SR. algorithm.
Here, we assume that computing f(x) and C(y, fx)) (see Def. 3 can be done in polynomial time. Regarding C, our assumption is met for commonly-used losses such as the absolute and squared error ones. In fact, computing losses of such form takes O(n) operations, i.e, the runtime is linear in the number of observations. Regarding f, our assumption met, e.g., for all functions that can be discovered by the SR algorithuns currently in SRBench (La Cava et al., 2021); a notable exception of practical interest are recursive functions taking exponential time to compute (see, e.g., Assoli et al. (2022) ). For non-recursive functions, f can be implemented directed acyelic graph, where nodes represent the functions from P, and edges represent compositions. To compute fx), it suffices to visit each node of the graph for each observation, thus requiring O( x n) operations, where L is the number of primitives in f. Fig. [ shows an example of such a graph, especially in the form of a tree, which is perhaps the most common way of encoding mathematical expressions in SR (see, e.g.; the SR algorithms benchmarked by La Cava et al. (2021)
We conclude this section with the following important definition.
Definition 4. Decision version of the SR problem (SR-Dec)
Given an SR instance and an e E R SR-Dec outputs YES if and only if:
1 a Ef E F: n
Σ c (y:, f(xi)) <c
(3)
Essentially, Def. 4 is the problem of deciding whether there exists function in the search space such that its empirical risk is smaller than chosen threshold
4 The result
We proceed directly by providing the main result of this paper.
Theorem 1. The SR problem is NP-hard.
Proof. Let us begin by stating that SR-Dec is in NP. Recall that the computations of f(x) and L(y, f(x)) take polynomial time (sce Sec. Of course, the check < c takes O(1) time. Thus, if is guessed by an oracle, then we can provide an answer to SR-Dec in polynomial time.
We proceed by considering the unbounded subset sum problem (USSP). USSP is a similar problem to the unbounded knapsack problem, where same item can be put in the knapsack an arbitrary number of times, and the weight of an item corresponds exactly to the profit gained by including that item in the knapsack. The decision version of USSP, USSP-Dec, is defined follows. Given 1,...,K k. (I: items), wy N (weight.
5

---

Published in Transactions on Machine Learning Research (10/2022)
of that item), and t e N (the target); USSP-Dec asks:
3m Σ mymy i?
(4)
where m, E No (multiplicity with which an item is picked). USSP-Dec is known to be NP-complete (Kellerer et al. 2004).
To prove that SR-Dec is NP-complete, we show that any instance of USSP-Dec can be reduced to some instance of SR-Dec in polynomial time. To this end, we will restrict SR-Dee as follows: (1) We pick the set of primitives P to be P = {+, (+..1...., raa, with d = k; (2) We set c = 1. Next, we craft. 'D to have a single observation (n = 1) and = k features. For the only observation in D (dropping the index for the observation number, since there is only one), we set 21 W1, = W2 ................................ Ix Wk: and y =tt.
In other words, we have set the search space F to contain only linear sums of the features in the data set DD i.e, functions of the form f(x) Σド= K, Importantly, m E No and Fa E N, meaning that the co-domain of any is No- Consequently, the smallest non-zero loss that can be achieved is 1. The only functions f that can achieve an error smaller than = 1 are those that interpolate the observation exactly, i.e., f(x) =
Then, the following holds:
(Eq. 13 with 1)
SfEF:LW.ff <
(5)
(C (v, f(x)) < 1 + f(x) 4) Ff yy
(6)
(Equivalence y = t due to D)
Bf F f(x) =t?
(7)
(Expanding F based on choice of P) B
Tymj my E No
f((x)= tt?
(8)
(Equivalence I, Wy. d = k due to D)
Wymj my E No
: f((x)= =t?
(9)
1=1
(Re-formulating in terms of m)
Em Σ
wymy = t?
(10)
J=1
In other words, there exist some instances of SR-Dec that can be re-formulated as USSP-Dec (cfr. Eqs. and (10 Now, since assembling P as stated above takes linear time in k picking = 1 takes (1) time, and constructing DD as stated above takes linear time in k, then any instance of USSB-Dec can be reduced to some instance of SR-Dec in polynomial time: SR-Dec is NP-complete.
We conclude the proof with a reductio ab absurdum. Let us assume that there exists an algorithm to compute an optimal f* for the SR problem (Def. 3) in polynomial time. An optimal f the one for which the loss is minimal, which means that using f Eq. 3 allows us to immediately answer SR-Dec. Since verifying that Cly, (x)) < takes polynomial time, we conclude that if the SR problem can be solved in polynomial time, then we can also solve SR-Dec in polynomial time. Therefore, the SR problem is Na-hard.
We remark that, in the proof of Theorem [ we construct P. so as not to contain Ri (nor any constant). Some readers might disagree with this quite broad definition of SR. In fact, some SR algorithms heavily rely on the presence of constants as well as on their optimization (e.g., FFX by McConaghy (201I) and FEAT by La Cava et al. (2018)). Not allowing for arbitrary constants be present in the functions of the scarch space might be seen as violation of the very definition of SR. In other words, some might think that P must contain R. We next show that SR remains NP-hard in this special case.
Corollary 1. The SR problem is Na-hard even when P must include R.
6

---

Published in Transactions on Machine Learning Research (10/2022)
Proof. We follow similar construction of the proof of Theorem | This time, set P to additionally contain R, Le., P (t, ,....., ) with d = k. This means that the function space F now contains functions of the form f(x = c 2ッj1 3ymy with mj E No and E I (sampled from R). As to D, we will now include two observations instead of single one. The first observation is set before, i.e., 11,1 = #2,1 W22-- ................................ Fk,1 = (d = k) and y1 = L. For the second observation, we set 2,2 0, = O,.... 2x, = 0 and now remains to determine how we should set 6.
.2....,,
0, ie., the value of all features and of the label are set to zero. It
Because of our construction of D, an I for which SR-Dee answered YES in the situation considered Theorem 1 i.c., with c= O, would still make SR-Dec answers YES for a with sufficiently small magnitude. Note that those functions interpolate D exactly if e 0. Thus, the magnilude of must. be such that lel" < e, with r the degree of the loss (Eq. (2)), because the corresponding empirical risk for those functions is 2/d" + lll"" = |cl", where the two summands on the left side of the equation are respective to the two observations in D.
Now the question becomes whether using c * 0 allows to answer YES to more functions than those that would interpolate when c c= 0. If that would be true, then we can no longer apply the strategy used in the proof of Theorem | to reduce from USSP-Dec. To be able to still use that strategy, we will now show that there exist instances of SR-Dee, in particular by picking a different e, such that even if P must include R, then only functions with 0 are candidates for YES answer. This would allow us to reduce once again from USSP-Dec (by appropriately picking c), because f(x) = 0 + Σj= m, Ei= Tymy. as in Eq. (8).
For 7 0 to allow SR-Dec to answer YES to more functions than those for when c = 0, e must contribute to lower the empirical risk.
For the second observation, e can only increase the risk, because the loss is:
-e(ett d
L(y2, f(x2)) 192 [(x2)] 0 c+ Σ my
= |c|--
(11)
For the first observation, however, using e * O can lower the respective loss and thus contribute to lowering the empirical error. In particular, consider that the first functions that are candidates to receive a YES answer thanks to + 0 are those that had a loss of 1 when c 0 (sce the proof of Theorem where we set c= 1). If we can have SR-Dec answer NO to these functions, then will necessarily answer NO also to all other functions that have loss larger than on the first observation when = O. We thus proceed by considering that the smallest, non-zero loss that can be obtained the first observation for c + 0, is |1 - -c".. This leads to the following empirical risk over our D:
Σ C(u,, f(x:)) = (11 -c" -10 -0-c))).
(12)
2
2
For any integer 2 0, the minimum is - - this is easy verify for r 0 and 1, while for Eq. (12) describes U-shaped curve that is symmetric around, and has minimum in c= }/- Therefore, it suffices to pick 2-;; even if is optimal, one cannot lower the empirical risk below 2 for any function whose loss (on the first observation) is not zero. In other words, imposing e = 2- ensures that SR-Dec will answer YES iff SR-Dec answers YES also for c= 0. This means that one can set = 0 and proceed with a reduction from USSP-Dec as in the proof Theorem [
Finally, we provide the following remark.
Remark. One can consider the structural risk ménimization selling whereby the following minimization is sought:
n E f(x)) (C(FF, with 0.
(13)
Then, SR-Dec can be restricted to automatically answer NO for any f that does not satisfy certain conditions, such as linearily. For example, one can pick P {+, X, Fd, R} to search in the space of arbilrary
7

---

Published in Transactions on Machine Learning Research (10/2022)
polgmomials, and pick C. such that C(f) 3o if degf k else for an arbitrary integer k 2 0. Combining this with the use of the 0-1 loss function enables to reduce SR-Dec from existing results in literature that consider linear classifiers (Feldman et al.) 2012)) or coding polynomials in finite fields (Guruswami & Vardy 2005).
5 Conclusion
Our main contribution here prove that symbolic regression (SR), i.e., the problem of discovering accurate model of data in the form of mathematical expression, is in fact NP-hard. In particular, we have provided formal definitions of what SR entails, and showed how the decision version of the unbounded subset sum problem can be reduced to a decision version of the SR problem. Except for the general definition of SR we considored, we have additionally shown that SR remains NP-hard oven when the set of primitives must contain distributions from which constants can be sampled, and provided a sketch of how an alternative proof can be constructed by using the 0-1 loss and previous results from the literature.
Having settled the matter the hardness of SR, hope that this note inspires further works lower and upper bounds of different SR variants. In fact, while we have shown that hardness holds in principle (by picking search space suitable for reduction from USSP-Dec), there might exist specific variants of SR (e-g., different search spaces, specific regularization terms, or specific type of data) that are more commonly encountered in practical applications. For such more specific variants, proving hardness or designing polynomial-timg algorithms would complement the SR status quo, which mostly focuses on heuristic algorithmic design. Cltimately, theoretical advances may greatly help advancing our knowledge of what is possible with SR.
References
117th US Congress. Algorithmic accountability act 2022. URL 117th- hathh-oogrresssoooo hicerres//ouuuee111/6660000 bi11/6580/
Arturs Backurs, Piotr Indyk, and Ludwig Schmidt. On the fine-grained complexity of empirical risk minimization: Kernel methods and noural networks. Advances in Newral Information Processing Systems, 30, 2017.
Avrim L Blum and Ronald Rivest. Training a 3-node neural network is NP-complete. Neural Networks, 5 (1):117-127, 1992.
Bogdan Burlacu, Gabriel Kronberger, and Michael Kommenda. Operon C++ an efficient genetic programming framework for symbolic regression. In Proceedings of the 2020 Genelic and Evolulionary Compulation Conference Companion, pp.
2020.
Tristan Cazenave. Monte-Carlo expression discovery. International Journal Artificial Intelligence Tools, 22(1):1250035, 2013.
A Cozad. Dala-and theory-driven lechniques for surrogale-bused oplimization. PhD thesis, Ph. D. thesis, Department of Chemical Engineering, Carnegie Mellon University, 2014.
Alison Cozad and Nikolaos V Sahinidis. A global MINLP approach to symbolic regression. Mathematical Programming, 170(1):97-119 2018.
Nichael Lynn Cramer. A representation for the adaptive generation of simple sequential programs. Proceedings of the First International Conference on Genetic Algorithms, pp. 183-187, 1985.
Stéphane d' Ascoli, Pierre-Alexandre Kamienny, Guillaume Lample, and François Charton. Deep symbolic regression for recurrent sequences. arXiv preprint arXiv:2201.0 00X::2200..04600, 2022.
Erik Derner, Jifí Kubalik, Nicola Ancona, and Robert Babuška. Constructing parsimonious analytic models for dynamic systems via symbolic regression. Applied Soft Compuling, 94:106432, 2020.
8

---

Published in Transactions on Machine Learning Research (10/2022)
Aniko Ekart and Sandor Z. Nemeth. Selection based the pareto nondomination criterion for controlling code growth in genetic programming. Genetic Programming and Evolvable Machines, 2(1):61-73, 2001.
European Commission. Artificial intelligence act, 2021. URL httpg://arrrfffccaaaaaelllleeeea
Brian C Falkenhainer and Ryszard Michalski. Integrating quantitative and qualitative discovery: The ABACUS system. Machine Learning, 144):367 401, 1986.
Vitaly Feldman, Venkatesan Guruswami, Prasad Raghavendra, and Yi Wu. Agnostic learning of monomials by halfspaces is hard. SIAM Journal on Computing, 41(6):1558-1590) 2012.
Bengt Fornberg and Julia Zuev. The runge phenomenon and spatially variable shape parameters in rbf interpolation. Computers & Mathematics with Applications, 54(3):379-93 398, 2007.
Donald Gerwin. Information processing, data inferences, and scientific generalization. Behavioral Science, 19 (5):314 325, 1974.
Venkatesan Guruswami and Alexander Vardy. Maaiimmm-lkeelhhon decoding of Reed-Solomon codes is NP-hard. IEEE Transactions on Information Theory, 51(7):2249-222 2256, 2005.
Daniel Hein, Steffen Udluft, and Thomas Runkler. Interpretable policies for reinforcement learning by genetic programming. Engineering Applications of Artificial Intelligence, 76:158-166 2018.
Alberto Hernandez, Adarsh Balasubramanian, Fenglin Yuan, Simon AM Mason, and Tim Mueller. Fast, accurate, and transferable many-body interatomic potentials by symbolic regression. npj Computational Materials, 5(1):1-11, 2019.
Joseph F Hicklin. Application of the genetic algorithm automatic progrom generation. PhD thesis, University of Idaho, 1986.
Xiyang Hu, Cynthia Rudin, and Margo Seltzer. Optimal sparse decision trees. Advances in Newral Information Processing Systems, 32, 2019.
Anna Jobin, Marcello lenca, and Effy Vayena. The global landscape of Al ethics guidelines. Nature Machine Intelligence, 1(9):389 399, 2019.
Pierre-Alexaudre Kamienny, Stéphane d' d'Ascoli, Guillaume Lample, and François Charton. End-Lo-end symbolic regression with transformers. arXiv preprint arXiv:2204 10532, 2022.
Lukas Kammerer, Gabriel Kronberger, Bogdan Burlacu, Stephan M. Winkler, Michael Kommenda, and Michael Affenzeller. Symbolic Regression by Exhaustive Search: Reducing the Search Space Using Syntactical Constraints and Efficient Semanlic Structure Deduplication, pp. 79 99. Springer International Publishing, 2020.
Hans Kellerer, Ulrich Pferschy, and David Pisinger. Introduction to NP-Completeness of knapsack problems. In Knapsack Problems, pp. 483-493. Springer, 2004.
John R Koza. Genetic programming: A paradigm for genetically breeding populations of computer programs to solve problems, volume 34. Stanford University, 1990.
John R Koza. Genetic programming as means for programming computers by natural selection. Statistics and Computing, 4(2):87-112, 1994.
John R. Koza et al. Hierarchical genetic algorithms operating on populations of computer programs. In International Joint Conference on Artificial Intelligence, volume 89, pp. 768 774, 1989.
Gabriel Kronberger, Michael Kommenda, Andreas Promberger, and Falk Nickel. Prodicting friction system performance with symbolic regression and genetic programming with factor variables. In Proceedings of the Genetic and Evolutionary Compulation Conference, pp. 1278 1278-1285, 1285, 2018.
9

---

Published in Transactions on Machine Learning Research (10/2022)
William La Cava, Tilak Raj Singh, James Taggart, Srinivas Suri, and Jason H. Moore. Learning concise representations for regression by evolving networks of trees. In International Conference on Learning Representations, 2018.
William La Cava, Patryk Orzechowski, Bogdan Burlacu, Fabricio Olivetti de Franca, Marco Virgolin, Ying Jin, Michael Kommenda, and Jason H Moore. Contemporary symbolic regression methods and their relative performance. In Advances in Newral Information Processing Systems Datasets and Benchmarks Track, 2021.
Pat Langley. Data-driven discovery of physical laws. Cognilive Science, 5(1):31 54, 1981.
Pablo Lemos, Niall Jeffrey, Miles Cranmer, Shirley Ho, and Peter Battaglia. Rediscovering orbital mechanics with machine learning. arXiv preprint Xii:2202000300 2022.
Haoran Li, Yang Weng, and Hanghang Tong. CoNSoLe: Convex neural symbolic learning. or preprinl arXiv:2206.002. 2022.
Qiang Lu, Jun Ren, and Zhiguang Wang. Using genetic programming with prior formula knowledge to solve symbolic regression problem. Computational Intelligence and Neuroscience, 2016.
Marcus Märtens and Dario Izzo. Symbolic regression for space applications: Disferentiable cartesian genelic programming powered by multi-objective memetic algorithms. arXiv preprint orXiv:2206.062) 2022.
Yoshitomo Matsubara, Naoya Chiba, Ryo Igarashi, Tatsunori Taniai, and Yoshilaka Ushiku. Rethinking symbolic regression datasets and benchmarks for scientific discovery. arXie preprint arXiv:2206. 10540, 2022.
Trent McConaghy. FFX: Fast, scalable, deterministic symbolic regression technology. In Genetic Programming Theory and Practice IX, pp. 235-260. Springer, 2011.
Robert I McKay, Nguyen Xuan Hoai, Peter Alexander Whigham, Yin Shan, and Michael O'neill. Grammar- based genetic programming: survey. Genetic Programming and Evolvable Machines, 11(3):365 396, 2010.
David J Montana. Strongly typed genetic programming. Evolutionary Computation,
1995.
Fabrício Olivetti de França. greedy search tree heuristic for symbolic regression. Information Sciences, 442-443:18 32, 2018. ISSN 0020-0255.
Michael O'Neill and Conor Ryan. Grammatical evolution. IEEE Transactions on Evolutionary Computation, 5(4):349 358, 2001.
Clemens Otte. Safe and interpretable machine learning: A methodological review. Compulational Intelligence in Intelligent Data Analysis, pp- 111-122, 2013.
Brenden K Petersen, Mikel Landajuela Larma, T Nathan Mundhenk, Claudio P Santiago, Soo K Kim, and Joanne T Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. arXie preprint arXin:1012.0487 2019.
Brenden K Petersen, Mikel Tandajuela Tarma, Terrell N Mundhenk, Claudio Prata Santiago, Soo Kyung Kim, and Joanne Tacry Kim. Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients. In International Conference on Learning Represenlations, 2020.
Ricardo Poli, William B Langdon, and Nicholas F McPhee. A Field Guide to Genetic Programmingg Lulu Press, 2008.
Daniel Rivero, Enrique Fernandez-Blanco and Alejandro Pazos. DoME: A deterministic technique for equation development and symbolic regression. Expert Systems with Applications, 198:116712, 2022.
Michael Schmidt and Hod Lipson. Distilling free-form natural laws from experimental data. Science, 324 (5923):81-88) 2009.
10

---

Published in Transactions on Machine Learning Research (10/2022)
Fangzheng Sun, Yang Liu, Jian-Xun Wang, and Hao Sun. Symbolic physics learner: Discovering governing equations via monte carlo tree search. arXiv preprint arXi::2266131334 2022.
Silviu-Marian Udrescu and Max Tegmark, AI Feynman: A physics-inspired method for symbolic regression. Science Advances, 6(16):eany6631 2020.
Vladimir Vapnik. The nature of slatistical learning theory. Springer, 1999.
Martin Vastl, Jonás Kulhánek, Jirí Kubalik, Erik Derner, and Robert Babuška. SymFormer: End-to-end
symbolic regression using transformer-based architecture. ar preprint arXiv:2205. 15764, 2022.
Sergiy Verstyuk and Michael Douglas. Machine learning the gravity equation for international trade. Available at SSRN 4053795, 2022.
Marco Virgolin, Andrea De Lorenzo, Eric Medvet, and Francesca Randone. Learning a formula of inter- pretability to learn interpretable formulas. International Conference on Paralled Problem Solving from Nature, pp. 93. Springer, 2020a.
Marco Virgolin, Ziyuan Wang, Tanja Alderliesten, and Peter A N Bosman. Machine learning for the prediction of pseudorealistio ic pediatric abdominal phantoms for radiation dose reconstruction. Journal of Medical Imaging, 7(4)0046501, 2020b.
Marco Virgolin, Tanja Alderliesten, Cees Witteveen, and Peter A N Bosman. Improving model-based genetic programming for symbolic regression of small expressions. Evolutionary Computation, 2922):211-237,
2021a.
Marco Virgolin, Andrea De Lorenzo, Francesca Randone, Eric Medvet, and Mattias Wahde. Model learning with personalized interpretability estimation (ML-PIE). In Proceedings of the Genetic and Evolutionary Computation Conference Companion, pp. 1355-1364, 2021b.
Tony Worm and Kenneth Chiu. Prioritized grammar enumeration: Symbolic regression by dynamic pro- gramming. In Proceedings of the Genetic and Evolutionary Computation Conference, pp. 1021-1028, 2013.
Ivan Zelinka, Zuzana Oplatkova, and Lars Nolle. Analytic programming symbolic regression by means of arbitrary evolutionary algorithms. International Journal of Simulation: Systems, Science and Technology, 6(9):44-56 56, 2005.
Hui Zou and Trevor Hastie. Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(2))301 320, 2005.
11


Reviews:
[{"review_id": "IjqDd9wwpm", "reviewer": "Reviewer_KHax", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n\n1. The paper is well-written. The discussion of related work is comprehensive. A strong case is built for the need to learn interpretable symbolic models that can be analyzed for safety. The paper also identifies papers where the NP-hardness of symbolic regression has been hinted or mentioned before. A detailed discussion on symbolic regression, SRbench and deep learning based SR is presented. \n\nWeaknesses\n\n1. The reviewer was a bit surprised that no prior work has proven symbolic regression to be NP-hard. It appears at the end of Page 3 and start of Page 4 that the SR problem definition (which would ideally differentiate between training on a set of data and generalizing to a validation set) is replaced with just the minimization of the loss function. At this point, prior work such as https://arxiv.org/abs/cs/0405005 become an instance of (1) and the result of say, \"Maximum-likelihood decoding of Reed-Solomon Codes is NP-hard\" (\"computing the degree k polynomial that disagrees with the minimum number of input point\") would imply the main result in the paper by suitably defining the L function. So, the reviewer is not convinced the main result in the paper is new. \n\n2. The reduction from unbounded subset sum problem to SR is interesting. But whether epsilon is part of the problem definition of SR or part of the solution is debatable. One can argue that the goal is to find f that minimizes epsilon, but the problem definition itself does not set the value of epsilon (something that is true in practice). Under this argument, it would not be okay to pick epsilon to be 0 when translating USSP-Dec to SR. \n\n3. See some specific requests on presentation in \"Requested Changes\". \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper studies the complexity of the problem of symbolic regression from data to learn a mathematical expression and shows that this problem is NP hard. While it has been widely believed that this problem is NP-hard (Lu et al 2016, Peterson et al 2019), the paper claims to be the first to provide a formal proof for the NP-hardness. \n\n", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "Definition 2 can be better worded: \"formed by compositions of elements of P and their compositions\" - composition is naturally treated as being recursive and unless, one wants to restrict composition to just two (or fixed number of) levels - there is no need to mention compositions twice. Page 4 anyway talks about restrictions on the space of f (in particular, rules out recursion) later. \n\nIn Definition 3, is P including the ephemeral random constant / distribution representing element R or is it just the set of functions and variables. If it is just the later (as the wording of Definition 3 suggests), the discussion of R before Definition 3 appears to be a distraction that can be eliminated from discussion. \n\n\n\n\n ", "broader_impact_concerns": "There are no broader impact concern for the paper. The authors are advised to add a broader impact statement building on their argument about the benefits of symbolic model in the introduction. But even without this statement, the broader impact of studying the complexity of symbolic regression is obvious. ", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "M7uDhTYRZC", "reviewer": "Reviewer_tATo", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "+ Well written, accessible even to non-expert readers.  The proofs are easy to follow.\n\n+ Good coverage of the SR literature.\n\n+ The mathematical setup and the proofs themselves seem sound to me.\n\n+ The contribution seems to be novel.\n\n- Missing discussion of significance and expected impact of this result.  The authors mention that the paper settles a well-known conjecture, which is true.  In a sense, the paper serves to retroactively justify the usage of genetic algorithm and other heuristic search procedures for this problem.  However, the significance of doing so is a bit unclear.  A more detailed discussion of the limits of this result and of its dependency on the choice of search space would have helped to give this result more depth.\n\n- Lacks connection to classical hardness results in statistical learning.  Specifically, considering that SR is cast essentially as an empirical risk minimization problem, I was expecting the authors to explore the connection between their results and the well-known hardness of ERM.  Could ERM be reduced to SR?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors prove that symbolic regression (SR) is NP-hard by reduction to the unbounded subset sum problem.  The reduction holds also for cases in which constants (and distributions) appear in the search space.", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "- Highlight the significance of your result by briefly discussing what consequences you expect it to have / what benefits you expect it to bring to research on SR.\n\n- Briefly discuss how the two theorems depend on the choice of search space.  The reduction seems to show that allowing for just sum, product, and equality is enough to trigger NP-hardness, but it doesn't say much for SR cases in which these operations are not allowed.  Please briefly clarify this point.\n\n- Briefly discuss links to classical hardness results in ML.\n\n- p 2: \"Early forms of GP where proposed\" -> were.\n\n- p 3: Search space of SR: The definition depends on the notion of composition between functions, which is defined, and on that of composition between functions and variables, which is not.  Please clarify.\n\n- p 4: \"i.e., distance\" -> a distance", "broader_impact_concerns": "This paper poses no ethical concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "mk5uGB2vIZ", "reviewer": "Reviewer_YHbt", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n\n* It is an important result in showing that SR is NP-hard\n* The paper is very well written and easy to follow\n* The reduction proof seems to be correct\n\n### Weaknesses\n\n* Two simplifications were assumed (i.e., no regularization and no recursive functions). However, I agree with the authors that this does not diminish the importance of their result.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Symbolic regression is one of the central techniques for interpretable and verifiable ML in many applications. So far, using evolutionary algorithms (EA) is one of the most successful approaches. However, with the black-box view of EAs, there comes the question of whether a polynomial algorithm could exist to find an optimal solution. The authors of the paper at hand show that by reduction of the unbounded subset sum problem, the symbolic regression problem is NP-hard and thus, it is unlikely that an optimal solution can be found efficiently. ", "claims_and_evidence_yn": "", "claims_explanation": "", "audience_yn": "", "audience_explanation": "", "requested_changes": "A bit nitpicking on the wording: \"this is already problematic for any SR algorithms\". I wouldn't say that this is problematic because SR algorithms are nevertheless fairly successful. Maybe let's say that it already poses a challenge to do this on the training set efficiently. \n\nIs it guaranteed that there is a unique solution to Equation 1? If not, I wonder whether it should be f^* \\in \\argmin since the task would not to find all solutions, but only one of them. \n", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The technical tools used in the proof do not seem especially complex
or novel, and the conclusion is also not surprising. That said, the
reviewers and I feel that a formal proof of a key fact about symbolic
regression is of inherent value. The paper is also well-written, and
the discussion of related work is quite substantive. The authors also
did a good job of revising the paper based on the reviewers' comments.
Given all this, I am recommending acceptance as is.
Decision: Accept (poster)