Об ошибках округления и способах борьбы с ними

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

В строке 11 складывается первая пара чисел. А в 12 строке вычисляется возможная ошибка округления. Хочу обратить внимание ^[6] на то, что с точки зрения ^[7] обычной математики ^[8], результат выполнения строки 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-битового числа, но и это не панацея, как мы убедились на вышеизложенном примере.

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