Работа с алгоритмом искусственной пчелиной семьи



Большинство методов в ИИ соотв. Машинное обучение требует решения постоянной задачи оптимизации. Обычно эта задача оптимизации является нелинейной и не может быть решена вручную. Существует огромное количество литературы о видах оптимизационных задач и возможных методах их эффективного решения.

Основное различие между этими методами можно сделать за счет использования производных. Так, например, градиентные методы, как следует из названия, требуют вычисления производных. То же самое относится и к методу Ньютона, который предназначен для поиска корней первой производной.

Эти методы, как известно, являются быстрыми, если они работают. Но также известно, что иногда они вообще не работают или испытывают трудности с попаданием в ловушку локальных минимумов.

В этой истории я хочу представить технику оптимизации, которая сама основана на ИИ. Как и большинство алгоритмов на основе ИИ, он недетерминирован.

Алгоритм

Он называется «Алгоритм искусственной пчелиной колонии» и следует схеме, наблюдаемой в природе в отношении пчел:

Фиксированное количество N «пчел» размещается случайным образом (равномерно) в домене. Каждая пчела запоминает текущую позицию, а алгоритм запоминает лучшую глобальную позицию.

Локальный поиск. На каждом этапе каждая пчела сначала пытается найти лучшую позицию, случайным образом выбирая другую пчелу и перемещаясь на случайное расстояние по направлению к своей коллеге или от нее. Он выполняет фактический ход только в том случае, если новая позиция будет лучше, чем текущая.

Фаза наблюдателя: когда все пчелы выполнили локальный поиск, они перемещаются в соответствии с категориальным распределением относительно значения пригодности. Подробно, каждой пчеле присваивается пригодность путем вычисления расстояния от ее значения до значения текущей «худшей» пчелы. Таким образом, чем выше это расстояние, тем лучше фитнес. Каждое значение пригодности, которое мы рассматриваем как категорию, взвешивается путем деления значения пригодности на сумму всех значений пригодности. Это обеспечивает распределение (категориальное) всех значений пригодности, и мы можем брать из него случайные выборки. Обратите внимание, чем выше приспособленность, тем больше вероятность того, что образец будет выбран из соответствующей категории. Делая это N раз, что на самом деле способствует N разу выбора некоторой пчелы, мы находим новое место для каждой пчелы. Таким образом, чем лучше пригодность местоположения, тем больше пчел было выделено для него.

Фаза разведки: после этого все пчелы, которые не улучшили свое текущее значение, повторив вышеуказанные ходы X раз, случайным образом перемещаются по всему домену, так же, как и при инициализации.

Вышеописанное повторяется несколько раз, чтобы улучшить найденный глобальный оптимум.

Кроме того, этот метод очень стабилен, работает для больших размерностей, требует, чтобы функция была только непрерывной, и очень прост в реализации.

Более того, он позволяет вносить множество пользовательских улучшений. Трудно оценить его скорость сходимости, но во многих приложениях было проверено, что он особенно хорошо работает при поиске глобального минимума, не попадая в ловушку локальных минимумов.

Таким образом, в случаях, когда требуется очень точное решение этого глобального минимума, этот метод лучше всего использовать, позволяя ему приблизиться к этой точке, а затем использовать самый крутой градиент или итерацию Ньютона для получения желаемой точности.

Более подробную информацию, включая изобретателя, можно найти по адресу:

Алгоритм создания искусственной пчелиной семьи. (2022, 8 июля). В Википедии. https://en.wikipedia.org/wiki/Artificial_bee_colony_algori

Выполнение

Мы рассмотрим возможную реализацию на Rust. Для тех, кому нужно немного подвести итоги в Rust, рекомендую прочитать моё краткое введение.

Чтобы получить случайные значения, мы помещаем крейт rand в наши зависимости (Cargo.toml):

[dependencies]
rand = { version = "0.8.5" }

На протяжении всей программы мы используем следующие типы:

use rand::{distributions::WeightedIndex,
 prelude::Distribution, rngs::ThreadRng, thread_rng, Rng};

pub struct Params {
    n_bees: usize,
    abandon_limit: usize,
    max_iter: usize,
}

type X = Vec<f64>;
type Target = dyn Fn(&X) -> f64;

struct Bee {
    position: X,
    value: f64,
    not_improved_since: usize,
    fitness: f64, // distance to worst
}

Итак, Bee описывается своими position, value, fitness и тем, как часто ему не удавалось улучшить свои value. X — это просто удобный заполнитель для домена, а Target — это тип функции, которую мы хотим свести к минимуму.

Далее мы перечислим некоторые простые служебные функции, которые необходимы:

fn compute_fitness_on_bees(bees: &mut Vec<Bee>) {
    let worst_value = bees[find_worst_bee_idx(bees)].value;
    for bee in bees.iter_mut() {
        bee.fitness = worst_value - bee.value;
    }
}

// returns the best position and value among all the bees
fn find_best_position_and_value(bees: &Vec<Bee>) -> (X, f64) {
    let best_bee = bees
        .iter()
        .min_by(|b1, b2| b1.value.partial_cmp(&b2.value).unwrap())
        .unwrap();
    (best_bee.position.clone(), best_bee.value)
}

// returns the index of the bee with highest (worst) value
fn find_worst_bee_idx(bees: &Vec<Bee>) -> usize {
    bees.iter()
        .enumerate()
        .max_by(|(_, b1), (_, b2)| b1.value.partial_cmp(&b2.value).unwrap())
        .unwrap()
        .0
}

// creates a random point within the n-dim. interval [a, b],
// that is a </= x </= b
fn create_random_position(a: &X, b: &X, rng: &mut ThreadRng) -> X {
    let mut res = vec![];
    for (u, v) in a.iter().zip(b.iter()) {
        res.push(rng.gen_range(*u..*v));
    }
    res
}

Наконец, мы реализуем каждую фазу одну за другой, как описано в алгоритме. Вы увидите, что все шаги просты для реализации:

// create n bees and position them randomly in [a,b]
fn init_bees(f: &Target, a: &X, b: &X, n_bees: usize) -> Vec<Bee> {
    let mut res = vec![];
    let mut rng = thread_rng();
    for _ in 0..n_bees {
        let position = create_random_position(a, b, &mut rng);
        res.push(Bee {
            fitness: 0.0,
            not_improved_since: 0,
            value: f(&position),
            position,
        });
    }
    res
}

// let each bee locally search by moving a random distance to or away from
// another random bee.
// At the same time, keep the globally best know value/position up to date.
fn do_local_search_phase(
    f: &Target,
    bees: &mut Vec<Bee>,
    best_position: &mut X,
    best_value: &mut f64,
    a: &X,
    b: &X,
) {
    let mut rng = thread_rng();
    for idx in 0..bees.len() {
        let mut new_position = bees.get(idx).unwrap().position.clone();
        let other_bee_idx = rng.gen_range(0..bees.len());
        let position_idx = rng.gen_range(0..new_position.len());
        let phi = rng.gen_range(-1.0..1.0);
        let x_i_other = bees.get(other_bee_idx).unwrap()
                   .position[position_idx];
        let x_i = new_position[position_idx];
        new_position[position_idx] = (x_i + phi * (x_i_other - x_i))
            .max(a[position_idx])
            .min(b[position_idx]);
        let new_value = f(&new_position);
        let old_value = bees.get(idx).unwrap().value;
        let bee = bees.get_mut(idx).unwrap();
        if new_value < old_value {
            bee.position = new_position;
            bee.value = new_value;
            bee.not_improved_since = 0;
        } else {
            bee.not_improved_since += 1;
        }
        if new_value < *best_value {
            *best_value = bee.value;
            *best_position = bee.position.clone();
        }
    }
}

// Re-distribute bees onto other locations and prioritize locations with
// higher fitness. This is done by laying a categorical distrubition over all
// the bees.
fn do_onlooker_phase(bees: &mut Vec<Bee>) {
    compute_fitness_on_bees(bees);
    let sum_fitness = bees.iter().map(|bee| bee.value).sum::<f64>();
    let weights = bees
        .iter()
        .map(|bee| bee.value / sum_fitness)
        .collect::<X>();
    let categorical_dist = WeightedIndex::new(&weights).unwrap();
    let mut rng = thread_rng();
    for idx in 0..bees.len() {
        let other_bee = bees.get(categorical_dist.sample(&mut rng)).unwrap();
        let new_position = other_bee.position.clone();
        let new_value = other_bee.value;
        bees[idx].position = new_position;
        bees[idx].value = new_value;
    }
}

// Abandon bees from positions if have not been improved for long.
// This avoids being trapped at local minima.
fn do_scout_phase(bees: &mut Vec<Bee>, 
abandon_limit: usize, a: &X, b: &X, f: &Target) {
    let mut rng = thread_rng();
    bees.iter_mut()
        .filter(|bee| bee.not_improved_since > abandon_limit)
        .for_each(|bee| {
            let new_position = create_random_position(a, b, &mut rng);
            let new_value = f(&new_position);
            bee.position = new_position;
            bee.value = new_value;
            bee.not_improved_since = 0;
        });
}

Теперь у нас есть все вместе, чтобы объединить эти фазы в итерацию:

pub fn optimize(f: &Target, a: &X, b: &X, params: Params) -> (X, f64) {
    let mut bees = init_bees(f, a, b, params.n_bees);
    let (mut best_position, mut best_value) = find_best_position_and_value(&bees);
    for _ in 0..params.max_iter {
        do_local_search_phase(f, &mut bees, &mut best_position, &mut best_value, a, b);
        do_onlooker_phase(&mut bees);
        do_scout_phase(&mut bees, params.abandon_limit, a, b, f);
    }
    (best_position, best_value)
}

В качестве тестовой функции я использую функцию Экли (см. здесь), которая представляет собой функцию с довольно большим количеством локальных минимумов:

fn main(){
   let ackley_fn = |x: &X| {
            -20.0 * f64::exp(
               -0.2 * f64::sqrt(0.5 * (x[0] * x[0] + x[1] * x[1]))
             ) - f64::exp(
               0.5 * (f64::cos(2.0 * PI * x[0]) + f64::cos(2.0 * PI * x[1]))
             ) + E + 20.0
        };

        println!("{}", ackley_fn(&vec![0.0, 0.0])); // this is the minimum

        let params = Params {
            n_bees: 200,
            abandon_limit: 20,
            max_iter: 500,
        };
        println!(
            "{:?}",
            optimize(&ackley_fn, &vec![-5.0, -5.0], &vec![5.0, 5.0], params)
        );
}

Результат:

([0.007315634561002593, 9.547793830446227e-5], 0.02211822639045735)

Вы можете найти весь код для игры здесь.

Спасибо за прочтение!