Continuation by Gaussian smoothing

The problem

A Rastrigin function tilted by a quadratic centered off its lattice, in 10 dimensions:

f(x) = Σᵢ (xᵢ − aᵢ)² + 10 (1 − cos 2πxᵢ),   a = (1.3, −0.7, 2.2, −1.6, 0.35, 3.25, −2.8, 0.7, −0.3, 1.8)

Each gene has a well near every integer, and the quadratic makes the well nearest aᵢ the deepest. The function is separable, so its global minimum is computed exactly, gene by gene: in each well around an integer k near aᵢ, where the term is convex (|x − k| ≤ 1/4), the root of its derivative 2(x − aᵢ) + 20π sin 2πx by bisection to the last bit, and of those the lowest. The minimum is f* = 0.82084153378296, at x ≈ (1.0015, −0.9985, 2.0010, −1.9980, 0.0018, 3.0013, −2.9990, 0.9985, −0.0015, 1.9990).

What makes it hard

A well near every integer of the box [−5, 5] in every gene: 11¹⁰, about 2.6 × 10¹⁰ local minima, each a trap for a local method. From xᵢ = −3, L-BFGS-B settles in the well near −3 of every gene.

Smoothing removes them. The function averaged over a Gaussian of standard deviation σ, E[f(x + σz)], has a closed form, since E[cos 2π(x + σz)] = e^(−2π²σ²) cos 2πx:

f_σ(x) = Σᵢ (xᵢ − aᵢ)² + σ² + 10 (1 − e^(−2π²σ²) cos 2πxᵢ)

Its second derivative is at least 2 − 40π² e^(−2π²σ²), so f_σ is convex for σ above 0.517, with its minimum near a; and f₀ = f. Graduated smoothing (Blake and Zisserman's graduated non-convexity) follows the minimum from the convex function to the exact one: σ = 0.6, 0.4, 0.3, 0.2, 0.1 and 0, each stage started where the last ended. As σ falls, the wells come back around the point, and it slides into the one nearest a, the deepest.

The stages move the point: the smoothed minimum lies between aᵢ and the deepest well, and closes in on the well as σ falls (0.38 from the global minimum at σ = 0.6, then 4.1e-2, 9.6e-3, 2.4e-3 and 4.4e-4); the last stage, on f itself, ends at it.

Representation

A Real genome of 10 genes in [−5, 5], starting from xᵢ = −3. The fitness is f_σ, to minimize, with its gradient, 2(xᵢ − aᵢ) + 20π e^(−2π²σ²) sin 2πxᵢ: Differentiable in Rust and gradient=True in Python. The stage's σ is a value the fitness function reads, an Arc<AtomicU64> in Rust and a variable in Python. The cosines, sines and exponential are genoxide's portable ones, and the terms are summed in the same order in both, so both versions take the same steps.

Algorithm

Continuation around Lbfgsb, through 6 stages. A stage ends when L-BFGS-B has converged (its projected gradient within 1e-10); then the on_stage closure sets the next σ, the point is evaluated again on the new function, and L-BFGS-B goes on from it. It drops its curvature pairs between stages by default: they describe the last stage's function. The run stops by itself (StopReason::Converged) after the last stage.

Then, as contrasts from the same start: σ = 0 at once, L-BFGS-B on f alone; and the stages with the curvature pairs kept (keep_pairs(true)).

Output

A row per stage: its σ, its rounds and evaluations (the first evaluation of a stage re-evaluates its start), the stage's own minimum value, and the distance to the global minimum (the largest difference of a gene) where it ended. The last stage ends within 2.7e-15 of the global minimum, at f* to all 14 printed digits.

The project page plays the three runs back: the distance to the global minimum and the stage's σ at every round.

Good results

The global minimum f* = 0.82084153378296, reached in 36 evaluations.

Run, from xᵢ = −3	Rounds (evaluations)	Result
6 stages of smoothing (this example)	30 (36)	f*, within 2.7e-15 of the global minimum
σ = 0 from the start	5 (6)	trapped at f = 146.38192099, 145.56 above f*, 6.0 from the global minimum
6 stages, L-BFGS-B's pairs kept	50 (56)	within 4.4e-9 of the global minimum

L-BFGS-B alone converges fast, to the wrong minimum: the well it starts in. The stages cost 30 evaluations more, and find the global one. Keeping the curvature pairs across stages costs 20 evaluations more and digits of accuracy: the pairs describe the last stage's function, whose curvature changes most where the wells come back, and mislead the first steps of the next.

Smoothing finds the global minimum here because the deepest well is the one nearest the minimum of the convex smoothing, a property of this function (a quadratic with a periodic ripple), not of smoothing in general: for functions without that structure, graduated smoothing leads to a good minimum, not always the best.

Reference: Blake, A. and Zisserman, A. (1987). Visual Reconstruction. MIT Press.

Known optimum: f* = 0.82084153378296, computed gene by gene to the last bit

Source: examples/continuation

Interactive run: tachsin.gr/projects/genoxide/examples/continuation

RustPython

cargo run --release --example continuation

//! Continuation: a tilted Rastrigin function in 10 dimensions, Σ (xᵢ − aᵢ)² + 10 (1 − cos 2πxᵢ),
//! minimized through 6 stages of its Gaussian smoothing, σ from 0.6 (convex) to 0 (the function
//! itself), each stage from the last, with L-BFGS-B.
//!
//! The function is separable, so its global minimum is computed exactly, gene by gene, and the
//! last stage ends at it. Then, as contrasts, σ = 0 from the same start, which a local method
//! can't take out of the nearest basin, and the stages with L-BFGS-B's curvature pairs kept.
//!
//! With `GENOXIDE_TRACE=<file>`, it also writes a trace of its runs for the plot on the example's
//! page, with `trace.rs`.
//!
//! ```text
//! cargo run --release --example continuation
//! ```

mod trace;

use genoxide::math;
use genoxide::prelude::*;
use std::f64::consts::{PI, TAU};
use std::sync::atomic::{AtomicU64, Ordering};
use std::sync::{Arc, Mutex};

// the centers of the quadratic, off the cosine's lattice, and the cosine's amplitude
const CENTERS: [f64; 10] = [1.3, -0.7, 2.2, -1.6, 0.35, 3.25, -2.8, 0.7, -0.3, 1.8];
const A: f64 = 10.0;
// the stages' smoothing: 0.6 is convex (above 0.517), 0 the function itself
const SIGMAS: [f64; 6] = [0.6, 0.4, 0.3, 0.2, 0.1, 0.0];
// every gene starts here
const START: f64 = -3.0;
// L-BFGS-B's largest projected gradient component at which a stage has converged
const TOLERANCE: f64 = 1e-10;

// what a run did: per stage, its σ, rounds, evaluations, best value and distance to the global
// minimum; and in all, its rounds, evaluations and last point
struct Run {
    stages: Vec<(f64, u64, u64, f64, f64)>,
    rounds: u64,
    evaluations: u64,
    point: Vec<f64>,
}

fn main() -> Result<()> {
    let exact = global_minimum();
    let minimum = value(&exact, 0.0);
    let mut trace = trace::Trace::from_env();

    println!(
        "A tilted Rastrigin function in {} dimensions, Σ (xᵢ − aᵢ)² + {A} (1 − cos 2πxᵢ), from \
         xᵢ = {START}",
        CENTERS.len()
    );
    println!(
        "L-BFGS-B through 6 stages of the function smoothed by a Gaussian of σ, each from the last"
    );
    let staged = run(&SIGMAS, false, &exact, &mut trace, 0)?;
    println!(
        " stage     σ  rounds  evaluations  value of the stage  distance to the global minimum"
    );
    for (index, &(sigma, rounds, evaluations, best, distance)) in staged.stages.iter().enumerate() {
        println!(
            "{:>6}  {sigma:>4.2}  {rounds:>6}  {evaluations:>11}  {best:>18.10}  {:>29}",
            index + 1,
            scientific(distance)
        );
    }
    let error = distance(&staged.point, &exact);
    println!(
        "converged after {} rounds and {} evaluations: f = {:.14}, within {} of the global \
         minimum",
        staged.rounds,
        staged.evaluations,
        value(&staged.point, 0.0),
        scientific(error)
    );
    println!("the global minimum, gene by gene by bisection: f* = {minimum:.14}");
    assert!(error < 1e-12);

    // the contrasts: σ = 0 from the same start, and the stages keeping the curvature pairs
    let cold = run(&SIGMAS[5..], false, &exact, &mut trace, 1)?;
    let trapped = value(&cold.point, 0.0);
    println!(
        "σ = 0 from the start: {} rounds and {} evaluations, trapped at f = {trapped:.8}, {:.8} \
         above f*, {} from the global minimum",
        cold.rounds,
        cold.evaluations,
        trapped - minimum,
        scientific(distance(&cold.point, &exact))
    );
    let paired = run(&SIGMAS, true, &exact, &mut trace, 2)?;
    println!(
        "the stages with L-BFGS-B's curvature pairs kept: {} rounds and {} evaluations, within {}",
        paired.rounds,
        paired.evaluations,
        scientific(distance(&paired.point, &exact))
    );
    trace.write();
    Ok(())
}

// L-BFGS-B through the stages of `sigmas` from the same start, keeping its pairs between them
// or not, until the last stage has converged; each round recorded in the trace as run `line`
fn run(
    sigmas: &[f64],
    keep_pairs: bool,
    exact: &[f64],
    trace: &mut trace::Trace,
    line: usize,
) -> Result<Run> {
    // the stage's σ, shared with the fitness function
    let sigma = Arc::new(AtomicU64::new(sigmas[0].to_bits()));
    let shared = Arc::clone(&sigma);
    let smoothed = Differentiable(move |x: &Reals, gradient: &mut [f64]| {
        let s = f64::from_bits(shared.load(Ordering::Relaxed));
        smoothed(x, s, gradient)
    });
    let lbfgsb = Lbfgsb::builder(Real::uniform(CENTERS.len(), -5.0..=5.0)?)
        .initial_genome(Reals::from(vec![START; CENTERS.len()]))
        .gradient_tolerance(TOLERANCE)
        .keep_pairs(keep_pairs)
        .minimize()
        .build()?;
    // each stage's distance to the global minimum when it ends
    let distances = Arc::new(Mutex::new(Vec::new()));
    let (ended, minimum) = (Arc::clone(&distances), exact.to_vec());
    let stages = sigmas.to_vec();
    let continuation = Continuation::builder(lbfgsb)
        .stages(sigmas.len())
        .on_stage(move |stage, _| {
            sigma.store(stages[stage].to_bits(), Ordering::Relaxed);
            Ok(())
        })
        .on_stage_finished(move |_, lbfgsb| {
            let point = lbfgsb.population()[0].genome();
            let mut ended = ended.lock().expect("not poisoned");
            ended.push(distance(point, &minimum));
        })
        .build()?;
    let mut engine = Engine::new(continuation, smoothed)
        .stop_when(Stop::evaluations(10_000))
        .control(|staged: &mut Continuation<Lbfgsb>, progress| {
            let point = staged.population()[0].genome();
            let s = sigmas[staged.stage()];
            trace.record(line, progress.generation(), distance(point, exact), s);
            Ok(())
        });
    let outcome = engine.run()?;
    assert_eq!(outcome.stop_reason(), StopReason::Converged);
    let distances = distances.lock().expect("not poisoned");
    let stages = engine
        .algorithm()
        .stages()
        .iter()
        .zip(distances.iter())
        .map(|(stage, &distance)| {
            let best = stage.best().score().expect("valid");
            let s = sigmas[stage.index()];
            (s, stage.generations(), stage.evaluations(), best, distance)
        })
        .collect();
    Ok(Run {
        stages,
        rounds: outcome.generations(),
        evaluations: outcome.evaluations(),
        point: engine.algorithm().population()[0].genome().to_vec(),
    })
}

// the function smoothed by a Gaussian of σ, E[f(x + σz)] for z standard normal, and its gradient
// into `gradient`. E[cos 2π(x + σz)] = e^(−2π²σ²) cos 2πx and E[(x + σz − a)²] = (x − a)² + σ²,
// so each gene's term is (x − a)² + σ² + A (1 − e^(−2π²σ²) cos 2πx), convex once
// 4π²A e^(−2π²σ²) < 2. With genoxide's portable cos, sin and exp, the same bits everywhere.
fn smoothed(x: &[f64], s: f64, gradient: &mut [f64]) -> f64 {
    let e = math::exp(-2.0 * PI * PI * s * s);
    let mut sum = 0.0;
    for ((g, &xi), &a) in gradient.iter_mut().zip(x).zip(&CENTERS) {
        let d = xi - a;
        sum += d * d + s * s + A * (1.0 - e * math::cos(TAU * xi));
        *g = 2.0 * d + 2.0 * PI * A * e * math::sin(TAU * xi);
    }
    sum
}

// the function smoothed by σ, without the gradient
fn value(x: &[f64], s: f64) -> f64 {
    smoothed(x, s, &mut vec![0.0; x.len()])
}

// the global minimum, gene by gene: in each basin around an integer k near the center a, where
// the term is convex (|x − k| ≤ 1/4), the root of its derivative 2(x − a) + 2πA sin 2πx by
// bisection to the last bit, and of those the lowest
fn global_minimum() -> Vec<f64> {
    let mut minimum = Vec::with_capacity(CENTERS.len());
    for &a in &CENTERS {
        let derivative = |x: f64| 2.0 * (x - a) + 2.0 * PI * A * math::sin(TAU * x);
        let term = |x: f64| (x - a) * (x - a) + A * (1.0 - math::cos(TAU * x));
        let mut best: Option<(f64, f64)> = None;
        for k in (a.floor() as i64 - 3)..=(a.ceil() as i64 + 3) {
            let (mut low, mut high) = (k as f64 - 0.25, k as f64 + 0.25);
            for _ in 0..100 {
                let middle = 0.5 * (low + high);
                if derivative(middle) > 0.0 {
                    high = middle;
                } else {
                    low = middle;
                }
            }
            let x = 0.5 * (low + high);
            if best.is_none_or(|(lowest, _)| term(x) < lowest) {
                best = Some((term(x), x));
            }
        }
        minimum.push(best.expect("a basin").1);
    }
    minimum
}

// the largest difference of a gene from the global minimum's
fn distance(x: &[f64], exact: &[f64]) -> f64 {
    x.iter()
        .zip(exact)
        .map(|(a, b)| (a - b).abs())
        .fold(0.0, f64::max)
}

// one digit after the point, e.g. 1.2e-7
fn scientific(value: f64) -> String {
    format!("{value:.1e}")
}

python examples/continuation/main.py

"""Continuation: a tilted Rastrigin function in 10 dimensions, Σ (xᵢ − aᵢ)² + 10 (1 − cos 2πxᵢ),
minimized through 6 stages of its Gaussian smoothing, σ from 0.6 (convex) to 0 (the function
itself), each stage from the last, with L-BFGS-B.

The function is separable, so its global minimum is computed exactly, gene by gene, and the last
stage ends at it. Then, as contrasts, σ = 0 from the same start, which a local method can't take
out of the nearest basin, and the stages with L-BFGS-B's curvature pairs kept.

With ``GENOXIDE_TRACE=<file>``, it also writes a trace of its runs for the plot on the example's
page, with trace.py.

    python examples/continuation/main.py
"""

import math

import numpy as np

import genoxide as gx

from trace import Trace

# the centers of the quadratic, off the cosine's lattice, and the cosine's amplitude
CENTERS = np.array([1.3, -0.7, 2.2, -1.6, 0.35, 3.25, -2.8, 0.7, -0.3, 1.8])
A = 10.0
# the stages' smoothing: 0.6 is convex (above 0.517), 0 the function itself
SIGMAS = (0.6, 0.4, 0.3, 0.2, 0.1, 0.0)
# every gene starts here
START = -3.0
# L-BFGS-B's largest projected gradient component at which a stage has converged
TOLERANCE = 1e-10
PI, TAU = math.pi, math.tau


def smoothed(x, s):
    """The function smoothed by a Gaussian of σ, E[f(x + σz)] for z standard normal, and its
    gradient. E[cos 2π(x + σz)] = e^(−2π²σ²) cos 2πx and E[(x + σz − a)²] = (x − a)² + σ², so each
    gene's term is (x − a)² + σ² + A (1 − e^(−2π²σ²) cos 2πx), convex once 4π²A e^(−2π²σ²) < 2.
    With genoxide's portable cos, sin and exp, and the terms summed one after the other, the same
    bits as the Rust example."""
    e = float(gx.math.exp(np.array([-2.0 * PI * PI * s * s]))[0])
    d = x - CENTERS
    terms = d * d + s * s + A * (1.0 - e * gx.math.cos(TAU * x))
    gradient = 2.0 * d + 2.0 * PI * A * e * gx.math.sin(TAU * x)
    # a cumulative sum adds one term after the other, as the Rust example's loop
    return float(np.cumsum(terms)[-1]), gradient


def value(x, s):
    """The function smoothed by σ, without the gradient."""
    return smoothed(x, s)[0]


def global_minimum():
    """The global minimum, gene by gene: in each basin around an integer k near the center a,
    where the term is convex (|x − k| ≤ 1/4), the root of its derivative 2(x − a) + 2πA sin 2πx by
    bisection to the last bit, and of those the lowest. A gene's basins are bisected side by side,
    as an array."""
    minimum = []
    for a in CENTERS:
        k = np.arange(math.floor(a) - 3, math.ceil(a) + 4, dtype=float)
        low, high = k - 0.25, k + 0.25
        for _ in range(100):
            middle = 0.5 * (low + high)
            rising = 2.0 * (middle - a) + 2.0 * PI * A * gx.math.sin(TAU * middle) > 0.0
            high = np.where(rising, middle, high)
            low = np.where(rising, low, middle)
        x = 0.5 * (low + high)
        terms = (x - a) * (x - a) + A * (1.0 - gx.math.cos(TAU * x))
        # the first of the lowest, as the Rust loop keeps it
        minimum.append(x[int(np.argmin(terms))])
    return np.array(minimum)


def distance(x, exact):
    """The largest difference of a gene from the global minimum's."""
    return float(np.max(np.abs(x - exact)))


def scientific(value):
    """One digit after the point, e.g. 1.2e-7."""
    mantissa, exponent = f"{value:.1e}".split("e")
    return f"{mantissa}e{int(exponent)}"


trace = Trace()


def run(sigmas, keep_pairs, exact, line):
    """L-BFGS-B through the stages of ``sigmas`` from the same start, keeping its pairs between
    them or not, until the last stage has converged; each round recorded in the trace as run
    ``line``."""
    stage = {"index": 0}
    ended = []
    last = {}

    def on_stage(index):
        stage["index"] = index

    def on_stage_finished(finished, point):
        ended.append(distance(point, exact))

    def record(running, progress):
        point = progress.population[0]
        trace.record(line, progress.generation, distance(point, exact), sigmas[stage["index"]])
        last["point"] = point

    lbfgsb = gx.Lbfgsb(
        gx.Real((-5.0, 5.0), length=len(CENTERS)),
        initial_genome=[START] * len(CENTERS),
        gradient_tolerance=TOLERANCE,
        keep_pairs=keep_pairs,
        objective="minimize",
    )
    continuation = gx.Continuation(
        lbfgsb, stages=len(sigmas), on_stage=on_stage, on_stage_finished=on_stage_finished
    )
    result = continuation.run(
        lambda x: smoothed(x, sigmas[stage["index"]]),
        gradient=True,
        evaluations=10_000,
        control=record,
    )
    assert result.stop_reason == "converged"
    stages = [
        (
            sigmas[finished.index],
            finished.generations,
            finished.evaluations,
            finished.best_fitness,
            ended[finished.index],
        )
        for finished in result.stages
    ]
    return stages, result, last["point"]


exact = global_minimum()
minimum = value(exact, 0.0)
print(
    f"A tilted Rastrigin function in {len(CENTERS)} dimensions, Σ (xᵢ − aᵢ)² + {A:g} (1 − cos "
    f"2πxᵢ), from xᵢ = {START:g}"
)
print("L-BFGS-B through 6 stages of the function smoothed by a Gaussian of σ, each from the last")
stages, staged, point = run(SIGMAS, False, exact, 0)
print(" stage     σ  rounds  evaluations  value of the stage  distance to the global minimum")
for index, (sigma, rounds, evaluations, best, gap) in enumerate(stages):
    print(
        f"{index + 1:>6}  {sigma:>4.2f}  {rounds:>6}  {evaluations:>11}  {best:>18.10f}  "
        f"{scientific(gap):>29}"
    )
error = distance(point, exact)
print(
    f"converged after {staged.generations} rounds and {staged.evaluations} evaluations: f = "
    f"{value(point, 0.0):.14f}, within {scientific(error)} of the global minimum"
)
print(f"the global minimum, gene by gene by bisection: f* = {minimum:.14f}")
assert error < 1e-12

# the contrasts: σ = 0 from the same start, and the stages keeping the curvature pairs
_, cold, point = run(SIGMAS[5:], False, exact, 1)
trapped = value(point, 0.0)
print(
    f"σ = 0 from the start: {cold.generations} rounds and {cold.evaluations} evaluations, trapped "
    f"at f = {trapped:.8f}, {trapped - minimum:.8f} above f*, "
    f"{scientific(distance(point, exact))} from the global minimum"
)
_, paired, point = run(SIGMAS, True, exact, 2)
print(
    f"the stages with L-BFGS-B's curvature pairs kept: {paired.generations} rounds and "
    f"{paired.evaluations} evaluations, within {scientific(distance(point, exact))}"
)
trace.write()

What it prints, from a seeded run:

A tilted Rastrigin function in 10 dimensions, Σ (xᵢ − aᵢ)² + 10 (1 − cos 2πxᵢ), from xᵢ = -3
L-BFGS-B through 6 stages of the function smoothed by a Gaussian of σ, each from the last
 stage     σ  rounds  evaluations  value of the stage  distance to the global minimum
     1  0.60       7            8      103.6083820481                         3.8e-1
     2  0.40       6            7       98.0869311788                         4.1e-2
     3  0.30       5            6       84.7785592764                         9.6e-3
     4  0.20       4            5       55.8118222263                         2.4e-3
     5  0.10       4            5       18.8330678713                         4.4e-4
     6  0.00       4            5        0.8208415338                        2.7e-15
converged after 30 rounds and 36 evaluations: f = 0.82084153378296, within 2.7e-15 of the global minimum
the global minimum, gene by gene by bisection: f* = 0.82084153378296
σ = 0 from the start: 5 rounds and 6 evaluations, trapped at f = 146.38192099, 145.56107945 above f*, 6.0e0 from the global minimum
the stages with L-BFGS-B's curvature pairs kept: 50 rounds and 56 evaluations, within 4.4e-9