# Symbolic Regression

7+ years ago I posted a question here regarding symbolic regression in Mathematica. At that time there did not seem to be much built-in support for this. I'm interested in multi-input models of the form y=f(x1,x2,x3,x4,x5) where I need to find the best fit "f" given a data set.

Would Mathematica accept this FindFormula[data,{x1,x2,x3,x4,x5}] ?

Here is the link to my original question: Does anyone have experience with performing symbolic regression using Mathematica?

Has any meaningful progress been made during the last 7 years in this regard in Mathematica ? I'm aware of TuringBot. Is that considered the state of the art in symbolic regression ?

Thanks for any knowledge on this subject.

Clarification and Concrete Example:

Here is a dataset consisting of 27 points. Each point has 3 independent variables, x1, x2 and x3 and 1 dependent variable, y.

The structure of this list is {{x1,x2,x3,y} . . .}

data = {{1000, 200, 0.00831, 34.5091}, {1000, 200, 0.01622, 21.5683}, {1000, 200, 0.02493, 12.8382}, {1000, 100, 0.00831, 5.1462}, {1000, 100, 0.01622, 4.3582}, {1000, 100, 0.02493, 3.4017}, {1000, 50, 0.00831, 0.6760}, {1000, 50, 0.01622, 0.6452}, {1000, 50, 0.02493, 0.5953}, {1100, 200, 0.00831, 37.9601}, {1100, 200, 0.01622, 23.7251}, {1100, 200, 0.02493, 14.1220}, {1100, 100, 0.00831, 5.6609}, {1100, 100, 0.01622, 4.7940}, {1100, 100, 0.02493, 3.7419}, {1100, 50, 0.00831, 0.7436}, {1100, 50, 0.01622, 0.7097}, {1100, 50, 0.02493, 0.6548}, {1200, 200, 0.00831, 41.4110}, {1200, 200, 0.01622, 25.8820}, {1200, 200, 0.02493, 15.4059}, {1200, 100, 0.00831, 6.1755}, {1200, 100, 0.01622, 5.2299}, {1200, 100, 0.02493, 4.0821}, {1200, 50, 0.00831, 0.8112}, {1200, 50, 0.01622, 0.7743}, {1200, 50, 0.02493, 0.7144}};

For datapoint #1 1000, 200 and 0.00831 combine algebraically in some unknown way to produce 34.5091 as an output.

y = f(x1,x2,x3)

or

34.5091 = f(1000,200,0.00831) etc.

The task is to find the simplest algebraic function f that explains the data. I'm not interested in an interpolating function, neural network or other non-algebraic result; an algebraic formula is desired.

Here is my attempt at this using Mathematica:

FindFormula[data, {x1, x2, x3}, 5, All]

And here is the Mathematica output:

FindFormula::wrgfmt: Argument {{1000.,200.,0.00831,34.5091},{1000.,200.,0.01622,21.5683},{1000.,200.,0.02493,12.8382},{1000.,100.,0.00831,5.1462},{1000.,100.,0.01622,4.3582},{1000.,100.,0.02493,3.4017},{1000.,50.,0.00831,0.676},{1000.,50.,0.01622,0.6452},<<12>>,{1200.,200.,0.02493,15.4059},{1200.,100.,0.00831,6.1755},{1200.,100.,0.01622,5.2299},{1200.,100.,0.02493,4.0821},{1200.,50.,0.00831,0.8112},{1200.,50.,0.01622,0.7743},{1200.,50.,0.02493,0.7144}} at position 1 does not have the right format. Data should be a numerical array of depth less or equal than 2.

From the help menus it appears that Mathematica can deal with just 1 independent variable.

Can this problem involving 3 independent variables be solved using Mathematica ?

If not, are there future plans for this capability ?

@ Mark Kotanchek Thank you for your efforts on my question but your answer is not turning out to be useful for us.

First of all, none of the 10 functions you supplied is equivalent to the theoretical function used to generate the 27 data points I supplied in my question.

Secondly, if I evaluate your first function at each of the 27 given points and compute the errors, I see that the maximum error your function produces is in excess of 172%, error = ((incorrect - correct)/correct)*100. This error magnitude (at a data point known to your software) far exceeds what would be considered acceptable and brings into question expected error magnitudes at unknown data points where such a function would need to be used and trusted. This error takes place at datapoint #7 {1000, 50, 0.00831} where the correct value is 0.67601 and your model prediction is 1.84302

Thirdly, any function found that truly describes the theoretical underpinnings of the data should be able to accurately extrapolate beyond the confines of the data set used to find the function. For example, the gravitational force model found in the referenced YouTube video should be able to return reasonable values when asked to extrapolate. When I use your first model to extrapolate a reasonable amount (no negative or zero valued inputs for example) errors exceeding 1000% can occur. I’d be happy to give you this 28th datapoint if you think it would be useful. But keep in mind, in my real world of engineering acquiring 27 accurate datapoints, equally spaced in a matrix is a luxury. Real-world scenarios (for us) are much messier, data starved, incomplete and irregular.

Could it be that your software requires more than 30 seconds of computation time to return a robust model ? Or do you think this is beyond the current state of the art and I should re-ask in another 7 years ?

Here is the error plot for your model #1. I see a kind of symmetry here, perhaps that holds clues for model improvement ?

Mark, here is the error plot for your 2nd model. Although it's greatly improved, the errors are still too large at the known points. I'm not aware of too many industries that would consider a 24% error acceptable. What your models are in competition with is the Mathematica function Interpolation that returns a data model with maximum error less than 10^-13, essentially zero. I realize that Interpolation does not return algebra and is not doing symbolic regression but error-prone algebra is of little value to those working on life-critical systems. Why don't you refrain from restricting yourself and utilize all available functions in your software and let it run not for 10 minutes but 10 hours (or more) so we can really see what it can do on this problem ? If you think 3 additional data points contained within the input space would help, tell me where in the volume you would like them. This would give your software a total of 30 to work with.

• It would be useful to link to your original questions. Jul 18 '21 at 19:47
• There is now FindFormula and FindDistribution. Jul 18 '21 at 19:56
• I have initiated closing this question due missing (1) links to original question (2) Issues regarding any answers to your original question (3) a specific example. Please edit your question. Jul 18 '21 at 20:10
• @Steve That would be helpful to include that in your question.
– JimB
Jul 18 '21 at 20:21
• A concrete example would be really helpful here. Mathematica input and some indication of desired output (if not obvious from the input). Jul 19 '21 at 5:09

There is a commercial solution from www.evolved-analytics.com which also includes all of the surrounding functionality to actually extract the value from the data.

The next version will be released in the next few days. Free trials are available as are academic versions.

DataModeler is 100% Mathematica and one of the versions is an add-on package so all of the 500+ functions and symbols are fully documented and available for custom workflows. For those not interested in such, the GUI is also 100% Mathematica (with a little linux under the hood to allow easy export to word or powerpoint documentation).

CAVEAT: I am the CTO and original developer so I am biased.

If you want to see how Symbolic Regression can be used from within Mathematica (Wolfram Language), you can check out my video here: https://www.youtube.com/watch?v=eMcBMvy_tlk. It's an EXTREMELY powerful product.

• I watched your video, and the scientist in me questions if the modeling success you show is a true generalized capability or is due to the relative simplicity of the function you chose to find. Newton's law of gravitation contains only 1 term and no exponents greater than 2; a function you would agree that is not that challenging in today's engineering world. It may be hard to find a needle in a haystack, but not so a battle tank. Do you care to take a crack at my 27 data point set ? Mark was unable to produce a useable model after at least 2 attempts using his tool. Jul 29 '21 at 18:04

I have been a user of DataModeler from Evolved Analytics for more than 10 years that is mentioned just above. It will be one of products that match your need. It has a function called SymbolicRegression that creates approximate functions with genetic programming for optimization. I have been using it for the analysis in the medical field.

I compared its ability of fitting to several datasets with other machine learning models, including Lasso, Ridge, Logistic, SVM, Random Forest, XGBoost and so on. Probably due to the high flexibility of symbolic regression, it achieved as much accuracy/AUC as or even better than these of other models.

• And even more trolling? Jul 19 '21 at 1:06
• @jagra Truth be told, Stephen had FindFormula implemented after playing with DataModeler. The package version is a commercial add-on for Mathematica targeted for serious analysis. The GUI hides the complexity of Mathematica and is generally preferred by the user community. Jul 19 '21 at 2:25
• @jagra Why are you accusing people of trolling? Did you not check the link provided to assess relevance to Mathematica? Jul 19 '21 at 4:48
• @berniethejet No. Not yet. I had submitted paper including those comparisons three times and they were not accepted mainly due to the small number of samples. Hopefully it will be accepted next time since I have increased the sample size. Jul 21 '21 at 23:34
• @berniethejet You, too!. Jul 24 '21 at 1:10

Just for fun, I devoted 30 seconds to modeling the data set ...

Each dot represents a model and the ones in red define the ParetoFront exploring the trade-off of model complexity vs accuracy. Somewhat arbitrarily, I decided to focus on those models having a complexity less than 100 and an R2 better than 0.98. The models along the ParetoFront in this region are:

From the contenders, I chose a diverse model set from which to form an ensemble,

The attraction of a model ensemble is that it defines a trustable model since the models will agree where constrained by data but will diverge (since we chose them for their diversity) when asked to venture into unknown regions of parameter space or if the system has undergone some sort of fundamental change.

For real-world analysis, understanding the relative importance of variables, identification of metavariables as well as the creation of trustable models is very important.

Real-world data here is defined as up to tens of thousands of variables and millions of data records. Of course, things get easier with smaller data sets. For the big data sets, more than 30 seconds of model search will be required.

Per JimB's comment, here is a graphic using the standard Mathematica machine learning techniques relative to a SymbolicRegression ModelEnsemble. The blue line represents the model prediction and the envelope the trust for a simple ball trajectory problem. The left side shows the prediction if we stop the available data at the halfway point and the right if we only have the extremal values and are looking to interpolate. (The red points are known whereas the yellow are the unknown reality.) The gray lines in the case of SymbolicRegression are the trajectories of the constituent models which comprise the ensemble.

• This is great for obtaining a compact prediction model but you've gone way overboard to claim that this provides a "trustable" model especially when if there's no option in your software to verify such a trust and even more so with the statement "things get easier with smaller data sets". That is wishful thinking.
– JimB
Jul 26 '21 at 19:28
• I understand and appreciate the great amount and quality of work. That doesn't mean one can skip the model verification process for any particular model and set of data. My concern is that folks that use your software might be led to believe that running your software is all that's necessary. Is there something in the software or your webpage that discusses model verification? Language is important. Using the term "describes" rather than "explains" is much more accurate.
– JimB
Jul 26 '21 at 20:32
• @JimB Obviously, I've left out a couple decades of development and industrial application in the above answer; however, you appear to have missed that we are choosing diverse models for the ensemble. Model diversity is achieved by choosing those models with the least correlated prediction residuals. The model search does get easier with fewer candidate variables (the haystack to be searched gets smaller) as well as with fewer data records (less time required for model assessment). This approach has successfully been deployed in chemical process industries for over 20 years. Jul 26 '21 at 20:36
• I have no doubt this approach has been successful for predictions much of the time but "finding the underlying model" certainly has happened much less frequently. Do you have a definition or reference for "Least correlated prediction residuals"? That term and "trustable models" seems to have much more of a marketing flavor. (Not to mention "diversity of models.")
– JimB
Jul 26 '21 at 20:44
• @JimB I am with you that all data and models are guilty until proven otherwise. One of the attractions of the model ensemble is that it makes obvious that there are lots of models which fall into the good-enough category. Even so, we always have the struggle of trying to explain that THE model doesn't exist and picking one and running with that is fraught with peril. As near as I can tell, the world isn't a 2nd order polynomial with no cross-terms. That said, our experience is that SymbolicRegression models do tend to degrade more gracefully than others when asked to extrapolate. Jul 26 '21 at 20:48