Combined method of RF and Symbolic regression models

[kyoungyi]

Hi Dr. Miles Cranmer,

Thank you for this wonderful tool. I am currently trying to get a mathematical operation (as an index) to estimate the impact of weather stress on agricultural yields by using your package "pysr", based on 40 years of county-monthly level data.

I have 9 features that were narrowed down from the Random Forest model (Wadeker et al., 2020, 2021) and I would like to get one complete equation to capture partial dependence (accumulated local effect) plots of RF for the final relationship.

Questions are as follows:

1. Selected features in a final equation (This is very important for me to understand!):

1.1) Compared to the important features (with interaction effects) from the RF model, the most important features do not show up in the final pysr equations (of multiple runs with different random states). I am confused by this inconsistency. Do you have any input on this?

1.2) In addition, the final equation often chooses a couple of features with two or three times of repetitions. Is there any way to force which features to be included in a final equation?

1.3) I assume that the final features absorb (or minor improvements to be selected) the effects of the remaining features that did not show up in the final equation? Or do I need to adjust the default parameters such as the ones for mutations or migration parts?

1.4) Due to the inconsistencies of selected features, I am not sure if a combination of RF+symbolic model could be good, in addition to the fact that pysr does not support plotting. Any inputs?

2. Parameters:

2.1) With changes in parameters such as populations and denoise, the running speed was incredibly slow (even on HPC) and it does not seem to stop with the timeout function as well.

2.2) FYI, I’ve constrained some parameters in a way like Lemos et al., 2022 and Wadeker et al., 2020 but left the others to default values. I have been testing with different parameters based on "PySRRegressor reference" page but with my limited understanding, it has not been easy. What are your strategies to best optimize parameters with different datasets?

• Training data structure: (8200, 9)
• Loss: Mean squared error (default)
• Optimizer iterations: 10
• Maxsize: 40
• Operations: [“plus”, “mult”, “sub”, “pow”, “div”, “square”, “cube”, “tan”, cos”, inv”]

Thank you again for your input and time in advance.
Also. if anyone has experience or recommendations, I will very much appreciate your comments!

[MilesCranmer]

1.1) Compared to the important features (with interaction effects) from the RF model, the most important features do not show up in the final pysr equations (of multiple runs with different random states). I am confused by this inconsistency. Do you have any input on this?

If they don't show up over multiple different runs, maybe those features aren't actually important? Computing feature importance with RF isn't a guarantee. Note that with the latest versions of PySR after 2022 (after Jay's paper was done) selecting important features isn't really important, because the crossover operation does this implicitly. So I would just give all the features you have and PySR should be okay to handle it.

1.2) In addition, the final equation often chooses a couple of features with two or three times of repetitions. Is there any way to force which features to be included in a final equation?

Not sure I follow. If a feature doesn't improve performance, it won't be included. PySR tries to find simple and accurate expressions; the simple aspect means fewer features.

You could write a custom loss function that requires features to be used (see #161 for related question) but not sure why you would want that.

1.3) I assume that the final features absorb (or minor improvements to be selected) the effects of the remaining features that did not show up in the final equation? Or do I need to adjust the default parameters such as the ones for mutations or migration parts?

Not sure I follow.

1.4) Due to the inconsistencies of selected features, I am not sure if a combination of RF+symbolic model could be good

I would just use PySR by itself and run it for long enough.

, in addition to the fact that pysr does not support plotting. Any inputs?

You can plot PySR outputs with matplotlib, like: plt.scatter(model.predict(X, i), y), where i is the row of the equation you wish to evaluate. See https://astroautomata.com/PySR/examples/#4-plotting-an-expression for an example.

2.1) With changes in parameters such as populations and denoise, the running speed was incredibly slow (even on HPC) and it does not seem to stop with the timeout function as well.

The denoise parameter uses a Gaussian process. GPs can be extremely slow for many points, as (I think) they scale with n^3. So it probably just didn't get to the search yet. Maybe use fewer points or just avoid using denoising?

2.2) FYI, I’ve constrained some parameters in a way like Lemos et al., 2022 and Wadeker et al., 2020 but left the others to default values. I have been testing with different parameters based on "PySRRegressor reference" page but with my limited understanding, it has not been easy. What are your strategies to best optimize parameters with different datasets?

My general tips would be to avoid using redundant operators, like how pow and square and cube are equivalent. The fewer operators the better. I usually do the following:

I run from IPython on the head node of a slurm cluster. Passing cluster_manager="slurm" will make PySR set up a run over the entire allocation. I set procs equal to the total number of cores over my entire allocation.

1. Use the default parameters.
2. Use only the operators I think it needs and no more.
3. Set niterations to some very large value, so it just runs for a week until my job finishes. If the equation looks good, I quit the job early.
4. Increase populations to 3*num_cores.
5. Set ncyclesperiteration to maybe 5000 or so, until the head node occupation is under 10%.
6. Set constraints and nested_constraints as strict as possible. These can help quite a bit with exploration.
7. Set maxdepth as strict as possible.
8. Set maxsize a bit larger than the final size you want. e.g., if you want a final equation of size 30, you might set this to 35, so that it has a bit of room to explore.
9. Set parsimony equal to about the minimum loss you would expect, divided by 5-10. e.g., if you expect the final equation to have a loss of 0.001, you might set parsimony=0.0001.
10. Set weight_optimize to some larger value, maybe 0.001. This is very important if ncyclesperiteration is large, so that optimization happens more frequently.
11. Set turbo to True. This may or not work, if there's an error just turn it off (some operators are not SIMD-capable). If it does work, it should give you a nice 20% speedup.

Since I am running in IPython, I can just hit "q" to stop the job, tweak the hyperparameters, and then start the search again.

Some things I try out to see if they help:

1. Play around with complexity_of_operators. Set operators you dislike (e.g., pow) to have a larger complexity.
2. Try setting adaptive_parsimony_scaling a bit larger, maybe up to 1000.
3. Sometimes I try using warmup_maxsize_by. This is useful if you find that the search finds a very complex equation very quickly, and then gets stuck. It basically forces it to start at the simpler equations and build up complexity slowly.
4. Play around with different losses:
   i. I typically try L2DistLoss() and L1DistLoss(). L1 loss is more robust to outliers compared to L2, so is often a good choice for a noisy dataset.
   ii. I might also provide the weights parameter to fit if there is some reasonable choice of weighting. For example, maybe I know the signal-to-noise of a particular row of y - I would set that SNR equal to the weights. Or, perhaps I do some sort of importance sampling, and weight the rows by importance.

Very rarely I might also try tuning the mutation weights, the crossover probability, or the optimization parameters. I never use denoise or select_k_features as I find they aren't very useful.

For large datasets I usually just randomly sample ~1000 points or so. In case all the points matter, I might use batching=True.

If I find the equations get very complex and I'm not sure if they are numerically precise, I might set precision=64.

Hope this helps!

[MilesCranmer]

Turned this into a docs page! https://astroautomata.com/PySR/tuning/

[kyoungyi]

This is amazing! Thank you so much for your swift and detailed answers. You covered most of my questions. I will follow your tips and try it again.

[kyoungyi]

Hi Miles, thank you again for your swift reply and thorough answers before.
I am still stuck with some problems and having a tough time figuring them out. Your help would be very appreciated again.

1. Poor predictive skills:
   Despite using the default or recommended parameters you provided earlier, I still get poor R^2 values (e.g., around 0.1) and low Pearson correlations between the predicted and true response values. I expected the R^2 to be similar or higher than the random forest model's R^2 (R^2 ~ 0.7) for the same training dataset. I am not sure what else I can adjust or how long I need to run the model (Time-Delayed Variables #275 ). You mentioned running it for a week but I noticed minor changes in the equations after a certain point (maybe after 70-80% progress). Do you have any suggestions for further testing?

2. Features selected by pysr:
   I let all my features (~15 features) into the model as you suggested, but only a couple of variables that I thought were less important (based on science/references/RF model) were repeatedly showing up in the final equations. For example, x1x1 +x1/x2 +sin(x1) + 0.02212x2 .... (while the others were not showing up much). Any thoughts on this? Any ways to diagnose errors here?

3. Formula prior:
   As a possible solution for question 2, I fixed a symbolic form of an expression with input variables ([Feature] add formula prior #285 ). For example, I wanted to have a form like f1(x1) + f2(x2,x3,x4,x5,x6,x7,x8,x9..) = y. However, it didn't follow the fixed form and f2 selects only one variable with many constants. I could have penalized the function to include all variables in f2 but I think the symbolic model should be able to get the important features and I am not familiar with julia language... Anyway, I just changed (2,3) and (1,4) to the corresponding index number of my variables in "return !t.constant && t.feature in (2, 3) and (1,4)". Would it be wrong? If not, what are the possible reasons for not following the fixed form? Also, is there a way to control the complexity of f1 and f2 separately?

Any further input or recommendations would be very appreciated!

[zhuyi-bjut]

@kyoungyi Hi! About Formula prior, I now have the same problem as you : pysr tends to build a formula with fewer variables and constants, thus ignoring some variables with physical meaning ( but that doesn 't mean it 's not important ). I wonder if you already have a solution ? Or related suggestions ? Thanks！

[kyoungyi]

Hi @zhuyi-bjut , I think this feature is still being developed. But I have tried to focus on a more refined dataset to reduce noise. Also, if you are using/have used another model, you can try to fit the symbolic regression with the predicted values of another model.