Алгоритм Кэхэна

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

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

Авторство алгоритма приписывают Уильяму Кэхэну^[2].

Алгоритм

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

function kahan_sum(input)
    var sum = 0.0
    var c = 0.0             // Сумма погрешностей.
    for i = 1 to len(input) do
        y = input[i] - c    // Пока всё хорошо: c - ноль.
        t = sum + y         // Увы, sum велика, y мало, так что младшие разряды y потеряны.
        c = (t - sum) - y   // (t - sum) восстанавливает старшие разряды y; вычитание y восстанавливает -(младшие разряды y)
        sum = t             // Алгебраически, c всегда должна равняться нулю. Берегитесь слишком оптимизирующих компиляторов!
                            // В следующий раз потерянные младшие разряды будут заново прибавлены к y.
    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 — 0.

  y = 3.14159 - 0                   y = input[i] - c
  t = 10000.0 + 3.14159
    = 10003.1                       Много разрядов потеряно!
  c = (10003.1 - 10000.0) - 3.14159 Это нужно вычислять как записано! 
    = 3.10000 - 3.14159             Восстановлена не вошедшая в t часть y , а не всё исходное y.
    = -0.0415900                    Завершающие нули показаны, потому что это шестиразрядная арифметика.
sum = 10003.1                       Таким образом не все разряды из input(i) включены в sum.

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

  y = 2.71828 - -0.0415900          Учитывается погрешность с предыдущего шага.
    = 2.75987                       Её порядок не слишком отличается от y. Большинство разрядов учтены.
  t = 10003.1 + 2.75987             Но только немногие разряды попадают в sum.
    = 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 уже имеет двойную точность, немногие системы могут предоставить учетверённую точность вычислений, и, если бы могли, ввод тоже мог бы иметь учетверённую точность.

Недостатки

К сожалению, алгоритм Кэхэна не гарантирует защиту от потери точности, связанной с наличием в сумме слагаемых с существенно различными порядками. Это может происходить, если последовательность суммирования выбрана неудачно. Чтобы убедиться в этом, достаточно в рассмотренном выше примере вместо числа 2,71828 взять −10000. На последнем шаге алгоритма имеем следующее:

  y = -10000.0 - -0.0415900 = -10000.0         Результат округляется до 6 значащих цифр
  t = 10003.1 + -10000.0 = 3.10000
  c = (3.10000 - 10003.1) - -10000.0 = -10000.0 + 10000.0 = 0
  sum = 3.10000

Точный результат: 3,14159, то есть произошла потеря точности.

Надо заметить, что если предварительно упорядочить слагаемые по убыванию их абсолютной величины, то в результате применения алгоритма получим результат 3,14159, где все знаки верны.

Примечания

↑ ¹ ² Higham, Nicholas J. (1993), "The accuracy of floating point summation", SIAM Journal on Scientific Computing, 14 (4): 783—799, doi:10.1137/0914050
↑ Kahan, William (January 1965), "Further remarks on reducing truncation errors", Communications of the ACM, 8 (1): 40, doi:10.1145/363707.363723

Ссылки

Evan Manning, Floating-point Summation Dr. Dobb’s Journal September, September 01, 1996

[Higham93-1] ¹ ² Higham, Nicholas J. (1993), "The accuracy of floating point summation", SIAM Journal on Scientific Computing, 14 (4): 783—799, doi:10.1137/0914050

[kahan65-2] Kahan, William (January 1965), "Further remarks on reducing truncation errors", Communications of the ACM, 8 (1): 40, doi:10.1145/363707.363723

[1]

[2]