Skip to content

Nguyen-1

The problem

Symbolic regression: given points (x, y), find a formula for y in x. Here the points come from

y = x³ + x² + x,

at 20 values of x drawn at random from [−1, 1]. It's the first of the twelve problems of Uy et al. (2011), known as the Nguyen problems: polynomials, sines and cosines, logarithms, a square root and four functions of two variables, all with the same function set. McDermott et al. (2012, Genetic programming needs better benchmarks, GECCO 2012: 791-798) list them among the problems most used to compare genetic programming methods, together with Koza's quartic.

The formula is found when it's recovered exactly: an error at the level of rounding, both on the 20 training points and on 100 test points drawn from [−1, 1] that the search never sees.

The problem's name, target, sampling and function set are as McDermott et al. (2012) restate them, not yet checked against the paper. Its points come from a fixed seed of genoxide's random stream: another library's 20 points differ, and so can how hard the problem is (the All twelve tab shows how much).

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

Little, for genetic programming: it's the smallest of the Nguyen polynomials, and a tree of 11 nodes, x·x·x + (x + x·x), is enough. The function set is Koza's: 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, without constants. As with Koza's quartic, sines, cosines and exponentials of x fit the smooth curve closely without being it, and a search can settle on such a near miss; here that's rare.

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::Nguyen1: 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 100 test points. gp::Gp sets Koza's limits and initialization: trees of depth at most 17 and at most 1024 nodes, and ramped half-and-half with depths 2 to 6, the initial population of each island by Gp::ramped_half_and_half (Koza's even division, 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. With linear scaling (the error of a + b × tree, a and b fitted by least squares), the search recovers the formula less often: in 13 of 20 runs I tried, against every one without. Scaling makes any tree with the right shape up to a shift and a stretch as good as the formula, and many trees of sines and exponentials come close to that shape. The trees are evaluated on all 20 points at once (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 is invalid.

Algorithm

The same as Koza's quartic: 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 (with the crossover points at function nodes 90% of the time, within the limits), and subtree mutation at a rate of 0.1. The run stops at an RMSE of 10⁻¹⁰ times the standard deviation of the 20 values, or after 200 generations.

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.

Good results

The optimum is the polynomial itself, recovered exactly. The run of output.txt found it after 4 generations and 16,671 evaluations, as

add(mul(x, mul(x, x)), add(x, mul(x, x))),

which is x·x·x + (x + x·x) = x³ + x² + x, in 11 nodes, with an RMSE of about 10⁻¹⁷ on the training and 10⁻¹⁶ on the test points, the rounding of a different order of operations.

Over 100 runs (the islands' seeds 100 s + 0 to 7 for s = 1 to 100, the first being output.txt's, with the same data), all 100 recovered the formula, after 3 generations and 13,500 evaluations in the median, and 7 generations at most.

Reference: Uy, N. Q., Hoai, N. X., O'Neill, M., McKay, R. I. and Galván-López, E. (2011). Semantically-based crossover in genetic programming: application to real-valued symbolic regression. Genetic Programming and Evolvable Machines 12(2): 91-119.

Known optimum: x³ + x² + x, exactly (an RMSE of 0, up to rounding)

Source: examples/nguyen_1

Interactive run: tachsin.gr/projects/genoxide/examples/nguyen-1

cargo run --release --example nguyen_1
//! Nguyen-1: find the formula x³ + x² + x from 20 points of it, by genetic programming.
//!
//! The first of Nguyen's twelve symbolic regression problems (Uy et al. 2011): 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 100 test points in [−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 nguyen_1
//! ```

mod trace;

use genoxide::gp::regression::problems::{Nguyen1, 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<()> {
    // Nguyen's function set and sampling: 20 training points uniform in [-1, 1], from a fixed
    // seed, and 100 test points from another; the RMSE of the tree itself, without linear scaling
    let problem = Nguyen1::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!("Nguyen-1: 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("nguyen_1", &problem, &regression, -1.0..=1.0);
    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 100 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/nguyen_1/main.py
"""Nguyen-1: find the formula x^3 + x^2 + x from 20 points of it, by genetic programming.

The first of Nguyen's twelve symbolic regression problems (Uy et al. 2011): 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 100 test points in [-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/nguyen_1/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: 8.76e-11,
    0.00e0."""
    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 [-1, 1], from a fixed seed,
# and 100 test points from another; the RMSE of the tree itself, without linear scaling
problem = gx.gp.regression.problems.Nguyen1()
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("Nguyen-1: 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("nguyen_1", problem, regression, (-1.0, 1.0))
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 100 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:

Nguyen-1: x^3 + x^2 + x from 20 points in [-1, 1]
8 islands of 500 trees, until the RMSE is at most 8.76e-11

Target after 4 generations and 16671 evaluations
RMSE on the 20 training points: 6.40e-18
RMSE on the 100 test points:    9.06e-17

the expression, 11 nodes of depth 3:
add(mul(x, mul(x, x)), add(x, mul(x, x)))