Допустим, у нас есть массив чисел с плавающей запятой, и мы хотим их суммировать. Можно наивно подумать, что их достаточно просто сложить, например, на Rust:
fn naive_sum(arr: &[f32]) -> f32 {
let mut out = 0.0;
for x in arr {
out += *x;
}
out
}Однако это запросто может привести к произвольно большой накопленной погрешности. Давайте проверим:
naive_sum(&vec![1.0; 1_000_000]) = 1000000.0
naive_sum(&vec![1.0; 10_000_000]) = 10000000.0
naive_sum(&vec![1.0; 100_000_000]) = 16777216.0
naive_sum(&vec![1.0; 1_000_000_000]) = 16777216.0Ой-ёй… Что произошло? Проблема в том .что следующее 32-битное число с плавающей запятой после 16777216 — это 16777218. Так что при вычислении 16777216 + 1, значение округляется до ближайшего числа с плавающей запятой, имеющей чётную мантиссу, то есть снова до 16777216. Мы зашли в тупик.
К счастью, есть более совершенные способы суммирования массива.
Попарное суммирование
Чуть более умный способ — это попарное суммирование. Вместо абсолютно линейной суммы с единственным накопителем он рекурсивно суммирует массив, разбивая его пополам, суммируя половины, а затем складывая суммы.
fn pairwise_sum(arr: &[f32]) -> f32 {
if arr.len() == 0 { return 0.0; }
if arr.len() == 1 { return arr[0]; }
let (first, second) = arr.split_at(arr.len() / 2);
pairwise_sum(first) + pairwise_sum(second)
}Этот способ более точен:
pairwise_sum(&vec![1.0; 1_000_000]) = 1000000.0
pairwise_sum(&vec![1.0; 10_000_000]) = 10000000.0
pairwise_sum(&vec![1.0; 100_000_000]) = 100000000.0
pairwise_sum(&vec![1.0; 1_000_000_000]) = 1000000000.0Однако это довольно медленно. Чтобы написать процедуру суммирования и максимально быстро работающую, и дающую достаточную точность, мы не должны выполнять рекурсию плоть до массива длиной 1, при этом тратится слишком много ресурсов на вызовы. Можно продолжать использовать наивное суммирование для малых размеров, а рекурсию применять только при больших размерах. При этом наихудший случай будет хуже на постоянный коэффициент, но, в свою очередь, сделает попарное суммирование почти таким же быстрым, как наивное.
Выбрав точку разделения как кратное 256, мы гарантируем, что базовый случай в рекурсии всегда будет равен точно 256 элементам, за исключением самого последнего блока. Это гарантирует использование самого оптимального редуцирования и всегда корректно прогнозирует условие цикла. Эта небольшая деталь повысила производительность в случае больших массивов на 40%!
fn block_pairwise_sum(arr: &[f32]) -> f32 {
if arr.len() > 256 {
let split = (arr.len() / 2).next_multiple_of(256);
let (first, second) = arr.split_at(split);
block_pairwise_sum(first) + block_pairwise_sum(second)
} else {
naive_sum(arr)
}
}Суммирование Кэхэна
Наихудший случай погрешности округления наивного сумирования масштабируется при суммировании 𝑛 элементов как 𝑂(𝑛𝜖), где 𝜖 — это машинный ноль вашего типа с плавающей запятой (здесь это 2−24). Показатели попарного суммирования улучшают этот результат до 𝑂((log𝑛)𝜖+𝑛𝜖2). Суммирование Кэхэна ещё больше улучшает показатель (до 𝑂(𝑛𝜖2)), полностью избавляя нас от члена 𝜖 и оставляя только 𝜖2 , который несущественен, если только вы не суммируете очень большое количество чисел.
Все эти границы масштабируются как ∑𝑖∣𝑥𝑖∣, то есть для суммирования Кэхена граница абсолютной погрешности в наихудшем случае всё равно остаётся квадратичной относительно 𝑛.
На практике показатели всех алгоритмов суммирования существенно лучше, чем границы их наихудших случаев, так как в большинстве ситуаций погрешности округляются не только строго вверх или вниз, а в среднем обнуляют друг друга.
pub fn kahan_sum(arr: &[f32]) -> f32 {
let mut sum = 0.0;
let mut c = 0.0;
for x in arr {
let y = *x - c;
let t = sum + y;
c = (t - sum) - y;
sum = t;
}
sum
}Суммирование Кэхэна заключается в хранении суммы в двух регистрах: настоящей общей суммы и небольшого значения коррекции погрешности 𝑐. Если бы мы использовали бесконечно точную арифметику, то 𝑐 всегда бы был нулём, но при плавающей запятой такого быть не может. Недостаток такого подхода заключается в том, что для добавления каждого числа к сумме теперь требуется не одна, а целых четыре операции.
Чтобы решить эту проблему, мы можем сделать нечто подобное тому, что проделали с попарным суммированием. Сначала можно накапливать блоки в суммы наивным образом, а потом комбинировать суммы блоков при помощи суммирования Кэхэна, чтобы снизить лишнюю трату ресурсов ценой точности:
pub fn block_kahan_sum(arr: &[f32]) -> f32 {
let mut sum = 0.0;
let mut c = 0.0;
for chunk in arr.chunks(256) {
let x = naive_sum(chunk);
let y = x - c;
let t = sum + y;
c = (t - sum) - y;
sum = t;
}
sum
}Точное суммирование
Я знаю по крайней мере два обобщённых способа получения корректно округлённой суммы последовательности чисел с плавающей запятой. Они логически вычисляют сумму с бесконечной точностью, а в конце округляют её до значения с плавающей запятой.
Первый способ основан на примитиве 2Sum, то есть на избавленном от погрешностей преобразовании из двух чисел 𝑥,𝑦 в 𝑠,𝑡 так, что 𝑥+𝑦=𝑠+𝑡, где 𝑡 — это малая погрешность. Многократно применяя это преобразование, пока погрешность не пропадёт, мы можем получить корректно округлённую сумму. Сложность заключается в отслеживании того, что и в каком порядке складывать, а для исчезновения всех членов в наихудшем случае требуется 𝑂(𝑛2) сложений. Именно это реализовано в math.fsum языка Python и в крейте Rust fsum, которые используют дополнительную память для хранения частичных сумм. Крейт accurate тоже реализует это при помощи преобразования на месте в i_fast_sum_in_place.
Другой способ заключается в хранении большого буфера целых чисел, по одному на экспоненту. В случае прибавления числа с плавающей запятой оно раскладывается на экспоненту и мантиссу, а мантисса прибавляется к соответствующему целому числу в буфере. Если происходит переполнение целого числа buf[i], то выполняется инкремент целого числа в buf[i + w], где w — ширина целого числа.
На самом деле, этот способ позволяет вычислять абсолютно точную сумму без какого-либо округления; по сути, это более свободное представление числа с фиксированной запятой, оптимизированное для накопления чисел с плавающей запятой. Этот способ требует времени 𝑂(𝑛), но использует большой, хотя и константный объём памяти (≈≈ 1 КБ в случае f32, ≈≈ 16 КБ в случае f64). Преимущество этого способа ещё и в том, что это онлайновый алгоритм — и прибавление числа к сумме, и получение текущего общего значения занимают амортизированное 𝑂(1).
Вариант этого способа реализован в крейте accurate как крейт OnlineExactSum, в котором для буфера используются числа с плавающей запятой, а не целые.
Раскрываем мощь компилятора
Кроме точности, у naive_sum есть ещё одна проблема. Компилятор Rust не может изменять порядок сложения чисел с плавающей запятой, потому что такое сложение не ассоциативно. Поэтому он не может автоматически векторизировать naive_sum, чтобы та использовала для вычисления суммы команды SIMD, и не может применять параллелизм на уровне команд.
Для решения этой проблемы в компиляторе Rust есть встроенные средства (intrinsics), выполняющие суммирование чисел с плавающей запятой, в то же время обеспечивая ассоциативность, например, std::intrinsics::fadd_fast. Однако такие команды невероятно опасны, потому что они предполагают, что и входные, и выходные значения являются конечными числами (не бесконечностями и не NaN), в противном случае они оказываются неопределённым поведением. Функционально это делает их неприменимыми, поскольку только в самых ограниченных ситуациях при вычислениях суммы мы знаем, что все входные числа являются конечными и что их сумма не может привести к переполнению.
Недавно я поделился своим разочарованием этими операторами с Беном Кимоком, и вместе мы предложили (а он реализовал) новое множество операторов: std::intrinsics::fadd_algebraic и ему подобные.
Я предложил, чтобы мы назвали эти операторы алгебраическими (algebraic), потому что они позволяют (теоретически) выполнять любые преобразования, допустимые в вещественной алгебре. Например, подстановки 𝑥−𝑥→0, 𝑐𝑥+𝑐𝑦→𝑐(𝑥+𝑦) или 𝑥6→(𝑥2)3. В общем случае эти операторы обрабатываются, как будто они выполняются с использованием вещественных чисел, и их можно отобразить на любое множество команд с плавающей запятой, которые будут эквивалентом исходного выражения при условии, что команды с плавающей запятой будут точными.
Стоит отметить, что вещественные числа не содержат NaN и бесконечностей, так что для валидности преобразований эти операторы предполагают их отсутствие, однако если такие значения встретятся, они не вызовут неопределённого поведения.
Кроме того, они позволяют генерировать команды fused multiply-add , так как в вещественной алгебре fma(𝑎,𝑏,𝑐)=𝑎𝑏+𝑐.
Благодаря этим новым командам можно очень легко сгенерировать автоматически векторизированную сумму:
#![allow(internal_features)]
#![feature(core_intrinsics)]
use std::intrinsics::fadd_algebraic;
fn naive_sum_autovec(arr: &[f32]) -> f32 {
let mut out = 0.0;
for x in arr {
out = fadd_algebraic(out, *x);
}
out
}Если выполнить компиляцию с -C target-cpu=broadwell , то мы увидим, что компилятор автоматически сгенерировал следующий короткий цикл с четырьмя накопителями и командами AVX2:
.LBB0_5:
vaddps ymm0, ymm0, ymmword ptr [rdi + 4*r8]
vaddps ymm1, ymm1, ymmword ptr [rdi + 4*r8 + 32]
vaddps ymm2, ymm2, ymmword ptr [rdi + 4*r8 + 64]
vaddps ymm3, ymm3, ymmword ptr [rdi + 4*r8 + 96]
add r8, 32
cmp rdx, r8
jne .LBB0_5Он обрабатывает 128 байтов данных с плавающей запятой (то есть 32 элемента) за 7 команд. Кроме того, все команды vaddps не зависят друг от друга, потому что накопление выполняется в разные регистры. Если проанализировать код при помощи uiCA, то можно увидеть. что по его оценке выполнение этого цикла занимает 4 такта, то есть обрабатывается по 32 байта на цикл. При частоте 4 ГГц это до 128 ГБ/с! Стоит отметить, что это гораздо больше, чем пропускная способность ОЗУ моей машины, так что этой скорости можно достичь, только суммируя данные, уже находящиеся в кэше.
Помня об этом, мы также может легко определить block_pairwise_sum_autovec и block_kahan_sum_autovec, заменив их вызовы к naive_sum на naive_sum_autovec.
Точность и скорость
Давайте попробуем сравнить разные способы суммирования. Попробуем в качестве бенчмарка просуммировать 100000 случайных float в интервале от -100000 до +100,000. Это 400 КБ данных, так что всё равно уместится в кэш моего AMD Threadripper 2950x.
Весь код выложен на Github. Компиляция выполнена с RUSTFLAGS=-C target-cpu=native и --release. Я получил следующие результаты:
Алгоритм | Производительность | Средняя абсолютная погрешность |
|---|---|---|
| 5.5 ГБ/с | 71.796 |
| 0.9 ГБ/с | 1.5528 |
| 1.4 ГБ/с | 0.2229 |
| 5.8 ГБ/с | 3.8597 |
| 5.9 ГБ/с | 4.2184 |
| 118.6 ГБ/с | 14.538 |
| 71.7 ГБ/с | 1.6132 |
| 98.0 ГБ/с | 1.2306 |
| 1.1 ГБ/с | 0.0015 |
| 1.9 ГБ/с | 0.0015 |
| 1.2 ГБ/с | 0.0000 |
Крейт accurate имеет ненулевую абсолютную погрешность потому, что пока он не реализует корректно округление до ближайшего, поэтому окончательный результат может быть смещён на одну единицу в последнем знаке.
Во-первых, хотелось бы отметить, что разница производительности между самым быстрым и самым медленным способом составляет более чем 100 раз. И это для суммирования массива! Это может быть не совсем честно, потому что самые медленные способы вычисляют нечто существенно более сложное, но всё равно существует разница в двадцать раз между достаточно разумной наивной реализацией и самым быстрым способом.
Мы обнаружили, что в общем случае способы _autovec, использующие fadd_algebraic, быстрее и точнее, чем те, которые используют обычное сложение с плавающей запятой. Они точнее по той же причине, по которой точнее попарная сумма: любое изменение порядка сложений лучше, поскольку используемая по умолчанию длинная цепочка сложений уже является наихудшим случаем точности точности в сумме.
Ограничившись Парето-оптимальными вариантами, мы получим следующие четыре реализации:
Алгоритм | Производительность | Средняя абсолютная погрешность |
|---|---|---|
| 118.6 ГБ/с | 14.538 |
| 98.0 ГБ/с | 1.2306 |
| 1.9 ГБ/с | 0.0015 |
| 1.2 ГБ/с | 0.0000 |
Отметим, что разницы в реализации могут влиять достаточно сильно; к тому же, вероятно, существуют ещё десятки способов компенсированного суммирования, которые я не сравнивал.
Думаю, в большинстве случаев здесь выигрывает block_kahan_autovec, он имеет хорошую точность (не происходит вырождение при больших входных значениях) и почти максимальную скорость. В большинстве сценариев использования дополнительная точность, получаемая благодаря корректному округлению сумм, не требуется, а она в 50-100 раз медленнее.
Разбив цикл на явный остаток плюс короткий цикл сумм по 256 элементов, мы можем ещё немного повысить производительность и избежать пары операций с плавающей запятой для последнего блока:
#![allow(internal_features)]
#![feature(core_intrinsics)]
use std::intrinsics::fadd_algebraic;
fn sum_block(arr: &[f32]) -> f32 {
arr.iter().fold(0.0, |x, y| fadd_algebraic(x, *y))
}
pub fn sum_orlp(arr: &[f32]) -> f32 {
let mut chunks = arr.chunks_exact(256);
let mut sum = 0.0;
let mut c = 0.0;
for chunk in &mut chunks {
let y = sum_block(chunk) - c;
let t = sum + y;
c = (t - sum) - y;
sum = t;
}
sum + (sum_block(chunks.remainder()) - c)
}Алгоритм | Производительность | Средняя абсолютная погрешность |
|---|---|---|
| 112.2 ГБ/с | 1.2306 |
Разумеется, можно попробовать менять 256 на другие значения; я выяснил, что при 128 производительность снижается примерно на 20%, а при 512 производительность особо не повышается, но страдает точность.
Заключение
Я считаю, что fadd_algebraic и схожие алгебраические intrinsics очень полезны в достижении высокоскоростных процедур работы с плавающей запятой, и что их нужно добавить и в другие языки. Глобальный -ffast-math недостаточно хорош, так как мы уже видели, что наилучшая реализация оказалась гибридом между автоматически оптимизированными для скорости вычислениями и вручную реализованными неассоциативными компенсированными операциями.
Кроме того, при использовании LLVM следует опасаться -ffast-math. Когда этот флаг установлен в LLVM, получение NaN или бесконечности становится неопределённым поведением. Понятия не имею, почему разработчики выбрали такую жёсткую позицию, делающую ненадёжной практически все использующие его программы. Если вы используете в качестве целевой платформы для своего языка LLVM, то избегайте флагов fast-math nnan и ninf.
