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

В обчислювальній математиці алгоритм Кехена (також відомий, як компенсаційне підсумовування[1]) — це алгоритм обчислення суми послідовності чисел з рухомою комою, який значно зменшує обчислювальну похибку у порівнянні з наївним підходом (простим послідовним підсумовуванням чисел з заокругленням результату на кожному кроці). Зменшення похибки досягається введенням додаткової змінної для збереження суми похибок.

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

Авторство алгоритму приписують Вільяму Кехену[3].

Алгоритм

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

function KahanSum(input)
    var sum = 0.0
    var c = 0.0                 // Сума похибок.
    for i = 1 to input.length 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             Відновлена частина y, яка на увійшла в t, а не все вхідне y.
  = -.0415900                     Завершуючі нулі показані тому, що це шестирозрядна арифметика.
sum = 10003.1                     Таким чином, не всі розряди з input[i] включені в sum.

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

y = 2.71828 - (-.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             В даному випадку забагато.
  = .040130                       Так чи інакше надлишок буде віднятий наступного разу.
sum = 10005.9                     Точний результат: 10005.85987, що коректно заокруглено до 6 знаків.

Таким чином, додавання відбувається з використанням двох змінних: sum зберігає суму, і c зберігає частину суми, яка не потрапила в sum на наступній ітерації.

Точність

Старанний аналіз похибок при компенсаційному підсумовуванні є необхідним для оцінки точності. Незважаючи на те, що він є більш точним, ніж наївне підсумовування, він все ще може давати великі відносні похибки для погано обумовлених сум.

Припустимо, що виникає потреба підсумувати ${\displaystyle n}$ значень ${\displaystyle x_{i}}$ для ${\displaystyle i=1,2,\ldots ,n}$.

Точна сума (розрахована з нескінченною точністю):

${\displaystyle S_{n}=\sum _{i=1}^{n}x_{i}}$.

З компенсаційним підсумовуванням вона набуває вигляду ${\displaystyle S_{n}+E_{n}}$, де похибка ${\displaystyle E_{n}}$ обмежена виразом[2]:

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

${\displaystyle \varepsilon }$ — машинний епсилон (наприклад, ε≈10^−16 для чисел з рухомою комою подвійної точності).

Як правило, нас цікавить величина відносної похибки ${\displaystyle |E_{n}|/|S_{n}|}$, яка обмежена виразом:

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

У виразі для оцінки відносної похибки дріб ${\displaystyle \sum |x_{i}|/\left|\sum {x_{i}}\right|}$ є числом обумовленості задачі підсумовування. По суті, число обумовленості виражає внутрішню чутливість задачі підсумовування до похибок, незалежно від способу розрахунку[4].

В мовах програмування

В принципі, достатньо агресивні оптимізувальні компілятори можуть перешкоджати реалізації алгоритму Кехена. Наприклад, якщо компілятор спрощує вирази у відповідності до правил асоціативності для дійсних чисел, він може «спростити» другий крок у послідовності t = sum + y; c = (t - sum) - y; до ((sum + y) - sum) - y;, а потім до c = 0;, виключаючи компенсацію похибок[5]. На практиці багато компіляторів не використовують правила асоціативності (які є лише наближеними в арифметиці чисел з рухомою комою) для оптимізації, за винятком явного використання опцій компілятора, які дозволяють «небезпечні» оптимізації[6][7][8][9], хоча компілятор Intel C++ дозволяє за замовчуванням оптимізації, які базуються на правилах асоціативності[10]. Оригінальна K&R версія мови програмування C дозволяла компілятору змінювати порядок операндів у виразах з рухомою комою у відповідності до правил асоціативності, але наступний стандарт ANSI C заборонив зміну порядку операндів для того, щоб зробити мову програмування C краще пристосованою для розробки додатків, які використовують чисельні методи (і більш подібною до мови програмування Fortran, яка забороняла зміну порядку операндів[11]), однак на практиці способи оптимізації залежать від опцій компіляторів.

В загальному, вбудовані функції для підсумовування в мовах програмування, як правило, не гарантують, що алгоритм частинного підсумовування буде використовувати хоча б алгоритм Кехена.

Примітки

1. Strictly, there exist other variants of compensated summation as well: see Higham, Nicholas (2002). Accuracy and Stability of Numerical Algorithms (2 ed). SIAM. с. 110—123.(англ.)
2. Higham, Nicholas J. (1993), The accuracy of floating point summation, SIAM Journal on Scientific Computing, 14 (4): 783—799, doi:10.1137/0914050(англ.)
3. Kahan, William (January 1965), Further remarks on reducing truncation errors, Communications of the ACM, 8 (1): 40, doi:10.1145/363707.363723(англ.)
4. L. N. Trefethen and D. Bau, Numerical Linear Algebra (SIAM: Philadelphia, 1997).(англ.)
5. Goldberg, David (March 1991), What every computer scientist should know about floating-point arithmetic (PDF), ACM Computing Surveys, 23 (1): 5—48, doi:10.1145/103162.103163, архів оригіналу (PDF) за 12 грудня 2019, процитовано 15 липня 2015(англ.)
6. GNU Compiler Collection manual, version 4.4.3: 3.10 Options That Control Optimization [Архівовано 3 березня 2016 у Wayback Machine.], -fassociative-math (Jan. 21, 2010).(англ.)
7. Compaq Fortran User Manual for Tru64 UNIX and Linux Alpha Systems [Архівовано 7 червня 2011 у Wayback Machine.], section 5.9.7 Arithmetic Reordering Optimizations (retrieved March 2010).(англ.)
8. Börje Lindh, Application Performance Optimization [Архівовано 2 червня 2010 у Wayback Machine.], Sun BluePrints OnLine (March 2002).(англ.)
9. Eric Fleegal, "Microsoft Visual C++ Floating-Point Optimization [Архівовано 15 липня 2015 у Wayback Machine.]", Microsoft Visual Studio Technical Articles (June 2004).(англ.)
10. Martyn J. Corden, "Consistency of floating-point results using the Intel compiler [Архівовано 15 липня 2015 у Wayback Machine.]," Intel technical report (Sep. 18, 2009).(англ.)
11. Tom Macdonald, "C for Numerical Computing", Journal of Supercomputing vol. 5, pp. 31–48 (1991).(англ.)