Back to the Formula -- LHC Edition

This paper uses symbolic regression to identify ML-informed observables for new physics searches. I think the paper is detailed, well-written and deserves publication after a few small questions are addressed. I do commend the authors on their choice of network naming.

I do not quite see how this approach will stand its ground when it comes to scalability and portability. The processes that the authors choose to investigate in this proof-of-principle analysis are well-motivated, and I do understand that something comprehensive seems to be in the making. Nonetheless, these processes are so simple at leading order in terms of involved mass scales and kinematics that it is hard to see the approach's relevance going forward. While it is good to see that known results are reproduced, a lot of it seems to rely on an educated choice of input parameters? Also, how much is this approach actually limited by ab initio choices? Given that the information extracted by the network is essentially the correlation of the matrix element, what is the selling point of this approach compared to e.g. the matrix element method?

Higher order corrections will introduce non-polynomial uncertainties and dependencies, is this a bottleneck for this approach?

In this study, the authors present a Machine Learning based symbolic regression method in order to extract analytic functions from LHC event shapes. The study is interesting and worthy of publication, but I would like to raise a few issues regarding the text.

The language of the study is highly casual and lacks various references or explanations of the concepts out of the knowledge of a general particle physics audience.

- For instance, just before equation 2 "Taylor expansion" has been used as a verb, i.e. "taylor the log likelihood".

- Neyman-Pearson lemma and Cramer-Rao bound requires citation. Ref [2] refers to a physics analysis paper rather than the original reference. Additionally, these methods are quite unfamiliar to the article's audience; the paper can benefit from a couple of lines of explanation.

- Which Delphes detector card is used has never been stated in the study.

- Parsimony has not been defined. It looks like a scaling constant, but it's unclear.

- Simulated annealing method has not been adequately cited.

- There are multiple references to a "score" measure with various score definitions in the paper. It would be clearer if authors could refer to function or method name or reference to an equation instead of repeatedly using the word "score".

- EWDim6 model has not been cited.

- Page 14 first paragraph reads, "...Eq 12 or Eq 28 do not accept enough more complex functions..." I believe there is a typo in enough/more where only one of them should be there.

- Delphes module internally uses FastJet. Although the antikt algorithm has been cited, FastJet has never been mentioned in the paper.

I also have one additional question: representation is one of the biggest problems of ML methods where even with the best optimization algorithms and statistically valid data, an NN may not be able to represent the given model. Since in regular classification or regression type of analyses, the change in the results can be measured by training the network multiple times and taking the mean of the results, which allows a better approximation of the final result. However, since this study aims to describe the data with an analytic function, for each training, the functional form should change. A similar averaging or ensembling method would be sufficient for the cases that this difference can be absorbed into the constants. But if each execution leads to a slightly different functional form, I can not imagine how to write the analytic function to be more representative of the model. I was wondering if the authors have faced such differences in the results; if they did, do they have a method to combine these results somehow?