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

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

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

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

• 1 Алгоритм
  • 1.1 Рабочий пример
• 2 Точность
• 3 Дальнейшие улучшения
• 4 Альтернативы
• 5 Возможная недействительность при оптимизации компилятора
• 6 Поддержка библиотеками
• 7 См. Также
• 8 Ссылки
• 9 Внешние ссылки

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

function KahanSum (input) var sum = 0.0 // Подготавливаем аккумулятор. var c = 0.0 // Текущая компенсация потерянных младших битов. for i = 1 от до 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 при новой попытке. return sum

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

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

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

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

y = 3,14159 - 0,00000 y = input [i] - c t = 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, что лучше, чем их отсутствие, но не так хорошо, как выполнение вычислений с удвоенной точностью ввода. Однако простое повышение точности вычислений в целом нецелесообразно; если входуже имеет двойную точность, некоторые системы обеспечивают четверную точность, а если да, то вводможет иметь четырехкратную точность.

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

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

$S\; n\; =\; \sum \; i\; =\; 1\; nxi\; \{\backslash \; displaystyle\; S\_\; \{n\}\; =\; \backslash \; sum\; \_\; \{i\; =\; 1\}\; ^\; \{n\}\; x\_\; \{i\}\}$(вычислено с бесконечным точности).

Вместо этого при компенсированном суммировании получается $S\; n\; +\; E\; n\; \{\backslash \; displaystyle\; S\_\; \{n\}\; +\; E\_\; \{n\}\}$, где ошибка $E\; n\; \{\; \backslash \; displaystyle\; E\_\; \{n\}\}$ограничено

$|\; E\; n\; |\; \le \; [2\; \epsilon \; +\; O\; (n\; \epsilon \; 2)]\; \sum \; i\; =\; 1\; n\; |\; х\; я\; |,\; \{\backslash \; displaystyle\; |\; E\_\; \{n\}\; |\; \backslash \; leq\; \{\backslash \; big\; [\}\; 2\; \backslash \; varepsilon\; +\; O\; (n\; \backslash \; varepsilon\; ^\; \{2\})\; \{\backslash \; big]\}\; \backslash \; sum\; \_\; \{i\; =\; 1\}\; ^\; \{n\}\; |\; x\_\; \{i\}\; |,\}$

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

$|\; E\; n\; |\; |\; S\; n\; |\; \le \; [2\; \epsilon \; +\; O\; (n\; \epsilon \; 2)]\; \sum \; i\; =\; 1\; n\; |\; х\; я\; |\; |\; \sum \; i\; =\; 1\; n\; x\; i\; |.\; \{\backslash \; displaystyle\; \{\backslash \; frac\; \{|\; E\_\; \{n\}\; |\}\; \{|\; S\_\; \{n\}\; |\}\}\; \backslash \; leq\; \{\backslash \; big\; [\}\; 2\; \backslash \; varepsilon\; +\; O\; (n\; \backslash \; varepsilon\; ^\; \{2\})\; \{\backslash \; big]\}\; \{\; \backslash \; frac\; \{\backslash \; sum\; \backslash \; limits\; \_\; \{i\; =\; 1\}\; ^\; \{n\}\; |\; x\_\; \{i\}\; |\}\; \{\backslash \; left\; |\; \backslash \; sum\; \backslash \; limits\; \_\; \{i\; =\; 1\}\; ^\; \{n\}\; x\_\; \{i\}\; \backslash \; right\; |\}\; \}.\}$

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

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

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

Точно так же Σ | x i | которое появляется в $E\; n\; \{\backslash \; displaystyle\; E\_\; \{n\}\}$выше, является границей наихудшего случая, которая возникает, только если все ошибки округления имеют один и тот же знак (и имеют максимально возможную величину). На практике более вероятно, что ошибки имеют случайный знак, и в этом случае члены Σ | x i | заменяются случайным блужданием, и в этом случае, даже для случайных входных данных с нулевым средним, ошибка $E\; n\; \{\backslash \; displaystyle\; E\_\; \{n\}\}$растет только как $O\; (\epsilon \; n)\; \{\backslash \; displaystyle\; O\; \backslash \; left\; (\backslash \; varepsilon\; \{\backslash \; sqrt\; \{n\}\}\; \backslash \; right)\}$(игнорируя член nε), с той же скоростью сумма $S\; n\; \{\backslash \; displaystyle\; S\_\; \{n\; \}\}$возрастает, отменяя факторы $n\; \{\backslash \; displaystyle\; \{\backslash \; sqrt\; \{n\}\}\}$при вычислении относительной ошибки. Таким образом, даже для асимптотически плохо обусловленных сумм относительная ошибка для компенсированного суммирования часто может быть намного меньше, чем можно было бы предположить при анализе наихудшего случая.

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

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

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

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

function KleinSum (input) var sum = 0.0 var cs = 0.0 var ccs = 0.0 для i = 1 от до input.length dovar t = sum + input [i] if | sum |>= | input [i] | затем c = (sum - t) + input [i] else c = (input [i] - t) + sum endif sum = tt = cs + c если | cs |>= | c | затем cc = (cs - t) + c else cc = (c - t) + cs endif cs = t ccs = ccs + cc next i return sum + cs + ccs

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

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

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

t = sum + y;

c = ( t - sum) - y; от

до

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

, а затем до

c = 0;

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

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

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

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

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

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

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