Category: Алгоритмы

19
Янв
2021

Как муравьи решают проблемы коммивояжёров

Для решения задачи коммивояжёра используются разные алгоритмы, один из них называется «муравьиным». Разбираемся, что он из себя представляет.
— Читать дальше «Как муравьи решают проблемы коммивояжёров»

15
Янв
2021

Механика «Тетриса» легла в основу алгоритма «идеального» заполнения гостиниц

Разработчики уже пытаются получить патент на своё «изобретение».
— Читать дальше «Механика «Тетриса» легла в основу алгоритма «идеального» заполнения гостиниц»

04
Янв
2021

🐍 Анимация градиентного спуска и ландшафта функции потерь на Python

Демонстрация работающих примеров визуализации ландшафта функции потерь и анимации процесса градиентного спуска.

Статья публикуется в переводе, автор оригинального текста Tobias Roeschl.

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

Ландшафт функции потерь сверточной нейронной сети с 56 слоями (VGG-56, <a href="https://arxiv.org/abs/1712.09913" target="_blank" rel="noopener noreferrer nofollow">источник</a>)
Ландшафт функции потерь сверточной нейронной сети с 56 слоями (VGG-56, источник)

Приведенное выше изображение демонстрирует крайне невыпуклый ландшафт функции потерь нейронной сети. Ландшафт функции потерь – это визуальное представление значений, которые принимает целевая функция на наших тренировочных данных. Поскольку наша цель – визуализация потерь в трех измерениях, мы должны выбрать два параметра, которые будут изменяться, а остальные параметры останутся неизменными. Однако стоит заметить, что существуют продвинутые техники (например, сокращение размерности или filter response normalization), которые можно использовать для аппроксимации ландшафта функции потерь в подпространстве параметров с более низкой размерностью. Трехмерное представление функции потерь нейронной сети VGG с 56 слоями изображено на рис. 1. Однако это выходит за пределы данной статьи.

Искусственная нейронная сеть, над которой мы будем работать, состоит из одного входного слоя (с 784 узлами), двух скрытых слоев (с 50 и 500 узлами соответственно) и выходного слоя (с 10 узлами). Мы будем везде использовать сигмоиду в качестве функции активации. Эта нейронная сеть не будет содержать никаких порогов. Тренировочные данные состоят из изображений размером 28*28 пикселей с рукописными цифрами от 0 до 9 из набора данных MNIST. Технически, мы могли бы выбрать любые два веса из 784*50+50*500+500*10 = 69.200 весов, используемых в нашей сети, но я произвольно решил выбрать веса w250,5 (2) и w251,5(2), соединяющие 250-й и 251-й узлы второго скрытого слоя с шестым выходным нейроном, соответственно. Шестой выходной нейрон предсказывает, можно ли на данном изображении увидеть цифру ‘5’. Следующий рисунок схематически изображает архитектуру нейронной сети, с которой мы собираемся работать. Для простоты и понятности некоторые связи между нейронами – и большинство подписей с весами – были намеренно опущены.

Рис. 2. Архитектура нейронной сети (создана автором в <a href="https://app.diagrams.net/" target="_blank" rel="noopener noreferrer nofollow">draw.io</a>)
Рис. 2. Архитектура нейронной сети (создана автором в draw.io)

На Python’е мы сначала импортируем набор данных MNIST. Поскольку рукописные цифры в этом наборе представлены в виде изображений с градациями серого цвета, мы можем нормализовать наши входные данные, преобразовав значения пикселей из диапазона 0-255 к диапазону 0-1. разделив эти значения на 255.

        # Импортируем библиотеки
import numpy as np
import gzip
from sklearn.preprocessing import OneHotEncoder
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import expit
import celluloid
from celluloid import Camera
from matplotlib import animation 

# Открываем файлы MNIST: 
def open_images(filename):
    with gzip.open(filename, "rb") as file:     
        data=file.read()
        return np.frombuffer(data,dtype=np.uint8, offset=16).reshape(-1,28,28).astype(np.float32) 

def open_labels(filename):
    with gzip.open(filename,"rb") as file:
        data = file.read()
        return np.frombuffer(data,dtype=np.uint8, offset=8).astype(np.float32) 
    
X_train=open_images("C:\\Users\\tobia\\train-images-idx3-ubyte.gz").reshape(-1,784).astype(np.float32) 
X_train=X_train/255 # rescale pixel values to 0-1

y_train=open_labels("C:\\Users\\tobia\\train-labels-idx1-ubyte.gz")
oh=OneHotEncoder(categories='auto') 
y_train_oh=oh.fit_transform(y_train.reshape(-1,1)).toarray() # one-hot-encoding of y-values
    

Чтобы создать ландшафты потерь, мы собираемся нарисовать трехмерный график потерь по отношению к двум упомянутым выше весам w250_5(2) и w251_5(2). Чтобы это сделать, мы должны определить функцию потерь среднеквадратичной ошибки по отношению к весам w_a и w_b. Потери нашей модели равны средней сумме квадратичных ошибок между предсказаниями модели и истинными значениями для каждого из 10 выходных нейронов. Для N примеров эта функция будет выглядеть так:

J(y,pred)=1N∑i=1N(∑k=110(yi,k−predi,k)2)

Здесь y и pred представляют, соответственно, матрицы актуальных и предсказанных значений y. Предсказанные значения рассчитываются прямым проходом входных данных по нейронной сети до выходного слоя. Выход каждого слоя служит входом для следующего слоя. Входная матрица умножается на матрицу весов соответствующего слоя. После этого применяется функция сигмоиды для получения выхода данного конкретного слоя. Эти матрицы весов инициализируются небольшими случайными числами с помощью генератора псевдо-случайных чисел numpy. Используя инициализацию (seed), мы можем обеспечить воспроизводимость данных. После этого мы заменяем два веса, которым разрешено изменяться, аргументами w_a и w_b. На Python мы можем реализовать функцию потерь следующим образом:

        hidden_0=50 # количество узлов первого скрытого слоя 
hidden_1=500 # количество узлов второго скрытого слоя 

# Зададим функцию потерь:
def costs(x,y,w_a,w_b,seed_):  
        np.random.seed(seed_) # задаем инициализатор генератора случайных чисел 
        w0=np.random.randn(hidden_0,784)  # матрица весов 1-го скрытого слоя
        w1=np.random.randn(hidden_1,hidden_0) # матрица весов 2-го скрытого слоя
        w2=np.random.randn(10,hidden_1) # матрица весов выходного слоя
        w2[5][250] = w_a # задаем значение веса w_250,5(2)
        w2[5][251] = w_b # задаем значение веса w_251,5(2)
        a0 = expit(w0 @ x.T)  # выход 1-го скрытого слоя 
        a1=expit(w1 @ a0)  # выход 2-го скрытого слоя 
        pred= expit(w2 @ a1) # выход последнего слоя 
        return np.mean(np.sum((y-pred)**2,axis=0)) # потери w.r.t. w_a и w_b
    

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

        # Зададим набор значений для ячеистой сети: 
m1s = np.linspace(-15, 17, 40)   
m2s = np.linspace(-15, 18, 40)  
M1, M2 = np.meshgrid(m1s, m2s) # создаем ячеистую сеть  

# Определим потери для каждой координаты ячеистой сети: 
zs_100 = np.array([costs(X_train[0:100],y_train_oh[0:100].T  
                               ,np.array([[mp1]]), np.array([[mp2]]),135)  
                        for mp1, mp2 in zip(np.ravel(M1), np.ravel(M2))])
Z_100 = zs_100.reshape(M1.shape) # Значения z для N=100

zs_10000 = np.array([costs(X_train[0:10000],y_train_oh[0:10000].T  
                               ,np.array([[mp1]]), np.array([[mp2]]),135)  
                       for mp1, mp2 in zip(np.ravel(M1), np.ravel(M2))])
Z_10000 = zs_10000.reshape(M1.shape) # Значения z для N=10,000


# Рисуем ландшафты функции потерь: 
fig = plt.figure(figsize=(10,7.5)) # создаем фигуру
ax0 = fig.add_subplot(121, projection='3d' )
ax1 = fig.add_subplot(122, projection='3d' )

fontsize_=20 # задаем размер шрифта для меток осей 
labelsize_=12 # задаем размер меток делений 

# Настраиваем дочерние рисунки (subplots): 
ax0.view_init(elev=30, azim=-20)
ax0.set_xlabel(r'$w_a$', fontsize=fontsize_, labelpad=9)
ax0.set_ylabel(r'$w_b$', fontsize=fontsize_, labelpad=-5)
ax0.set_zlabel("costs", fontsize=fontsize_, labelpad=-30)
ax0.tick_params(axis='x', pad=5, which='major', labelsize=labelsize_)
ax0.tick_params(axis='y', pad=-5, which='major', labelsize=labelsize_)
ax0.tick_params(axis='z', pad=5, which='major', labelsize=labelsize_)
ax0.set_title('N:100',y=0.85,fontsize=15) # задаем заголовок дочернего рисунка

ax1.view_init(elev=30, azim=-30)
ax1.set_xlabel(r'$w_a$', fontsize=fontsize_, labelpad=9)
ax1.set_ylabel(r'$w_b$', fontsize=fontsize_, labelpad=-5)
ax1.set_zlabel("costs", fontsize=fontsize_, labelpad=-30)
ax1.tick_params(axis='y', pad=-5, which='major', labelsize=labelsize_)
ax1.tick_params(axis='x', pad=5, which='major', labelsize=labelsize_)
ax1.tick_params(axis='z', pad=5, which='major', labelsize=labelsize_)
ax1.set_title('N:10,000',y=0.85,fontsize=15)

# Поверхностные графики потерь (= ландшафты функции потерь):  
ax0.plot_surface(M1, M2, Z_100, cmap='terrain', #поверхностный график
                             antialiased=True,cstride=1,rstride=1, alpha=0.75)
ax1.plot_surface(M1, M2, Z_10000, cmap='terrain', #поверхностный график
                             antialiased=True,cstride=1,rstride=1, alpha=0.75)
plt.tight_layout()
plt.show()
    
Рис. 3. Ландшафты функции потерь для разных размеров тренировочного набора (изображение автора)
Рис. 3. Ландшафты функции потерь для разных размеров тренировочного набора (изображение автора)

Рисунок 3 изображает два примера ландшафтов функции потерь для одних и тех же весов (w250_5(2) и w251_5(2)) и одинаковых начальных значений всех весов. Левый ландшафт был построен на первых 100 изображениях из набора данных MNIST, а правый – на первых 10.000 изображениях. При ближайшем рассмотрении левого графика там можно увидеть типичные особенности невыпуклых поверхностей: локальные минимумы, плато, хребты (иногда называемые “седловыми точками”) и “глобальный” минимум. Однако термин “минимум” следует использовать осторожно, поскольку мы видели лишь небольшой набор значений, и не проводили анализ производной.

Рис. 4 (создан автором с помощью <a href="https://app.diagrams.net/" target="_blank" rel="noopener noreferrer nofollow">draw.io</a>)
Рис. 4 (создан автором с помощью draw.io)

Градиентный спуск

Эти “географические барьеры” – полная противоположность выпуклым и гладким ландшафтам функции потерь, которые можно увидеть для линейной регрессии или логистической регрессии. Считается, что эти “барьеры” замедляют или даже препятствуют достижению глобального минимума методом градиентного спуска, и, следовательно, оказывают негативное влияние на качество модели. Чтобы исследовать это явление, я решил анимировать градиентный спуск для данного ландшафта функции потерь из трех различных “стартовых позиций”. Градиентный спуск, в целом, заключается в обновлении параметров модели (в данном случае – весов) в соответствии со следующим уравнением:

we+1=we−αΔJ(we)

Здесь Delta J представляет градиент функции потерь, w – веса всей модели, e представляет текущую эпоху обучения, а aplha – скорость обучения.

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

∂J∂wij(2)=∂J∂outi(2)∂outi(2)∂ini(2)∂ini(2)∂wij(2)

Здесь wij определяется как вес между j-м узлом предыдущего слоя и i-м узлом текущего слоя, которым в нашем случае является выходной слой. Вход i-го нейрона выходного слоя обозначен просто как ini(2), что эквивалентно сумме активаций предыдущего слоя, умноженных на соответствующие веса связей, ведущих к этому узлу. Выход i-го нейрона выходного слоя обозначен как outi(2), что соответствует сигмоиде(ini(2)). Решая приведенное выше уравнение, мы получаем:

∂J∂wij(2)=(outi(2)−targeti)⋅outi(2)⋅(1−outi(2))⋅outj(1)

Здесь outj(1) соответствует активации j-го узла предыдущего слоя, соединенного с i-м нейроном выходного слоя связью с весом wij. Переменная targeti соответствует целевому выводу каждого из 10 нейронов выходного слоя. Возвращаясь к Рис. 2, outj(1) будет соответствовать активации h250 или h251, в зависимости от веса, по которому мы хотим рассчитать частную производную. Прекрасное объяснение обратного распространения (backpropagation), включая математические детали дифференцирования, можно найти здесь.

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

        # Сохраняем значения потерь и весов в списках: 
weights_2_5_250=[] 
weights_2_5_251=[] 
costs=[] 

seed_= 135 # для инициализации генератора случайных чисел
N=100 # размер выборки 

# Задаем нейронную сеть: 
class NeuralNetwork(object):
    def __init__(self, lr=0.01):
        self.lr=lr
        np.random.seed(seed_) # Инициализируем генератор случайных чисел
        # Инициализируем матрицы весов: 
        self.w0=np.random.randn(hidden_0,784)  
        self.w1=np.random.randn(hidden_1,hidden_0)
        self.w2=np.random.randn(10,hidden_1)
        self.w2[5][250] = start_a # задать стартовое значение для w_a
        self.w2[5][251] = start_b # задать стартовое значение для w_b
    
    def train(self, X,y):
        a0 = expit(self.w0 @ X.T)  
        a1=expit(self.w1 @ a0)  
        pred= expit(self.w2 @ a1)
        # Частные производные функции потерь w.r.t. весов выходного слоя: 
        dw2= (pred - y.T)*pred*(1-pred)  @ a1.T / len(X)   # ... среднее по выборке
        # Обновляем веса: 
        self.w2[5][250]=self.w2[5][250] - self.lr * dw2[5][250] 
        self.w2[5][251]=self.w2[5][251] - self.lr * dw2[5][251] 
        costs.append(self.cost(pred,y)) # дописываем значения потерь в список
    
    def cost(self, pred, y):
        return np.mean(np.sum((y.T-pred)**2,axis=0))
    
# Начальные значения w_a/w_b: 
starting_points = [  (-9,15),(-10.1,15),(-11,15)] 

for j in starting_points:
    start_a,start_b=j
    model=NeuralNetwork(10) # установить скорость обучения в 10
    for i in range(10000):  # 10,000 эпох           
        model.train(X_train[0:N], y_train_oh[0:N]) 
        weights_2_5_250.append(model.w2[5][250]) # дописываем значения весов в список
        weights_2_5_251.append(model.w2[5][251]) # дописываем значения весов в список

# Create sublists of costs and weight values for each starting point: 
costs = np.split(np.array(costs),3) 
weights_2_5_250 = np.split(np.array(weights_2_5_250),3)
weights_2_5_251 = np.split(np.array(weights_2_5_251),3)
    

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

        fig = plt.figure(figsize=(10,10)) # создаем фигуру
ax = fig.add_subplot(111,projection='3d' ) 
line_style=["dashed", "dashdot", "dotted"] # стили линий
fontsize_=27 # задаем размер шрифта для меток осей 
labelsize_=17 # задаем размер шрифта для меток делений
ax.view_init(elev=30, azim=-10)
ax.set_xlabel(r'$w_a$', fontsize=fontsize_, labelpad=17)
ax.set_ylabel(r'$w_b$', fontsize=fontsize_, labelpad=5)
ax.set_zlabel("costs", fontsize=fontsize_, labelpad=-35)
ax.tick_params(axis='x', pad=12, which='major', labelsize=labelsize_)
ax.tick_params(axis='y', pad=0, which='major', labelsize=labelsize_)
ax.tick_params(axis='z', pad=8, which='major', labelsize=labelsize_)
ax.set_zlim(4.75,4.802) # задаем диапазон значений по z в графике

# Определяем, график каких эпох рисовать:
p1=list(np.arange(0,200,20))
p2=list(np.arange(200,9000,100))
points_=p1+p2

camera=Camera(fig) # создаем объект Camera
for i in points_:
    # Рисуем три траектории градиентного спуска...
    #... каждая начинается из своей начальной точки
    #... и имеет собственный уникальный стиль линии:
    for j in range(3): 
        ax.plot(weights_2_5_250[j][0:i],weights_2_5_251[j][0:i],costs[j][0:i],
                linestyle=line_style[j],linewidth=2,
                color="black", label=str(i))
        ax.scatter(weights_2_5_250[j][i],weights_2_5_251[j][i],costs[j][i],
                   marker='o', s=15**2,
               color="black", alpha=1.0)
    # Surface plot (= loss landscape):
    ax.plot_surface(M1, M2, Z_100, cmap='terrain', 
                             antialiased=True,cstride=1,rstride=1, alpha=0.75)
    ax.legend([f'epochs: {i}'], loc=(0.25, 0.8),fontsize=17) # позиция надписи
    plt.tight_layout() 
    camera.snap() # сделать снимок после каждой итерации
    
animation = camera.animate(interval = 5, # интервал между фреймами в миллисекундах
                          repeat = False,
                          repeat_delay = 0)
animation.save('gd_1.gif', writer = 'imagemagick', dpi=100)  # сохраним анимацию  
    
Рис. 5. Траектории градиентного спуска (создано автором)
Рис. 5. Траектории градиентного спуска (создано автором)

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

        fig = plt.figure(figsize=(10,10)) # создаем фигуру
ax0=fig.add_subplot(2, 1, 1) 
ax1=fig.add_subplot(2, 1, 2) 

# Задаем параметры дочерних рисунков (subplots): 
ax0.set_xlabel(r'$w_a$', fontsize=25, labelpad=0)
ax0.set_ylabel(r'$w_b$', fontsize=25, labelpad=-20)
ax0.tick_params(axis='both', which='major', labelsize=17)
ax1.set_xlabel("epochs", fontsize=22, labelpad=5)
ax1.set_ylabel("costs", fontsize=25, labelpad=7)
ax1.tick_params(axis='both', which='major', labelsize=17)

contours_=21 # задаем количество контурных линий
points_=np.arange(0,9000,100) # определяем, какие эпохи рисовать

camera = Camera(fig) # создаем объект Camera
for i in points_:
    cf=ax0.contour(M1, M2, Z_100,contours_, colors='black', # контурный график
                     linestyles='dashed', linewidths=1)
    ax0.contourf(M1, M2, Z_100, alpha=0.85,cmap='terrain') # контурные графики с заливкой 
    
    for j in range(3):
        ax0.scatter(weights_2_5_250[j][i],weights_2_5_251[j][i],marker='o', s=13**2,
               color="black", alpha=1.0)
        ax0.plot(weights_2_5_250[j][0:i],weights_2_5_251[j][0:i],
                linestyle=line_style[j],linewidth=2,
                color="black", label=str(i))
        
        ax1.plot(costs[j][0:i], color="black", linestyle=line_style[j])
    plt.tight_layout()
    camera.snap()
    
animation = camera.animate(interval = 5,
                          repeat = True, repeat_delay = 0)  # создаем анимацию 
animation.save('gd_2.gif', writer = 'imagemagick')  # сохраняем анимацию в gif
    
Рис.6. Траектории градиентного спуска в 2D (создано автором)
Рис.6. Траектории градиентного спуска в 2D (создано автором)

Обе анимации иллюстрируют, что при невыпуклом ландшафте функции потерь градиентный спуск может “зависнуть” на локальном минимуме, седловой точке или плато. Для преодоления некоторых из этих проблем было изобретено множество вариантов градиентного спуска (ADAGRAD, Adam и т.д.) Однако я хочу разъяснить, что не все ландшафты функции потерь сильно невыпуклые в заданном диапазоне значений w_a и w_b. Выпуклость функции потерь зависит, среди прочего, от количества скрытых слоев – у глубоких нейронных сетей ландшафт функции потерь обычно сильно невыпуклый.

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

Рис. 7. N=500, w200-30(1), w200-31(1). Создано автором (<a href="https://gist.github.com/dbc60bb7c6ef0b8ae3f8ceab8c334b57.git" target="_blank" rel="noopener noreferrer nofollow">код</a>)
Рис. 7. N=500, w200-30(1), w200-31(1). Создано автором (код)
Рис. 8 N=1000, w5_5(1), w5_6(1). Создано автором (<a href="https://gist.github.com/4f87869ae3e94eb6331fb9e213e9f343.git" target="_blank" rel="noopener noreferrer nofollow">код</a>).
Рис. 8 N=1000, w5_5(1), w5_6(1). Создано автором (код).

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

Полный блокнот с исходными текстами визуализации можно найти в моем GitHub’е.

Ссылки

База рукописных цифр MNIST (Янн ЛеКун, Коринна Кортес и Крис Буржс).

  1. Ли Хао и др. “Визуализация ландшафтов функции потерь нейронных сетей”. Достижения нейронных систем обработки информации, 2018.
  2. https://machinelearningmastery.com/how-to-normalize-center-and-standardize-images-with-the-imagedatagenerator-in-keras/
  3. https://machinelearningmastery.com/why-training-a-neural-network-is-hard/
  4. https://mattmazur.com/2015/03/17/a-step-by-step-backpropagation-example/
  5. Мэттью Стэйб, Сашанк Дж. Редди, Сэтьен Кэйл, Санджив Кумар, Суврит Сра. “Избежание седловых точек адаптивными градиентными методами” (2019).
  6. Ян Дофин и др. “Обнаружение и уничтожение проблемы седловых точек при невыпуклой оптимизации в пространстве с большим количеством измерений”. NIPS (2014).
  7. А. Хороманска, М. Хенафф, М. Матье, Г.Б. Аруз и Я. ЛеКун. “Ландшафты функции потерь нейронных сетей со множеством слоев”. Журнал исследований машинного обучения, 38, 192-204.

Приложение

Демонстрационное изображение (создано автором)
Демонстрационное изображение (создано автором)
25
Дек
2020

Как работает лифт в небоскребах? Алгоритмы + задачи с собеседований

Алгоритм по которому работает лифт в высотном здании должен учитывать множество факторов. Он сложнее чем у обычного лифта.
— Читать дальше «Как работает лифт в небоскребах? Алгоритмы + задачи с собеседований»

08
Дек
2020

На GitHub опубликовали алгоритм для расшифровки «пикселизированных» изображений

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

30
Ноя
2020

Хакатон Code Battle Online: Tetris

Игровой мини-хакатон для начинающих IT-специалистов с базовыми навыками в C#, Java, JavaScript, Python или C++.
— Читать дальше «Хакатон Code Battle Online: Tetris»

24
Ноя
2020

Стоит прочитать: обзор книги Кормена и Лейзерсона «Алгоритмы. Построение и анализ»

В книге охватывается основной спектр современных алгоритмов: сортировки, графовые алгоритмы, динамическое программирование и тому подобное.
— Читать дальше «Стоит прочитать: обзор книги Кормена и Лейзерсона «Алгоритмы. Построение и анализ»»

03
Ноя
2020

🐍Сложность алгоритмов и операций на примере Python

Определить вычислительную сложность отдельных операций просто, но как вычислить сложность целой функции? Попробуем ответить на этот вопрос в небольшой статье.

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

Такой тип анализа называется статическим, так как не требует запуска программы, – в противовес динамическому, или эмпирическому, анализу, при котором проводятся измерения параметров работающего кода.

Классы сложности операций в Python

Самое базовое понятие статического анализа – O(1). Этот класс сложности имеют операции, которые выполняются за константное время, например, создание переменной или сложение небольших чисел.

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

Параметр N, который вы встретите далее, – это размер структуры данных – len(data_structure).

Рассмотрим основные операции некоторых структур данных в Python.

Списки (lists)

Операция Пример Сложность Примечания
Получение элемента l[i] O(1)
Сохранение элемента l[i] = 0 O(1)
Размер списка len(l) O(1)
Добавление элемента в конец списка l.append(5) O(1)
Удаление последнего элемента (pop) l.pop() O(1) То же, что и l.pop(-1)
Очищение списка l.clear() O(1) То же, что и l = []
Получение среза l[a:b] O(b-a) l[1:5] => O(1), l[:] => O(len(l) – 0) = O(N)
Расширение l.extend(...) O(len(…)) Зависит от размера расширения
Создание list(...) O(len(…)) Зависит от размера итерируемой структуры (…)
Сравнение списков (==, !=) l1 == l2 O(N)
Вставка l[a:b] = ... O(N)
Удаление элемента (del) del l[i] O(N) Зависит от i. O(N) – в худшем случае
Проверка наличия x in/not in l O(N) Линейный поиск в списке
Копирование l.copy() O(N) То же, что и l[:]
Удаление значения (remove) l.remove(...) O(N)
Удаление элемента (pop) l.pop(i) O(N) O(N-i). Для l.pop(0) => O(N)
Получение минимального/максимального значения min(l)/max(l) O(N) Линейный поиск в списке
Разворачивание списка l.reverse() O(N)
Перебор for v in l: O(N) В худшем случае, без прерывания цикла (return/break)
Сортировка l.sort() O(N Log N)
Умножение k*l O(k N) 5*l => O(N), len(l)*l => O(N2)

Кортежи (tuples) поддерживают все операции, которые не изменяют структуру данных – и они имеют такие же классы сложности, как у списков.

Множества (sets)

Операция Пример Сложность Примечания
Размер множества len(s) O(1)
Добавление элемента s.add(5) O(1)
Проверка наличия значения x in/not in s O(1) Для списков и кортежей => O(N)
Удаление значения (remove) s.remove(..) O(1) Для списков и кортежей => O(N)
Удаление значения (discard) s.discard(..) O(1)
Удаление значения (pop) s.pop() O(1) Удаляемое значение выбирается “рандомно”
Очищение множества s.clear() O(1) То же, что и s = set()
Создание set(...) O(len(…)) Зависит от размера итерируемой структуры (…)
Сравнение множеств (==, !=) s != t O(len(s)) То же, что и len(t)
Сравнение множеств (<=/<) s <= t O(len(s)) issubset
Сравнение множеств (>=/>) s >= t O(len(t)) issuperset s <= t == t >= s
Объединение (union) s | t O(len(s)+len(t))
Пересечение (intersection) s & t O(len(s)+len(t))
Разность (difference) s – t O(len(s)+len(t))
Симметричная разность s ^ t O(len(s)+len(t))
Перебор множества for v in s: O(N) В худшем случае, без прерывания цикла (return/break)
Копирование s.copy() O(N)

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

Неизменяемые множества (frozen sets) поддерживают все операции, которые не изменяют структуру данных – и они имеют такие же классы сложности, как у обычных множеств.

Словари (dict и defaultdict)

Операция Пример Сложность Примечания
Получение элемента d[k] O(1)
Сохранение элемента d[k] = v O(1)
Размер словаря len(d) O(1)
Удаление элемента (del) del d[k] O(1)
get/setdefault d.get(k) O(1)
Удаление (pop) d.pop(k) O(1)
Удаление (popitem) d.popitem() O(1) Удаляемое значение выбирается “рандомно”
Очищение словаря d.clear() O(1) То же, что и s = {} или s = dict()
Получение ключей d.keys() O(1) То же для d.values()
Создание словаря dict(...) O(len(…))
Перебор элементов for k in d: O(N) Для всех типов: keys, values, items. В худшем случае, без прерывания цикла

Как видно, большинство операций со словарями имеют сложность O(1).

Тип defaultdict поддерживает все эти операции с теми же классами сложности. Таким образом, вызов конструктора в том случае, если значения не найдены в defaultdict, имеет сложность O(1).

Тонкости анализа

Обратите внимание, что операция for i in range(...) имеет сложность O(len(...)). Для for i in range(1, 10) она равна O(1).

Если len(alist) – это N, тогда for i in range(len(alist)) будет иметь сложность O(N), так как цикл выполняется N раз.

При этом for i in range(len(alist)/2) также имеет сложность O(N). В этом случае цикл выполняется N/2 раз, и мы можем отбросить константу 1/2. При увеличении размера списка вдвое выполняемая работа также удваивается.

Точно так же for i in range(len(alist)/1000000) имеет сложность O(N). Это важно понять, так как нас интересует, что происходит, когда N стремится к бесконечности.

***

При сравнении двух списков на равенство, класс сложности должен быть O(N), как указано в таблице выше. Однако в реальности это значение нужно умножить на O==(…), где O==(…) это класс сложности для операции сравнения (==) двух значений в списке. Если мы работаем с целыми числами (int), то сложность сравнения будет равна O(1), если со строками (string), то в худшем случае мы получим O(len(string)).

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

Составные классы сложности

Разобравшись со сложностью отдельных операций мы переходим к их комбинированию.

Закон сложения для O-нотации

O(f(n))+O(g(n))=O(f(n)+g(n))

При сложении двух классов сложности складываются функции этих классов. В конечном счете, O(f(n) + g(n)) приводит нас к большему из двух исходных классов – так как меньший член выражения просто отбрасывается.

Таким образом,

O(N)+O(log⁡N)=O(N+log⁡N)=O(N)

так какN растет быстрее, чем log N:

limx→∞log⁡NN=0

Это правило позволяет вычислить класс сложности для последовательности операций. Например, сначала мы выполняем выражение со сложностью O(f(n)), а следом за ним – O(g(n)). Тогда выполнение обоих выражений (одно за другим) будет иметь сложность O(f(n)) + O(g(n)), то есть O(f(n) + g(n)).

Пример

Вызов функции f(...) имеет сложность O(N), а вызов g(...)O (N * log N). Вызываем эти функции друг за другом:

        f(...)
g(...)
    

Сложность этого действия будет равна:

O(N)+O(Nlog⁡N)=O(N+Nlog⁡N)=O(Nlog⁡N)

Если мы вызовем функцию f(...) дважды:

        f(...)
f(...)
    

то получим:

O(N)+O(N)=O(N+N)=O(2N)=O(N)

Константа 2 в вычислениях отбрасывается.

Условия

Отдельно разберем исполнение условий (if). Сначала выполняется само условие, а затем один из блоков (if или else).

        if test:
  block 1
else:
  block 2
    

Предположим, что вычисление условия test имеет сложность O(T), блока block 1O(B1), а блока block 2O(B2).

Тогда сложность всего кода будет равна:

O(T)+max(O(B1),O(B2))

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

Подставим реальные значения:

  • testO(N),
  • block 1O(N2),
  • block 2O(N).

и вычислим сложность кода:

O(N)+max(O(N2),O(N))=
O(N)+O(N2)=O(N+N2)=O(N2)

Если бы операция test имела класс сложности O(N3), то общая сложность кода составила бы O(N3).

O(N3)+max(O(N2),O(N))=
O(N3)+O(N2)=O(N3+N2)=O(N3)

Фактически, общий класс сложности для if-условий можно записать еще проще:

O(T)+O(B1)+O(B2)

Для первого примера в этом случае получим:

O(N)+O(N2)+O(N)=O(N2)

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

Закон умножения для O-нотации

O(f(n))∗O(g(n))=O(f(n)∗g(n))

Если мы повторяем O(N) раз некоторый процесс с классом сложности O(f(N)), то результирующий класс сложности будет равен:

O(N)×O(f(N))=O(N×f(N))

Предположим, некоторая функция f(...) имеет класс сложности O(N2). Выполним ее в цикле N раз:

        for i in range(N):
    f(...)
    

Сложность этого кода будет равна:

O(N)×O(N2)=O(N×N2)=O(N3)

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

Статический анализ на практике

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

Для всех трех примеров размер списка обозначим как N, а сложность операции сравнения элементов примем за O(1).

Алгоритм 1

Список уникален, если каждое его значение не встречается ни в одном последующем индексе. Для значения alist[i] последующим фрагментом списка будет срез alist[i+1:].

        def is_unique1 (alist : [int]) -> bool:
  for i in range(len(alist)):             # 1
    if alist[i] in alist[i+1:]:           # 2
      return False                        # 3
  return True                             # 4
    

Определим сложность для каждой строки метода:

  1. O(N) – для каждого индекса. Создание объекта range требует выполнения трех последовательных операций: вычисления аргументов, передачу их в __init__ и выполнение тела __init__. Две последние имеют класс сложности O(1). Сложность len(alist) также O(1), поэтому общая сложность выражения range(len(alist)) – O(1) + O(1) + O(1) = O(1).
  2. O(N) – получение индекса + сложение + создание среза + проверка in: O(1) + O(1) + O(N) + O(N) = O(N)
  3. O(1) – в худшем случае никогда не выполняется, можно проигнорировать.
  4. O(1) – в худшем случае всегда выполняется.

Таким образом, класс сложности целой функции is_unique1 равен:

O(N)×O(N)+O(1)=O(N2)

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

***

Возможно, вы хотели написать так:

O(N)×(O(N)+O(1))+O(1)

ведь выражение if cостоит из вычисления самого условия (O(N)) и блока return False (O(1)). Но в худшем случае этот блок никогда не будет выполнен и цикл продолжится, поэтому мы не включаем его в формулу. Но даже если добавить его, то ничего не изменится, так как O(N) + O(1) – это по-прежнему O(N).

Кроме того, в худшем случае вычисляемый срез списка – это alist[1:]. Его сложность – O(N-1) = O(N). Однако, когда i == len(alist) этот срез будет пуст. Средний срез содержит N/2 значений, что по-прежнему даем нам сложность O(N).

***

Вместо срезов можно использовать вложенный цикл. Сложность функции при этом не изменится.

        def is_unique1 (alist : [int]) ->bool:
  for i in range(len(alist)):          # O(N)
    for j in range(i+1, len(alist)):   # O(N)
      if alist[i] == alist[j]          # O(1)
        return False                   # O(1)
 return True                           # O(1)
    

Класс сложности целой функции тот же:

O(N)×O(N)×O(1)+O(1)=O(N2)

Алгоритм 2

Список уникален, если после его сортировки рядом не находятся одинаковые значения.

Сначала мы копируем список, чтобы не изменять исходную структуру сортировкой – функции обычно не должны влиять на входящие параметры.

        def is_unique2 (alist : [int]) -> bool:
  copy = list(alist)             # 1
  copy.sort()                    # 2
  for i in range(len(alist)-1):  # 3
    if copy[i] == copy[i+1]:     # 4
      return False               # 5
  return True                    # 6
    

Сложность по строчкам:

  1. O(N).
  2. O(N log N) – для быстрой сортировки.
  3. O(N) – на самом деле N-1, но это то же самое. Операции получения размера списка и вычитания имеют сложность O(1).
  4. O(1) – сложение, две операции получения элемента по индексу и сравнение – все со сложностью O(1).
  5. O(1) – в худшем случае никогда не выполняется.
  6. O(1) – в худшем случае всегда выполняется.

Общий класс сложности функции is_unique2:

O(N)+O(N×log⁡N)+O(N)×O(1)+O(1)=
O(N+Nlog⁡N+O(N×1)+1)=
O(N+Nlog⁡N+N+1)=O(Nlog⁡N+2N+1)=O(Nlog⁡N)

Сложность этой реализации меньше, чем is_unique1. Для больших значений N is_unique2 будет выполняться значительно быстрее.

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

Класс сложности O(N log N) хуже, чем O(N), но уже значительно лучше, чем O(N2).

Фактически, в код метода можно внести одно упрощение:

        # было
copy = list(alist)    # O(N)
copy.sort()           # O(N log N)

# стало
copy = sorted(alist)  # O(N log N)
    

Функция sorted создает список с теми же значениями и возвращает его отсортированную версию. Поэтому нам не требуется явно создавать копию.

Это изменение ускорит выполнение кода, но не повлияет на его сложность, так как O(N + N log N) = O(N log N). Ускорение – это всегда хорошо, но намного важнее найти алгоритм с минимальной сложностью.

Нужно отметить, что is_unique2 работает только в том случае, если все значения в списке сравнимы между собой (для сортировки используется оператор сравнения <). Если список заполнен одновременно строками и числами, ничего не получится. В то же время is_unique1 использует только оператор ==, который может сравнивать значения разных типов без выбрасывания исключений.

Алгоритм 3

Список уникален, если при превращении в множество (set) его размер не изменяется.

        def is_unique3 (alist : [int]) -> bool:
  aset = set(alist)               # O(N)
  return len(aset) == len(alist)  # O(1)
    

Рассчитать класс сложности для всей функции очень просто:

O(N)+O(1)=O(N+1)=O(N)

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

Тело функции можно записать в одну строчку:

        return len(set(alist)) == len(alist)
    

Сложность при этом не изменится.

В отличие от is_unique2, эта реализация может работать и со смешанными списками (числа и строки). В то же время требуется, чтобы все значения были хешируемыми/неизменяемыми. Например, is_unique3 не будет работать для списка списков.

***

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

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

Кроме того, при анализе важно учитывать дополнительные ограничения реализации (например, типы данных).

Приоритетные очереди

Рассмотрим еще один важный пример зависимости вычислительной сложности от реализации алгоритма – приоритетные очереди.

Этот тип данных поддерживает две операции:

  • добавление значения;
  • извлечение значения с самым высоким приоритетом из оставшихся.

Разные реализации приоритетных очередей имеют разные классы сложности этих операций:

add remove
Реализация 1 O(1) O(N)
Реализация 2 O(N) O(1)
Реализация 3 O(log N) O(log N)

Реализация 1

Новое значение добавляется в конец списка (list) или в начала связного списка (linked list) – O(1).

Чтобы найти приоритетное значение для удаления осуществляется поиск по списку – O(N).

Легко добавить, но трудно удалить.

Реализация 2

Для нового значения определяется правильное место – для этого перебираются элементы списка (или связного списка) – O(N).

Зато самое приоритетное значение всегда находится в конце списка (или в начале связного списка) – O(1).

Легко удалить, но трудно добавить.

Реализация 3

Используется двоичная куча, которая позволяет реализовать обе операции со “средней” сложностью O(log N), что для больших значений ближе к O(1), чем к O(N).

Теперь проведем анализ этих реализаций на реальных задачах.

Сортировка

Чтобы отсортировать N значений с помощью приоритетной очереди, нужно сначала добавить все значения в очередь, а затем удалить их все.

Реализация 1

N×O(1)+N×O(N)=O(N)+O(N2)=O(N2)

Реализация 2

N×O(N)+N×O(1)=O(N2)+O(N)=O(N2)

Реализация 3

N×O(log⁡N)+N×O(log⁡N)=
O(Nlog⁡N)+O(Nlog⁡N)=O(Nlog⁡N)

Примечание: N * O(…) – это то же самое, что и O(N) * O(…).

Первая и вторая реализации выполняют одну операцию быстро, а другую медленно. В итоге общий класс сложности определяет самая медленная операция.

Третья реализация оказывается быстрее всех, так как ее среднее время O(N log N) ближе к O(1), чем к O(N).

10 максимальных значений

Теперь найдем с помощью приоритетной очереди 10 максимальных значений. Для этого в очередь помещаются N элементов, а извлекаются только 10.

Реализация 1

N×O(1)+10×O(N)=O(N)+O(N)=O(N)

Реализация 2

N×O(N)+10×O(1)=O(N2)+O(1)=O(N2)

Реализация 3

N×O(log⁡N)+10×O(log⁡N)=
O(Nlog⁡N)+O(log⁡N)=O(Nlog⁡N)

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

***

Мораль этого примера заключается в том, что иногда не существует “самой лучшей реализации”. В зависимости от задачи оптимальными могут быть разные варианты.

А чтобы найти лучшее решение, важно понимать, что такое вычислительная сложность и как ее определить 🙂

09
Сен
2020

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

Графы стали мощным средством моделирования и представ…

18
Авг
2020

Какой самый сложный алгоритм вы использовали в своей работе — рассказывают эксперты

Помимо привычных задач порой попадаются алгоритмы, которые надолго запоминаешь. Спросили у экспертов о таких алгоритмах.
— Читать дальше «Какой самый сложный алгоритм вы использовали в своей работе — рассказывают эксперты»

21
Май
2020

Зачем программисту изучать алгоритмы

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

14
Апр
2020

Когда антивирус бессилен: распознавание вредоносных программ методами компьютерного зрения

Небольшая заметка про отображение бинарных файлов в виде картинок и применение средств computer vision для обнаружения троянов и вирусов.

Любой софт,
предназначенный специально д…

04
Мар
2020

Хороший, плохой, злой: как Яндекс использует нейросети для борьбы со спамом и матом

Разработчики Яндекса автоматизировали борьбу со спамом и матом в своих сервисах. Рассказываем, какие инструменты они для этого использовали.

Любой сайт, где пользователи могут са…

24
Фев
2020

Как шифровать информацию с помощью роста кристаллов

Учёные из Университета Глазго показали, что истинно случайные числа для криптографии можно создавать, наблюдая за ростом кристаллов. Рассказываем, как это работает.

Эпоха автомат…

04
Фев
2020

Python и динамическое программирование на примере задачи о рюкзаке

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

Момент настал – вы переезжаете в Сан-Франциско, чтобы стать известным дата сайентистом. Будете работать в крупной компании и скоро прославитесь. Но не всё сразу. Поначалу придётся ютиться в жилище, много меньшем, чем нынешний дом. Нужно решить, что взять с собой. Вы аналитик и намерены решить задачу эффективно. Обсудим алгоритм?

Описание проблемы

Представим, что на данный момент вы живёте в довольно большой квартире площадью 100 м2, а вещи занимают 50 м2. Вы знаете, что площадь нового жилья составит 40 м2 и хотите, чтобы вещи занимали не более 20 м2. Вы любите, когда на полу есть свободное место, да и не всё можно поставить впритык друг к другу. Нужно составить список вещей, определить площадь, занимаемую каждой вещью и составить рейтинг ценности предметов. Конечная цель – максимизировать материальную ценность того, что разместится на 20 квадратных метрах.

Перечисляем предметы

Вы составили список вещей и выразили занимаемую площадь в квадратных дециметрах (1 дм2 = 0.01 м2). Каждому предмету сопоставили ценность по шкале от 0 до 100. Получилась сводка данных о 29 предметах, которую можно оформить как словарь вида {'название предмета':(площадь, ценность) }:

            stuffdict = {'couch_s':(300,75), 
             'couch_b':(500,80), 
             'bed':(400,100), 
             'closet':(200,50), 
             'bed_s':(200,40), 
             'desk':(200,70), 
             'table':(300,80),
             'tv_table':(200,30),
             'armchair':(100,30),
             'bookshelf':(200,60), 
             'cabinet':(150,20),
             'game_table':(150,30),
             'hammock':(250,45),
             'diner_table_with_chairs':(250,70),
             'stools':(150,30),
             'mirror':(100,20),
             'instrument':(300,70),
             'plant_1':(25,10),
             'plant_2':(30,20),
             'plant_3':(45,25),
             'sideboard':(175,30),
             'chest_of_drawers':(25,40),
             'guest_bed':(250,40),
             'standing_lamp':(20,30), 
             'garbage_can':(30,35), 
             'bar_with_stools':(200,40), 
             'bike_stand':(100,80),
             'chest':(150,25),
             'heater':(100,25)
            }
        

Готово! Теперь видно, что кровать ('bed') занимает 400 дм2, а её ценность составляет 100, игровой стол ('game_table') занимает 150 дм2 квадратных метра и имеет ценность 30.

Максимазация ценности

Как максимизировать суммарную ценность объектов так, чтобы суммарная площадь не превышала 2000 дм2? Можно попробовать перебрать все возможные комбинации и рассчитать суммарную ценность для каждой из комбинаций. Решение получится, если выбрать максимальную суммарную ценность для 2000 дм2. Но вот незадача: и комбинаторика, и теория множеств утверждают, что 29 предметов дают 2²⁹ возможных комбинаций выбора предметов.

То есть более 536 миллионов. Похоже, такой перебор займёт некоторое время. Нельзя ли быстрее?

Стоит подумать о других вариантах. Что если начать с наиболее ценного предмета, затем следующего за ним по ценности, и так до тех пор, пока не заполнятся 20 квадратных метров? Такой алгоритм явно быстрее. Но даст ли он оптимальное решение? Есть сомнения. Как быстро решить задачу и найти лучшее решение?

Примечание: описанные выше случаи соответствуют полному перебору (метод «грубой силы», англ. brute force) и жадному (англ. greedy) алгоритму.

***

Динамическое программирование

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

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

Этапы решения задачи с помощью динамического программирования

  1. Создаём словарь со списком элементов и их параметрами (площадь, ценность).
  2. Создаём списки значений площади и ценности.
  3. Используем списки для построения таблиц мемоизации.
  4. Получаем элементы из последней строки таблицы. Последняя строка таблицы мемоизации содержит оптимальное решение. Возможно, существует несколько оптимальных решений. Например, когда есть два предмета с одинаковой площадью и стоимостью или ряд предметов имеют суммарные площадь и ценность, как у другого предмета.

Рассмотрим, как реализовать этот план на практике. Первый шаг мы предусмотрительно выполнили ранее, перейдём ко второму.

Создаём списки значений площади и ценности

Разделяем списки значений исходного словаря, например, так:

            def get_area_and_value(stuffdict):
    area = [stuffdict[item][0] for item in stuffdict]
    value = [stuffdict[item][1] for item in stuffdict]        
    return area, value
        

Используем списки для мемоизации

Пусть n – общее число предметов, A – их максимально допустимая суммарная площадь в новом жилище (2000 дм2). Составим таблицу из n + 1 строк и A + 1 столбцов. Строки пронумеруем индексом i, столбцы – индексом a (чтобы помнить что в столбцах площадь, area). То есть в качестве номеров столбцов мы рассматриваем дискретные значения площади, отличающиеся друг от друга на 1.

            def get_memtable(stuffdict, A=2000):
      area, value = get_area_and_value(stuffdict)
      n = len(value) # находим размеры таблицы
      
      # создаём таблицу из нулевых значений
      V = [[0 for a in range(A+1)] for i in range(n+1)]

      for i in range(n+1):
            for a in range(A+1):
                  # базовый случай
                  if i == 0 or a == 0:
                        V[i][a] = 0

                  # если площадь предмета меньше площади столбца,
                  # максимизируем значение суммарной ценности
                  elif area[i-1] <= a:
                        V[i][a] = max(value[i-1] + V[i-1][a-area[i-1]], V[i-1][a])

                  # если площадь предмета больше площади столбца,
                  # забираем значение ячейки из предыдущей строки
                  else:
                        V[i][a] = V[i-1][a]       
      return V, area, value
        

Вначале создаётся и заполняется таблица в виде вложенных списков. Нулевую строку и нулевой столбец заполняем нулями. Это базовый случай: когда площадь или количество элементов равны нулю, значение ячейки равно нулю. Далее таблица значений V будет заполняться строка за строкой слева направо и сверху вниз:


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

            V[i][a] = max(value[i-1] + V[i-1][a-area[i-1], V[i-1][a])
        

То есть выбираем максимальное из двух значений:

  1. Сумма ценности текущего предмета value[i-1] и величины элемента из предыдущей строки i-1 с площадью, меньшей на величину площади текущего предмета area[i-1]. Обратите внимание: нужно помнить, что элементы в таблице отличаются по нумерации на единицу из-за нулевой строки.
  2. Значение элемента предыдущей строки с той же площадью, то есть из того же столбца, что текущая ячейка. То же значение устанавливается в случае, если площадь текущей ячейки меньше, чем площадь текущего элемента (см. блок else)

За счёт такого подхода при одной и той же суммарной площади элементов происходит максимизация суммарной ценности.

Забираем нужные элементы из последней строки таблицы

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

            def get_selected_items_list(stuffdict, A=2000):
      V, area, value = get_memtable(stuffdict)
      n = len(value)
      res = V[n][A]      # начинаем с последнего элемента таблицы
      a = A              # начальная площадь - максимальная
      items_list = []    # список площадей и ценностей
    
      for i in range(n, 0, -1):  # идём в обратном порядке
            if res <= 0:  # условие прерывания - собрали "рюкзак" 
                  break
            if res == V[i-1][a]:  # ничего не делаем, двигаемся дальше
                  continue
            else:
                  # "забираем" предмет
                  items_list.append((area[i-1], value[i-1]))
                  res -= value[i-1]   # отнимаем значение ценности от общей
                  a -= area[i-1]  # отнимаем площадь от общей
            
      selected_stuff = []

      # находим ключи исходного словаря - названия предметов
      for search in items_list:
            for key, value in stuffdict.items():
                  if value == search:
                        selected_stuff.append(key)
            
      return selected_stuff
        

Итак, мы нашли список:

            >>> stuff = get_selected_items_list(stuffdict)
>>> print(stuff)
['bike_stand', 'garbage_can', 'standing_lamp', 'chest_of_drawers',
'plant_3', 'plant_2', 'diner_table_with_chairs', 'bookshelf',
'armchair', 'table', 'desk', 'bed', 'couch_s']
        

Проверим суммарные площадь и ценность собранных предметов:

            >>> totarea = sum([stuffdict[item][0] for item in stuff])
>>> totvalue = sum([stuffdict[item][1] for item in stuff])
>>> print(totarea)
2000
>>> print(totvalue)
715
        

Получилось! Собираем вещи и отправляемся в путь (звучит песня Mamas & Paps – San Francisco).

***

Бонус: Тепловая карта таблицы

Тепловая карта ниже отображает рассмотренную выше таблицу. По оси абсцисс отложена доступная площадь, по оси ординат – список предметов. Цветовая шкала соответствует суммарной ценности. Можно видеть, как значения возрастают по мере продвижения по таблице снизу вверх.


Для построения использовалась библиотека matplotlib:

            def plot_memtable(V, stuffdict):
    plt.figure(figsize=(30,15))
    item_list = list(stuffdict.keys())
    item_list.insert(0, 'empty')
    sns.heatmap(V, yticklabels=item_list)
    plt.yticks(size=25)
    plt.xlabel('Area', size=25)
    plt.ylabel('Added item', size=25)
    plt.title('Value for Area with Set of Items', size=30)
    plt.show()
        

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

15
Янв
2020

Как программисту создать картинку без Фотошопа

Нужны уникальные картинки, но рисовать — слишком муторно и сложно? Пусть это сделает алгоритм. Рассказываем про generative art — искусство программистов.
— Читать дальше «Как программисту создать картинку без Фотошопа»

16
Дек
2019

6 алгоритмов сортировки в гифках

Освежите знания алгоритмов сортировки с помощью нашей визуализированной подборки. Примеры кода на С#.

Сортировка – важный аспект программирования. Отсортированные данные позволяют быстрее и эффективнее рабо…

17
Окт
2019

OpenAI научила роборуку собирать кубик Рубика на весу

Со сложными конфигурациями из 26 поворотов она справляется в 20 % случаев. Если надо сделать 15 поворотов, количество успешных сборок возрастает до 60 %.
— Читать дальше «OpenAI научила роборуку собирать кубик Рубика на весу»

12
Окт
2019

6 шагов по созданию проектов машинного обучения

Статья расскажет, как приступить к созданию проекта с машинным обучением. Какие данные необходимо собирать, как правильно моделировать и развёртывать.
— Читать дальше «6 шагов по созданию проектов машинного обучения»

11
Окт
2019

Алгоритмы и структуры данных на C++: деревья отрезков

Рассмотрена древовидная структура данных для эффективной обработки запросов к массивам по диапазону. Затронуты построение дерева, получение значения на отрезке, изменение элемента и изменение на отрезке.

13
Сен
2019

Организаторы Международной математической олимпиады готовят конкурс на создание алгоритмов

Алгоритмам дадут столько же времени, сколько и участникам-людям, и решать они будут классические олимпиадные математические задачи.
— Читать дальше «Организаторы Международной математической олимпиады готовят конкурс на создание алгоритмов»

12
Сен
2019

ИИ научился управлять 3D-фигуркой в соответствии с текстовыми указаниями

В будущем такой алгоритм позволит управлять роботами или персонажами в видеоиграх простым языком: пойди туда, сделай то, принеси это.
— Читать дальше «ИИ научился управлять 3D-фигуркой в соответствии с текстовыми указаниями»