PySR: Symbolic Regression That Returns an Equation Instead of Weights

The gap PySR fills between a fitted model and a formula

A gradient boosted model or a neural network will give you predictions. It will not give you a formula you can write on a whiteboard, differentiate by hand, or check against a physical law. PySR targets that second outcome. The README defines symbolic regression as a machine learning task where the goal is to find an interpretable symbolic expression that optimizes some objective. The word interpretable is doing real work there: the output is an expression tree built from operators you supply, not a set of coefficients attached to fixed basis functions. The intended user is someone who already suspects a relationship exists in the data and wants the search to propose its algebraic form. The README also points at symbolic distillation of neural networks, citing 2006.11287, where a trained network is converted into an analytic equation. That is a second audience: people who have a working black box and want a readable approximation of it. The README is explicit that symbolic regression works best on low-dimensional datasets, which is a scoping statement rather than a marketing line.

A Python front end over a Julia search engine

PySR is not a single-language library. The README states that PySR is developed alongside the Julia library SymbolicRegression.jl, which forms the search engine, and that fitting launches a Julia process which does a multithreaded search for equations. The Python class is therefore a controller: it holds configuration, serializes data and operators across the language boundary, and receives back a population of candidate expressions. The algorithm details live in the PySR paper, arXiv:2305.01582, not in the README, so anyone who needs to reason about selection pressure or mutation rates has to go to the paper. What the README does show is the shape of the search: with niterations=40 the documentation describes the run as containing hundreds of thousands of mutations and equation evaluations. The result is not one equation but a Pareto front, printed as a table with pick, score, equation, loss and complexity columns. You choose a row, and model.predict(X, 3) returns predictions from the third equation in that table. That table, not a single fitted object, is the artifact you take away.

Installation pulls in Julia on first import

The install line is short: pip install pysr, or conda install -c conda-forge pysr. The README notes that Julia dependencies will be installed at first import, which means the expensive and failure-prone part of setup happens later, inside a Python process, not at the moment you run pip. The README documents one concrete failure: a hard crash at import with a message like GLIBCXX_... not found, caused by another Python dependency loading an incorrect libstdc++ library. The stated fix is to prepend the Julia libstdc++ directory to LD_LIBRARY_PATH, with the example path $HOME/.julia/juliaup/julia-1.10.0+0.x64.linux.gnu/lib/julia/ and the caveat that this likely differs on your system. For environments where you cannot modify the host, the repository ships a Dockerfile and an Apptainer.def. The Apptainer route is specifically for clusters without root access, built with apptainer build --notest pysr.sif Apptainer.def and launched with apptainer run pysr.sif. If your cluster already has a module system that pins libstdc++, the Docker or Apptainer path is the one to take rather than patching shell profiles.

Configuring the search: operators, loss, and size limits

The quickstart constructs a PySRRegressor with four settings that matter. maxsize=20 caps expression size, which bounds the search. niterations=40 sets how long the search runs, and the inline comment in the README says to increase it for better results, so 40 is a demonstration value, not a recommended default. binary_operators=["+", "*"] and unary_operators=["cos", "exp", "sin", "inv(x) = 1/x"] define the grammar the search can build from. The custom operator is written in Julia syntax, and it needs a matching Python definition passed through extra_sympy_mappings={"inv": lambda x: 1 / x}, otherwise SymPy cannot manipulate the resulting expression. The loss is also Julia syntax: elementwise_loss="loss(prediction, target) = (prediction - target)^2". That split is the sharpest practical constraint in the whole API. Every operator and every loss you add has to be expressed twice, once for the Julia search and once for the Python-side symbolic layer, and the two must agree. Adding a new operator is therefore a two-file change, not a one-line change.

Where PySR is the wrong tool

The README states plainly that symbolic regression works best on low-dimensional datasets. The demo uses five features, and the printed equations mostly involve x0 and x3, which is what you would expect when the true relation only touches two of them. Scale the feature count and the operator combinations grow, so the same niterations budget explores a shrinking fraction of the space. If your problem is genuinely high-dimensional and you want a predictive model, a gradient boosted tree will be cheaper and more accurate. There is a second limitation that is easy to miss: the output is a Pareto front, and picking a row is a judgement call. The score column is a heuristic for that choice, not an error bar. Nothing in the README establishes that the selected equation is the true generating function, and an expression that fits 100 points well can still be the wrong functional form. A third constraint is determinism. The README describes an evolutionary search with hundreds of thousands of mutations, and it does not document a seed parameter or reproducibility guarantee, so two runs on the same data may surface different equations. If you need an auditable, repeatable fit, that matters.

How PySR differs from sparse regression libraries such as PySINDy

The nearest alternative in spirit is sparse identification of nonlinear dynamics, commonly implemented as PySINDy. The two solve different problems. PySINDy takes a fixed library of candidate terms that you write down in advance, then fits sparse coefficients over that library. The functional form is your input. PySR searches over the space of functional forms itself, combining operators through mutation and selection, so the form is the output. That difference has consequences. PySINDy gives you a convex or near-convex fitting problem with a coefficient vector you can inspect, which makes it fast and reproducible. PySR gives you a combinatorial search that can propose a term you would never have put in the library, at the cost of runtime and run-to-run variation. If you already know the candidate terms from domain theory, PySINDy is the more direct route. If you do not, PySR is doing work that a fixed library cannot do. The choice is about whether the uncertainty is in the coefficients or in the form.

Maintenance, upgrades, and the licence

PySR is not archived, the default branch is master, and the release cadence is fast: v2.2.0 on 2026-09-02, v2.2.1 the same day, and v2.3.0 on 2026-09-07. Three releases in under a week is a sign of active development, and also a sign that pinning a version is wise for anything you need to reproduce. The coupling to SymbolicRegression.jl means a PySR upgrade can move the Julia side with it, and the first import after an upgrade may reinstall Julia dependencies. Plan for that in CI: cache the Julia depot directory, or use the Docker or Apptainer image so the toolchain is fixed at build time. The licence is Apache-2.0, a permissive licence that includes an explicit patent grant and requires you to retain notices and state significant changes. That is a summary, not legal advice; check the LICENSE file for the terms that bind you. If you publish results, the README asks you to cite arXiv:2305.01582.

Who should adopt PySR, and what to check first

Adopt PySR if your dataset is low-dimensional, you can name a plausible operator set, and the deliverable is a formula rather than a score. The workflow is short: fit, read the equations_ table, pick a row by complexity and loss, then call model.predict with that index. Do not adopt it as a drop-in replacement for a general regressor, and do not adopt it if you cannot state your operators up front, because the grammar you pass in bounds what the search can ever return. Before you commit, run the first import on the machine that will actually execute the fit, since that is where the Julia install and the GLIBCXX failure surface. Then check that the expression you select uses the features you expect; in the README's own demo the fitted terms concentrate on x0 and x3, which is the behaviour you want to confirm rather than assume. Treat the complexity column as a budget line, since maxsize caps it and a larger cap costs search time.

Editorial conclusion

Adopt PySR when you have a low-dimensional dataset and need a closed-form expression you can read, differentiate or hand to a domain expert, and when you can accept a Julia toolchain appearing at first import. Do not adopt it if you need a fitted model whose only job is prediction, or if your feature count is large enough that the combinatorial search space stops being tractable. Before committing, verify on your own machine that the first import succeeds and that the equations_ table contains an expression whose complexity you are willing to defend, because that table is the actual deliverable.

Official sources

1. astroautomata/PySR on GitHub
2. License: Apache-2.0
3. Project website
4. README
5. Releases