Метод кросс-энтропии: простейшая эвристика для сложнейших задач / forpes.ru

Главная
Метод кросс-энтропии: простейшая эвристика для сложнейших задач

Метод кросс-энтропии: простейшая эвристика для сложнейших задач +14

27.06.2025 10:01

uchitel 6 2200 Источник

Сколько эвристик вы знаете? Муравьи, отжиг, генетика, рой частиц, пчелы, светлячки, кукушки, гуси, совы, летучие мыши, осьминоги, дельфины, киты, шимпанзе, гориллы, львы, слоны, гравитация, электромагнетизм, вода, музыка… количество эвристик уже давно перевалило за полсотни. Все они вдохновлены природными процессами и явлениями, что делает их наглядными и понятными.

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

Самый минимальный теоретический минимум

Метод кросс-энтропии (далее CEM — cross-entropy method) вдохновлен лишь двумя статистическими принципами: минимизацией кросс-энтропии и сэмплированием по важности. Надо сразу отметить, что теория, которая стоит за этими принципами, строится на довольно зубодробительной математике. Однако, концептуально все достаточно просто. Хотя и немного контринтуитивно.

Кросс-энтропия

Если совсем по простому, то кросс-энтропия позволяет оценить разницу между двумя распределениями. Есть реальность, есть модель реальности, а кросс-энтропия — это просто мера, которая позволяет оценить насколько второе соответствует первому.

Давайте сделаем визуальный пример. Допустим, у нас есть выборка из «истинного» распределения вероятностей $inline$ , тогда благодаря минимизации кросс-энтропиии мы можем найти такое распределение $inline$ , которое будет максимально похоже на $inline$ .

Код для отрисовки

import numpy as np
from scipy.stats import norm, beta, uniform, bernoulli
from scipy.special import softmax
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

# Оценка кросс-энтропии методом Монте-Карло:
# Истинное (целевое) распределение p
a_p, b_p = 7, 3

# Модельное распределение q_1
a_q_1, b_q_1 = 3, 5

# Модельное распределение q_2
a_q_2, b_q_2 = 4, 1

# Модельное распределение q_3
a_q_3, b_q_3 = 6.5, 2.5

# Генерируем выборку из p
N = 10000
samples_p = beta.rvs(a=a_p, b=b_p, size=N)

# Вычисляем логарифм плотности q_1 в точках из выборки p
log_q_1 = beta.logpdf(samples_p, a=a_q_1, b=b_q_1)
log_q_2 = beta.logpdf(samples_p, a=a_q_2, b=b_q_2)
log_q_3 = beta.logpdf(samples_p, a=a_q_3, b=b_q_3)

# Оценка перекрёстной энтропии
h_p_q_1 = -np.mean(log_q_1)
h_p_q_2 = -np.mean(log_q_2)
h_p_q_3 = -np.mean(log_q_3)

print("Перекрёстная энтропия H(p, q):", h_p_q_1)

x = np.linspace(0, 1, 300)
p = beta.pdf(x, a=a_p, b=b_p)
q_1 = beta.pdf(x, a=a_q_1, b=b_q_1)
q_2 = beta.pdf(x, a=a_q_2, b=b_q_2)
q_3 = beta.pdf(x, a=a_q_3, b=b_q_3)

plt.figure(figsize=(10, 5), dpi=150)
plt.plot(x, p, label='p')
plt.plot(x, q_1, label=r'$q_{1}, H(p, q_{1})=$' +f'{h_p_q_1:.3}')
plt.plot(x, q_2, label=r'$q_{2}, H(p, q_{2})=$' +f'{h_p_q_2:.3}')
plt.plot(x, q_3, label=r'$q_{3}, H(p, q_{3})=$' +f'{h_p_q_3:.3}')

plt.legend()

plt.title(r'''Оценка перекрестной энтропии методом Монте-Карло
между $p$ и тремя распределениями: $q_{1}$, $q_{2}$, $q_{3}$''')
plt.show()

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

Сэмплирование по важности

Сэмплирование по важности — это метод снижения дисперсии, который можно использовать в методе Монте-Карло. Идея состоит в том, что при моделировании разные значения входных случайных величин оказывают разное влияние на оцениваемый параметр. То есть для нас некоторые значения оказываются «важнее», чем другие. Ну, а суть метода заключается в выборе таких распределений, которые «поощряют» важные значения.

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

Код для отрисовки

from numpy.random import randint

# "Длина" симуляции
N = 100
# Количество симуляций
n_iter = 3000

# Списки с результатами оценок:
data_uniform_gen = []
data_modified_gen = []
for _ in range(n_iter):
    # Получаем выборку из "равномерного" генератора
    samples_uniform_gen = randint(1, 100, size=N)
    # Оцениваем вероятность
    p_a_uniform_gen = np.cumsum(samples_uniform_gen == 1) / np.cumsum(samples_uniform_gen < 7)
    # Сохраняем результаты
    data_uniform_gen.append(p_a_uniform_gen)

    # Теперь все то же самое для "модифицированного" генератора
    samples_modified_gen = np.array([randint(1, 31) if np.random.rand() < 2/3 else randint(31, 100) for _ in range(N)])
    p_a_modified_gen = np.cumsum(samples_modified_gen == 1) / np.cumsum(samples_modified_gen < 7)
    data_modified_gen.append(p_a_modified_gen)

# Вычисляем отклонения
std_uniform_gen = np.nanstd(np.vstack(data_uniform_gen), axis=0)
std_modified_gen = np.nanstd(np.vstack(data_modified_gen), axis=0)

plt.figure(figsize=(10, 5), dpi=150)
plt.plot(std_uniform_gen, label='"Равномерный" генератор')
plt.plot(std_modified_gen, label='"Модифицированный" генератор')

plt.legend()

plt.ylabel(r'$\sigma$', rotation=90)
plt.xlabel('N')
plt.title("Отклонения оценок вероятности для двух генераторов")

plt.show()

Как видим, отклонение оценок у модифицированного генератора значительно меньше. Однако, возникает вполне закономерный вопрос, а почему бы нам сразу не подбрасывать игральную кость? А еще лучше, сразу вычислить в уме вероятность выпадения 1? Зачем такие танцы с бубном, если мы знаем, что речь идет об игральном кубике?

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

Другой причиной выделения «интересных» подмножеств может быть оптимизация некоторой целевой функции в виде черного ящика. Вероятность отыскания оптимального значения с помощью случайного перебора может быть невероятно мала, но если мы видим, что в одной области количество субоптимальных решений больше, чем в другой, то почему бы не сузить область поиска?

Сэмплирование по важности — это просто метод, который говорит нам, что для моделирования редких событий мы можем использовать более подходящие распределения. А благодаря кросс-энтропии — мы можем убедиться в том, что выбираемые распределения, действительно, являются наиболее подходящими.

Как все это работает?

Допустим у нас есть вот такая целевая функция:

$Z(x) = \mathrm{sinc}(0.5x) - 0.01x^{2} + 0.2x.$

Выглядит вот так:

Код для отрисовки

# Целевая функция
def objective(x):
    return np.sinc(0.5*x) - 0.01*x**2 + 0.2*x

x = np.linspace(-1.5, 20, 500)
z = objective(x)

plt.figure(figsize=(5, 3), dpi=150)
plt.plot(x, z)
plt.xlabel('x')
plt.ylabel('Z(x)')
plt.title(r'$Z(x) = \mathrm{sinc} (0.5 x) - 0.01 x^{2} + 0.2 x$')
plt.show()

Давайте найдем ее максимум на интервале $inline$ методом кросс-энтропии.

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

Шаги алгоритма следующие:

Инициализация параметров распределения. Например, непрерывное равномерное распределение $inline$ .
Генерация выборки из текущего распределения.
Оценка качества решений для всех значений из выборки.
Отбор подмножества лучших решений — «элиты».
Обновление параметров распределения на основе лучших решений.
Проверка условий остановки.

Давайте все это закодим.

Код и отрисовка

max_iterations = 10     # максимальное число итераций
sample_size = 50        # размер выборки
elite_size = 15         # количество лучших экземпляров выборки для обновления распределения

loc, scale = -5, 25     # Задаем начальные параметры распределения

elite_data = []

for i in range(max_iterations):
    # генерируем выборку
    samples = uniform.rvs(loc=loc, scale=scale, size=sample_size)
    # Оцениваем качество решений
    scores = objective(samples)
    # Отбор лучших решений
    elite_idx = np.argsort(scores)[-elite_size:]
    elite_samples = samples[elite_idx]
    # Сохраним элитные выборки что бы потом на них глянуть
    elite_data.append(elite_samples)
    #print(scores[-1])
    # Обновление параметров
    loc, scale = uniform.fit(elite_samples)
    # Смотрим что получается
    #print(f"Итерация {i+1:2d}: лучшее z(x) = {scores[-1]:.3f} при x = {elite_samples[-1]:.3f}")

plt.figure(figsize=(5, 3), dpi=150)

# Визуализируем
y_coords = np.linspace(0, 1, len(elite_data))
for i in range(len(elite_data)):
    plt.plot(elite_data[i], [y_coords[i]] * elite_size, 'b|' )

plt.plot(x, z)
plt.xlabel('x')
plt.ylabel('Z(x)')
plt.title(r'$Z(x) = \mathrm{sinc} (0.5 x) - 0.01 x^{2} + 0.2 x$')
plt.show()

Проще простого. Но где кросс-энтропия? И почему для обновления параметров мы используем максимизацию правдоподобия?

Дело вот в чем. Мы ввели параметрическое семейство распределений $inline$ , где $inline$ — это наша переменная, а $inline$ — это параметры. В нашем случае это непрерывное равномерное распределение. Далее, у нас есть набор «хороших» решений для целевой функции — элита, которые соответствуют некоторому подмножеству сгенерированной выборки: $\mathrm{elite} = [x_{1}, ..., x_{k}]$ . Формально CEM можно рассматривать как минимизацию кросс-энтропии между эмпирическим распределением элиты и нашей моделью:

$H(p_{\mathrm{elite}}(x), \; p(x, v)) = - \int p_{\mathrm{elite}}(x) \log p(x, v) dx.$

Минимизация этой величины равносильна максимизации правдоподобия, но только не по всем экземплярам выборки, а по ее элите. Опять же, чисто формально, если $p_{\mathrm{elite}}(x)$ и $inline$ принадлежат одному параметрическому семейству, то задача обновления параметров может быть сформулирована как:

$v^{t+1} = \underset{v}{\mathrm{argmax}} \sum_{x \in \mathrm{elite}} \log p(x, v),$

а это не что инное, как максимизация правдоподобия элиты в рамках модели $inline$ .

Но кого вообще интересуют эти формальности? Мы используем Монте-Карло, т.е. случайный перебор, у нас генерируется случайная выборка, в которой мы выделяем элиту. С помощью все того же метода Монте-Карло мы можем посчитать кросс-энтропию между эмпирическим распределением элиты и любым другим эмпирическим распределением, которое взбредет нам в голову. И в этом нет ничего «противозаконного», мы просто знаем, что есть метод сэплирования по важности, который поощряет нас за то, что мы выбираем распределения, которые поощряют сэмплирование важных для решения задачи решений.

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

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

Улучшайзинг алгоритма

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

Давайте возьмем задачу о рюкзаке 0-1 (все как обычно, но каждый предмет в единственном экземпляре), а данные вот отсюда:

Вместимость рюкзака: 6404180
Веса предметов: [382745, 799601, 909247, 729069, 467902, 44328, 34610, 698150, 823460, 903959, 853665, 551830, 610856, 670702, 488960, 951111, 323046, 446298, 931161, 31385, 496951, 264724, 224916, 169684]
Ценности предметов: [825594, 1677009, 1676628, 1523970, 943972, 97426, 69666, 1296457, 1679693, 1902996, 1844992, 1049289, 1252836, 1319836, 953277, 2067538, 675367, 853655, 1826027, 65731, 901489, 577243, 466257, 369261]
Оптимальное решение: [1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1]

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

Код cem_v1

# --- Входные данные ---
# Веса
weights = np.array([382745, 799601, 909247, 729069, 467902, 44328,
                    34610, 698150, 823460, 903959, 853665, 551830,
                    610856, 670702, 488960, 951111, 323046, 446298,
                    931161, 31385, 496951, 264724, 224916, 169684])
# Ценности
values = np.array([825594, 1677009, 1676628, 1523970, 943972, 97426,
                   69666, 1296457, 1679693, 1902996, 1844992, 1049289,
                   1252836, 1319836, 953277, 2067538, 675367, 853655,
                   1826027, 65731, 901489, 577243, 466257, 369261])

# Максимальная вместимость рюкзака
capacity = 6404180

def cem_v1(weights, values, capacity, printer):
    """
    Выполняет поиск оптимального решения задачи о рюкзаке 0-1. Решением
    является бинарный вектор в котором 0 - оставленный предмет,
    а 1 - выбранный предмет. Но данная функция возвращает только общие
    ценность и вес выбранных предметов, которые соответствуют найденному решению.
    Так же возвращается количество итераций, выполненных в процессе поиска.
    
    Args:
    weights (array) - веса предметов;
    values (array) - ценности предметов;
    capacity (int) - вместимость рюкзака;
    printer (bool) - если True, то выводит результат на экран и ничего не возвращает,
                     если False, то результат возвращается без вывода на экран.

    Returns:
    best_score - общая ценность выбранных предметов;
    best_weight - общий вес выбранных предметов;
    num_of_iters - количество выполненных итераций.
    """
    # --- Параметры CEM ---
    sample_size = 100       # количество сгенерированных решений за итерацию
    elite_size = 10         # размер элиты

    # --- Вспомогательные переменные ---
    flag = False            # флаг остановки алгоритма
    min_delta_p = 0.001     # минимальная дельта вероятности для остановки
    num_of_iters = 0        # переменная для подсчета количества итераций
    n_items = len(weights)  # Общее количество предметов
    
    # --- Инициализация ---
    # Вероятности того, что предмет будет взят (сначала все равны)
    p = np.full(n_items, 0.5)
    
    # --- Целевая функция с штрафом за перегрузку ---
    def objective(x):
        total_weight = np.sum(x * weights)
        if total_weight > capacity:
            return -1000  # штраф за перегрузку
        return np.sum(x * values)
    
    while not flag:
        
        # Генерация бинарных образцов. Генерируем с помощью binomial, что
        # удобно, если переменные не 0-1. Но в большинстве комбинаторных
        # задач лучше использовать choice.
        samples = np.random.binomial(n=1, p=p, size=(sample_size, n_items))
    
        # Вычисление оценок
        scores = np.array([objective(x) for x in samples])
    
        # Отбор элиты
        elite_indices = np.argsort(scores)[-elite_size:]  # максимизация
        elite_samples = samples[elite_indices]
    
        # Обновление вероятностей
        p_new = elite_samples.mean(axis=0)
    
        # Проверяем критерий остановки алгоритма
        flag = np.all(np.abs(p_new - p) < min_delta_p)
    
        p = p_new
        num_of_iters += 1
    
    # Результат
    best_vector = elite_samples[-1]    # Лучшее найденное решение
    best_score = scores[elite_indices[-1]]
    best_weight = np.sum(best_vector * weights)
    
    if printer == True:
        print(f"Number of iterations {num_of_iters}: Best score = {best_score}, Best weight = {best_weight}")
    else:
        return best_score, best_weight, num_of_iters

cem_v1(weights=weights, values=values, capacity=capacity, printer=True)

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

Код для отрисовки

# Вычисляем общие ценность и вес для оптимального решения
x_star = np.array([1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1])
total_value_star = sum(x_star * values)
total_weight_star = sum(x_star * weights)

# Соберем данные о 3000 результатах работы алгоритма функции cem_v1()
data_on_solutions_cem_v1 = np.array([
    cem_v1(weights=weights,
           values=values,
           capacity=capacity,
           printer=False)
    for _ in range(3000)
])

fig, ax = plt.subplots(1, 3, figsize=(12, 4), dpi=150)
sns.histplot(data_on_solutions_cem_v1[:, 0], ax=ax[0])
ax[0].axvline(total_value_star, c='k')
ax[0].set_title('Общая ценность')

sns.histplot(data_on_solutions_cem_v1[:, 1], ax=ax[1])
ax[1].axvline(total_weight_star, c='k')
ax[1].set_ylabel('')
ax[1].set_title('Общий вес')

sns.histplot(data_on_solutions_cem_v1[:, 2], discrete=True, ax=ax[2])
ax[2].set_ylabel('')
ax[2].set_title('Количество итераций');

Первое, что следует попробовать для улучшения результатов — это настройка объемов выборки и элиты. Но это как-то слишком просто и неинтересно.

Штрафование

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

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

Код для отрисовки

p = np.full(len(weights), 0.5)
samples = np.random.binomial(n=1, p=p, size=(5000, len(weights)))
total_weight_data = np.array([np.sum(x * weights) for x in samples])
total_value_data = np.array([np.sum(x * values) for x in samples])

plt.figure(figsize=(10, 5), dpi=150)
plt.scatter(total_weight_data, total_value_data, s=1, alpha=0.3)
plt.axvline(total_weight_star, c='r')
plt.axhline(total_value_star, c='r')

plt.xlabel('Вес выбранных предметов');
plt.ylabel('Ценность выбранных предметов')
plt.title('Зависимость общей ценности от общего веса выбранных предметов')

plt.show()

Очевидно, что иногда, особенно на первых итерациях, нам необходимо, чтобы очень хорошие, но недопустимые решения все же попадали в элиту. Следовательно награждать всех «негодяев» одинаковым штрафом — не лучшая идея. Лучше всего для этой цели использовать какую-нибудь более-менее адекватную штрафную функцию. Конкретно в нашем случае (по крайней мере на данной стадии реализации алгоритма) это не даст никакого прироста эффективности, но вообще, функциональное штрафование иногда может оказаться очень полезной фичей.

Взвешивание

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

def p_temper_softmax(x, T, len_x):
    """
    Оценивает вероятность появления единиц в каждой позиции вектора решений
    на основе данных из "элитной" выборки с помощью функции softmax с температурой.

    Args:
    x - выборка с элитными векторами;
    T - температура;
    len_x - объем элитной выборки.

    Returns:
    Массив с вероятностями появления единицы в соответствующих позициях
    вектора решений.
    """
    
    ones = x.sum(axis=0)
    x = np.vstack((len_x - ones, ones))
    return softmax(x/ T, axis=0)[-1]

Скорее кодим cem_v2, делаем температуру побольше и смотрим, что получается.

Код cem_v2 и отрисовки

def cem_v2(weights, values, capacity, printer):
    """
    Выполняет поиск оптимального решения задачи о рюкзаке 0-1. Решением
    является бинарный вектор в котором 0 - оставленный предмет,
    а 1 - выбранный предмет. Но данная функция возвращает только общие
    ценность и вес выбранных предметов, которые соответствуют найденному решению.
    Так же возвращается количество итераций, выполненных в процессе поиска.
    
    Args:
    weights (array) - веса предметов;
    values (array) - ценности предметов;
    capacity (int) - вместимость рюкзака;
    printer (bool) - если True, то выводит результат на экран и ничего не возвращает,
                     если False, то результат возвращается без вывода на экран.

    Returns:
    best_score - общая ценность выбранных предметов;
    best_weight - общий вес выбранных предметов;
    num_of_iters - количество выполненных итераций.
    """
    
    # --- Параметры CEM ---
    sample_size = 100       # количество сгенерированных решений за итерацию
    elite_size = 10         # размер элиты

    # --- Вспомогательные переменные ---
    flag = False            # флаг остановки алгоритма
    min_delta_p = 0.001     # минимальная дельта вероятности для остановки
    num_of_iters = 0        # переменная для подсчета количества итераций
    n_items = len(weights)  # Общее количество предметов
    
    # --- Инициализация ---
    # Вероятности того, что предмет будет взят (сначала все равны)
    p = np.full(n_items, 0.5)
    
    # --- Целевая функция с штрафом за перегрузку ---
    def objective(x):
        total_weight = np.sum(x * weights)
        if total_weight > capacity:
            return -1000  # штраф за перегрузку
        return np.sum(x * values)
    
    while not flag:
        
        # Генерация бинарных образцов
        samples = np.random.binomial(n=1, p=p, size=(sample_size, n_items))
    
        # Вычисление оценок
        scores = np.array([objective(x) for x in samples])
    
        # Отбор элиты
        elite_indices = np.argsort(scores)[-elite_size:]  # максимизация
        elite_samples = samples[elite_indices]
    
        # Обновление вероятностей
        p_new = p_temper_softmax(x=elite_samples, T=4, len_x=elite_size)
    
        # Проверяем критерий остановки алгоритма
        flag = np.all(np.abs(p_new - p) < min_delta_p)
    
        p = p_new
        num_of_iters += 1
        if num_of_iters > 50:
            break
    
    # Результат
    best_vector = elite_samples[-1]    # Лучшее найденное решение
    best_score = scores[elite_indices[-1]]
    best_weight = np.sum(best_vector * weights)
    
    if printer == True:
        print(f"Number of iterations {num_of_iters}: Best score = {best_score}, Best weight = {best_weight}")
    else:
        return best_score, best_weight, num_of_iters

# Соберем данные о 3000 результатах работы алгоритма функции cem_v1()
data_on_solutions_cem_v2 = np.array([
    cem_v2(weights=weights,
           values=values,
           capacity=capacity,
           printer=False)
    for _ in range(3000)
])

fig, ax = plt.subplots(1, 3, figsize=(12, 4), dpi=150)

sns.histplot(data_on_solutions_cem_v1[:, 0], ax=ax[0], element="step", fill=False, label='v1')
sns.histplot(data_on_solutions_cem_v2[:, 0], ax=ax[0], element="step", fill=False, label='v2')
ax[0].axvline(total_value_star, c='k')
ax[0].legend()
ax[0].set_title('Общая ценность')

sns.histplot(data_on_solutions_cem_v1[:, 1], ax=ax[1], element="step", fill=False, label='v1')
sns.histplot(data_on_solutions_cem_v2[:, 1], ax=ax[1], element="step", fill=False, label='v2')
ax[1].axvline(total_weight_star, c='k')
ax[1].set_ylabel('')
ax[1].legend()
ax[1].set_title('Общий вес')

sns.histplot(data_on_solutions_cem_v1[:, 2], discrete=True, ax=ax[2], element="step", fill=False, label='v1')
sns.histplot(data_on_solutions_cem_v2[:, 2], discrete=True, ax=ax[2], element="step", fill=False, label='v2')
ax[2].set_ylabel('')
ax[2].legend()
ax[2].set_title('Количество итераций');

Что ж, общая ценность немного подтянула хвост куда надо, а значит мы движемся в верном направлении. Правда количество итераций увеличилось — но кому это интересно. Двигаемся дальше!

Взвешивание по целевой функции

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

def p_weighted_score(elite_scores, elite_samples, size_samples, T):
    """
    Оценивает вероятности на основе взвешивания элитных решений по
    соответствующим значениям целевой функции.

    Args:
    elite_scores - значения целевой функции на соответствующие элите;
    elite_samples - элита;
    size_samples - размер элиты;
    T - температура.

    Returns:
    Массив с вероятностями появления единицы в соответствующих позициях
    вектора решений.
    
    """
    # Нормализуем значения целевой функции и  делим на температуру
    mu = np.mean(elite_scores)
    sigma = np.std(elite_scores) + 1    # прибавляем 1 что бы при sigma = 0 не получались значения nan
    elite_scores = (elite_scores - mu) / (T * sigma)

    # Вычисляем веса
    weights_softmax = softmax(elite_scores).reshape(size_samples, -1)

    # Оцениваем вероятность с помощью взвешенного среднего
    p = np.sum(weights_softmax * elite_samples, axis=0)
    p = np.clip(p, 0, 1)    # избавляемся от нанометрических выходов из [0, 1]
    return p

Кодим cem_v3 и с нетерпением ждем картинки:

Код cem_v3 и отрисовка

def cem_v3(weights, values, capacity, printer):
    """
    Выполняет поиск оптимального решения задачи о рюкзаке 0-1. Решением
    является бинарный вектор в котором 0 - оставленный предмет,
    а 1 - выбранный предмет. Но данная функция возвращает только общие
    ценность и вес выбранных предметов, которые соответствуют найденному решению.
    Так же возвращается количество итераций, выполненных в процессе поиска.
    
    Args:
    weights (array) - веса предметов;
    values (array) - ценности предметов;
    capacity (int) - вместимость рюкзака;
    printer (bool) - если True, то выводит результат на экран и ничего не возвращает,
                     если False, то результат возвращается без вывода на экран.

    Returns:
    best_score - общая ценность выбранных предметов;
    best_weight - общий вес выбранных предметов;
    num_of_iters - количество выполненных итераций.
    """
    # --- Параметры CEM ---
    sample_size = 100       # количество сгенерированных решений за итерацию
    elite_size = 10         # размер элиты

    # --- Вспомогательные переменные ---
    flag = False            # флаг остановки алгоритма
    min_delta_p = 0.001     # минимальная дельта вероятности для остановки
    num_of_iters = 0        # переменная для подсчета количества итераций
    n_items = len(weights)  # Общее количество предметов
    
    # --- Инициализация ---
    # Вероятности того, что предмет будет взят (сначала все равны)
    p = np.full(n_items, 0.5)
    
    # --- Целевая функция с штрафом за перегрузку ---
    def objective(x):
        total_weight = np.sum(x * weights)
        # Если общий вес больше вместимости, то вычисляем размер
        # штрафа с помощью линейной функции
        if total_weight > capacity:
            return 11000000  # штраф за перегрузку
        return np.sum(x * values)
    
    while not flag:
        
        # Генерация бинарных образцов
        samples = np.random.binomial(n=1, p=p, size=(sample_size, n_items))
    
        # Вычисление оценок
        scores = np.array([objective(x) for x in samples])
    
        # Отбор элиты
        elite_indices = np.argsort(scores)[-elite_size:]  # максимизация
        elite_scores = scores[elite_indices]
        elite_samples = samples[elite_indices]
    
        # Обновление вероятностей
        p_new = p_weighted_score(
            elite_scores=elite_scores,
            elite_samples=elite_samples,
            size_samples=elite_size,
            T=0.8
        )
        #print(elite_samples)
    
        # Проверяем критерий остановки алгоритма
        flag = np.all(np.abs(p_new - p) < min_delta_p)
    
        p = p_new
        num_of_iters += 1
    
    # Результат
    best_vector = elite_samples[-1]    # Лучшее найденное решение
    best_score = scores[elite_indices[-1]]
    best_weight = np.sum(best_vector * weights)
    
    if printer == True:
        print(f"Number of iterations {num_of_iters}: Best score = {best_score}, Best weight = {best_weight}")
    else:
        return best_score, best_weight, num_of_iters

# Соберем данные о 3000 результатах работы алгоритма функции cem_v1()
data_on_solutions_cem_v3 = np.array([
    cem_v3(weights=weights,
           values=values,
           capacity=capacity,
           printer=False)
    for _ in range(3000)
])

fig, ax = plt.subplots(1, 3, figsize=(12, 4), dpi=150)

sns.histplot(data_on_solutions_cem_v1[:, 0], ax=ax[0], element="step", fill=False, label='v1')
sns.histplot(data_on_solutions_cem_v2[:, 0], ax=ax[0], element="step", fill=False, label='v2')
sns.histplot(data_on_solutions_cem_v3[:, 0], ax=ax[0], element="step", fill=False, label='v3')
ax[0].axvline(total_value_star, c='k')
ax[0].legend()
ax[0].set_title('Общая ценность')

sns.histplot(data_on_solutions_cem_v1[:, 1], ax=ax[1], element="step", fill=False, label='v1')
sns.histplot(data_on_solutions_cem_v2[:, 1], ax=ax[1], element="step", fill=False, label='v2')
sns.histplot(data_on_solutions_cem_v3[:, 1], ax=ax[1], element="step", fill=False, label='v3')
ax[1].axvline(total_weight_star, c='k')
ax[1].set_ylabel('')
ax[1].legend()
ax[1].set_title('Общий вес')

sns.histplot(data_on_solutions_cem_v1[:, 2], discrete=True, ax=ax[2], element="step", fill=False, label='v1')
sns.histplot(data_on_solutions_cem_v2[:, 2], discrete=True, ax=ax[2], element="step", fill=False, label='v2')
sns.histplot(data_on_solutions_cem_v3[:, 2], discrete=True, ax=ax[2], element="step", fill=False, label='v3')
ax[2].set_ylabel('')
ax[2].legend()
ax[2].set_title('Количество итераций');

Как видим фича сработала на ура!

Банальное увеличение объема выборки и элиты

Теперь можно со спокойной душой ради любопытства взять и увеличить в cem_v3 объем выборки и элиты в 2 раза. Это, конечно, нельзя назвать какой-то прям серьезной модификацией, но мы закодим отдельно cem_v4 и посмотрим что получится:

Код cem_v4 и отрисовка

def cem_v4(weights, values, capacity, printer):
    """
    Выполняет поиск оптимального решения задачи о рюкзаке 0-1. Решением
    является бинарный вектор в котором 0 - оставленный предмет,
    а 1 - выбранный предмет. Но данная функция возвращает только общие
    ценность и вес выбранных предметов, которые соответствуют найденному решению.
    Так же возвращается количество итераций, выполненных в процессе поиска.
    
    Args:
    weights (array) - веса предметов;
    values (array) - ценности предметов;
    capacity (int) - вместимость рюкзака;
    printer (bool) - если True, то выводит результат на экран и ничего не возвращает,
                     если False, то результат возвращается без вывода на экран.

    Returns:
    best_score - общая ценность выбранных предметов;
    best_weight - общий вес выбранных предметов;
    num_of_iters - количество выполненных итераций.
    """
    # --- Параметры CEM ---
    sample_size = 200      # количество сгенерированных решений за итерацию
    elite_size = 20        # размер элиты

    # --- Вспомогательные переменные ---
    flag = False            # флаг остановки алгоритма
    min_delta_p = 0.001     # минимальная дельта вероятности для остановки
    num_of_iters = 0        # переменная для подсчета количества итераций
    n_items = len(weights)  # Общее количество предметов
    
    # --- Инициализация ---
    # Вероятности того, что предмет будет взят (сначала все равны)
    p = np.full(n_items, 0.5)
    
    # --- Целевая функция с штрафом за перегрузку ---
    def objective(x):
        total_weight = np.sum(x * weights)
        # Если общий вес больше вместимости, то вычисляем размер
        # штрафа с помощью линейной функции
        if total_weight > capacity:
            return 11000000  # штраф за перегрузку
        return np.sum(x * values)
    
    while not flag:
        
        # Генерация бинарных образцов
        samples = np.random.binomial(n=1, p=p, size=(sample_size, n_items))
    
        # Вычисление оценок
        scores = np.array([objective(x) for x in samples])
    
        # Отбор элиты
        elite_indices = np.argsort(scores)[-elite_size:]  # максимизация
        elite_scores = scores[elite_indices]
        elite_samples = samples[elite_indices]
    
        # Обновление вероятностей
        p_new = p_weighted_score(
            elite_scores=elite_scores,
            elite_samples=elite_samples,
            size_samples=elite_size,
            T=0.8
        )
        #print(p_new)
    
        # Проверяем критерий остановки алгоритма
        flag = np.all(np.abs(p_new - p) < min_delta_p)
    
        p = p_new
        num_of_iters += 1
    
    # Результат
    best_vector = elite_samples[-1]    # Лучшее найденное решение
    best_score = scores[elite_indices[-1]]
    best_weight = np.sum(best_vector * weights)
    
    if printer == True:
        print(f"Number of iterations {num_of_iters}: Best score = {best_score}, Best weight = {best_weight}")
    else:
        return best_score, best_weight, num_of_iters

# Соберем данные о 3000 результатах работы алгоритма функции cem_v1()
data_on_solutions_cem_v4 = np.array([
    cem_v4(weights=weights,
           values=values,
           capacity=capacity,
           printer=False)
    for _ in range(3000)
])

fig, ax = plt.subplots(1, 3, figsize=(12, 4), dpi=150)

sns.histplot(data_on_solutions_cem_v1[:, 0], ax=ax[0], element="step", fill=False, label='v1')
sns.histplot(data_on_solutions_cem_v2[:, 0], ax=ax[0], element="step", fill=False, label='v2')
sns.histplot(data_on_solutions_cem_v3[:, 0], ax=ax[0], element="step", fill=False, label='v3')
sns.histplot(data_on_solutions_cem_v4[:, 0], ax=ax[0], element="step", fill=False, label='v4')
ax[0].axvline(total_value_star, c='k')
ax[0].legend()
ax[0].set_title('Общая ценность')

sns.histplot(data_on_solutions_cem_v1[:, 1], ax=ax[1], element="step", fill=False, label='v1')
sns.histplot(data_on_solutions_cem_v2[:, 1], ax=ax[1], element="step", fill=False, label='v2')
sns.histplot(data_on_solutions_cem_v3[:, 1], ax=ax[1], element="step", fill=False, label='v3')
sns.histplot(data_on_solutions_cem_v4[:, 1], ax=ax[1], element="step", fill=False, label='v4')
ax[1].axvline(total_weight_star, c='k')
ax[1].set_ylabel('')
ax[1].legend()
ax[1].set_title('Общий вес')

sns.histplot(data_on_solutions_cem_v1[:, 2], discrete=True, ax=ax[2], element="step", fill=False, label='v1')
sns.histplot(data_on_solutions_cem_v2[:, 2], discrete=True, ax=ax[2], element="step", fill=False, label='v2')
sns.histplot(data_on_solutions_cem_v3[:, 2], discrete=True, ax=ax[2], element="step", fill=False, label='v3')
sns.histplot(data_on_solutions_cem_v4[:, 2], discrete=True, ax=ax[2], element="step", fill=False, label='v4')
ax[2].set_ylabel('')
ax[2].legend()
ax[2].set_title('Количество итераций');

Так может быть надо было с самого начала просто бахнуть по 1000 и 100 для объемов выборки и элиты? Если сделать так в cem_v3, то гарантия получения оптимума будет практически 100%-ой.

С одной стороны да, увеличение объемов в соответствии с возможностями сжигаемого железа — это первое, что надо делать, особенно если для заморочек нет времени. С другой стороны соотношение объемов, как и все вышеприведенные фичи, тоже по разному сказывается на результате. Поэтому, если для заморочек время все-таки есть, то лучше заморочиться.

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

В общем, все ограничевается воображением, а это очень круто, особенно если учесть простоту реализации.

Стохастическое программирование

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

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

В заключение

В общем-то на этом пока все. Жму F5 ожидая ваших комментариев.

Кстати, если вам интересна тема оптимизации, то обязательно подписывайтесь на (спокойно, не мой) канал NoML Community в телеге, там много полезной информации по OR (и не только), а еще там проводят интересные семинары.

Комментарии (6)

avshkol
27.06.2025 10:56
#28493778
Сколько эвристик вы знаете? Муравьи, отжиг, генетика, рой частиц, пчелы, светлячки, кукушки, гуси, совы, летучие мыши, осьминоги, дельфины, киты, шимпанзе, гориллы, львы, слоны, гравитация, электромагнетизм, вода, музыка…
После первой же фразы полез в llm за разъяснениями
Примеры указанных эвристик и краткое описание некоторых из них:

Муравьиный алгоритм — имитирует поведение муравьёв, которые метят путь феромонами, что позволяет находить кратчайший путь к цели. Используется для оптимизации маршрутов, например, задачи коммивояжёра1 2 3 5.

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

Генетические алгоритмы — основаны на принципах естественного отбора и генетики, используют операции скрещивания и мутации для поиска решений (не раскрыты в результатах, но известны как класс эвристик).

Рой частиц (Particle Swarm Optimization) — имитация поведения стай птиц или косяков рыб для поиска оптимальных решений (не в результатах, но популярная эвристика).

Пчелиные алгоритмы — моделируют поиск пищи пчёлами, применяются для оптимизации задач (не в результатах).

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

Гравитация, электромагнетизм, вода, музыка — эвристики, вдохновлённые физическими и природными процессами (например, гравитационные алгоритмы, алгоритмы, основанные на электромагнитных взаимодействиях, или музыкальные эвристики) — применяются для поиска решений в сложных задачах (не раскрыты в результатах).

В целом, эвристики — это методы и алгоритмы, которые упрощают поиск решений в сложных задачах, часто имитируя поведение живых организмов или природных явлений6 8.
1. uchitel Автор
  27.06.2025 10:56
  #28494086
  Не хотелось раздувать текст. Сделал ставку на то, что если слова "эвристика" и "муравьи" будут рядом, то в целом должно быть понятно о чем речь. Ну и на ваше любопытство я тоже ставил :)
  
  Что касается эл-эл-эмок, то уже больше года не могу понять какой контент теперь нужен как мне, так и читателям. После долгих раздумий пришел к выводу, что надо что-то вроде контента-путеводителя по темам. Теперь нет никакого смысла в доскональных и академически-строгих разъяснениях чего бы то ни было. Как считаете?
  1. avshkol
    27.06.2025 10:56
    #28495004
    Досконально разъяснения несложной темы очень заходят. Недавно была про электрообогреватели великолепная статья
    
    uchitel Автор
    27.06.2025 10:56
    #28495120
    Хороший науч-поп на Хабре всегда заходит.
    
    Хотя лично мне больше нравится, что статьи выскакивают в поиске, то есть это не только приятное времяпрепровождение, но еще и польза. А уж если попал из поиска и статья оказалась не только полезной, но еще и интересной, то от этого вдвойне приятнее.

VAE
27.06.2025 10:56
#28493842
>Улучшайзинг алгоритма
Это, конечно, круто...
1. uchitel Автор
  27.06.2025 10:56
  #28494112
  Но ведь нормальный улучшотрон получился :)

Метод кросс-энтропии: простейшая эвристика для сложнейших задач +14

Самый минимальный теоретический минимум

Кросс-энтропия

Сэмплирование по важности

Как все это работает?

Улучшайзинг алгоритма

Штрафование

Взвешивание

Взвешивание по целевой функции

Банальное увеличение объема выборки и элиты

Стохастическое программирование

В заключение

Комментарии (6)

avshkol

uchitel Автор

avshkol

uchitel Автор

VAE

uchitel Автор