Evolution strategies and genetic programming: from the one-fifth rule to Kepler’s law

The (1+1)-ES and the one-fifth rule, populations and self-adaptation on the Rastrigin function, CMA-ES, programs as trees, symbolic regression that finds Kepler’s third law in planetary data, and program bloat.

14 minutes Matěj Žalmánek 2026-10-15

Series AI basics Part 28 of 52

Full series syllabus 52 of 83 topics published

What you will take away

How evolution strategies tune real numbers and how they set their own step size.
Why a population beats a single individual on rugged terrain.
How genetic programming grows whole formulas — and why they tend to bloat.

Useful to know first: Genetic algorithms: population, crossover, mutation and selection on a backpack for a hike

In the 1960s two students at the Technical University of Berlin, Ingo Rechenberg and Hans-Paul Schwefel, stood in a wind tunnel over a jointed plate. They wanted the shape with the least drag, but they had neither the equations nor the computers to work it out. So they rolled dice: they turned the joints at random, measured the drag, and if it was lower, they kept the new shape. That experiment gave rise to the evolution strategy — evolution that tunes numbers instead of ones and zeros. A quarter of a century later John Koza showed that evolution can grow whole programs as well.

This follows on from the article Genetic algorithms. There an individual was a string of bits; here it will be a vector of real numbers and a program tree.

1. Evolving numbers: the (1+1)-ES

Algorithm — (1+1)-ES

Input: a function to be minimized f : ℝⁿ → ℝ (n ≥ 1), a starting point x₀ ∈ ℝⁿ, a step size σ > 0 and a number of generations K ∈ ℕ. Output: a point x ∈ ℝⁿ and its value f(x). An individual is a point x ∈ ℝⁿ.

Set x ← x₀.
For k = 1, …, K repeat:
1. draw z ~ N(0, I) (a vector of n independent normally distributed numbers, definition below) and create the offspring y = x + σ · z;
2. if f(y) ≤ f(x), set x ← y; otherwise keep x.
After K generations stop and return x and f(x). (In practice one often stops earlier, as soon as f(x) drops below a chosen threshold.)

So the value f(x) never increases. In this basic form σ is fixed; how to change it is described by the one-fifth rule below.

In plain words

One parent, one child. The child is the parent with a little random noise in every number; if it is better, it takes the parent’s place. The step size σ says how big that noise is — whether it makes small adjustments or big jumps.

y = x + σ · z, z ~ N(0, I)

y: offspring — the new candidate we are going to evaluate
x: parent — the best point found so far; the offspring is created in its neighbourhood, because good points tend to lie close to other good points
σ: step size — a single number that sets how far from the parent we search; the noise z has variance 1, only σ gives it the scale of the problem
z ~ N(0, I): direction and length of the attempt — mean 0 means no direction is favoured in advance; the identity matrix I means all coordinates change independently and equally strongly

Mutation in an evolution strategy: the parent plus a random step whose length is controlled by σ.

Definition — the normal distribution N(0, 1) and N(0, I)

Let Z be a real random variable. Z has the standard normal distribution N(0, 1) if and only if it has the density φ : ℝ → (0, ∞),

φ(z) = 1/√2π · exp(−z²/2), i.e. P(a ≤ Z ≤ b) = ∫_a^b φ(z) dz for all a ≤ b.

Then EZ = 0 and var Z = 1. About 68% of values lie in the interval [−1, 1] and 95% in [−1.96, 1.96].

Let n ≥ 1 and let I ∈ ℝ^n×n be the identity matrix. A random vector z = (z₁, …, z_n) has the distribution N(0, I) if and only if z₁, …, z_n are independent and each has the distribution N(0, 1). Then E‖z‖² = ∑_i=1ⁿ Ez_i² = n, so the length ‖z‖ is typically around √n (about 3.2 for n = 10).

Our program makes normal numbers with the Box–Muller method: if u₁, u₂ are independent and uniformly distributed on (0, 1), then

z = √−2 ln u₁ · cos(2πu₂)

has the distribution N(0, 1).

Let us try it on the simplest possible function, the sphere f(x) = x₁² + … + x₁₀² in ten dimensions. The minimum is obviously at zero; we start at the point (3, …, 3), where f = 90. How should we choose σ? A large value can get downhill quickly, but near the minimum it keeps overshooting. A small one is precise, but takes forever to reach the bottom.

2. The one-fifth rule

Definition — the one-fifth rule (Rechenberg, 1973)

Let a (1+1)-ES run with step size σ > 0. A mutation is successful if and only if the offspring is not worse than the parent, f(y) ≤ f(x). Let p_s ∈ [0, 1] be the share of successful mutations (over the most recent stretch of the run). The one-fifth rule changes σ as follows:

if p_s > 1/5, increase σ;
if p_s < 1/5, decrease σ;
if p_s = 1/5, keep σ.

The form used in this article changes σ after every mutation:

σ ← σ · exp(1/3) after a success, σ ← σ · exp(−1/12) after a failure.

If a mutation succeeds with probability p, ln σ changes per step by p · 1/3 − (1 − p) · 1/12 on average, and that is zero if and only if p = 1/5: 1/5 · 1/3 = 4/5 · 1/12. So the step size stays the same on average exactly when one fifth of mutations succeed.

In plain words

If almost every attempt succeeds, you are being needlessly careful — step it up. If almost none does, you are jumping too far — scale it back. Rechenberg worked out that progress is fastest when roughly every fifth attempt succeeds. The algorithm thus sets its own step size.

Both branches of the rule can be written as one formula; s = 1 after a success and 0 after a failure:

σ ← σ · exp( 5/12 · ( s − 1/5 ) )

σ: step size — it changes by multiplication, not addition: it stays positive and the step can shrink by many orders of magnitude (from 1 down to thousandths), as closeness to the minimum requires
5/12: adaptation rate — after a success exp(5/12 · 4/5) = exp(1/3) ≈ 1.40, after a failure exp(−5/12 · 1/5) = exp(−1/12) ≈ 0.92; a larger number would react faster, but σ would fluctuate more
s: outcome of the last attempt (1 = the offspring was not worse) — the only information the rule needs
1/5: target success rate — when exactly one fifth of attempts succeed, the mean of s − 1/5 is zero and σ stays the same on average

After N attempts with S successes the step is multiplied by exp(5/12 · (S − N/5)): out of ten attempts, two successes leave the step unchanged, four enlarge it 2.30 times, one shrinks it to 0.66.

(1+1)-ES on the sphere in ten dimensions fixed σ = 1σ = 0.1σ = 0.011/5 rule

The one-fifth rule gets below 10⁻⁶ in 587 evaluations — and the line in the chart falls steadily, because the step shrinks as fast as the minimum gets closer. With a fixed σ = 1 the algorithm soon hits a “floor” around one, below which it cannot get with such large jumps. With σ = 0.01 it goes down, but slowly: even after 5,000 evaluations it is not below 10⁻⁴.

import math, random

def es_1_plus_1(f, x, sigma=1.0, kroku=1000):
    fx = f(x)
    for _ in range(kroku):
        y = [xi + sigma * random.gauss(0, 1) for xi in x]
        fy = f(y)
        if fy <= fx:
            x, fx = y, fy
            sigma *= math.exp(1 / 3)     # success: step up
        else:
            sigma *= math.exp(-1 / 12)   # failure: scale back
    return x, fx

Example — the first twelve steps of the run in the chart

Start x = (3, …, 3), f(x) = 90, σ = 1. The σ column is the step size after the update:

step	`f`(`y`)	success?	`f`(`x`)	`σ`
1	115.31	no	90.00	0.920
2	95.78	no	90.00	0.846
3	123.94	no	90.00	0.779
4	107.86	no	90.00	0.717
5	101.54	no	90.00	0.659
6	113.11	no	90.00	0.607
7	94.03	no	90.00	0.558
8	95.09	no	90.00	0.513
9	86.27	yes	86.27	0.717
10	114.52	no	86.27	0.659
11	97.72	no	86.27	0.607
12	102.44	no	86.27	0.558

after 8 failures (s = 0): σ = 1 · [exp(5/12 · (0 − 1/5))]⁸ = exp(−8/12) ≈ 0.513

step 9 succeeds (86.27 < 90, s = 1): σ = 0.513 · exp(5/12 · (1 − 1/5)) = 0.513 · exp(1/3) ≈ 0.513 · 1.396 ≈ 0.717

(Colours as in the formula of the one-fifth rule above.)

With σ = 1 a step in ten dimensions has a length of about √10 ≈ 3.2 and the first eight offspring are worse than the start, so the rule pulls the step in. Successes then come more often and σ settles where roughly every fifth attempt succeeds — over the whole run 99 of 587 attempts succeeded, i.e. 17%. As the point approaches the minimum, the step shrinks with it: at the end of the run, after 580 evaluations, σ ≈ 0.00055.

3. Populations and self-adaptation

The sphere has a single minimum. Real problems often have lots of them — and there one individual is not enough. A typical test function is the Rastrigin function: f(x) = 10n + ∑ (x_i² − 10 cos 2πx_i). It is a sphere into which someone has carved a regular grid of dimples; the global minimum 0 is at the origin, surrounded by hundreds of local ones.

f(x) = 10n + ∑_i=1ⁿ ( x_i² − 10 cos 2πx_i )

10n: upward shift — at the origin the cosine subtracts at most 10 per coordinate, the constant adds it back so that the global minimum has a value of exactly 0
x_i²: sphere — the overall slope towards the origin; without it all dimples would be equally deep and no minimum would be global
10 cos 2πx_i: ripples with period 1 and amplitude 10 — they dig a dimple at every integer; the amplitude is much larger than the slope of the sphere between neighbouring dimples, so local search gets stuck in the nearest dimple

The Rastrigin function: a smooth sphere with a regularly rippled surface.

Example — the start and the neighbouring dimples (n = 2)

start (4.2, −3.8): x₁² = 17.64, 10 cos(2π · 4.2) = 10 cos(0.4π) ≈ 3.09

x₂² = 14.44, 10 cos(2π · (−3.8)) = 10 cos(0.4π) ≈ 3.09, 10n = 20

f = 20 + (17.64 − 3.09) + (14.44 − 3.09) ≈ 45.90 (colours as in the formula above)

bottom of the dimples next to the origin, at (±1, 0) and (0, ±1): f ≈ 0.99 (at the point 0.995, not exactly 1 — the sphere pulls it towards the origin)

bottom of the diagonal dimples at (±1, ±1): f ≈ 1.99

The value 1.49 in generation 8 (see below) therefore means a population in a neighbouring dimple, not yet at its bottom.

Definition — (μ, λ)-ES and (μ + λ)-ES with self-adaptation

Let f : ℝⁿ → ℝ be the function to be minimized and μ, λ ∈ ℕ. An individual is a pair (x, σ) ∈ ℝⁿ × (0, ∞) — a point and its own step size. A population is a μ-tuple of individuals (parents). Let τ > 0 be the learning rate of the step, usually τ ≈ 1/√2n. One generation:

Self-adaptation and mutation: for j = 1, …, λ pick a parent (x, σ) uniformly at random, independently draw ξ_j ~ N(0, 1) and z_j ~ N(0, I) and create the offspring (y_j, σ′_j), where σ′_j = σ · exp(τ · ξ_j), y_j = x + σ′_j · z_j. So the offspring first mutates the step and only then moves its position with the new step.
Selection:
- in the (μ, λ)-ES (requires λ > μ) the new population is formed by the μ offspring with the smallest values f(y_j) and the parents die;
- in the (μ + λ)-ES the new population is formed by the μ individuals with the smallest value of f out of all μ + λ parents and offspring together.

Generations are repeated until the chosen number of generations (or evaluations) is used up.

In plain words

The step size is inherited like any other trait. Whoever has a well-tuned step has better children — and so their step wins out too. Nobody has to count successes; evolution tunes itself. The (μ, λ) variant, in which parents do not survive, also cannot be blocked by a single individual that got lucky: every generation starts afresh.

σ′ = σ · exp(τ · N(0, 1)), τ = 1/√2n

σ′: the offspring’s step — the position is moved with it straight away, so the offspring’s fitness reveals whether the step was good
σ: the parent’s step — inherited; the offspring starts from what worked for the parent
τ: learning rate of the step — how much the step may change between generations; it shrinks as the number of dimensions n grows, because a single fitness value then carries less information about whether the step was good
N(0, 1): random change — symmetric around zero, so growing and shrinking are equally likely; the direction is decided only by selection

The step changes by multiplication with an exponential: it stays positive and doubling is as likely as halving.

Example — mutating the step for n = 2

τ = 1/√2 · 2 = 1/2

N(0, 1) = +1: σ′ = σ · exp(0.5 · 1) ≈ 1.65 σ; N(0, 1) = −1: σ′ = σ · exp(0.5 · (−1)) ≈ 0.61 σ

(Colours as in the formula above.)

A typical offspring therefore tries a step from roughly 40% smaller to 65% larger than its parent’s. In our run the best individual’s step fell in 15 generations from 2 to 0.011 — almost two hundred times.

import math, random

def es_mi_lambda(f, start, mi=15, lam=100, sigma=2.0, generaci=80):
    n = len(start)
    tau = 1 / math.sqrt(2 * n)
    rodice = [(start, sigma)] * mi
    for _ in range(generaci):
        deti = []
        for _ in range(lam):
            x, s = random.choice(rodice)                    # random parent
            s = s * math.exp(tau * random.gauss(0, 1))      # first mutate the step
            y = [xi + s * random.gauss(0, 1) for xi in x]   # then the position with it
            deti.append((y, s))
        deti.sort(key=lambda d: f(d[0]))
        rodice = deti[:mi]          # (μ, λ): only the best offspring survive
    return rodice[0]

(15, 100)-ES on the Rastrigin function; the cross is the global minimum

The population starts around the point (4.2, −3.8), where f ≈ 45.9. At first it is spread over several dimples at once and contracts towards the deeper ones. In generation 8 the best individual is still in a neighbouring dimple (f ≈ 1.49), but the step size has already fallen from 2 to 0.27. In generation 15 the population is at the global minimum and the step has shrunk to a hundredth so that it can fine-tune it.

Out of a hundred runs with different seeds, the (15, 100)-ES found the global minimum in 32 cases. For comparison, the (1+1)-ES with the one-fifth rule and the same budget of 8,000 evaluations: in just one. On spheres a single individual is fastest; on rugged terrain the population wins — and even it often ends up in one of the neighbouring dimples.

4. CMA-ES

Self-adaptation tunes only the length of the step. What if the minimum lies in a narrow slanted valley, where it pays to take long steps in one direction and short ones in another? CMA-ES (Hansen and Ostermeier, 2001) learns the whole covariance matrix of the mutation from successful steps — the shape of the cloud from which it draws offspring. During the search the cloud rotates and stretches along the valley. For optimizing continuous functions whose gradient we do not know, CMA-ES is today one of the most reliable methods.

Definition — mutation with a covariance matrix

Let n ≥ 1 and let C ∈ ℝ^n×n be a symmetric (C = C^T) positive definite matrix (v^TCv > 0 for all v ∈ ℝⁿ, v ≠ 0). A random vector z has the distribution N(0, C) if and only if it has the same distribution as

Au, where u ~ N(0, I) and A ∈ ℝ^n×n is any matrix with AA^T = C

(the distribution does not depend on the choice of A; such an A always exists, e.g. the Cholesky decomposition). Then Ez = 0 and E[z_iz_j] = C_ij for all i, j = 1, …, n; so C is the covariance matrix of the vector z.

Let m ∈ ℝⁿ be the centre of the population and σ > 0 the step size. Mutation with the covariance matrix C creates the offspring

y = m + σ · z, z ~ N(0, C).

Points with the same density lie on an ellipsoid; its axes are the eigenvectors of C and the lengths of the axes are proportional to the square roots of the eigenvalues. For C = I it is a sphere and we get an ordinary evolution strategy.

In simplified form (the real CMA-ES adds so-called evolution paths and its own control of σ, see N. Hansen, The CMA Evolution Strategy: A Tutorial, 2016) the matrix moves towards the successful steps after every generation:

C ← (1 − c) C + c ∑_i=1^μ w_i z_iz_i^T

C: covariance matrix of the mutation — the shape of the cloud from which offspring are drawn
1 − c: memory — the old matrix moves only slightly towards the new information, because one generation is a small and random sample
c: learning rate of the matrix, a number between 0 and 1 — how much weight the steps of the last generation get
∑ w_i: weighted average over the μ best offspring; the weights are positive, sum to 1 and a better offspring has a larger one — we learn only from steps that paid off
z_iz_i^T: the successful step z_i = y_i − m/σ as a matrix — stretches the cloud exactly in its direction

Directions in which things went well get a larger variance in the next generation.

Example — one successful step along the diagonal

C = I (a circular cloud), a single successful step z = (1, 1), so μ = 1 and w₁ = 1, c = 0.2 (colours as in the formula above):

zz^T = [[1, 1], [1, 1]]

C ← 0.8 · [[1, 0], [0, 1]] + 0.2 · 1 · [[1, 1], [1, 1]] = [[1, 0.2], [0.2, 1]]

eigenvalues: 1 + 0.2 = 1.2 in the direction (1, 1) and 1 − 0.2 = 0.8 in the direction (1, −1)

The cloud has stretched along the successful direction and narrowed across it — after a few such generations it lies along the valley.

5. Genetic programming

Definition — a program as a tree

Let 𝓕 be a finite set of functions, where each h ∈ 𝓕 has an arity ar(h) ∈ ℕ, ar(h) ≥ 1 (the number of arguments), and 𝓣 a finite set of terminals (variables and constants). The set of programs over 𝓕 and 𝓣 is the smallest set of rooted trees for which:

every terminal t ∈ 𝓣 is a program — a tree with the single node t;
if h ∈ 𝓕, k = ar(h) and g₁, …, g_k are programs, then the tree h(g₁, …, g_k) with root h and subtrees g₁, …, g_k in this order is a program too.

So the internal nodes are functions with as many children as their arity, and the leaves are terminals. The size of a program |g| is the number of its nodes. For given values of the variables, the value of a program is defined recursively from the leaves to the root: the value of a terminal is its value, and the value of h(g₁, …, g_k) is h applied to the values of g₁, …, g_k.

Closure: let D be the set of values (here ℝ). The set 𝓕 is closed if and only if every h ∈ 𝓕 is defined for all inputs, i.e. it is a function h : D^ar(h) → D; then every program has a value. That is why protected division is used, a ÷ b = a/b for |b| ≥ 10⁻⁹ and a ÷ b = 1 otherwise, and the square root of the absolute value, √|a|.

In plain words

The formula a/5 · √a is a tree: multiplication at the top, division and square root below it, a and a five at the bottom. Genetic programming crosses programs by snipping a branch off one and grafting a branch from the other in its place. Mutation replaces a branch with a newly grown random one. Protected operations make sure no program crashes; it can only produce nonsense — and selection weeds that out.

John Koza of Stanford University described the method in 1992 in the book Genetic Programming and demonstrated it on dozens of problems, from controlling a robot to finding formulas.

Example — evaluating the program tree for the Earth

The tree of the program a/5 · √a has multiplication at the root, the division a ÷ 5 on the left and √a on the right. From the leaves to the root, a = 149.6 (million km):

left branch: a/5 = 149.6/5 = 29.92

right branch: √a = √149.6 ≈ 12.231

root: 29.92 × 12.231 ≈ 365.95 days (actual 365.26, deviation +0.19%)

def gp(populace, generaci, p_krizeni=0.9, max_hloubka=7):
    for _ in range(generaci):
        dalsi = [nejlepsi(populace)]                  # elitism: the single best one
        while len(dalsi) < len(populace):
            if random.random() < p_krizeni:           # crossover by swapping subtrees
                a, b = turnaj(populace), turnaj(populace)
                dite = nahrad(a, nahodny_uzel(a), podstrom(b, nahodny_uzel(b)))
            else:                                     # mutation: branch → new random one
                a = turnaj(populace)
                dite = nahrad(a, nahodny_uzel(a), nahodny_strom(hloubka=2))
            if hloubka(dite) <= max_hloubka:         # throw away children that are too deep
                dalsi.append(dite)
        populace = dalsi
    return nejlepsi(populace)        # smallest error; on a tie the smaller tree

6. Symbolic regression: Kepler’s law

Johannes Kepler searched for ten years for the relationship between a planet’s distance from the Sun and its orbital period; he published the third law in 1619. Let us give genetic programming the same task. It has eight planets at its disposal — the semi-major axis a in millions of kilometres and the orbital period T in days —, the operations +, −, ×, ÷ and √, the variable a and the digits 1 to 9. The error of a program is the mean relative deviation from the actual orbital periods. The population has 300 programs and 40 generations.

Definition — symbolic regression

Let (x₁, y₁), …, (x_N, y_N) ∈ X × ℝ be data (X a set of inputs, e.g. ℝ^d), 𝓟 the set of programs over chosen functions 𝓕 and terminals 𝓣 (definition above; the variables correspond to the components of the input) and ℓ : ℝ × ℝ → [0, ∞] a loss function. Symbolic regression is the optimization problem of finding

g* ∈ argmin_{g ∈ 𝓟} ∑_i=1^N ℓ(g(x_i), y_i),

where g(x_i) is the value of the program g for the input x_i. Unlike ordinary regression, the shape of the model with unknown parameters is not given in advance — the expression itself is searched for, built only from the given building blocks.

In plain words

Linear regression from the article on statistics looks for the best straight line. Symbolic regression does not look for numbers in a known formula but for the formula itself — and a human can read the result.

Precisely, the error of a program g is the mean of the absolute differences of logarithms:

E(g) = 1/8 ∑_i=1⁸ | ln g(a_i) − ln T_i |

1/8 ∑: average over eight planets — the error does not depend on how much data we have and can be read as the “typical deviation per planet”
| … |: absolute value — a prediction 1% too large is as bad as one 1% too small, and errors of different signs do not cancel out
ln g(a_i) − ln T_i: the difference of logarithms is ln g(a_i)/T_i, i.e. the relative deviation: for small deviations ln(1 + r) ≈ r. Orbital periods differ almost a thousandfold (88 days for Mercury, 60,190 for Neptune); with an ordinary difference only Neptune would decide the error

The mean relative deviation, in which Mercury has the same weight as Neptune. A program with a negative or infinite result gets an error of ∞.

The best program of the initial, random population deviates from the actual orbital periods by 32% on average. By the third generation it has a deviation of 8% and the right shape a^1.5, just with a constant of 0.22. In the fifteenth generation this comes along:

Best program in generation 15 ((a ÷ 5) × √(a))

T = a · √a/5 — that is, T proportional to a^3/2, or T² proportional to a³. That is Kepler’s third law. The mean deviation is 0.30% and the constant 1/5 is very close to the true value 0.1996 for these units. Out of twenty runs with different seeds, eighteen found the exponent 1.5.

Where does 0.1996 come from? Newton’s law of gravitation gives, for a circular orbit (and for an ellipse with semi-major axis a, when the planet’s mass is negligible compared with the Sun’s), the form of the law including the constant:

T = 2π/√GM_☉ · a^3/2

T: orbital period — the time of one revolution
2π: one revolution is an angle of 2π; the formula comes from the angular velocity ω as T = 2π/ω
√GM_☉: the strength of the Sun’s gravity (G the gravitational constant, M_☉ the mass of the Sun) — stronger attraction means a faster orbit, hence it is in the denominator
a^3/2: distance — gravity weakens with a² and the path is longer in proportion to a; equating the centripetal and gravitational forces, ω²a = GM_☉/a², gives ω ∝ a^−3/2

Kepler’s third law derived from Newtonian gravity; genetic programming found the shape a^3/2 and the constant from just eight pairs of numbers.

Example — the constant in millions of km and days

The Sun’s gravitational parameter GM_☉ = 1.327124 · 10²⁰ m³/s² (the nominal value of the International Astronomical Union, Resolution B3 of 2015). The unit of length is 10⁹ m, so a^3/2 in metres is (10⁹)^3/2 = √10²⁷ times larger. The constant in seconds is 2π/√GM_☉ · √10²⁷ (colours as in the formula above):

√10²⁷/1.327124 · 10²⁰ = √7.535 · 10⁶ ≈ 2,745.0 s

2π · 2,745.0 s ≈ 17,247 s

17,247 s/86,400 s/day ≈ 0.1996 days

The program with the constant 1/5 = 0.2 therefore differs by 0.19% — exactly the deviation we got for the Earth above.

7. Bloat

But evolution did not stop at the beautiful formula. By generation 40 it had “improved” the deviation to 0.23% with this:

Best program in generation 40 (((a + ((((2 × 8) − a) ÷ 5) ÷ a)) ÷ 5) × √(a))

Instead of six nodes the tree has sixteen, and the added branches do nothing but fine-tune the constant and work around rounding in the data. Physically it is a worse formula — overfitted to eight planets. Meanwhile the average program size in the population grew from 9.8 nodes at the start to 23.8 in the twentieth generation.

Definition — bloat and parsimony pressure

Let P₀, P₁, … be the populations of programs in successive generations, |g| the size of a program g (the number of nodes of the tree) and F(g) ≥ 0 its error on the data (smaller is better). Denote by s̄_t the average size of the programs in P_t.

Bloat is the phenomenon in which s̄_t grows markedly with increasing t without a corresponding decrease in the error F — programs gain nodes that do not lower the error at all, or only negligibly. (It is an observed phenomenon, not a notion with a single precise boundary between growth that is still “corresponding” and growth that no longer is.)

Parsimony pressure is a modification of selection that favours smaller programs. Let c > 0. Two forms are used:

size penalty: selection compares programs by the adjusted error F′(g) = F(g) + c · |g|;
lexicographic comparison: a program g is better than h if and only if F(g) < F(h), or F(g) = F(h) ∧ |g| < |h| — size decides only when the errors are tied.

In plain words

Programs accumulate code that does nothing or makes a tiny correction — much as stretches without a function pile up in the genomes of living organisms. A large tree is also harder for crossover to break, so it multiplies more easily. Our algorithm prefers the smaller tree at equal error, but when a larger tree lowers the error by a hair, it wins. The same dilemma as with model overfitting: accuracy on the data, or simplicity that generalizes?

F′(g) = F(g) + c · |g|

F′(g): the adjusted error used for selection — smaller is better
F(g): the program’s error on the data — still the main thing we want to reduce
c: the cost of one node — says how much the error must improve for an extra node to pay off; it is chosen so that it decides only for negligible improvements
|g|: program size, the number of nodes in the tree — a measure of complexity that can be computed without understanding the program

An extra node has to pay for itself by lowering the error by at least c.

Example — how large a c would have saved Kepler’s law

Tree from generation 15: error 0.2986%, 6 nodes. Tree from generation 40: error 0.2312%, 16 nodes. The smaller tree wins when

0.002986 + 6c < 0.002312 + 16c

10c > 0.000674 ⇒ c > 0.000674/10 ≈ 6.7 · 10⁻⁵

(Colours as in the formula above: both sides of the inequality are adjusted errors F′.)

Ten extra nodes improved the error by less than a hundredth of a percent. With a cost above 0.007% per node, the tree from generation 40 would have lost in selection to the readable law from generation 15.

8. Evolution strategies and programs today

Evolution strategies instead of reinforcement learning. In 2017 OpenAI showed that a simple evolution strategy can train neural networks to control simulated robots and play Atari games about as well as reinforcement learning — and it can be spread across thousands of processors that exchange only a few numbers.
Discovering laws from data. In 2009 Schmidt and Lipson had symbolic regression derive laws of motion and conservation laws from measurements of pendulums. Today’s libraries such as PySR are used in physics and astronomy, where a readable formula is worth more than a black box.
Evolving machine learning algorithms. The AutoML-Zero project (Google, 2020) let evolution assemble learning algorithms from elementary mathematical operations and in doing so rediscovered, for example, backpropagation.
Programs with language models. The FunSearch and AlphaEvolve systems from the previous article are genetic programming in which mutation is done by a language model instead of a random branch.

9. Summary

Evolution strategies mutate real numbers with normal noise of step size σ.
The one-fifth rule and self-adaptation set σ automatically; CMA-ES also learns the shape of the mutations.
On spheres one individual is enough; on rugged terrain a population helps: on the Rastrigin function 32 successes out of 100 against one.
Genetic programming grows programs as trees; symbolic regression found Kepler’s third law in data about the planets.
Programs tend to bloat and overfit — parsimony pressure helps against that.

Want to get your hands on it? Have an AI agent write a symbolic regression and give it your own data — say, the period of a pendulum versus the length of the string — and watch whether it finds √L.

Do you see it differently? Write to support@coden-app.com.

Evolution strategies and genetic programming: from the one-fifth rule to Kepler’s law

1. Evolving numbers: the (1+1)-ES

2. The one-fifth rule

3. Populations and self-adaptation

4. CMA-ES

5. Genetic programming

6. Symbolic regression: Kepler’s law

7. Bloat

8. Evolution strategies and programs today

9. Summary

An Agent on your computer, driven from your phone