Этот пост был первоначально опубликован на https://sahnimanas.github.io 26 августа 2019 г.

На моем не слишком потрепанном процессоре ноутбука я могу запускать самые распространенные модели CNN за (максимум) 10–100 миллисекунд с такими библиотеками, как TensorFlow. В 2019 году даже смартфон может запускать «тяжелые» модели CNN (например, ResNet) менее чем за полсекунды. Представьте себе мое удивление, когда я рассчитал свою простую реализацию сверточного слоя и обнаружил, что для одного слоя требуется более 2 секунды!

Неудивительно, что современные библиотеки глубокого обучения имеют высокооптимизированные реализации большинства операций производственного уровня. Но что же такое черная магия, которую используют эти библиотеки, а мы, простые смертные, не используем? Как они могут повысить производительность в 100 раз? Что именно нужно делать для «оптимизации» или ускорения работы нейронной сети? Это вопросы, которые я часто задаю (и получаю), когда говорю о высокопроизводительных / эффективных DNN.

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

Многое из того, что я здесь освещаю, взято из основополагающей статьи Анатомия высокопроизводительного умножения матриц Гото и др., Которая легла в основу алгоритмов, используемых в библиотеках линейной алгебры, таких как OpenBLAS. ; и эти полезные уроки от Dr. Стефан Хаджис и Крис Роуз .

Наивные свертки

«Преждевременная оптимизация - корень всех зол». - Дональд Кнут

Прежде чем приступить к оптимизации, давайте сначала разберемся с нашими базовыми показателями и узкими местами. Вот наивная свертка numpy / for-loop, которую вы бы написали в Deep Learning 101:

''' Convolve `input` with `kernel` to generate `output` 
    input.shape = [input_channels, input_height, input_width]
    kernel.shape = [num_filters, input_channels, kernel_height, kernel_width] 
    output.shape = [num_filters, output_height, output_width]
'''
for filter in 0..num_filters
  for channel in 0..input_channels
    for out_h in 0..output_height
      for out_w in 0..output_width
        for k_h in 0..kernel_height
          for k_w in 0..kernel_width
            output[filter, channel, out_h, out_h] += 
                kernel[filter, channel, k_h, k_w] 
                * input[channel, out_h + k_h, out_w + k_w]

Ой. Это 6 вложенных for циклов для одной конверсии (7, если вы перебираете пакеты из нескольких входных данных). И мы еще не рассматриваем шаг, расширение или какие-либо другие параметры. Если здесь я вставлю размеры, скажем, для первого уровня MobileNet и запущу его на простом языке C, это займет колоссальные 22 секунды! При наиболее агрессивных оптимизациях компилятора, таких как -O3 или -Ofast, время сокращается до 2,2 секунды. Но это все еще ужасно медленно только для первого слоя.

Что, если я запустил тот же слой, используя, скажем, Caffe? На том же ПК это заняло всего 18 мс. Это более чем в 100 раз ускорение! На моем процессоре вся сеть работает примерно за 100 мсек.

В чем заключается узкое место и где мы должны начать оптимизацию?

Самый внутренний цикл выполняет 2 операции с плавающей запятой (умножение и сложение), и для размеров, которые я использовал, он был выполнен примерно 85,16 миллиона раз; то есть эта свертка требует 170 миллионов операций с плавающей запятой (MFLOP). Пиковая производительность моего процессора, по данным Intel, составляет 80 миллиардов FLOP в секунду, что означает, что он может теоретически завершить работу за 0,002 секунды. Ясно, что мы далеки от этого, и также ясно, что возможности обработки сырых данных здесь более чем достаточно.

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

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

Некоторые предпосылки

FLOP / с

Нашей метрикой для «производительности» или скорости будет пропускная способность, измеряемая в количестве FL всплесков P oint O операций на S секунда или FLOP / s. Большая операция с большим количеством операций с плавающей запятой, естественно, будет выполняться медленнее, поэтому скорость FLOP / s позволяет более согласованно сравнивать производительность.

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

  • есть 2 физических ядра
  • каждое ядро ​​имеет частоту 2,5 ГГц или 2,5 × 10 циклов ЦП в секунду
  • в каждом цикле он может обрабатывать 32 FLOP (с использованием AVX и FMA - подробнее об этом чуть позже)

Это дает пиковую производительность 2 x (2,5 x 10⁹ циклов в секунду) x (32 FLOPs / цикл) = 160 GFLOP / с. Это теоретический пик моего процессора. Точно так же для одного ядра это число составляет 80 GFLOP / s.

Порядки хранения и основные строки

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

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

А как насчет порядка самих размеров? Обычно для 4-мерных тензоров, таких как CNN, вы услышите о порядках хранения, таких как NCHW, NHWC и т. Д. Я буду предполагать NCHW на протяжении всей этой статьи - если у меня есть N блоков C-каналов изображений HxW, то все изображения с одинаковые N смежны; внутри этого блока все пиксели одного и того же канала C являются смежными и так далее.

Галогенид

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

Это упрощает эксперименты с различными оптимизациями за счет разделения алгоритма (что вычислять) и расписания (как / когда вычислять). Мы можем оставить алгоритм фиксированным и поэкспериментировать с разными расписаниями.

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

От свертки к GEMM

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

Можем ли мы вместо этого превратить это в проблему, которую легче решить? Возможно матричное умножение?

Умножение матриц, или matmul, или Ge neralized M atrix M ultiplication (GEMM), лежит в основе глубокого обучения ». Он используется в полносвязных слоях, RNN и т. Д., А также может использоваться для реализации сверток.

В конце концов, Conv - это скалярный продукт фильтра с входными патчами. Если мы разместим фильтр в 2D-матрице, а входные патчи - в другой, то умножение этих двух матриц приведет к вычислению одного и того же скалярного произведения. А матричное умножение - в отличие от CNN - тщательно изучалось и оптимизировалось в течение нескольких десятилетий, что является важной проблемой во многих научных областях.

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

Рассмотрим эту нормальную прямую свертку 3x3:

Ниже представлена ​​та же операция, реализованная в виде умножения матриц. Правая матрица является результатом im2col - она ​​должна быть построена путем копирования пикселей из исходного изображения. Левая матрица имеет веса сверток, которые уже сохранены в памяти таким образом.

Обратите внимание, что матричное произведение дает нам результат конверсии напрямую - нет необходимости в дополнительном «преобразовании» обратно в исходную форму.

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

С помощью im2col мы теперь преобразовали операцию свертки в умножение матриц. Теперь мы можем подключить более универсальные и популярные библиотеки линейной алгебры, такие как OpenBLAS, Eigen и т. Д., Чтобы позаботиться об эффективном вычислении этого matmul, опираясь на десятилетия оптимизации и тщательной тонкой настройки.

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

Ускорение GEMM

Наивный

В оставшейся части этого поста я предполагаю, что GEMM выполняется как C + = A * B, где C имеет размер (MxN), A - (MxK), а B - (KxN).

Как и раньше, сначала давайте примерим простое умножение матриц из учебника:

for i in 0..M:
  for j in 0..N: 
    for k in 0..K:
      C[i, j] += A[i, k] * B[k, j]

Или в галиде:

Halide::Buffer C, A, B;
Halide::Var x, y; 
C(x,y) += A(k, x)*B(y, k); // loop bounds, dims, etc. are taken care of automatically

Самая внутренняя строка выполняет 2 операции с плавающей запятой (умножение и сложение) и выполняется M * N * K раз, поэтому количество FLOP для этого GEMM составляет 2MNK.

Давайте измерим его производительность для различных размеров матриц:

Мы едва достигаем 10% максимальной производительности! Хотя мы рассмотрим способы ускорения вычислений, повторяющаяся тема будет заключаться в том, что недостаточно просто вычислить данные быстро, если мы не можем получить данные. быстро.

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

Кеширование

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

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

Первое, на что вы должны обратить внимание, - это схема, по которой мы получаем доступ к нашим данным. Мы перемещаемся по строкам на A и по столбцам на B.

Их хранилище также является строковым, поэтому, как только мы находим A[i, k], следующий элемент в строке, A[i, k+1] уже кэшируется. Прохладный. Но посмотрите, что происходит с B:

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

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

Давайте изменим порядок циклов с i,j,k на i,k,j:

for i in 0..M:
  for k in 0..K:
    for j in 0..N:
      C[i,j] += A[i,k]*B[k,j]

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

Это простое изменение простого переупорядочивания циклов дает значительное ускорение:

Черепица

Чтобы еще больше улучшить переупорядочение, необходимо рассмотреть еще одну проблему с кешированием.

Для каждой строки A мы перебирали весь столбец B. С каждым шагом по B мы загружаем некоторые из его новых столбцов и удаляем некоторые старые из кеша. Когда мы переходим к следующему ряду A, мы начинаем все заново с первых столбцов. Мы постоянно добавляем и удаляем одни и те же данные из кеша или обрабатываем их.

Если бы все наши данные умещались в кеше, этого бы не произошло. Если бы мы работали с матрицами меньшего размера, они могли бы счастливо жить вместе, не подвергаясь многократному выселению. К счастью для нас, мы можем разбить матричное умножение на подматрицы. Чтобы вычислить небольшой фрагмент rxc C, нам нужно всего r строк из A и c столбцов Б. Разобьем C на плитки размером 6x16:

C(x, y) += A(k, y) * B(x, k);
C.update() 
 .tile(x, y, xo, yo, xi, yi, 6, 16) 
/* in pseudocode: 
for xo in 0..N/16: 
  for yo in 0..M/6: 
    for yi in 6: 
      for xi in 0..16: 
        for k in 0..K: 
          C(...) = ... 
*/

Мы разделили размеры x, y на внешний xo,yo и внутренний xi,yi. Мы потратим наши усилия на оптимизацию микроядра для меньшего блока 6x16 (отмеченного xi,yi) и запустим это микроядро по всем блокам (повторяется xo,yo).

Векторизация и FMA

Большинство современных процессоров поддерживают SIMD или S одну I инструкцию M Ultimate D ata. . Как следует из названия, SIMD можно использовать для выполнения одной и той же операции / инструкции (например, сложение, умножение и т. Д.) Над несколькими значениями одновременно в одном и том же цикле ЦП. Если мы можем запускать инструкции SIMD, скажем, для 4 точек данных за раз, это сразу дает четырехкратное ускорение.

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

Второй "прием", который мы использовали при вычислении пиковых значений FLOP, был FMA - F использовал M ultiply- A dd. . Хотя умножение и сложение считаются двумя отдельными операциями с плавающей запятой, они настолько распространены, что доступны специальные аппаратные блоки, которые объединяют 2 и выполняют их как одну инструкцию. Этим обычно занимается компилятор.

На процессорах Intel мы можем использовать SIMD (называемые AVX и SSE) для обработки до 8 чисел с плавающей запятой в одной инструкции. Оптимизация компилятора часто сама по себе определяет возможности векторизации, но на всякий случай мы возьмем все в свои руки.

C.update()
 .tile(x, y, xo, yo, xi, yi, 6, 16) 
 .reorder(xi, yi, k, xo, yo)
 .vectorize(xi, 8)
/*
in pseudocode: 
for xo in 0..N/16:
  for yo in 0..M/6:
    for k in 0..K:
      for yi in 6:
        for vectorized xi in 0..16: (in batches of 8)
          C(...) = ... 
*/

Резьба

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

C.update()
.tile(x, y, xo, yo, xi, yi, 6, 16) .reorder(xi, yi, k, xo, yo) .vectorize(xi, 8)
.parallel(yo)
/* in pseudocode:
for xo in 0..N/16 in steps of 16: 
  for parallel yo in steps of 6: 
    for k in 0..K: 
      for yi in 6:
      for vectorized xi in 0..16 in batches of 8:
        C(...) = ...
*/

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

Разворачивание

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

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

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

C.update()
 .tile(x, y, xo, yo, xi, yi, 6, 16)
 .reorder(xi, yi, k, xo, yo)
 .vectorize(xi, 8)
 .unroll(xi)
 .unroll(yi) 
/* in pseudocode: 
for xo in 0..N/16: 
  for parallel yo:
    for k in 0..K:
      C(xi:xi+8, yi+0)
      C(xi:xi+8, yi+1) 
      ...
      C(xi:xi+8, yi+5) 
      C(xi+8:xi+16, yi+0) 
      C(xi+8:xi+16, yi+1) 
      ... 
      C(xi+8:xi+16, yi+5)
*/

Теперь мы можем работать со скоростью до 60 GFLOP / с.

Собираем вместе

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

Вот расписание от Halide, которое было более тщательно оптимизировано:

matrix_mul(x, y) += A(k, y) * B(x, k);
out(x, y) = matrix_mul(x, y);
out.tile(x, y, xi, yi, 24, 32)
 .fuse(x, y, xy).parallel(xy) 
 .split(yi, yi, yii, 4)
 .vectorize(xi, 8)
 .unroll(xi)
 .unroll(yii);
matrix_mul.compute_at(out, yi)
 .vectorize(x, 8)
 .unroll(y); 
matrix_mul.update(0)
 .reorder(x, y, k)
 .vectorize(x, 8)
 .unroll(x)
 .unroll(y)
 .unroll(k, 2);

Таким образом, он делает следующее:

  • Разбивает out на плитки размером 32x24. Далее каждая плитка разбивается на части размером 8x24.
  • Вычисляет matmul 8x24 во временном буфере (matrix_mul) с аналогичным переупорядочиванием, векторизацией и разворачиванием.
  • Копирует временный matrix_mul обратно в out с векторизацией, развертыванием и т. Д.
  • Распараллеливает этот процесс по всем плиткам 32x24.

Наконец, мы можем достичь скорости более 120 GFLOP, что довольно близко к пиковой производительности в 160 GFLOP, и соответствуем таким библиотекам производственного уровня, как OpenBLAS. С аналогичным точно настроенным кодом для im2col, за которым следует GEMM, та же свертка теперь выполняется за ~ 20 мс.

Обратите внимание, что этот подход im2col + GEMM является популярным универсальным методом, используемым в большинстве библиотек глубокого обучения. Однако специализация является ключевой - для некоторых часто используемых размеров, различных архитектур (GPU) и различных параметров работы (таких как расширение, группировка и т. Д.) Эти библиотеки могут снова иметь более специализированные реализации, использующие аналогичные трюки и предположения, специфичные для этих случаев. .

Эти микроядра создаются в результате итеративного процесса, основанного на пробах и ошибках. Программисты часто только интуитивно понимают, что должно / не должно работать хорошо, и / или должны думать об объяснениях, основанных на результатах. Похоже, идеально подходит для исследований в области глубокого обучения, не так ли?

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

Являясь независимой редакцией, Heartbeat спонсируется и публикуется Comet, платформой MLOps, которая позволяет специалистам по данным и группам машинного обучения отслеживать, сравнивать, объяснять и оптимизировать свои эксперименты. Мы платим участникам и не продаем рекламу.

Если вы хотите внести свой вклад, отправляйтесь на наш призыв к участникам. Вы также можете подписаться на наши еженедельные информационные бюллетени (Deep Learning Weekly и Comet Newsletter), присоединиться к нам в » «Slack и подписаться на Comet в Twitter и LinkedIn для получения ресурсов, событий и гораздо больше, что поможет вам быстрее и лучше строить лучшие модели машинного обучения.