Koza's quartic
The problem
Symbolic regression: given points (x, y), find a formula for y in x, not just its coefficients but its shape. Here the points come from the quartic polynomial
y = x⁴ + x³ + x² + x,
at 20 values of x drawn at random from [−1, 1]. Koza (1992) used it to introduce genetic programming for symbolic regression, and it became the problem most papers tested on: McDermott et al. (2012, Genetic programming needs better benchmarks, GECCO 2012: 791-798) found it the most used of all, and argued that it's a toy that says little about a method. It stays the classic first problem, as here.
The formula is found when it's recovered exactly: an error at the level of rounding, both on the 20 training points and on 101 test points spread evenly over [−1, 1] that the search never sees.
The Python version builds the same islands with gx.gp and evaluates the trees in Rust
(problem.regression(linear_scaling=False)): it prints the same output and writes the same
trace.
What makes it hard
Nothing tells the search what shape the formula has. It builds formulas from Koza's function set, addition, subtraction, multiplication, protected division (1 where the denominator is 0), sine, cosine, the exponential and a protected logarithm (ln |a|, and 0 at 0), and the variable x. There are no constants: a formula needs x / x for a 1.
Many formulas come close without being it. Sines, cosines and exponentials of x fit a smooth curve on [−1, 1] to within a few thousandths, and a search can settle there: the error is small but not zero, and the formula grows as it adds corrections. Only the polynomial itself, in some arrangement, has an error at rounding level.
Representation
A tree of a genetic program (gp::Tree): the functions at the inner nodes, x at the leaves, stored
as a flat array in prefix order. The problem is gp::regression::problems::Koza1: its primitives
are gp::regression::Math values in a gp::PrimitiveSet of one type (add, sub, mul, the
protected division pdiv, sin, cos, exp and the protected logarithm plog), and its data
the 20 training and 101 test points. gp::Gp sets Koza's limits and initialization: trees of depth at most 17 (the root at depth
0) and at most 1024 nodes, and ramped half-and-half with depths 2 to 6. The initial population
of each island is Gp::ramped_half_and_half: Koza's even division, the same number of trees of each
depth, half by the full method and half by grow, without duplicates.
The fitness is the root mean squared error on the 20 points, minimized: gp::regression::Regression
without linear scaling, so the tree itself has to fit the points. The trees are evaluated on
all 20 points at once, a column of values per node (Tree::evaluate_columns), with genoxide's
math functions, so the values are the same bits on every platform. A tree with a value that isn't
finite (an exponential that overflows) is invalid.
Algorithm
Eight genetic algorithms of 500 trees each, as islands in a ring that pass their two best trees on every 10 generations. Each has tournaments of 7, subtree crossover at a rate of 0.9 (Koza's, with the crossover points at function nodes 90% of the time, and within the limits), and subtree mutation at a rate of 0.1 (a random node's subtree replaced by one grown to a depth of at most 4). The run stops at an RMSE of 10⁻¹⁰ times the standard deviation of the 20 values, or after 200 generations.
Islands, because a single population often settles on a near miss: with one population of 500 (Koza's size) and the same settings, 9 of the 20 seeds I tried found the formula in 200 generations. Separate islands settle on different near misses, and the formula found on one spreads to the others.
Output
The first two lines give the problem and the setting. Then come how the run stopped, after how many
generations and evaluations, the RMSE of the best tree on the training and the test points, and the
tree itself, with its size and depth, written as genoxide writes trees (Tree::display): each
function by its name with its arguments in parentheses.
The project page plays this run back.
Good results
The optimum is the quartic itself, recovered exactly. The run of output.txt found it after 4
generations and 17,061 evaluations, as
add(add(x, mul(x, x)), mul(add(x, mul(x, x)), mul(x, x))),
which is (x + x²) + (x + x²)·x² = x⁴ + x³ + x² + x, in 15 nodes. Its RMSE, about 10⁻¹⁶ on the training and the test points, is the rounding of a different order of operations.
Over seeds 1 to 100 (the islands' seeds, with the same data), 96 runs recovered the formula, after 7 generations and 27,000 evaluations in the median, and 22 generations at most. The other 4 ended after 200 generations on near misses of 131 to 385 nodes.
Known optimum: x⁴ + x³ + x² + x, exactly (an RMSE of 0, up to rounding)
Source: examples/koza_quartic
Interactive run: tachsin.gr/projects/genoxide/examples/koza-quartic
cargo run --release --example koza_quartic
//! Koza's quartic: find the formula x⁴ + x³ + x² + x from 20 points of it, by genetic
//! programming.
//!
//! The classic first problem of genetic programming (Koza 1992): trees of Koza's functions (+,
//! −, ×, protected division, sin, cos, exp and a protected logarithm) and the variable x, evolved
//! by a genetic algorithm with subtree crossover and mutation, fitted to the root mean squared
//! error on the points. The run stops at exact recovery: an error at the level of rounding, on
//! the 20 training points and on 101 test points across [−1, 1], with the expression printed.
//!
//! With `GENOXIDE_TRACE=<file>`, it also writes a trace of its run for the plot on the example's
//! page, with `trace.rs`.
//!
//! ```text
//! cargo run --release --example koza_quartic
//! ```
mod trace;
use genoxide::gp::regression::problems::{Koza1, RegressionProblem};
use genoxide::gp::{Gp, SubtreeCrossover, SubtreeMutation, Tree};
use genoxide::prelude::*;
// the islands, their trees, and the generations between migrations
const ISLANDS: u64 = 8;
const POPULATION: usize = 500;
const INTERVAL: u64 = 10;
fn main() -> Result<()> {
// Koza's function set and sampling: 20 training points uniform in [-1, 1], from a fixed seed,
// and 101 test points evenly spaced; the RMSE of the tree itself, without linear scaling
let problem = Koza1::new();
let regression = problem.regression().clone().linear_scaling(false);
let dataset = problem.dataset();
let (training, test) = (dataset.training(), dataset.test().expect("a test sample"));
// exact recovery: an error of at most 1e-10 of the values' standard deviation
let tolerance = 1e-10 * training.deviation();
println!("Koza's quartic x^4 + x^3 + x^2 + x from 20 points in [-1, 1]");
println!(
"{ISLANDS} islands of {POPULATION} trees, until the RMSE is at most {tolerance:.2e}\n"
);
// Koza's limits and initialization: depth 17, ramped half-and-half of depths 2 to 6
let gp = Gp::builder(problem.primitives().clone()).build()?;
let set = gp.primitives().clone();
let islands = (0..ISLANDS)
.map(|island| {
let seed = 100 + island;
// Koza's even division among the depths and methods, without duplicates
let initial =
gp.ramped_half_and_half(POPULATION, &mut StreamRng::seed_from_u64(seed))?;
Ga::builder(gp.clone())
.population_size(POPULATION)
.initial_genomes(initial)
.select(Tournament::new(7)?)
.crossover(SubtreeCrossover::new())
.mutate(SubtreeMutation::new())
.crossover_rate(0.9)
.mutation_rate(0.1)
.minimize()
.seed(seed)
.build()
})
.collect::<Result<Vec<_>>>()?;
let islands = Islands::builder(islands)
.interval(INTERVAL)
.migrants(2)
.build()?;
// with GENOXIDE_TRACE=<file>, a trace of the run for the plot on the example's page
let mut trace = trace::Trace::from_env(&problem, ®ression);
let outcome = Engine::new(islands, regression.clone())
.stop_when(Stop::target(tolerance).or(Stop::generations(200)))
.on_generation(|snapshot| trace.record(snapshot))
.run()?;
let best: &Tree = outcome.best_genome();
let rmse = |sample| regression.error(best, sample).unwrap_or(f64::NAN);
println!(
"{:?} after {} generations and {} evaluations",
outcome.stop_reason(),
outcome.generations(),
outcome.evaluations()
);
println!("RMSE on the 20 training points: {:.2e}", rmse(training));
println!("RMSE on the 101 test points: {:.2e}", rmse(test));
println!(
"\nthe expression, {} nodes of depth {}:\n{}",
best.len(),
best.depth(&set),
best.display(&set)
);
trace.write();
Ok(())
}
python examples/koza_quartic/main.py
"""Koza's quartic: find the formula x^4 + x^3 + x^2 + x from 20 points of it, by genetic
programming.
The classic first problem of genetic programming (Koza 1992): trees of Koza's functions (+, -,
x, protected division, sin, cos, exp and a protected logarithm) and the variable x, evolved by a
genetic algorithm with subtree crossover and mutation, fitted to the root mean squared error on
the points. The run stops at exact recovery: an error at the level of rounding, on the 20
training points and on 101 test points across [-1, 1], with the expression printed.
The trees and their errors are computed in Rust (``gx.gp``), with genoxide's portable math, so
the run is the Rust example's, to the bit, on every platform.
With ``GENOXIDE_TRACE=<file>``, it also writes a trace of its run for the plot on the example's
page, with trace.py.
python examples/koza_quartic/main.py
"""
import genoxide as gx
from trace import Trace
# the islands, their trees, and the generations between migrations
ISLANDS = 8
POPULATION = 500
INTERVAL = 10
def scientific(value):
"""``value`` with 2 decimals and an exponent, as Rust's ``{:.2e}`` writes it: 1.27e-10,
0.00e0."""
if value != value:
return "NaN"
mantissa, exponent = f"{value:.2e}".split("e")
return f"{mantissa}e{int(exponent)}"
# Koza's function set and sampling: 20 training points uniform in [-1, 1], from a fixed seed,
# and 101 test points evenly spaced; the RMSE of the tree itself, without linear scaling
problem = gx.gp.regression.problems.Koza1()
regression = problem.regression(linear_scaling=False)
dataset = problem.dataset()
training, test = dataset.training, dataset.test
# exact recovery: an error of at most 1e-10 of the values' standard deviation
tolerance = 1e-10 * training.deviation()
print("Koza's quartic x^4 + x^3 + x^2 + x from 20 points in [-1, 1]")
print(f"{ISLANDS} islands of {POPULATION} trees, until the RMSE is at most {scientific(tolerance)}")
print()
# Koza's limits and initialization: depth 17, ramped half-and-half of depths 2 to 6
gp = gx.gp.Gp(problem.primitives())
islands = []
for island in range(ISLANDS):
seed = 100 + island
islands.append(
gx.Ga(
gp,
population_size=POPULATION,
# Koza's even division among the depths and methods, without duplicates
initial_genomes=gp.ramped_half_and_half(POPULATION, seed),
select=gx.Tournament(7),
crossover=gx.gp.SubtreeCrossover(),
mutation=gx.gp.SubtreeMutation(),
crossover_rate=0.9,
mutation_rate=0.1,
objective="minimize",
seed=seed,
)
)
islands = gx.Islands(islands, interval=INTERVAL, migrants=2)
# with GENOXIDE_TRACE=<file>, a trace of the run for the plot on the example's page
trace = Trace(problem, regression)
result = islands.run(
regression, target=tolerance, generations=200, on_generation=trace.on_generation
)
best = result.best_genome
def rmse(sample):
error = regression.error(best, sample)
return float("nan") if error is None else error
print(
f"{result.stop_reason.capitalize()} after {result.generations} generations and "
f"{result.evaluations} evaluations"
)
print(f"RMSE on the 20 training points: {scientific(rmse(training))}")
print(f"RMSE on the 101 test points: {scientific(rmse(test))}")
print()
print(f"the expression, {len(best)} nodes of depth {best.depth}:")
print(best)
trace.write()
What it prints, from a seeded run:
Koza's quartic x^4 + x^3 + x^2 + x from 20 points in [-1, 1]
8 islands of 500 trees, until the RMSE is at most 1.27e-10
Target after 4 generations and 17061 evaluations
RMSE on the 20 training points: 1.12e-16
RMSE on the 101 test points: 1.30e-16
the expression, 15 nodes of depth 4:
add(add(x, mul(x, x)), mul(add(x, mul(x, x)), mul(x, x)))