Алгоритм суммирования Кахана

В численном анализе , то алгоритм суммирования Кагана , также известный как компенсированным суммирование , ^[1] значительно уменьшает числовую ошибку в общей сложности, полученного путем добавления последовательности конечнозонной точности чисел с плавающей точкой , по сравнению с очевидным подходом. Для этого используется отдельная компенсация хода (переменная для накопления небольших ошибок).

В частности, простое последовательное суммирование n чисел имеет ошибку наихудшего случая, которая растет пропорционально n , и среднеквадратичную ошибку, которая растет, как и для случайных входных данных (ошибки округления образуют случайное блуждание ). ^[2] При скомпенсированном суммировании граница ошибки наихудшего случая фактически не зависит от n , поэтому большое количество значений может быть суммировано с ошибкой, которая зависит только от точности с плавающей запятой . ^[2]${\displaystyle {\sqrt {n}}}$

Алгоритм приписывается Кэхэн . ^[3] Аналогичными более ранними методами являются, например, линейный алгоритм Брезенхэма , отслеживающий накопленную ошибку в целочисленных операциях (хотя впервые задокументированную примерно в то же время ^[4] ) и дельта-сигма модуляцию ^[5]

Алгоритм [ править ]

В псевдокоде алгоритм будет;

function KahanSum (input) var sum = 0.0 // Подготавливаем аккумулятор. var c = 0.0 // Текущая компенсация потерянных младших битов. for i = 1 to input.length do // Вход массива имеет индексированные элементы input [1] для input [input.length]. var y = input [i] - c // c равно нулю в первый раз. var t = sum + y // Увы, сумма большая, y маленькая, поэтому младшие цифры y теряются. c = (t - sum) - y // (t - sum) отменяет старшую часть y ; вычитание y восстанавливает отрицательное значение (младшая часть y ) sum = t // Алгебраически c всегда должно быть равно нулю. Остерегайтесь чрезмерно агрессивных оптимизирующих компиляторов! next i // В следующий раз потерянная низкая часть будет добавлена ​​к y при новой попытке. сумма возврата

Пример работы [ править ]

Этот пример будет дан в десятичном формате. Компьютеры обычно используют двоичную арифметику, но принцип тот же. Предположим, мы используем шестизначную десятичную арифметику с плавающей запятой, sumдостигли значения 10000.0, а следующие два значения input[i]- 3.14159 и 2.71828. Точный результат - 10005,85987, который округляется до 10005,9. При простом суммировании каждое входящее значение будет выровнено sum, и многие младшие цифры будут потеряны (путем усечения или округления). Первый результат после округления будет 10003,1. Второй результат будет 10005,81828 до округления и 10005,8 после округления. Это не так.

Однако при скомпенсированном суммировании мы получаем правильный округленный результат 10005,9.

Предположим, что cимеет начальное значение ноль.

 y = 3,14159 - 0,00000 y = вход [i] - c т = 10000,0 + 3,14159 = 10003.14159 Но сохраняются только шесть цифр. = 10003.1 Многие цифры потеряны! c = (10003.1 - 10000.0) - 3.14159 Это должно быть вычислено так, как написано! = 3,10000 - 3,14159 Восстановленная ассимилированная часть y по сравнению с исходным полным y . = -0.0415900 Конечные нули показаны, потому что это шестизначная арифметика.sum = 10003.1 Таким образом, несколько цифр из ввода (i ) соответствуют цифрам суммы .

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

 y = 2.71828 - (-0.0415900) Учитывается недостаток предыдущего этапа. = 2.75987 Его размер похож на y : большинство цифр сходятся. t = 10003,1 + 2,75987 Но немногие встречаются цифры суммы . = 10005.85987 И результат округляется = 10005,9 до шести цифр. c = (10005.9 - 10003.1) - 2.75987 Это извлекает все, что вошло. = 2,80000 - 2,75987 В данном случае многовато. = 0,040130 Но неважно, в следующий раз избыток будет вычтен.sum = 10005.9 Точный результат: 10005.85987, округлено до 6 цифр.

Таким образом, суммирование выполняется с помощью двух аккумуляторов: sumудерживает сумму и cнакапливает части, не ассимилированные sum, чтобы подтолкнуть младшую часть sumв следующий раз. Таким образом, суммирование продолжается с cвводом «защитных цифр» , что лучше, чем их отсутствие, но не так хорошо, как выполнение вычислений с удвоенной точностью ввода. Однако простое повышение точности вычислений в целом нецелесообразно; если inputон уже имеет двойную точность, немногие системы обеспечивают четырехкратную точность , а если да, то inputмогут иметь четырехкратную точность.

Точность [ править ]

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

Предположим, что суммируется n значений x _i для i = 1, ..., n . Точная сумма составляет

${\displaystyle S_{n}=\sum _{i=1}^{n}x_{i}}$ (вычислено с бесконечной точностью).

Вместо этого при скомпенсированном суммировании получается , где ошибка ограничена ^[2]${\displaystyle S_{n}+E_{n}}$${\displaystyle E_{n}}$

${\displaystyle |E_{n}|\leq {\big [}2\varepsilon +O(n\varepsilon ^{2}){\big ]}\sum _{i=1}^{n}|x_{i}|,}$

где ε - машинная точность используемой арифметики (например, ε  ≈ 10 ⁻¹⁶ для стандартной двойной точности с плавающей запятой IEEE ). Обычно интересующей величиной является относительная погрешность , которая поэтому ограничена сверху величиной${\displaystyle |E_{n}|/|S_{n}|}$

${\displaystyle {\frac {|E_{n}|}{|S_{n}|}}\leq {\big [}2\varepsilon +O(n\varepsilon ^{2}){\big ]}{\frac {\sum \limits _{i=1}^{n}|x_{i}|}{\left|\sum \limits _{i=1}^{n}x_{i}\right|}}.}$

В выражении для оценки относительной погрешности дробь Σ | x _i | / | Σ x _i | - число обусловленности задачи суммирования. По сути, число обусловленности представляет внутреннюю чувствительность задачи суммирования к ошибкам, независимо от того, как она вычисляется. ^[6] Граница относительной погрешности каждого ( обратно стабильного ) метода суммирования с помощью фиксированного алгоритма с фиксированной точностью (т.е. не тех, которые используют арифметику произвольной точности или алгоритмов, требования к памяти и времени которых изменяются в зависимости от данных), пропорциональна это номер условия. ^[2] Anплохо обусловленная задача суммирования - это проблема, в которой это отношение велико, и в этом случае даже скомпенсированное суммирование может иметь большую относительную ошибку. Например, если слагаемые x _i являются некоррелированными случайными числами с нулевым средним значением, сумма представляет собой случайное блуждание , и число обусловленности будет расти пропорционально . С другой стороны, для случайных входов с ненулевым средним асимптоты числа обусловленности к конечной константе as . Если все входные данные неотрицательны , то номер условия равен 1.${\displaystyle {\sqrt {n}}}$${\displaystyle n\to \infty }$

При заданном числе обусловленности относительная погрешность скомпенсированного суммирования фактически не зависит от n . В принципе, существует O ( nε ² ), который линейно растет с n , но на практике этот член фактически равен нулю: поскольку конечный результат округляется до точности ε , член nε ² округляется до нуля, если n не равно примерно 1. / ε или больше. ^[2] При двойной точности это соответствует примерно 10 ¹⁶ n , что намного больше, чем у большинства сумм. Таким образом, при фиксированном числе обусловленности ошибки скомпенсированного суммирования фактически равны O ( ε), независимо от  n .

Для сравнения, относительная погрешность наивного суммирования (простое последовательное сложение чисел с округлением на каждом шаге) возрастает при умножении на число условия. ^[2] Эта ошибка наихудшего случая редко наблюдается на практике, потому что она возникает только в том случае, если все ошибки округления имеют одинаковое направление. На практике гораздо более вероятно, что ошибки округления имеют случайный знак с нулевым средним значением, так что они образуют случайное блуждание; в этом случае наивное суммирование имеет среднеквадратичную относительную ошибку, которая увеличивается при умножении на число обусловленности. ^[7] Однако это все еще намного хуже, чем компенсированное суммирование. Однако, если суммирование может быть выполнено с удвоенной точностью, то ε заменяется на${\displaystyle O(\varepsilon n)}$${\displaystyle O\left(\varepsilon {\sqrt {n}}\right)}$ε ² , и наивное суммирование имеет ошибку наихудшего случая, сравнимую с членом O ( nε ² ) в компенсированном суммировании с исходной точностью.

Точно так же Σ | х _я | указанное выше является оценкой наихудшего случая, которая возникает только в том случае, если все ошибки округления имеют один и тот же знак (и имеют максимально возможную величину). ^[2] На практике более вероятно, что ошибки имеют случайный знак, и в этом случае члены в Σ | х _я | заменяются случайным блужданием, и в этом случае, даже для случайных входных данных с нулевым средним, ошибка растет только как (игнорируя член nε ² ), с той же скоростью растет сумма , отменяя${\displaystyle E_{n}}$${\displaystyle E_{n}}$${\displaystyle O\left(\varepsilon {\sqrt {n}}\right)}$${\displaystyle S_{n}}$${\displaystyle {\sqrt {n}}}$коэффициенты при вычислении относительной ошибки. Таким образом, даже для асимптотически плохо обусловленных сумм относительная ошибка для скомпенсированного суммирования часто может быть намного меньше, чем можно было бы предположить при анализе наихудшего случая.

Дальнейшие улучшения [ править ]

Ноймайер ^[8] представил улучшенную версию алгоритма Кахана, которую он называет «улучшенным алгоритмом Кахана – Бабушки», который также охватывает случай, когда следующий добавляемый член больше по абсолютной величине, чем текущая сумма, эффективно меняя роль о том, что велико, а что мало. В псевдокоде алгоритм выглядит так:

function NeumaierSum (input) var sum = 0.0 var c = 0.0 // Текущая компенсация потерянных младших битов. для i = 1 для input.length do  var t = sum + input [i] if | sum | > = | input [i] | then c + = (sum - t) + input [i] // Если сумма больше, младшие цифры input [i] теряются. else c + = (input [i] - t) + sum // Остальные младшие цифры суммы теряются. endif сумма = т следующий я return sum + c // Коррекция применяется только один раз в самом конце.

Для многих последовательностей чисел оба алгоритма согласуются, но простой пример Петерса ^[9] показывает, как они могут различаться. Для суммирования с двойной точностью алгоритм Кахана дает 0,0, тогда как алгоритм Ноймайера дает правильное значение 2,0.${\displaystyle [1.0,+10^{100},1.0,-10^{100}]}$

Также возможны модификации более высокого порядка с большей точностью. Например, вариант, предложенный Кляйном ^[10], который он назвал «итерационным алгоритмом Кахана – Бабушки» второго порядка. В псевдокоде алгоритм выглядит так:

функция KleinSum (вход) var sum = 0,0 var cs = 0,0 var ccs = 0,0 var c var cc для i = 1 для input.length do  var t = sum + input [i] if | sum | > = | input [i] | тогда c = (сумма - t) + ввод [i] еще c = (input [i] - t) + sum endif сумма = т т = cs + c если | cs | > = | c | тогда cc = (cs - t) + c еще cc = (c - t) + cs endif cs = t ccs = ccs + cc следующий я сумма возврата + cs + ccs

Альтернативы [ править ]

Хотя алгоритм Кахана обеспечивает рост ошибки при суммировании n чисел, только немного худший рост может быть достигнут путем попарного суммирования : один рекурсивно делит набор чисел на две половины, суммирует каждую половину, а затем складывает две суммы. ^[2] Это имеет то преимущество, что требует того же количества арифметических операций, что и простое суммирование (в отличие от алгоритма Кахана, который требует в четыре раза больше арифметических операций и имеет задержку в четыре раза больше простого суммирования) и может быть вычислен параллельно. Базовым случаем рекурсии в принципе может быть сумма только одного (или нуля) чисел, но для амортизации${\displaystyle O(1)}$${\displaystyle O(\log n)}$накладные расходы на рекурсию, обычно используется более крупный базовый случай. Эквивалент попарного суммирования используется во многих алгоритмах быстрого преобразования Фурье (БПФ) и отвечает за логарифмический рост ошибок округления в этих БПФ. ^[11] На практике при ошибках округления случайных знаков среднеквадратические ошибки попарного суммирования фактически растут как . ^[7]${\displaystyle O\left({\sqrt {\log n}}\right)}$

Другой альтернативой является использование арифметики произвольной точности , которая в принципе не требует округления вообще, что требует гораздо больших вычислительных затрат. Способ получения точно округленных сумм с использованием произвольной точности - адаптивное расширение с использованием нескольких компонентов с плавающей запятой. Это минимизирует вычислительные затраты в типичных случаях, когда не требуется высокая точность. ^{[12] [9]} Другой метод, использующий только целочисленную арифметику, но с большим накопителем, был описан Кирхнером и Кулишем ; ^[13] аппаратная реализация была описана Мюллером, Рубом и Рюллингом. ^[14]

Возможная недействительность при оптимизации компилятора [ править ]

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

t = sum + y;

c = (t - sum) - y;

к

c = ((sum + y) - sum) - y;

а затем в

c = 0;

таким образом устраняется компенсация ошибок. ^[15] На практике многие компиляторы не используют правила ассоциативности (которые являются приблизительными в арифметике с плавающей запятой) в упрощениях, если это явно не указано в параметрах компилятора, включающих «небезопасные» оптимизации, ^{[16] [17] [18 ] [19],} хотя компилятор Intel C ++ является одним из примеров, который по умолчанию разрешает преобразования на основе ассоциативности. ^[20] Исходная версия языка программирования C K&R C позволяла компилятору переупорядочивать выражения с плавающей запятой в соответствии с правилами ассоциативности вещественной арифметики, но последующий ANSI Cстандарт запрещает переупорядочивание, чтобы сделать C более подходящим для числовых приложений (и более похожим на Fortran , который также запрещает переупорядочивание), ^[21] хотя на практике параметры компилятора могут повторно включить переупорядочение, как упоминалось выше.

Поддержка библиотеками [ править ]

В общем, встроенные функции «суммы» в компьютерных языках обычно не дают никаких гарантий того, что будет использован конкретный алгоритм суммирования, не говоря уже о суммировании Кахана. ^{[ необходима цитата ]} Стандарт BLAS для подпрограмм линейной алгебры явно избегает предписания какого-либо определенного порядка вычислений операций по соображениям производительности ^[22], а реализации BLAS обычно не используют суммирование Кахана.

Стандартная библиотека компьютерного языка Python определяет функцию fsum для точно округленного суммирования с использованием алгоритма Шевчука ^[9] для отслеживания нескольких частичных сумм.

В языке Julia реализация sumфункции по умолчанию выполняет попарное суммирование для высокой точности с хорошей производительностью ^[23], но внешняя библиотека предоставляет реализацию варианта Ноймайера, названного sum_kbnдля случаев, когда требуется более высокая точность. ^[24]

В C # языке, HPCsharp NuGet пакет реализует Ньюмаьер вариант , и попарно суммирования : как в качестве скалярной, данных параллельно с использованием- SIMD инструкций процессора, и параллельно многоядерность. ^[25]

См. Также [ править ]

• Алгоритмы расчета дисперсии , включающие устойчивое суммирование

Ссылки [ править ]

Внешние ссылки [ править ]

• Суммирование с плавающей запятой, журнал доктора Добба, сентябрь 1996 г.