Skip to content

Two spirals

The problem

Two spirals wind around the origin three times, one inside the other: each point of the first has its mirror image through the origin on the second. A classifier gets the 194 points, 97 on each spiral, and has to tell which spiral each is on. Lang and Witbrock (1988) posed it as a benchmark for training neural networks, and it has been one since, for learning rules and for neuroevolution alike.

The points are Lang and Witbrock's: point i of the first spiral, for i from 0 to 96, is at the angle i π / 16 and the radius 6.5 (104 − i) / 104, at (r sin θ, r cos θ), and the second spiral's point is at (−x, −y). Here they're scaled by 1/6.5, to [−1, 1]², and computed with genoxide's portable math::sin and math::cos, so they're the same bits on every platform.

The Python version computes them with gx.math.sin and gx.math.cos, and runs the network in Rust with gx.nn.Mlp.forward: it prints the same output and writes the same trace.

What makes it hard

The two classes are interleaved all the way around, so no simple boundary separates them: the boundary a network learns has to be a spiral too, and fit between points that are close together near the center. Lang and Witbrock needed a network with shortcut connections, and a single hidden layer of the usual kind trained by backpropagation often fails to learn it.

For an evolutionary method there's a second difficulty: the network needs thousands of weights. CMA-ES, genoxide's method for continuous problems of up to a few hundred genes, keeps a covariance matrix of the genes, 2,545² = 6.5 million entries here, and learns it too slowly at this size.

Representation

A multilayer perceptron (nn::Mlp) of 2 inputs, two hidden layers of 48 tanh units and a tanh output, with biases: 2,545 weights, a Real genome in [−3, 3] each (Mlp::representation). Its output, in (−1, 1), is positive for the first spiral and negative for the second.

The fitness is the mean squared error of the outputs to the targets 1 and −1, minimized. The run stops when the generation's best network classifies all 194 points, or after 1,000,000 evaluations.

Algorithm

OpenAI's evolution strategy (OpenEs, Salimans et al. 2017), the baseline of neuroevolution at scale. It keeps no population, only the mean of a search distribution, a network, and moves it along an estimate of the gradient of the fitness:

  • Each generation, it draws 50 perturbations ε, standard normal, and asks for the networks mean + σ ε and mean − σ ε (mirrored sampling, Brockhoff et al. 2010): 100 samples.
  • It replaces their fitness by their ranks, spread evenly over [−0.5, 0.5], the best at 0.5.
  • The gradient estimate is the sum of the perturbations weighted by their samples' ranks, and Adam (Kingma and Ba, 2015) moves the mean along it.
  • Its cost per sample is linear in the number of weights, without a covariance matrix.

The settings: σ 0.01 of each weight's range (0.06), Adam's learning rate 0.01 of the range, the mean evaluated each generation too (evaluate_mean, often the best network), and small initial weights, uniform in [−0.25, 0.25]: large ones saturate the tanh units. The samples are drawn in parallel (parallel_breeding), a random stream per pair, and evaluated in parallel: the results are the same on any number of threads.

Output

The first line gives the evaluations and generations until a network classified all the points, and the second its mean squared error. Then comes a map of the network's decision over [−1, 1]², a character per cell: # where it says the first spiral, . where it says the second, and the points as A (first spiral) and B (second). Every A lies in the # region and every B in the . region, along two spiral bands.

The project page plays the run back: the network's decision over the plane as it learns.

Good results

The goal is every point classified. The run of output.txt reaches it after 135,744 evaluations, in 1,343 generations.

Over seeds 1 to 20, all 20 runs classified every point, after 165,413 evaluations on average (a median of 165,186, from 111,100 to 236,340) and 1,637 generations, about 35 seconds on 4 cores.

For contrast, sep-CMA-ES (Cmaes with cmaes::Covariance::Diagonal, CMA-ES's variant for thousands of genes, with its default population of 4 + ⌊3 ln 2545⌋ = 27) from the same initial weights, with the same fitness and stop, and an initial step of 0.01 or 0.05 of the range, three seeds each: 5 of the 6 runs classified every point, after 251,829 to 435,375 evaluations; the sixth stopped at its budget of 1,080,000 evaluations with 193 points of 194. It works, at this size, but needs half again to twice the evaluations of OpenEs and learns only one variance per weight; full CMA-ES isn't practical here (its matrix alone has 6.5 million entries, updated every generation).

Backpropagation, which Lang and Witbrock used, computes the gradient exactly from one pass over the points; OpenEs estimates it from about 100 evaluations per generation, each a pass over the points. That is the price of not needing gradients, and what lets the same method evolve controllers and other networks whose fitness has none, such as the pole-balancing pages'.

Reference: Lang, K. J. and Witbrock, M. J. (1988). Learning to tell two spirals apart. Proceedings of the 1988 Connectionist Models Summer School: 52-59, who posed the problem (no DOI). The method, to which the link points: Salimans, T., Ho, J., Chen, X., Sidor, S. and Sutskever, I. (2017). Evolution strategies as a scalable alternative to reinforcement learning. arXiv:1703.03864.

Known optimum: All 194 points classified

Source: examples/two_spirals

Interactive run: tachsin.gr/projects/genoxide/examples/two-spirals

cargo run --release --example two_spirals

//! Two spirals: evolve the 2,545 weights of a neural network that tells two interleaved spirals
//! apart, by OpenAI's evolution strategy.
//!
//! Lang and Witbrock's (1988) benchmark: 194 points on two spirals that wind three times around
//! the origin, 97 each, one the mirror image of the other through the origin. A network with two
//! hidden layers of 48 tanh units (2,545 weights with the biases) outputs a value in [−1, 1] for
//! a point; its sign is the spiral. The fitness is the mean squared error to the targets 1 and −1,
//! minimized by `OpenEs` from small random weights, until the network classifies every point.
//!
//! 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 two_spirals
//! ```

mod trace;

use genoxide::math::{cos, sin};
use genoxide::nn::{Activation, Mlp};
use genoxide::prelude::*;
use std::sync::Arc;
use std::sync::atomic::{AtomicBool, Ordering};

// the points of each spiral
const PER_SPIRAL: u32 = 97;
// the largest radius, to which the coordinates are scaled: the points lie in [−1, 1]²
const RADIUS: f64 = 6.5;

// Lang and Witbrock's points, scaled to [−1, 1]², with the targets 1 and −1: point i of the first
// spiral at the angle i π / 16 and the radius 6.5 (104 − i) / 104, and its mirror image through
// the origin on the second
pub fn points() -> Vec<([f64; 2], f64)> {
    let mut points = Vec::new();
    for i in 0..PER_SPIRAL {
        let angle = f64::from(i) * std::f64::consts::PI / 16.0;
        let radius = RADIUS * f64::from(104 - i) / 104.0;
        let (x, y) = (radius * sin(angle) / RADIUS, radius * cos(angle) / RADIUS);
        points.push(([x, y], 1.0));
        points.push(([-x, -y], -1.0));
    }
    points
}

// 2 inputs, two hidden layers of 48 tanh units and a tanh output, with biases
pub fn network() -> Result<Mlp> {
    Ok(Mlp::new([2, 48, 48, 1], Activation::Tanh)?.output_activation(Activation::Tanh))
}

// the points on the right side of 0
fn classified(network: &Mlp, weights: &[f64], points: &[([f64; 2], f64)]) -> usize {
    let mut network = network.with(weights).expect("the network's weights");
    let mut output = [0.0];
    let right = |(input, target): &&([f64; 2], f64)| {
        network.forward(input, &mut output);
        output[0] * target > 0.0
    };
    points.iter().filter(right).count()
}

fn main() -> Result<()> {
    let points = points();
    let mlp = network()?;
    // the mean squared error to the targets
    let error = |weights: &Reals| -> Option<f64> {
        let mut network = mlp.with(weights).ok()?;
        let mut output = [0.0];
        let squares = points.iter().map(|(input, target)| {
            network.forward(input, &mut output);
            (output[0] - target) * (output[0] - target)
        });
        Some(squares.sum::<f64>() / points.len() as f64)
    };
    // small random weights to start from: large ones saturate the tanh units
    let initial = Real::uniform(mlp.parameters(), -0.25..=0.25)?
        .random_genome(&mut StreamRng::seed_from_u64(1));
    let open_es = OpenEs::builder(mlp.representation(-3.0..=3.0)?)
        .population_size(100)
        .sigma(0.01)
        .optimizer(open_es::Optimizer::adam(0.01))
        .evaluate_mean(true)
        .initial_mean(initial)
        .parallel_breeding(true)
        .minimize()
        .seed(1)
        .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(&mlp);
    let solved = Arc::new(AtomicBool::new(false));
    let stop = Arc::clone(&solved);
    let mut solution = None;
    let outcome = Engine::new(open_es, error)
        .parallel(true)
        .stop_when(Stop::custom(move |_| stop.load(Ordering::Relaxed)))
        .stop_when(Stop::evaluations(1_000_000))
        .on_generation(|snapshot| trace.record(snapshot))
        // the generation's best network, if it classifies every point
        .on_generation(|snapshot| {
            let best = snapshot.population().best(Objective::Minimize);
            let Some(best) = best else { return };
            if classified(&mlp, best.genome(), &points) == points.len() {
                let progress = snapshot.progress();
                let error = best.fitness().and_then(Fitness::score);
                solution = Some((best.genome().clone(), error, progress.evaluations()));
                solved.store(true, Ordering::Relaxed);
            }
        })
        .run()?;

    let Some((weights, error, evaluations)) = solution else {
        let best = outcome.best_genome();
        println!(
            "not solved after {} evaluations: {} of {} points classified",
            outcome.evaluations(),
            classified(&mlp, best, &points),
            points.len()
        );
        return Ok(());
    };
    println!(
        "all {} points classified after {evaluations} evaluations in {} generations, by a \
         network of {} weights",
        points.len(),
        outcome.generations(),
        mlp.parameters()
    );
    println!("mean squared error {:.6}", error.unwrap_or(f64::NAN));
    println!();
    // the network's decision over [−1, 1]², a character per cell: # for the first spiral's side
    // and . for the second's, the points as A and B
    let mut network = mlp.with(&weights)?;
    let (columns, rows) = (61, 31);
    let mut map = vec![vec![' '; columns]; rows];
    let mut output = [0.0];
    for (row, line) in map.iter_mut().enumerate() {
        let y = 1.0 - 2.0 * row as f64 / (rows - 1) as f64;
        for (column, cell) in line.iter_mut().enumerate() {
            let x = -1.0 + 2.0 * column as f64 / (columns - 1) as f64;
            network.forward(&[x, y], &mut output);
            *cell = if output[0] > 0.0 { '#' } else { '.' };
        }
    }
    for ([x, y], target) in &points {
        let column = ((x + 1.0) / 2.0 * (columns - 1) as f64).round() as usize;
        let row = ((1.0 - y) / 2.0 * (rows - 1) as f64).round() as usize;
        map[row][column] = if *target > 0.0 { 'A' } else { 'B' };
    }
    for line in map {
        println!("{}", line.into_iter().collect::<String>());
    }
    trace.write();
    Ok(())
}

python examples/two_spirals/main.py

"""Two spirals: evolve the 2,545 weights of a neural network that tells two interleaved spirals
apart, by OpenAI's evolution strategy.

Lang and Witbrock's (1988) benchmark: 194 points on two spirals that wind three times around the
origin, 97 each, one the mirror image of the other through the origin. A network with two hidden
layers of 48 tanh units (2,545 weights with the biases) outputs a value in [-1, 1] for a point;
its sign is the spiral. The fitness is the mean squared error to the targets 1 and -1, minimized
by OpenEs from small random weights, until the network classifies every point.

The network runs in Rust (``gx.nn.Mlp.forward``), and the points use genoxide's portable sine and
cosine (``gx.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/two_spirals/main.py
"""

import math

import numpy as np

import genoxide as gx

from trace import Trace

# the points of each spiral
PER_SPIRAL = 97
# the largest radius, to which the coordinates are scaled: the points lie in [-1, 1]²
RADIUS = 6.5


def spiral_points():
    """Lang and Witbrock's points, scaled to [-1, 1]², a row each, and their targets 1 and -1:
    point i of the first spiral at the angle i π / 16 and the radius 6.5 (104 - i) / 104, and its
    mirror image through the origin on the second."""
    points, targets = [], []
    for i in range(PER_SPIRAL):
        angle = i * math.pi / 16.0
        radius = RADIUS * (104 - i) / 104.0
        x, y = radius * gx.math.sin(angle) / RADIUS, radius * gx.math.cos(angle) / RADIUS
        points += [[x, y], [-x, -y]]
        targets += [1.0, -1.0]
    return np.array(points), np.array(targets)


def round_half_up(value):
    """``value``, 0 or more, rounded to the nearest whole number, halves up, as Rust rounds it
    (Python's ``round`` rounds halves to even)."""
    floor = math.floor(value)
    return floor + 1 if value - floor >= 0.5 else floor


points, targets = spiral_points()
# 2 inputs, two hidden layers of 48 tanh units and a tanh output, with biases
mlp = gx.nn.Mlp([2, 48, 48, 1], "tanh", output_activation="tanh")


def classified(weights):
    """The points on the right side of 0."""
    outputs = mlp.forward(weights, points)[:, 0]
    return int(np.count_nonzero(outputs * targets > 0.0))


def error(weights):
    """The mean squared error to the targets, the squares added in order, as Rust adds them."""
    differences = mlp.forward(weights, points)[:, 0] - targets
    return float(np.cumsum(differences * differences)[-1]) / len(points)


# small random weights to start from: large ones saturate the tanh units
initial = gx.Real((-0.25, 0.25), length=mlp.parameters).random_genome(1)
open_es = gx.OpenEs(
    mlp.representation((-3.0, 3.0)),
    population_size=100,
    sigma=0.01,
    optimizer=gx.Adam(0.01),
    evaluate_mean=True,
    initial_mean=initial,
    parallel_breeding=True,
    objective="minimize",
    seed=1,
)
# with GENOXIDE_TRACE=<file>, a trace of the run for the plot on the example's page
trace = Trace(mlp, points, targets)
solution = None


def on_generation(progress):
    """Records the generation, and stops the run once the generation's best network classifies
    every point."""
    global solution
    trace.record(progress)
    # the generation's best, the first of equals; an invalid score is NaN
    best = int(np.nanargmin(progress.scores))
    weights = progress.population[best]
    if classified(weights) == len(points):
        solution = (weights, float(progress.scores[best]), progress.evaluations)
        return False
    return True


result = open_es.run(error, evaluations=1_000_000, parallel=True, on_generation=on_generation)

if solution is None:
    print(
        f"not solved after {result.evaluations} evaluations: "
        f"{classified(result.best_genome)} of {len(points)} points classified"
    )
else:
    weights, mean_squared_error, evaluations = solution
    print(
        f"all {len(points)} points classified after {evaluations} evaluations in "
        f"{result.generations} generations, by a network of {mlp.parameters} weights"
    )
    print(f"mean squared error {mean_squared_error:.6f}")
    print()
    # the network's decision over [-1, 1]², a character per cell: # for the first spiral's side
    # and . for the second's, the points as A and B
    columns, rows = 61, 31
    cells = np.array(
        [
            [-1.0 + 2.0 * column / (columns - 1), 1.0 - 2.0 * row / (rows - 1)]
            for row in range(rows)
            for column in range(columns)
        ]
    )
    outputs = mlp.forward(weights, cells)[:, 0].reshape(rows, columns)
    grid = [["#" if output > 0.0 else "." for output in line] for line in outputs]
    for (x, y), target in zip(points, targets):
        column = round_half_up((x + 1.0) / 2.0 * (columns - 1))
        row = round_half_up((1.0 - y) / 2.0 * (rows - 1))
        grid[row][column] = "A" if target > 0.0 else "B"
    for line in grid:
        print("".join(line))
    trace.write()

What it prints, from a seeded run:

all 194 points classified after 135744 evaluations in 1343 generations, by a network of 2545 weights
mean squared error 0.064620

............................##A#####A########################
......................#........##########A###################
.........................B....B....##########################
....................B..............B....######A##############
...............B............####.......B......###############
................######A###A###A###A######...B....#A##########
...........B######A########........###A##......B...##########
...........###A#####....B..B..B..B...####A#.........##A######
........B..######....B..............B..#####A##...B..########
...........#A##...B....##A##A#A#A#....B..########....########
.....B..######.B.......A##########A#.....B.###A##...B.##A####
......###A#........A#A.......B...###AA####B..###A#.....######
...B..#####..B...#A#....B.BB..BBBB....A####B.#####......#A###
.....###A#..B..####...BB########..BB...A########A#...B..#####
.....#####.....##A#..B.###AAAAA....B..#A##..B.####.......####
..B.###A##..B..#A##..B.##AA###....BB###A###.B..#A#...B...#A##
....#######....#A###.B..#A###.BBBBB####A###B...###........###
....###A###.B..######B...AA##...#####AA##.....##A#..B....####
...B.######....##A###.B..##AAAA##AA#A##...B..##A#........A###
......#####.B...##A###.BB######A#......B.B.#####...B....#####
....B...A##...B...#A###...B..........B..#####A......###A#####
...................###A#....B.B.B.#B######A#....B...#########
........##A###..B.....##A##############A#..........#A########
......B..########..B...####A##A##A##A##.......B...###########
..........###A#######.B.....########......B....##A###########
..........B..###A#####....B...B###B...B.....#################
................#####A##........#.....#######A###############
..............B...#######A##############A####################
.........................#####A####A#########################
...................B.......####.#############################
........................B.....B..############################