Discovering Equations with Genetic Programming
ML for Science and Engineering · Lecture 12
In SINDy, we chose a fixed library of candidate terms:
This works when we have a good guess about what terms might appear.
But what if the true equation is:
No finite polynomial library will contain this term. We need to search over functional forms themselves.
Given data $(x, y)$, find a mathematical expression $f$ such that $y \approx f(x)$.
Every mathematical expression has a natural tree structure:
Internal nodes = operators ($+, \sin$)
Leaves = variables ($x, y$) or constants
Click an expression to see its tree representation:
Count the expression trees you could write with $k$ binary operators, $u$ unary ones and $v$ leaves. Move the sliders.
Each extra node multiplies the count, so a computer a thousand times faster moves the wall by about two nodes.
Inspired by biological evolution: a population of candidate solutions evolves over generations.
Tournament selection: pick $k$ individuals at random, keep the best one.
Other strategies: roulette wheel (probability proportional to fitness), rank selection, elitism (always keep the best).
Pick a random subtree in each parent, swap them to create offspring:
Selection and crossover only recombine what the population already has. Mutation is where new material comes from. Press a button.
population = [random_tree() for _ in range(pop_size)]
for generation in range(max_generations):
# Evaluate fitness
for tree in population:
tree.fitness = mse(tree, data) + alpha * tree.size
# Selection + reproduction
offspring = []
while len(offspring) < pop_size:
p1 = tournament_select(population, k=5)
p2 = tournament_select(population, k=5)
c1, c2 = crossover(p1, p2)
c1 = maybe_mutate(c1, rate=0.1)
c2 = maybe_mutate(c2, rate=0.1)
offspring.extend([c1, c2])
# Elitism: keep best from previous generation
population = elites(population, n=5) + offspring
Target: discover $f(x)$ from noisy data. Watch the population evolve.
2,902 planets outside the solar system, every one of them measured: how long its year is, how far it orbits, and how heavy its star is.
NASA Exoplanet Archive, confirmed planets with all three quantities measured. The CSV is on this site. Discovered 1992 to 2026, mostly by transit and radial velocity.
Year against orbital distance. Every dot is a real planet around a real star.
Take the logarithm of both axes and the cloud becomes a line. A line in log space is a power law: $P = C\,a^{\alpha}$.
Least squares on 2,902 planets gives $\alpha = 1.457$, with $R^2 = 0.9896$.
Kepler's third law says the exponent is exactly $3/2$. We got 1.457. Close, and wrong. A search that stops here reports a law that does not exist.
The scatter around the line is not noise in the measurements. It is a variable we left out.
Residual scatter, before and after the star's mass enters the expression.
Fit again on the same 2,902 planets: $\alpha = 1.497$, $\beta = -0.479$. Kepler to three digits, and the scatter falls by 68 per cent.
Median error across every planet: 0.19 per cent. Newton could have written it, and did.
Fitting a curve is not finding a law. The search has to range over which variables matter, not only over exponents.
Every point on this surface is one candidate expression $P = C\,a^{\alpha} M^{\beta}$. Brighter is a better fit to all 2,902 planets.
The search does not know any physics. It tries exponents, measures the error, and walks downhill.
Press Search.
It lands on $(3/2, -1/2)$ because that is where the data is, not because anyone told it about gravity.
Kepler had Tycho Brahe's tables and took roughly a decade. This is the same inference on more planets, in the time it takes to let go of the mouse.
Every candidate is a dot: how big it is against how wrong it is, on the 2,902 planets. Click one.
The front is the set of candidates that nothing smaller beats. Everything above it is dominated: another expression is both simpler and more accurate.
The interesting place is the knee, where the error falls off a cliff and then stops improving. Here the knee is at size 8, and the expression sitting on it is Kepler's third law.
Notice the size 13 candidate: it is more complicated than the law and worse. Complexity that buys nothing is how a search fools you.
from pysr import PySRRegressor
import numpy as np
X = np.random.randn(100, 2) # two input variables
y = X[:, 0] * np.sin(X[:, 1]) # true relationship
model = PySRRegressor(
niterations=100,
binary_operators=["+", "-", "*", "/"],
unary_operators=["sin", "cos", "exp"],
populations=20, # parallel populations
maxsize=25, # max expression complexity
)
model.fit(X, y)
print(model) # Pareto front of equations
PySR outputs a Pareto front of equations ranked by complexity vs. loss. Uses Julia backend for speed.
| Method | Search space | Interpretable? | Extrapolates? |
|---|---|---|---|
| Neural network | Weights of fixed architecture | No | Poorly |
| SINDy | Coefficients of fixed library | Yes | If correct |
| Symbolic regression | All expressions | Yes | If correct |
| Physics-informed NN | NN + known PDE constraint | Partially | Somewhat |