Skip to content

Shekel 5

The problem

Shekel's function puts a well at each of m points aᵢ, to minimize:

f(x) = −Σᵢ₌₁ᵐ 1 / ((x − aᵢ)ᵀ(x − aᵢ) + cᵢ),   each xⱼ in [0, 10]

in 4 dimensions. Each term is −1/cᵢ at its center and falls off with the squared distance from it, so well i is 1/cᵢ deep and about √cᵢ wide. Here m = 5, with the first five rows of Dixon and Szegö's table (Shekel 7 and Shekel 10 take more of them):

i aᵢ cᵢ depth 1/cᵢ
1 4, 4, 4, 4 0.1 10
2 1, 1, 1, 1 0.2 5
3 8, 8, 8, 8 0.2 5
4 6, 6, 6, 6 0.4 2.5
5 3, 7, 3, 7 0.4 2.5

The function is Shekel's (1971), with the constants that Dixon and Szegö (1978) tabulate; they call it SQRIN5. Neither is online: genoxide takes the constants from Yao, Liu and Lin's (1999, table XIV) reprint.

The minimum is near a₁ = (4, 4, 4, 4), but not at it: the other wells pull it aside a little, and f(4, 4, 4, 4) = −10.153196, 4e-6 higher. genoxide's problems::Shekel5 gives the minimum as −10.153199679058227 at (4.000037152819676, 4.00013327659156, 4.000037152819676, 4.00013327659156): the point where the gradient is 0, computed to 40 digits by Newton's method. Later papers quote −10.1532 from Dixon and Szegö. It's the best known minimum, not proven global. Jamil and Yang (2013) put the minimum at (4, 4, 4, 4), with a value that is neither the minimum nor the value there.

What makes it hard

Each well is a local minimum, at nearly the value of its own term: −10.1532 at a₁, −5.1008 at a₃, −5.0552 at a₂, −2.6829 at a₄ and −2.6305 at a₅. The wells are narrow, and the rest of the box is a low plateau: at the center, (5, 5, 5, 5), f is −0.58, and in 0.1% of the box only is f below −1. The deepest well is also the narrowest (√0.1 = 0.32 against 0.45 to 0.63). Local searches from 2,000 random points end in a₁'s well 41% of the time, in a₄'s 39%, and in the other three 20%.

A method that learns where good points are from the points it has seen is misled by the first well it finds: a well of depth 2.5 or 5 is far better than the plateau, and pulls the search in before the deepest well has been sampled.

Representation

A Real genome of 4 genes, each in [0, 10]: the point x itself. The fitness is f(x), to minimize. The function, its bounds and its best known minimum are genoxide's problems::Shekel5.

Algorithm

Particle swarm optimization (Kennedy and Eberhart, 1995, Proceedings of ICNN'95: 1942-1948), with 80 particles and Clerc and Kennedy's constriction coefficients (2002, IEEE Transactions on Evolutionary Computation 6(1): 58-73), genoxide's defaults. Each particle moves towards its own best point and a best point of the swarm. Two topologies, each from seeds 1 to 30:

  • global: every particle follows the best point of the whole swarm. As soon as one particle finds a well, the whole swarm heads there;
  • ring: each particle follows the best of itself and its two neighbors on a ring. A good point spreads one neighbor per step, so parts of the swarm keep exploring other wells for longer.

Each run stops once its value is within 1e-6 of the best known minimum, or after 25,000 evaluations. A run is counted in the well whose center is nearest its best point.

The swarm is larger than the usual 20 to 50 particles, and the budget larger than most runs need: with 40 particles and 10,000 evaluations, the ring didn't reach the minimum in 2% to 3% of the runs (seeds 1 to 1,000), and with 80 particles and 25,000 evaluations it reaches it in all of them.

Output

The first line gives the best known minimum, the seeds, the budget and the swarm's size. Then a row per well that some run ended in, with how many runs of each topology ended there; then how many runs came within 1e-6 of the best known minimum, and the median number of evaluations they needed. In Python, run evaluates the function in Rust, so both versions print the same table.

The page's plot shows the 30 runs of the swarm with the global topology, each at its best point so far, and a curve of the best and the median run's distance above the best known minimum, on a logarithmic axis. The points are drawn at their (x₁, x₂), over the function on the plane x₃ = x₁, x₄ = x₂, which holds every center aᵢ of Shekel 5: a run in well i is drawn at its center.

The project page plays this run back.

Good results

A good result reaches −10.1532 in every run. The ring topology does: all 30 runs come within 1e-6 of it, after a median of 11,760 evaluations. With the global topology, for contrast, 14 runs end in the deepest well and 13 of them come within 1e-6 of its minimum, after a median of 7,040 evaluations; the other 16 end in the other four wells, 8 of them in a₃'s at −5.10. The ring is slower where both succeed, since good points spread more slowly.

Over seeds 1 to 1,000, the ring topology reaches the minimum in every run, after at most 17,680 evaluations, and the global topology in 47% of the runs. With 40 particles and 10,000 evaluations, they reach it in 97% and 39%.

Reference: Shekel, J. (1971). Test functions for multimodal search techniques. Proceedings of the 5th Annual Princeton Conference on Information Sciences and Systems, Princeton University. Constants as tabulated in Dixon, L. C. W. and Szegö, G. P. (1978). The global optimisation problem: an introduction. In Towards Global Optimisation 2, North-Holland: 1-15.

Known optimum: −10.15320 near (4, 4, 4, 4) (best known)

Source: examples/shekel5

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

cargo run --release --example shekel5

//! Shekel 5: minimize Shekel's function with 5 wells in 4 dimensions with particle swarms, from
//! 30 seeds, with a global and a ring topology.
//!
//! The function has a well at each of 5 points, the deepest at (4, 4, 4, 4). A swarm whose
//! particles all follow the best position found so far (the global topology) can gather in
//! another well before a particle falls into the deepest; a ring topology spreads good positions
//! slowly and keeps exploring longer. The table counts the runs that end in each well. The
//! function, its bounds and its best known minimum come from genoxide's `problems::Shekel5`.
//!
//! 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 shekel5
//! ```

mod trace;

use genoxide::prelude::*;
use genoxide::problems::{Problem, Shekel5};

// the wells: Shekel's first points aᵢ
const WELLS: [[f64; 4]; 5] = [
    [4.0, 4.0, 4.0, 4.0],
    [1.0, 1.0, 1.0, 1.0],
    [8.0, 8.0, 8.0, 8.0],
    [6.0, 6.0, 6.0, 6.0],
    [3.0, 7.0, 3.0, 7.0],
];
const SEEDS: u64 = 30;
const BUDGET: u64 = 25_000;
const PARTICLES: usize = 80;
// a run stops once its error to the best known minimum is at most this
const ERROR: f64 = 1e-6;

fn main() -> Result<()> {
    let problem = Shekel5;
    let minimum = problem.optimum().expect("known").value();
    let target = minimum + ERROR;
    let topologies = [pso::Topology::Global, pso::Topology::Ring { neighbors: 1 }];
    // per topology: the runs that end in each well, those that reach the target, and their
    // evaluations
    let mut wells = [[0u64; WELLS.len()]; 2];
    let mut reached = [0u64; 2];
    let mut evaluations = [Vec::new(), Vec::new()];
    // with GENOXIDE_TRACE=<file>, a trace of the runs with the global topology for the plot on
    // the example's page
    let mut trace = trace::Trace::from_env();
    for (t, topology) in topologies.into_iter().enumerate() {
        for seed in 1..=SEEDS {
            let swarm = Pso::builder(problem.representation())
                .population_size(PARTICLES)
                .topology(topology)
                .minimize()
                .seed(seed)
                .build()?;
            let outcome = Engine::new(swarm, problem)
                .stop_when(Stop::target(target).or(Stop::evaluations(BUDGET)))
                .on_generation(|snapshot| {
                    if t == 0 {
                        trace.record(snapshot);
                    }
                })
                .run()?;
            wells[t][nearest(outcome.best_genome())] += 1;
            if outcome.stop_reason() == StopReason::Target {
                reached[t] += 1;
                evaluations[t].push(outcome.evaluations());
            }
        }
    }

    println!(
        "Shekel 5: best known minimum {minimum:.5}, {SEEDS} seeds, {BUDGET} evaluations at most \
         per run, {PARTICLES} particles"
    );
    println!("runs ending in the well at  PSO, global  PSO, ring");
    for (i, well) in WELLS.iter().enumerate() {
        if wells[0][i] + wells[1][i] == 0 {
            continue;
        }
        let at = format!("a{} = ({})", i + 1, coordinates(well));
        println!("{at:<26}  {:>11}  {:>9}", wells[0][i], wells[1][i]);
    }
    println!(
        "{:<26}  {:>11}  {:>9}",
        "error below 1e-6", reached[0], reached[1]
    );
    let [global, ring] = evaluations.map(|mut evaluations| median(&mut evaluations));
    println!("{:<26}  {global:>11}  {ring:>9}", "evaluations (median)");
    trace.write();
    Ok(())
}

// the index of the well nearest `x`
fn nearest(x: &Reals) -> usize {
    let distance = |well: &[f64; 4]| -> f64 { (0..4).map(|j| (x[j] - well[j]).powi(2)).sum() };
    (0..WELLS.len())
        .min_by(|&a, &b| distance(&WELLS[a]).total_cmp(&distance(&WELLS[b])))
        .expect("wells")
}

// a well's coordinates, e.g. 3, 7, 3, 7 or 7, 3.6, 7, 3.6
fn coordinates(well: &[f64; 4]) -> String {
    let text: Vec<String> = well.iter().map(|x| x.to_string()).collect();
    text.join(", ")
}

// the median of the evaluations of the runs that reach the target, rounded down; 0 without any
fn median(evaluations: &mut [u64]) -> u64 {
    evaluations.sort_unstable();
    let middle = evaluations.len() / 2;
    match evaluations.len() {
        0 => 0,
        n if n % 2 == 1 => evaluations[middle],
        _ => (evaluations[middle - 1] + evaluations[middle]) / 2,
    }
}

python examples/shekel5/main.py

"""Shekel 5: minimize Shekel's function with 5 wells in 4 dimensions with particle swarms, from 30
seeds, with a global and a ring topology.

The function has a well at each of 5 points, the deepest at (4, 4, 4, 4). A swarm whose particles
all follow the best position found so far (the global topology) can gather in another well before
a particle falls into the deepest; a ring topology spreads good positions slowly and keeps
exploring longer. The table counts the runs that end in each well. The function, its bounds and
its best known minimum come from genoxide's problems.Shekel5, which run evaluates in Rust.

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/shekel5/main.py
"""

import genoxide as gx

from trace import Trace

# the wells: Shekel's first points aᵢ
WELLS = [
    (4.0, 4.0, 4.0, 4.0),
    (1.0, 1.0, 1.0, 1.0),
    (8.0, 8.0, 8.0, 8.0),
    (6.0, 6.0, 6.0, 6.0),
    (3.0, 7.0, 3.0, 7.0),
]
SEEDS = 30
BUDGET = 25_000
PARTICLES = 80
# a run stops once its error to the best known minimum is at most this
ERROR = 1e-6


def nearest(x):
    """The index of the well nearest ``x``."""
    distances = [sum((x[j] - well[j]) ** 2 for j in range(4)) for well in WELLS]
    return distances.index(min(distances))


def coordinates(well):
    """A well's coordinates, e.g. 3, 7, 3, 7 or 7, 3.6, 7, 3.6."""
    return ", ".join(f"{x:g}" for x in well)


def median(evaluations):
    """The median of the evaluations of the runs that reach the target, rounded down; 0 without
    any."""
    evaluations = sorted(evaluations)
    middle = len(evaluations) // 2
    if not evaluations:
        return 0
    if len(evaluations) % 2:
        return evaluations[middle]
    return (evaluations[middle - 1] + evaluations[middle]) // 2


problem = gx.problems.Shekel5()
minimum = problem.optimum.value
target = minimum + ERROR
# per topology: the runs that end in each well, those that reach the target, and their evaluations
wells = [[0] * len(WELLS), [0] * len(WELLS)]
reached = [0, 0]
evaluations = [[], []]
# with GENOXIDE_TRACE=<file>, a trace of the runs with the global topology for the plot on the
# example's page
trace = Trace([problem.genome.bounds] * problem.genome.length, problem.optimum)
for t, ring in enumerate([None, 1]):
    for seed in range(1, SEEDS + 1):
        swarm = gx.Pso(
            problem.genome,
            population_size=PARTICLES,
            ring=ring,
            objective="minimize",
            seed=seed,
        )
        on_generation = trace.on_generation if t == 0 else None
        result = swarm.run(
            problem, target=target, evaluations=BUDGET, on_generation=on_generation
        )
        wells[t][nearest(result.best_genome)] += 1
        if result.stop_reason == "target":
            reached[t] += 1
            evaluations[t].append(result.evaluations)

print(
    f"Shekel 5: best known minimum {minimum:.5f}, {SEEDS} seeds, {BUDGET} evaluations at most "
    f"per run, {PARTICLES} particles"
)
print("runs ending in the well at  PSO, global  PSO, ring")
for i, well in enumerate(WELLS):
    if wells[0][i] + wells[1][i] == 0:
        continue
    at = f"a{i + 1} = ({coordinates(well)})"
    print(f"{at:<26}  {wells[0][i]:>11}  {wells[1][i]:>9}")
print(f"{'error below 1e-6':<26}  {reached[0]:>11}  {reached[1]:>9}")
print(f"{'evaluations (median)':<26}  {median(evaluations[0]):>11}  {median(evaluations[1]):>9}")
trace.write()

What it prints, from a seeded run:

Shekel 5: best known minimum -10.15320, 30 seeds, 25000 evaluations at most per run, 80 particles
runs ending in the well at  PSO, global  PSO, ring
a1 = (4, 4, 4, 4)                    14         30
a2 = (1, 1, 1, 1)                     3          0
a3 = (8, 8, 8, 8)                     8          0
a4 = (6, 6, 6, 6)                     4          0
a5 = (3, 7, 3, 7)                     1          0
error below 1e-6                     13         30
evaluations (median)               7040      11760