Accuracy against size
The problem
Symbolic regression of Nguyen-7 (Uy et al. 2011): find a formula for y in x from the points of
y = ln(x + 1) + ln(x² + 1),
at 20 values of x drawn at random from [0, 2], with Koza's function set (addition, subtraction, multiplication, protected division, sine, cosine, the exponential and a protected logarithm, and x). The formula is a tree of 13 nodes, plog(x + x/x) + plog(x·x + x/x), the 1s built from x / x. But genetic programming rarely finds it: the All twelve tab of the Nguyen problems shows the usual search ending on trees of hundreds of nodes that fit the points to three or four decimals.
That's the common case in symbolic regression on real data, where there is no exact formula to find: every error can be lowered a little more by a larger tree. The question is then not which formula, but which trade-off. Minimizing the error and the size together, as two objectives, answers it in one run: the Pareto front, the smallest error found for each size, from the single variable x to trees of 50 nodes, where the small ones are formulas a person can read.
The Python version runs the same NSGA-II with gx.gp, both objectives evaluated in Rust
(gx.gp.WithSize(regression)): it prints the same output and writes the same trace.
What makes it hard
A search on the error alone grows its trees (bloat): each small improvement adds nodes, and a larger tree fits the training points better whether or not it's closer to the formula. The second objective stops that, but it pulls the other way: a population that keeps small trees has fewer large ones to recombine, and the most accurate end of the front is less accurate than a single-objective search would reach.
Representation
Trees of gp::regression::problems::Nguyen7's primitives (add, sub, mul, the protected
division pdiv, sin, cos, exp, the protected logarithm plog and x), with Koza's limits
and initialization in gp::Gp: depth at most 17, at most 1024 nodes, ramped half-and-half with
depths 2 to 6.
The two objectives, both minimized: the root mean squared error on the 20 points after linear
scaling (gp::regression::Regression: the error of a + b × tree, with a and b fitted by least
squares, so a tree only needs the shape), and the number of nodes. A tree with a value that isn't
finite at a point is invalid.
Algorithm
NSGA-II (Nsga2, Deb et al. 2002) on the trees, with a population of 1000, subtree crossover at a
rate of 0.9 and subtree mutation at a rate of 0.1, for 200 generations. NSGA-II keeps the trees of
the best non-dominated fronts, the most spread-out first, and removes duplicates, so the population
stays spread over the sizes instead of converging on one.
Output
The first two lines give the problem and the setting. Then comes the Pareto front, one line per point, from the smallest tree to the most accurate: its size, its RMSE on the training points and on 100 test points from [0, 2] that the search never saw, and, up to 25 nodes, the formula with its scaling. Trees that differ only in the order of their arguments have the same error and size: the front lists one of them. Two points that look alike (sizes 1 and 3, x and plog(exp(x))) differ in the last digits of their errors, by rounding.
The project page plays the run back: the front of each generation, as the logarithm of the error against the size.
Good results
This run's front has 27 points, from x alone (size 1, RMSE 0.022) to a tree of 50 nodes (RMSE 1.1 × 10⁻⁴). Along it:
- 7 nodes cut the error of a straight line to a third: −0.8709 + 1.2597 × (cos(cos(cos(cos(x)))) + x), RMSE 0.0074;
- 13 nodes cut it tenfold, to 0.0020;
- 25 nodes reach 4.2 × 10⁻⁴ on the training and 6.7 × 10⁻⁴ on the test points.
Beyond about 25 nodes, the training error keeps falling, to 1.1 × 10⁻⁴, but the test error rises, to 2 × 10⁻³: the larger trees fit the 20 points more closely and the formula less. The front shows where that starts, which a search on the error alone can't.
The exact formula, 13 nodes, isn't on the front. In 12 more runs I tried, seeds 2 to 4 with these settings and 9 with a population of 2000, 500 generations or a mutation rate of 0.3, the best error stayed between 6 × 10⁻⁵ and 5 × 10⁻⁴, without exact recovery.
Known optimum: The formula, ln(x + 1) + ln(x² + 1), is a tree of 13 nodes; this run's front doesn't reach it, and shows the best error found for each size instead
Source: examples/accuracy_and_size
Interactive run: tachsin.gr/projects/genoxide/examples/accuracy-and-size
cargo run --release --example accuracy_and_size
//! Accuracy against size: the trade-off between a formula's error and its size on Nguyen-7,
//! ln(x + 1) + ln(x² + 1), by NSGA-II on trees.
//!
//! Genetic programming rarely recovers Nguyen-7 exactly: the search settles on trees of hundreds
//! of nodes that fit the points to a few decimals. Minimizing the error and the size together, as
//! two objectives, gives instead the whole trade-off at once, the Pareto front: for each size, the
//! smallest error found, from a constant to large and accurate trees, with the small ones readable.
//!
//! 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 accuracy_and_size
//! ```
mod trace;
use genoxide::Objective::Minimize;
use genoxide::engine::FitnessFunction;
use genoxide::gp::regression::problems::{Nguyen7, RegressionProblem};
use genoxide::gp::{Gp, SubtreeCrossover, SubtreeMutation, Tree};
use genoxide::prelude::*;
const POPULATION: usize = 1000;
const GENERATIONS: u64 = 200;
fn main() -> Result<()> {
// Nguyen's function set and sampling: 20 training points uniform in [0, 2], from a fixed seed,
// and 100 test points from another; the RMSE after linear scaling
let problem = Nguyen7::new();
let regression = problem.regression().clone();
let test = problem.dataset().test().expect("a test sample");
println!("Nguyen-7: ln(x + 1) + ln(x^2 + 1) from 20 points in [0, 2]");
println!(
"NSGA-II, {POPULATION} trees for {GENERATIONS} generations, minimizing the RMSE and the size\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 initial = gp.ramped_half_and_half(POPULATION, &mut StreamRng::seed_from_u64(1))?;
let nsga2 = Nsga2::builder(gp, [Minimize, Minimize])
.population_size(POPULATION)
.initial_genomes(initial)
.crossover(SubtreeCrossover::new())
.mutate(SubtreeMutation::new())
.crossover_rate(0.9)
.mutation_rate(0.1)
.seed(1)
.build()?;
// the two objectives: the training RMSE, and the number of nodes; a tree whose value isn't
// finite at a point is invalid
let objectives = |tree: &Tree| -> Option<[f64; 2]> {
let error = regression.evaluate(tree)?;
Some([error, tree.len() as f64])
};
// with GENOXIDE_TRACE=<file>, a trace of the run for the plot on the example's page
let mut trace = trace::Trace::from_env();
let outcome = MultiEngine::new(nsga2, objectives)
.stop_when(Stop::generations(GENERATIONS))
.parallel(true)
.on_generation(|snapshot| trace.record(snapshot))
.run()?;
// the front, from the smallest tree to the most accurate: one tree per point, since trees that
// differ only in the order of their arguments have the same error and size
let mut front: Vec<(&Tree, f64)> = outcome
.front()
.iter()
.filter_map(|individual| Some((individual.genome(), individual.fitness()?.values()?[0])))
.collect();
front.sort_by(|(a, error_a), (b, error_b)| {
a.len().cmp(&b.len()).then(error_a.total_cmp(error_b))
});
front.dedup_by(|(a, error_a), (b, error_b)| a.len() == b.len() && error_a == error_b);
println!(
"the Pareto front after {} evaluations: {} points",
outcome.evaluations(),
front.len()
);
println!(
"{:>5} {:>9} {:>9} a + b * (expression), up to 25 nodes",
"size", "RMSE", "test RMSE"
);
for (tree, error) in front {
let test_error = regression.error(tree, test).unwrap_or(f64::NAN);
let expression = match regression.scaling(tree) {
Some(scaling) if tree.len() <= 25 => format!(
"{:.4} {} {:.4} * ({})",
scaling.intercept,
if scaling.slope < 0.0 { '-' } else { '+' },
scaling.slope.abs(),
tree.display(regression.primitives())
),
_ => String::new(),
};
println!(
"{:>5} {error:>9.2e} {test_error:>9.2e} {expression}",
tree.len()
);
}
trace.write();
Ok(())
}
python examples/accuracy_and_size/main.py
"""Accuracy against size: the trade-off between a formula's error and its size on Nguyen-7,
ln(x + 1) + ln(x^2 + 1), by NSGA-II on trees.
Genetic programming rarely recovers Nguyen-7 exactly: the search settles on trees of hundreds of
nodes that fit the points to a few decimals. Minimizing the error and the size together, as two
objectives, gives instead the whole trade-off at once, the Pareto front: for each size, the
smallest error found, from a constant to large and accurate trees, with the small ones readable.
Both objectives are computed in Rust (``gx.gp.WithSize``), in parallel and 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/accuracy_and_size/main.py
"""
import genoxide as gx
from trace import Trace
POPULATION = 1000
GENERATIONS = 200
def scientific(value):
"""``value`` with 2 decimals and an exponent, as Rust's ``{:.2e}`` writes it: 2.15e-2."""
if value != value:
return "NaN"
mantissa, exponent = f"{value:.2e}".split("e")
return f"{mantissa}e{int(exponent)}"
# Nguyen's function set and sampling: 20 training points uniform in [0, 2], from a fixed seed,
# and 100 test points from another; the RMSE after linear scaling
problem = gx.gp.regression.problems.Nguyen7()
regression = problem.regression()
test = problem.dataset().test
print("Nguyen-7: ln(x + 1) + ln(x^2 + 1) from 20 points in [0, 2]")
print(
f"NSGA-II, {POPULATION} trees for {GENERATIONS} generations, minimizing the RMSE and the size"
)
print()
# Koza's limits and initialization: depth 17, ramped half-and-half of depths 2 to 6
gp = gx.gp.Gp(problem.primitives())
nsga2 = gx.Nsga2(
gp,
objectives=["minimize", "minimize"],
population_size=POPULATION,
initial_genomes=gp.ramped_half_and_half(POPULATION, 1),
crossover=gx.gp.SubtreeCrossover(),
mutation=gx.gp.SubtreeMutation(),
crossover_rate=0.9,
mutation_rate=0.1,
seed=1,
)
# the two objectives: the training RMSE, and the number of nodes; a tree whose value isn't finite
# at a point is invalid
objectives = gx.gp.WithSize(regression)
# with GENOXIDE_TRACE=<file>, a trace of the run for the plot on the example's page
trace = Trace()
result = nsga2.run(
objectives, generations=GENERATIONS, parallel=True, on_generation=trace.on_generation
)
# the front, from the smallest tree to the most accurate: one tree per point, since trees that
# differ only in the order of their arguments have the same error and size
valid = [
(tree, float(values[0]))
for tree, values in zip(result.front_genomes, result.front_objectives)
if values[0] == values[0]
]
front = []
for tree, error in sorted(valid, key=lambda pair: (len(pair[0]), pair[1])):
if not front or (len(front[-1][0]), front[-1][1]) != (len(tree), error):
front.append((tree, error))
print(f"the Pareto front after {result.evaluations} evaluations: {len(front)} points")
print(f"{'size':>5} {'RMSE':>9} {'test RMSE':>9} a + b * (expression), up to 25 nodes")
for tree, error in front:
test_error = regression.error(tree, test)
scaling = regression.scaling(tree)
expression = ""
if scaling is not None and len(tree) <= 25:
intercept, slope = scaling
sign = "-" if slope < 0.0 else "+"
expression = f"{intercept:.4f} {sign} {abs(slope):.4f} * ({tree})"
test_text = scientific(float("nan") if test_error is None else test_error)
print(f"{len(tree):>5} {scientific(error):>9} {test_text:>9} {expression}")
trace.write()
What it prints, from a seeded run:
Nguyen-7: ln(x + 1) + ln(x^2 + 1) from 20 points in [0, 2]
NSGA-II, 1000 trees for 200 generations, minimizing the RMSE and the size
the Pareto front after 201000 evaluations: 27 points
size RMSE test RMSE a + b * (expression), up to 25 nodes
1 2.15e-2 2.52e-2 -0.0667 + 1.4324 * (x)
3 2.15e-2 2.52e-2 -0.0667 + 1.4324 * (plog(exp(x)))
5 1.98e-2 2.47e-2 1.8355 + 0.6416 * (sub(x, exp(cos(x))))
6 1.36e-2 1.64e-2 0.0280 + 0.8741 * (mul(x, add(cos(x), x)))
7 7.44e-3 1.10e-2 -0.8709 + 1.2597 * (add(cos(cos(cos(cos(x)))), x))
8 6.63e-3 8.61e-3 -0.0338 - 1.2885 * (sub(sin(sin(x)), add(x, sin(x))))
9 6.35e-3 9.97e-3 0.3603 + 1.2770 * (add(x, plog(cos(sin(cos(cos(cos(x))))))))
10 2.39e-3 2.42e-2 1.0392 - 1.2709 * (sub(cos(sin(cos(cos(exp(sin(plog(x))))))), x))
13 1.95e-3 2.71e-3 -0.8051 + 1.5876 * (add(x, sin(pdiv(cos(cos(x)), add(pdiv(x, sin(x)), x)))))
14 1.73e-3 1.80e-3 -0.7760 + 1.5764 * (add(x, sin(sin(pdiv(cos(cos(x)), add(pdiv(x, sin(x)), x))))))
16 1.10e-3 1.27e-3 1.0307 + 0.5021 * (add(add(sub(x, exp(pdiv(sin(sin(x)), x))), cos(sin(cos(sin(x))))), x))
17 1.07e-3 1.49e-3 -0.6796 + 1.4806 * (add(pdiv(cos(cos(x)), add(pdiv(x, mul(cos(cos(pdiv(x, x))), x)), x)), x))
18 1.03e-3 1.14e-3 -0.6613 + 1.4773 * (add(x, sin(pdiv(cos(cos(x)), add(pdiv(x, mul(cos(cos(pdiv(x, x))), x)), x)))))
20 9.01e-4 1.49e-3 -0.6836 + 1.4878 * (add(pdiv(cos(cos(x)), add(pdiv(x, mul(cos(cos(cos(plog(cos(sin(sin(x))))))), x)), x)), x))
21 4.87e-4 1.23e-3 -0.6814 + 1.4846 * (add(pdiv(cos(cos(x)), add(pdiv(x, mul(cos(cos(cos(plog(cos(sin(sin(sin(x)))))))), x)), x)), x))
22 4.74e-4 1.16e-3 -0.6776 + 1.4828 * (add(pdiv(cos(cos(x)), add(pdiv(x, mul(cos(cos(cos(plog(cos(sin(cos(cos(sin(x))))))))), x)), x)), x))
24 4.68e-4 1.12e-3 -0.6772 + 1.4828 * (add(pdiv(cos(cos(x)), add(pdiv(x, mul(cos(cos(cos(mul(plog(cos(cos(cos(x)))), cos(sin(x)))))), x)), x)), x))
25 4.23e-4 6.70e-4 -0.6572 + 1.4752 * (add(sin(pdiv(cos(cos(x)), add(pdiv(x, mul(cos(cos(cos(mul(plog(cos(cos(x))), mul(cos(x), x))))), x)), x))), x))
26 4.23e-4 6.65e-4
30 4.19e-4 1.20e-3
31 3.96e-4 1.20e-3
42 3.11e-4 2.20e-3
44 2.92e-4 2.15e-3
45 2.80e-4 2.09e-3
48 2.78e-4 2.10e-3
49 2.67e-4 2.25e-3
50 1.09e-4 2.12e-3