Привет, Хабр! Эта статья посвящена тестированию точности библиотек математических функций (libm). Мы рассмотрим, где эти библиотеки используются, почему они должны быть не только высокопроизводительными, но и высокоточными, а также поймем, откуда в корректных, на первый взгляд, вычислениях берутся ошибки и как их избежать. Узнаем, как устроено большинство тестов в стандартных математических библиотеках и почему они не всегда работают. И наконец, ответим на вопрос, как одним тестом полностью покрыть код математической функции. Без воды, регистрации и громоздких формул!
Сейчас, наверное, нет человека, который ни разу не видел аббревиатуры AI/ML, AR/VR, CV, HPC или не слышал об искусственном интеллекте, машинном обучении, дополненной или виртуальной реальности, компьютерном зрении. Но немногие задумываются, за счет чего такие приложения работают настолько быстро, «в режиме реального времени».
Причина в том, что ниже уровня воды у этих айсбергов находится высокопроизводительный математический бэкенд с низкоуровневыми оптимизациями под конкретное железо и процессоры. Математические библиотеки такого класса, как правило, можно разделить на три большие группы.
BLAS / LAPACK. Реализации основных операций линейной алгебры (сложение, вычитание, умножение) над плотными матрицами и векторами, а также решатели систем линейных уравнений с плотной матрицей и построение матричных разложений.
Sparse BLAS / Solvers. Реализации тех же операций, но уже над разреженными матрицами. Поскольку такие матрицы содержат очень большое число нулей, для них разработан ряд специальных форматов хранения, позволяющих не тратить память компьютера на многочисленные нулевые элементы. С учетом этих особенностей и выбранного формата хранения изменяются и алгоритмы. «Разреженность» здесь велика. В задачах глубокого обучения матрицы весов содержат от 50 до 90% нулевых элементов, а в приложениях математического моделирования доля нулей в матрицах систем линейных алгебраических уравнений часто может быть выше 99%.
О том, как разрабатываемые матричные расширения RISC-V могут помочь в оптимизации библиотек разреженной линейной алгебры, читайте в блоге YADRO.
Libm (стандартная математическая библиотека). Это реализации элементарных математических функций: тригонометрических, гиперболических, степенных, логарифмических, экспоненциальных и других. Именно такие функции мы рассмотрим далее.
Где используется libm
Популярные сейчас нейронные сети содержат не только нейроны, но и функции их активации: экспоненты, логарифмы, гиперболические тангенсы и т. д. В графических приложениях и компьютерном зрении часто требуются преобразования координат точек. Для этого координаты умножают на матрицу перехода между системами координат, состоящую из тригонометрических функций — синусов и косинусов. При моделировании физико-технических процессов на каждом расчетном шаге по времени нужно обновлять правую часть решаемой системы линейных алгебраических уравнений. В зависимости от математической модели, это требует вычисления тех или иных элементарных математических функций.
При прогнозировании погоды, например для расчета правой части системы, необходимо вычислить экспоненту. Если точность вычислений недостаточно высока, то с каждым шагом моделирования точность прогноза падает: мы можем терять тучку за тучкой, пока не исчезнет целый грозовой фронт. В итоге из-за ошибок в злосчастной экспоненте нам часто обещают солнце, пока небо стремительно затягивается и начинается ливень.
Откуда берутся эти ошибки? Конечно же, со временем года они никак не связаны. Ошибки возникают, потому что точность любых вычислений конечна и зависит от множества факторов. А ошибки аппроксимации математических функций нужно правильно минимизировать. Эти утверждения пока звучат очень обтекаемо, но далее мы с ними разберемся.
Чтобы подчеркнуть актуальность темы, отмечу, что реализация всех функций libm — это мастхэв в компиляторах языков программирования. Некоторые платформы даже содержат собственные реализации — например, свой libm есть в Bionic из AOSP (Android). Наиболее критичные для отдельных приложений элементарные математические функции реализованы также и в высокопроизводительном математическом бэкенде, используемом фреймворками. TensorFlow использует реализации из библиотеки XNNPACK, а OpenCV — из библиотеки Eigen.
В чем проблема с точностью
Распространенный способ тестирования точности libm — это сверка результата вычислений с эталонными значениями для предопределенного набора входных данных. Такой набор обычно состоит из нескольких десятков тысяч точек. Но если мы говорим о 32-битных числах с плавающей точкой (float), то их уже более 4 миллиардов. То есть при сверке проверяется порядка 0,001% всех точек. Вероятность того, что мы ошибаемся где-то между ними, мягко говоря, ненулевая.
Посмотрим более внимательно на стандартную библиотеку glibc, содержащую реализацию функций libm. Например, проверим точность реализации гиперболического тангенса (tanh), что часто встречается в активационных функциях нейронных сетей.
Возьмем примерно 37 500 точек, 64-битных чисел (double), вычислим в них tanh при помощи glibc и построим гистограмму точности: для каждого бита отметим, в скольких точках в нем содержится ошибка в вычисленном значении. Получаем, что уже в 52-м бите возникает ошибка в 450 точках, а начиная с 53-го бита ошибка есть во всех остальных. Это явно не соответствует стандартам точности libm.
Таким образом, даже в glibc есть ошибки из-за несовершенства системы тестирования точности.
Где мы теряем точность?
Создавая алгоритм вычисления математической функции, мы выполняем следующие шаги:
Вспоминаем свойства функции: четность, периодичность, масштабируемость. С их помощью мы можем выбрать отрезки, на которых необходимо рассчитывать значения функции через ее аппроксимацию, и отрезки, на которых можно получать значения функции по уже вычисленным, просто из свойств функции.
Подбираем алгоритм приведения аргумента к выбранному диапазону. Помните, как в школе делали преобразования графиков элементарных функций — сдвиги вдоль координатных осей, растяжение? О них и речь.
Смотрим, какие числа можно посчитать заранее, чтобы не тратить на это время при каждом вычислении значения функции. Такие обязательно найдутся, и мы их сводим в таблицу.
Аппроксимируем функцию на необходимых отрезках (мы их определили в шаге 1). На каждом таком отрезке мы строим полиномиальную аппроксимацию с заданной точностью.
Восстанавливаем результат: приводим аргумент в обратном порядке по сравнению с шагом 2.
Если внимательно посмотреть на каждый пункт, понимаешь, что, скорее всего, ошибки появляются на четвертом шаге — аппроксимации (приближении) функции.
Как можно аппроксимировать функцию, например, синуса? Вспомним простейший способ — ряд Тейлора. Думаю, многие его помнят из курса матанализа. Если мы аппроксимируем синус в окрестности нуля только одним членом ряда Тейлора, то это будет приближение обычной линейной функцией. Она хорошо аппроксимирует синус только в очень небольшой окрестности нуля.
Понятно, что нужно взять больше членов. Построим, например, ряд Тейлора до пятого порядка. На графике это красная линия F(х). Видно, что синус она приближает явно лучше, чем просто линейная функция, но с некоторого значения аргумента мы невооруженным глазом все равно видим расхождение красной и зеленой линий, то есть ошибку аппроксимации.
Спасет ли нас здесь большее число слагаемых? Можем, например, построить ряд Тейлора до седьмого члена, до девятого, одиннадцатого и так далее… Скорость расчета при этом, конечно же, будет снижаться, ведь нам придется совершать большее число арифметических операций. Но это еще полбеды: как ни парадоксально, численная устойчивость от этого тоже пострадает и точность не будет повышаться так стремительно, как это предсказывает теория. Почему это происходит?
Сделаем небольшое отступление на тему численной устойчивости.
Точность не бесконечна. Мы пытаемся отобразить бесконечное множество действительных чисел в числа, которые можно представить только конечным числом бит. Вот, например, как будет выглядеть число 0,1 в одинарной точности, то есть при использовании 32-битного формата с плавающей точкой:
float(0.1) = 0.1000000015Ошибка накапливается после каждого арифметического действия. Это вызвано округлением каждого результата. Чем больше членов ряда Тейлора мы берем, тем больше действий и округлений совершаем.
Числа с плавающей точкой неравномерно заполняют числовую ось. Почему это так, мы поймем чуть позже, а пока просто посмотрим на картинку и немного подумаем, к чему это может приводить.
Если очень маленькие числа складываются с очень большими, погрешности также получаются немаленькие, ведь чем больше числа, тем реже они расположены. Это примерно как пытаться сложить маленький протон и огромное Солнце: частица просто потеряется на фоне светила. А вот если сначала сложить миллиард протонов, а потом уже результат прибавлять к чему-то покрупнее, то все будет не так плачевно.
Законы арифметики не работают. Это следствие предыдущего пункта. У нас нет ни ассоциативности, ни дистрибутивности.
От перемены мест слагаемых сумма может поменяться очень значительно, как в примере с Солнцем и протонами. Если необходимо вычислить сумму элементов массива чисел сильно разного порядка, лучше отсортировать их в порядке возрастания, чтобы сначала складывать маленькие числа. Это уменьшит итоговую ошибку.
Какую аппроксимацию выбрать
Мы поняли, что увеличение числа членов ряда Тейлора нас не спасет. Давайте подумаем, что было особенно плохо в ранее построенной аппроксимации пятого порядка. Для наглядности построим график ошибки аппроксимации, то есть модуля разности нашего синуса и суммы ряда Тейлора, которым мы пытались его приблизить:
Что мы видим? Вблизи нуля ошибка мала, но на границах отрезка она начинает резко возрастать. Какой бы отрезок мы ни взяли, ошибка будет выглядеть именно так. Поэтому ряд Тейлора совсем не подходит для высокоточных реализаций функций из libm.
Что лучше использовать? Не буду углубляться в теорию, а сразу раскрою оптимальный способ — минимаксная аппроксимация. Один полином от другого может отличаться только коэффициентами, и при минимаксной аппроксимации коэффициенты выбираются так, чтобы минимизировать максимальную абсолютную погрешность на всем интервале аппроксимации. Для этого нам даже не нужно ставить оптимизационные задачи и решать их. Алгоритм расчета коэффициентов минимаксной аппроксимации давно придуман — большое спасибо Чебышёву, Паде, Ремезу и другим классикам теории аппроксимации.
Не буду описывать алгоритм здесь, чтобы статья не превратилась в конспект лекции по методам вычислений, но приведу график ошибки минимаксной аппроксимации синуса пятого порядка:
Видно, что по сравнению с аппроксимацией рядом Тейлора того же порядка, наблюдается совершенно другая картина: нет резких увеличений ошибки на краях отрезка, на всем интервале аппроксимации мы видим «лепестки» на графике ошибки. Но главное — абсолютная ошибка на рассматриваемом интервале при использовании минимаксной аппроксимации оказывается в 10 тысяч раз меньше, чем в случае приближения рядом Тейлора того же порядка. За это ее и любят. Если в библиотеках используются другие способы аппроксимации, то, вероятно, возникнут проблемы с точностью.
Мы как-то реализовали гиперболический тангенс при помощи минимаксной аппроксимации, и, в отличие от glibc, у нас не возникло ошибок ни в 52-м бите, ни в единой из тех 37,5 тысяч точек. Реализация tanh на основе минимаксной аппроксимации еще и работает на 70% быстрее, чем неточный аналог из glibc. Это связано с тем, что мы можем взять меньше членов аппроксимирующего полинома, не потеряв при этом в точности.
Итак, мы выяснили, как лучше аппроксимировать функции, чтобы реализации были высокоточными. Теперь перейдем к не менее важной теме: как мы эту точность проверяем. Чтобы понять это, прежде всего нужно вспомнить, что же представляют собой числа с плавающей точкой.
Как выглядят числа с плавающей точкой?
Для начала рассмотрим числа одинарной точности (float, FP32). Для хранения каждого такого числа используется 32 бита: 1 бит отвечает за кодирование знака числа (sign), 8 бит отводится на порядок (exponent) и 23 — на мантиссу (fraction). Пронумеруем биты с нуля справа налево (от мантиссы к знаку) и обозначим их значения (0 или 1) как bi, где i — номер бита. Тогда — знак,
— порядок,
— мантисса в двоичной системе счисления, что подчеркивается нижним индексом 2. Значение всего кодируемого числа в десятичной системе счисления рассчитывается по формуле:
Число 127 в формуле — это нулевое смещение (bias), которое рассчитывается по числу бит, отводящихся на порядок:
Переведем в десятичную систему счисления те части формулы, которые выше записаны в двоичной:
Теперь вместо двоичных чисел формула содержит суммы битовых значений, умноженных на 2 в некоторой степени. Для бит, составляющих мантиссу, эти степени отрицательны, то есть вся дробная часть числа заключена именно здесь.
Поясню формулу на примере. При помощи онлайн-калькулятора посмотрим, как в 32-битном представлении выглядит число 12,5:
Чекбоксы в строке Binary обозначают биты. Если в чекбоксе стоит галочка, значит, этот бит равен единице, иначе — нулю. Знак (–1)0 = 1 — это положительное число. Порядок 130 за вычетом нулевого смещения 127 дает нам 3, а 23 = 8. В мантиссе закодировано десятичное число 0,5625, можете проверить это сами по формуле выше. В итоге получаем: 8 × (1 + 0,5625) = 12,5.
Итак, за дробную часть числа отвечает мантисса. Как же она изменится, если вместо чисел одинарной точности мы будем использовать числа двойной точности (double, FP64)? Вместо 32 бит у нас уже 64 бита, и на мантиссу из них отводится вместо 23 целых 52 бита! Число бит, которые отводятся на мантиссу, увеличивается более чем в 2 раза, и формула для расчета десятичного значения числа меняется соответствующим образом:
Из формул расчета десятичных значений закодированных чисел видно, что в десятичной системе счисления они распределены по числовой оси неравномерно.
Как же измерить точность?
Числа расположены неравномерно, значит, привычная глазу линейка для оценки расстояний между значениями не работает. У нас есть результат вычисления математической функции, число с плавающей точкой, и есть «эталон» — это ожидаемый результат в квазибесконечной точности. Но как понять, насколько велика погрешность вычисления, расстояние между ними?
Для этого достаточно договориться о единице измерения. Расстояние между соседними числами обозначается как 1 ульп (ulp — unit in the last place). Относительно него и будем оценивать погрешность вычисления математической функции. Поделим расстояние от результата до эталона на то, что является одним ульпом — то есть на расстояние от эталона до соседнего числа той же точности. Стандарт libm требует, чтобы ошибка не превышала 0,5 ульпа с учетом округления.
Мы договорились о единицах измерения. Но остался еще один вопрос: с чем же мы сравниваем результаты?
Откуда брать эталон в квазибесконечной точности
Для повседневных расчетов мы часто пользуемся калькулятором. Наверняка вы замечали, что у калькуляторов на разных устройствах есть несколько режимов: обычный, инженерный, научный и так далее. Можно предположить, что и для расчета эталона с какой-то невероятной точностью существует некий продвинутый калькулятор. И такие действительно есть — это системы компьютерной алгебры, прикладные программы для символьных вычислений и числовых операций произвольной точности.
Ученые особенно любят Maple или Scilab, инженеры — Mathcad или Matlab, а разработчики — Sollya, поскольку эта библиотека имеет удобный C-интерфейс и ее можно вызывать прямо из тестов libm. Для кратного повышения точности внутри Sollya используется надежная библиотека GNU MPFR (Multiple Precision Floating-Point Reliable Library) с менее удобным интерфейсом. Она позволяет выбирать степень повышения точности относительно ширины типа данных. Для наших целей достаточно восьмикратной точности.
Как правильно тестировать точность libm
Прежде всего нужно не повторять чужих ошибок и генерировать случайные исходные данные вместо использования постоянного датасета. Также необходимо обеспечить равномерное, с точностью до устройства числа, распределение чисел. Так можно получить стопроцентное покрытие всех веток кода математической функции за тест, а в пределе, при условии регулярных ночных запусков, через некоторое время можно будет протестировать всю числовую ось. В виде псевдокода для 32-битных чисел это будет выглядеть следующим образом:
for (sign = 0; sign < 2; sign++)
{
for (exponent = e_min; exponent < e_max; exponent += e_step)
{
for (fraction = f_min; fraction< f_max; fraction += f_step)
{
value = (sign << 31) | (exponent << 23) | fraction
TEST(value)
}
}
}Для точки value необходимо запустить тестируемую реализацию математической функции, получить результат и передать его в Sollya, где вычисляется эталон и разность между ним и переданным значением. Если эта разность меньше половины расстояния между соседними числами рассматриваемой точности, 0,5 ульпа, то результат вычисления функции признается точным, а тест — пройденным.
Удобно использовать графическое представление результатов тестирования: для каждой функции строить графики рассчитанных значений, абсолютной ошибки, относительной ошибки и гистограмму распределения ошибок в точках по битам.
Графики выше соответствуют реализации экспоненты, удовлетворяющей стандартам точности для float. Мы воспользовались свойством масштабируемости экспоненты, и с некоторой периодичностью аппроксимация повторяется: мы явно это видим на графике относительной ошибки.
Но что, если мы неверно построили аппроксимацию? Такой пример представлен на графиках ниже: мы сразу замечаем, на каких отрезках аппроксимация испорчена, видим порядок этих ошибок и биты, содержащие ошибки. В результате можно достаточно быстро определить, в каких конкретно точках мы ошиблись, на какие ветки кода стоит обратить внимание и где необходимо перестроить аппроксимацию.
Подведем итоги
Мы разобрались, как можно построить систему тестирования точности реализаций математических функций без ограничений по фреймворкам так, чтобы в пределе обеспечить стопроцентное покрытие значений тестами. Узнали, какие аппроксимации позволяют строить высокопроизводительные реализации функций libm, удовлетворяющие стандартам точности (0,5 ульпа): на графике ниже показано, какое ускорение позволяют получить минимаксы по сравнению с реализациями из glibc и AOSP Bionic.
Думаю, даже тем, кто не был знаком с методами вычислений, стало понятно, что численная неустойчивость — штука коварная. Числа с плавающей точкой распределены неравномерно, законы арифметики нарушаются, и важен даже порядок вычислений. Не забывайте об этом при работе с математическими функциями и не доверяйте чересчур сильно прогнозам погоды, если не уверены, какую реализацию экспоненты они используют.
Использование численно неустойчивых алгоритмов может приводить к падениям миллионов тестов. Подробности — в статье моего коллеги Андрея Соколова.