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

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

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

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

Алгоритм

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

function KahanSum(input)
    // Prepare the accumulator.
    var sum = 0.0
    // A running compensation for lost low-order bits.
    var c = 0.0
    // The array input has elements indexed input[1] to input[input.length].
    for i = 1 to input.length do
        // c is zero the first time around.
        var y = input[i] - c
        // Alas, sum is big, y small, so low-order digits of y are lost.         
        var t = sum + y
        // (t - sum) cancels the high-order part of y;
        // subtracting y recovers negative (low part of y)
        c = (t - sum) - y
        // Algebraically, c should always be zero. Beware
        // overly-aggressive optimizing compilers!
        sum = t
    // Next time around, the lost low part will be added to y in a fresh attempt.
    next i

    return sum

Этот алгоритм также можно переписать для использования алгоритма Fast2Sum :

function KahanSum2(input)
    // Prepare the accumulator.
    var sum = 0.0
    // A running compensation for lost low-order bits.
    var c = 0.0
    // The array input has elements indexed 
    for i = 1 to input.length do
        // c is zero the first time around.
        var y = input[i] + c
        // sum + c is an approximation to the exact sum.
        (sum,c) = Fast2Sum(sum,y)
    // Next time around, the lost low part will be added to y in a fresh attempt.
    next i
    return sum

Рабочий пример

Алгоритм не требует какого-либо конкретного выбора системы счисления , только для арифметики «нормализации сумм с плавающей запятой перед округлением или усечением». Компьютеры обычно используют двоичную арифметику, но, чтобы пример было легче читать, он будет дан в десятичном формате. Предположим, мы используем шестизначную десятичную арифметику с плавающей запятой , 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 = input[i] - c
  t = 10000.0 + 3.14159             t = sum + y
    = 10003.14159                   Normalization done, next round off to six digits.
    = 10003.1                       Few digits from input[i] met those of sum. Many digits have been lost!
  c = (10003.1 - 10000.0) - 3.14159 c = (t - sum) - y  (Note: Parenthesis must be evaluated first!)
    = 3.10000 - 3.14159             The assimilated part of y minus the original full y.
    = -0.0415900                    Because c is close to zero, normalization retains many digits after the floating point. 
sum = 10003.1                       sum = t

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

  y = 2.71828 - (-0.0415900)        Most digits meet, since c is of a size similar to y.
    = 2.75987                       The shortfall (low-order digits lost) of previous iteration successfully reinstated.
  t = 10003.1 + 2.75987             But still only few meet the digits of sum.
    = 10005.85987                   Normalization done, next round to six digits.
    = 10005.9                       Again, many digits have been lost, but c helped nudge the round-off.
  c = (10005.9 - 10003.1) - 2.75987 Estimate the accumulated error, based on the adjusted y.
    = 2.80000 - 2.75987             As expected, the low-order parts can be retained in c with no or minor round-off effects.
    = 0.0401300                     In this iteration, t was a bit too high, the excess will be subtracted off in next iteration.
sum = 10005.9                       Exact result is 10005.85987, sum is correct, rounded to 6 digits.

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

Точность

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

Предположим, что кто-то суммирует $n$ ценности $x_{i}$, для $i=1,\,\ldots ,\,n$. Точная сумма

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

При компенсированном суммировании вместо этого получаем $S_{n}+E_{n}$, где ошибка $E_{n}$ ограничен

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

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

${\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|}}.$

В выражении для границы относительной ошибки дробь $\Sigma |x_{i}|/|\Sigma x_{i}|$ – число обусловленности задачи суммирования. По сути, число обусловленности представляет собой внутреннюю чувствительность задачи суммирования к ошибкам, независимо от того, как оно вычисляется. Относительная граница ошибки для каждого ( обратно стабильного ) метода суммирования с использованием фиксированного алгоритма с фиксированной точностью (т.е. не тех, которые используют арифметику произвольной точности или алгоритмов, требования к памяти и времени которых меняются в зависимости от данных), пропорциональна этому числу обусловленности. . Плохо обусловленная задача суммирования — это задача, в которой это отношение велико, и в этом случае даже компенсированное суммирование может иметь большую относительную ошибку. Например, если слагаемые $x_{i}$ являются некоррелированными случайными числами с нулевым средним значением, сумма представляет собой случайное блуждание , а число обусловленности будет расти пропорционально ${\sqrt {n}}$. С другой стороны, для случайных входных данных с ненулевым значением число обусловленности асимптотируется до конечной константы как $n\to \infty$. Если все входные данные неотрицательны , то номер условия равен 1.

Учитывая число обусловленности, относительная ошибка компенсированного суммирования фактически не зависит от $n$. В принципе, есть $O(n\varepsilon ^{2})$ который растет линейно с $n$, но на практике этот член фактически равен нулю: поскольку конечный результат округляется до точности $\varepsilon$, $n\varepsilon ^{2}$ срок округляется до нуля, если только $n$ примерно $1/\varepsilon$ или больше. В двойной точности это соответствует $n$ примерно $10^{16}$, намного больше, чем большинство сумм. Итак, при фиксированном числе обусловленности ошибки компенсированного суммирования фактически равны $O(\varepsilon )$, независимо от $n$.

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

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

Дальнейшие улучшения

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

function KahanBabushkaNeumaierSum(input)
    var sum = 0.0
    var c = 0.0                       // A running compensation for lost low-order bits.

    for i = 1 to input.length do
        var t = sum + input[i]
        if |sum| >= |input[i]| then
            c += (sum - t) + input[i] // If sum is bigger, low-order digits of input[i] are lost.
        else
            c += (input[i] - t) + sum // Else low-order digits of sum are lost.
        endif
        sum = t
    next i

    return sum + c                    // Correction only applied once in the very end.

Это усовершенствование аналогично замене Fast2Sum на 2Sum в алгоритме Кахана, переписанном с помощью Fast2Sum.

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

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

function KahanBabushkaKleinSum(input)
    var sum = 0.0
    var cs  = 0.0
    var ccs = 0.0
    var c   = 0.0
    var cc  = 0.0

    for i = 1 to input.length do
        var t = sum + input[i]
        if |sum| >= |input[i]| then
            c = (sum - t) + input[i]
        else
            c = (input[i] - t) + sum
        endif
        sum = t
        t = cs + c
        if |cs| >= |c| then
            cc = (cs - t) + c
        else
            cc = (c - t) + cs
        endif
        cs = t
        ccs = ccs + cc
    end loop 

    return sum + cs + ccs

Альтернативы

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

Другая альтернатива — использовать арифметику произвольной точности , которая в принципе вообще не требует округления и требует гораздо больших вычислительных усилий. Способ правильного округления сумм с произвольной точностью состоит в адаптивном расширении с использованием нескольких компонентов с плавающей запятой. Это сведет к минимуму вычислительные затраты в обычных случаях, когда высокая точность не требуется. Другой метод, использующий только целочисленную арифметику, но большой аккумулятор, был описан Киршнером и Кулишом ; аппаратная реализация была описана Мюллером, Рюбом и Рюллингом.

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

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

t = sum + y;

c = (t - sum) - y;

к

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

а затем

c = 0;

тем самым устраняя компенсацию ошибки. На практике многие компиляторы не используют правила ассоциативности (которые являются лишь приблизительными в арифметике с плавающей запятой) при упрощениях, если только это явно не указано в параметрах компилятора, разрешающих «небезопасные» оптимизации. хотя компилятор Intel C++ является одним из примеров, который по умолчанию допускает преобразования на основе ассоциативности. Первоначальная K&R C версия языка программирования C позволяла компилятору переупорядочивать выражения с плавающей запятой в соответствии с правилами ассоциативности вещественной арифметики, но последующий стандарт ANSI C запрещал переупорядочение, чтобы сделать C более подходящим для числовых приложений ( и больше похоже на Фортран , который тоже запрещает переупорядочение), хотя на практике параметры компилятора могут снова включить переупорядочение, как упоминалось выше.

Портативный способ локально запретить такую ​​оптимизацию — разбить одну из строк исходной формулировки на два утверждения и сделать два промежуточных продукта нестабильными :

function KahanSum(input)
    var sum = 0.0
    var c = 0.0

    for i = 1 to input.length do
        var y = input[i] - c
        volatile var t = sum + y
        volatile var z = t - sum
        c = z - y
        sum = t
    next i

    return sum

Поддержка библиотеками

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

Стандартная библиотека компьютерного языка Python определяет функцию fsum для точного суммирования. Начиная с Python 3.12, встроенная функция sum() использует суммирование Ноймайера.

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

На C# языке пакет Nuget HPCsharp реализует вариант Ноймайера и попарное суммирование : как скалярное, так и параллельное по данным с использованием инструкций процессора SIMD , и параллельное многоядерное суммирование.

См. также

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

Внешние ссылки

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