milescranmer/pysr/pysr

Using PySR Effectively

PySR evolves expression trees through the Julia backend SymbolicRegression.jl, returning readable equations on an accuracy/complexity Pareto front. This guide distills the documentation and several hundred real user threads. Full docs: https://pysr.ai/

Quick start

# pip install pysr  (Julia is downloaded automatically on first import; no separate install)
import numpy as np
from pysr import PySRRegressor
X = 2 * np.random.randn(100, 5)
y = 2 * np.cos(X[:, 3]) + X[:, 0] ** 2 - 2
model = PySRRegressor(
    binary_operators=["+", "-", "*", "/"],
    unary_operators=["cos"],
    niterations=100,
)
model.fit(X, y)
print(model)            # Pareto front: complexity, loss, equation
model.predict(X)        # uses the auto-selected "best" equation
model.sympy()           # SymPy expression
model.latex()           # LaTeX string

model.equations_ is a pandas DataFrame containing complexity, loss, score, equation, sympy_format, lambda_format. Select any front row by index in predict, sympy, latex, jax, pytorch; without an index, methods use the auto-selected equation. Indexing lets you inspect and export a specific Pareto candidate instead of only the default row.

Keep the process alive

The first .fit() pays Julia startup/JIT, roughly 1-2 minutes; subsequent same-process fits start almost instantly. Repeated fresh processes repay startup and compilation for every experiment, a common agent failure.

• Keep a long-lived IPython/Jupyter kernel or REPL. New PySRRegressor instances are cheap; Julia runtime state lives per-process.
• Scripts work. Iterate parameters inside one process instead of repeatedly relaunching python script.py. Process reuse is an optimization rather than a requirement.
• Reuse an environment containing pysr; fresh environments re-resolve and precompile Julia packages.
• Expect “Compiling Julia backend...” initially. Avoid treating normal compilation as a hang.

Candidate equations with guesses

Inject knowledge through modifiable/discardable guesses or mandatory template structure. Combine both using component/parameter guesses. Constructor strings require Julia syntax, input variable names, and enabled operators:

model = PySRRegressor(
    operators={2: ["+", "*"]},
    guesses=["x0 + 0.5 * x1", "x0 * (x1 + 1.0)"],
)

Single output: candidate list. Multiple outputs: lists in target-column order, guesses=[["x0 + x1"], ["x0 * x1"]]. Templates: dictionaries keyed by component/parameter names; #1, #2, etc. identify component arguments.

Guesses enter populations throughout search. should_optimize_constants=True optimizes constants/template parameters before insertion. fraction_replaced_guesses, default 0.001, controls the cycle-end replacement fraction.

Guesses are pushed into the hall of fame before generation one. Inspect accepted seeds at their complexity in model.equations_ or immediate hall_of_fame.csv; confirm insertion and optimized constants. Oversize seeds produce only Julia warning Guess expression '...' has complexity 23 > maxsize (10), using default x1, x2, ... names regardless of supplied names. They remain unrecorded but still migrate; descendants become recordable after mutations satisfy maxsize. Violations of constraints, nested_constraints, or maxdepth disappear silently. Lower-loss candidates replace seeds at the same complexity: only one expression survives per complexity.

Keep constraints compatible and set maxsize slightly above the largest seed; use the warning's complexity. Cost rises steeply, PySR warns above 40, and several-times-seed sizes can exhaust the budget without completing one generation despite appearing alive. Use warmup_maxsize_by against early bloat.

Set guesses on the estimator, never .fit(). Continue with model.set_params(guesses=[...], warm_start=True), replacing guesses while keeping data representation/search space fixed.

Recommended workflow

1. Subsample: usually 1,000-5,000 representative rows suffice even from millions. Many features, heavy noise, or rare regimes benefit from more. PySR 2.0 automatically batches larger datasets; batching=False forces full-data search evaluations. Batching still reevaluates hall-of-fame candidates fully. Fewer rows proportionally accelerate search.
2. Minimize operators: select domain-plausible operators. Redundant pow/square/cube or -/neg enlarge search. Polynomials: ["+", "-", "*"], omit / and ^.
3. Use defaults: populations, parsimony, mutation weights, and ncycles_per_iteration were tuned by large-scale search in 2024-2025. Avoid older forum/paper recipes.
4. Debug briefly: check operators, loss, and sensible equations in a few-minute run before one long search.
5. Budget the final run: evolution lacks conventional convergence and can jump families after hours. Set large niterations; control wall-clock with timeout_in_seconds, total evaluations with max_evals, or stopping with early_stop_condition="stop_if(loss, complexity) = loss < 1e-6 && complexity < 10" (Julia (loss, complexity) function). max_evals supports compute-matched comparisons; niterations * populations * population_size does not count evaluations. Monitor partial outputs/<run_id>/hall_of_fame.csv; stop gracefully with q then Enter in IPython.
6. Inspect the entire front: validate candidates and report tradeoffs.

Choosing an equation

model_selection="best" maximizes score within 1.5x minimum loss; "accuracy" minimizes loss. Score is negative log-loss slope per complexity unit; a large score marks a sharp accuracy gain for little added complexity. These heuristics control default predict.

Print model.equations_; inspect sharp loss drops, held-out model.predict(X_test, index=i) when overfitting is plausible, and limits x -> 0/inf against domain expectations. Present 2-3 candidates with tradeoffs.

Losses and weights

Default: MSE.

• MAE (mean absolute error): elementwise_loss="L1DistLoss()" works with or without weights. For unweighted data, the equivalent custom loss is elementwise_loss="my_loss(x, y) = abs(x - y)". Useful for outliers, heavy-tailed noise, or costs proportional to absolute error; targets the conditional median rather than the mean.
• Uncertainty: model.fit(X, y, weights=1/sigma**2); built-ins apply weights automatically. Custom: elementwise_loss="myloss(x, y, w) = w * abs(x - y)^2".
• Wide target range: largest values dominate MSE; use log-space loss:

elementwise_loss = """function loss_fnc(prediction, target)
    scatter_loss = abs(log((abs(prediction)+1e-20) / (abs(target)+1e-20)))
    sign_loss = 10 * (sign(prediction) - sign(target))^2
    return scatter_loss + sign_loss
end"""

• Relative/percentage error: divide by target; prediction denominators reward infinite predictions. Guard near-zero targets.
• LossFunctions.jl strings work: "HuberLoss(1.0)", "LPDistLoss{3}()".
• Binary classification: +1/-1 targets, margin loss "L2MarginLoss()"; apply sigmoid yourself for probabilities.
• Asymptotes/boundaries, including y -> 0 as x -> inf, exact boundary values, and known limits: add synthetic points with very large weights, usually preferable to custom objectives. Strict enforcement: estimate limits numerically inside loss_function, add graded penalties.
• Keep losses deterministic for caching and non-negative; loss_scale="linear" permits negatives, including log-likelihoods.

elementwise_loss receives scalar (prediction, target) or (prediction, target, weight); never sum/broadcast. Whole-vector/tree objectives require loss_function.

Scaling and features

• Normalization is optional: constants sample near N(0,1), mutate multiplicatively; extreme scales slow search. Prefer natural units, avoiding nuisance constants and hidden physical meaning. Rescale visibly struggling searches.
• Roughly 10 features or fewer need no special handling.
• Tens: increase maxsize and rows; search selects features reasonably well, automatic batching handles volume. Requiring all 30+ features needs maxsize well above 100.
• Above ~50: engineer domain aggregates or template decompositions. Fields/images/sequences/graphs need meaningful features rather than pixel columns, or appropriately biased neural networks followed by symbolic component distillation, arXiv:2006.11287. Fallback: gradient-boosting preselection, select_k_features=k.
• Omitted variables failed to justify complexity through accuracy; force inclusion through an absence-penalizing custom loss.