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

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

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

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

Нетрудно убедиться, что в каждой колонке исходных данных числа одни и те же, только идут в разном порядке. Числа красивые, круглые, то есть степени двойки. Многозначное число - это 2 в степени 25. И также несложно выполнить необходимые сложения в уме. Сумма пяти чисел в каждой колонке равна единице. А вот программное решение не так однозначно.

Эту задачу очень удобно решать с помощью векторных вычислений, доступных во всех современных процессорах. Мы используем систему команд AVX512, доступную в процессорах компаний Intel и AMD. Далее мы предполагаем, что Вам знакомы азы векторного программирования, и вместо них сосредоточимся на сути происходящего.

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

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

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

Любопытно, что пятая инструкция загрузки вполне разумно уехала вниз. Процессоры Intel и AMD оборудованы разными конвейерами внеочередных вычислений, и даже разные модели от одного производителя имеют различные исполнительные блоки, но порядок выполнения нашей учебной программы предсказать несложно. После первых двух команд загрузки данных (строки 4 и 5), процессор просмотрит список последующих команд. И, не дожидаясь выполнения команд загрузки в строках 6 и 7, начнет выполнение сложения в строке 8. По мере того, как будут загружены данные в строках 6 и 7, не дожидаясь завершения сложения в строке 8, будет выполнено сложение в строке 10. На этом наши возможности параллелизации исчерпаны, и остальные команды будут выполняться последовательно. Но правую ветку дерева операций удалось исполнить параллельно левой.

Давайте попробуем выполнить этот код и напишем небольшую тестовую программу с использованием Google Test.

Прогоним тест и получим неожиданный результат: во всех четырех столбцах получились разные суммы.

Корень проблемы заключается в особенности двоичного представления чисел с плавающей точкой. Стандарт языка C не оговаривает, как именно числа хранятся в памяти. Поэтому в начале 1980-х разработчики сопроцессора для работы с плавающей точкой Intel 8087 были свободны в своем творчестве. Продукт этого творчества превратился в стандарт представления плавающих чисел IEEE 754 и расползся по различным архитектурам. В спецификации расширения команд AVX512 уже указано двоичное представление именно в формате IEEE 754. Вот как по этому стандарту выглядит длинное число 2 в степени 25 в памяти:

Самый старший бит - это знак числа, 0 для положительных чисел. Затем, согласно стандарту, следует показатель степени двойки этого числа плюс 127. 25 плюс 127 равно 152. А дальше хранится мантисса числа. Как видим, в данном случае, состоящая из одних нулей. Дело в том, что, по стандарту IEEE 754, у нормального числа с плавающей точкой мантисса всегда имеет вид 1.xxxxxx. И, поскольку начинается мантисса всегда с единицы, эта единица не хранится в памяти, а подразумевается. Когда число загружается из памяти или из регистра в арифметико - логическое устройство (АЛУ), эта единица оказывается в нужном месте мантиссы. А также, мантисса расширяется на несколько бит. Чаще всего на восемь.

Давайте посмотрим, как проходит процедура сложения. Итак, 2 в 25 степени мы складываем с двойкой.

Загрузим оба слагаемых в АЛУ.

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

После этого можно произвести операцию сложения.

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

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

Теперь давайте посмотрим, что происходит при сложении отрицательного числа -(2²⁵) c двойкой.

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

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

Затем к мантиссе прибавляется единица.

Теперь мы готовы произвести сложение.

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

После нормализации, мы выгружаем число обратно в регистр, опять же, отбрасывая дополнительные разряды мантиссы.

Но при этом результат получается точным, округления не произошло. В данном случае дополнительные биты представления чисел в АЛУ спасли ситуацию.

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

Как бороться с ошибками округления? В этом нам поможет алгоритм Кагана. Давайте посмотрим, как он выглядит.

В строке 11 складывается первая пара чисел. А в 12 строке вычисляется возможная ошибка округления. Хочу обратить внимание на то, что с точки зрения обычной математики, результат выполнения строки 12 тождественно равен нулю. Иногда ретивый компилятор может догадаться об этом и оптимизировать ваш код, просто убрав вычисление в 12 строке. Следите, чтобы он этого не сделал, и проверяйте cгенерированный машинный код. Затем в строках 14 и 15 складывается вторая пара чисел из правой ветки дерева операций. В строках 17, 18 и 19 происходит очередное сложение и перевычисление накопившейся ошибки. Наконец, в строке 21 происходит заключительное сложение, и результат корректируется путем прибавления обеих накопленных ошибок. Запускаем Google Test этого кода и видим, что теперь результат совпадает с ожидаемым.

Мы победили погрешность вычислений, утроив количество необходимых арифметических операций: 13 против 4 исходных. Не удивительно, что в погоне за максимальной производительностью, измеряемой в количестве "полезных" операций в секунду, ни один из известных мне ускорителей не занимается коррекцией по Кагану.

Быстрому увеличению погрешности способствует и сокращение разрядности исходных данных, и неудачный выбор формата "коротких" чисел. Стандарт IEEE 754 определяет для 16-битового числа с плавающей точкой (FP16) 5 бит экспоненты и 10 бит мантиссы. Но в алгоритмах машинного обучения широко применяют формат BrainFloat16, полученный из 32-битового числа простым отсечением младших 16 бит мантиссы. Переводить из одного формата в другой просто, вот только в числе осталось по прежнему 8 бит экспоненты и всего 7 бит мантиссы. А мы видели, что именно длина мантиссы определяет склонность к возникновению ошибок округления. Конечно, несколько помогает то, что сумма нарастающим итогом часто хранится в виде 32-битового числа, но и это не панацея, как мы убедились на вышеизложенном примере.

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