Алгоритм Кехена
В обчислювальній математиці алгоритм Кехена (також відомий, як компенсаційне підсумовування[1]) — це алгоритм обчислення суми послідовності чисел з рухомою комою, який значно зменшує обчислювальну похибку у порівнянні з наївним підходом (простим послідовним підсумовуванням чисел з заокругленням результату на кожному кроці). Зменшення похибки досягається введенням додаткової змінної для збереження суми похибок.
Зокрема, наївне підсумовування n чисел в найгіршому випадку дає похибку, яка росте пропорційно 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 на наступній ітерації.
Точність
Старанний аналіз похибок при компенсаційному підсумовуванні є необхідним для оцінки точності. Незважаючи на те, що він є більш точним, ніж наївне підсумовування, він все ще може давати великі відносні похибки для погано обумовлених сум.
Припустимо, що виникає потреба підсумувати значень для .
Точна сума (розрахована з нескінченною точністю):
.
З компенсаційним підсумовуванням вона набуває вигляду , де похибка обмежена виразом[2]:
,
— машинний епсилон (наприклад, ε≈10−16 для чисел з рухомою комою подвійної точності).
Як правило, нас цікавить величина відносної похибки , яка обмежена виразом:
.
У виразі для оцінки відносної похибки дріб є числом обумовленості задачі підсумовування. По суті, число обумовленості виражає внутрішню чутливість задачі підсумовування до похибок, незалежно від способу розрахунку[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]), однак на практиці способи оптимізації залежать від опцій компіляторів.
В загальному, вбудовані функції для підсумовування в мовах програмування, як правило, не гарантують, що алгоритм частинного підсумовування буде використовувати хоча б алгоритм Кехена.
Примітки
- 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.(англ.)
- 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.(англ.)
- L. N. Trefethen and D. Bau, Numerical Linear Algebra (SIAM: Philadelphia, 1997).(англ.)
- 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.(англ.)
- GNU Compiler Collection manual, version 4.4.3: 3.10 Options That Control Optimization, -fassociative-math (Jan. 21, 2010).(англ.)
- 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).(англ.)
- Börje Lindh, Application Performance Optimization, Sun BluePrints OnLine (March 2002).(англ.)
- Eric Fleegal, "Microsoft Visual C++ Floating-Point Optimization", Microsoft Visual Studio Technical Articles (June 2004).(англ.)
- Martyn J. Corden, "Consistency of floating-point results using the Intel compiler," Intel technical report (Sep. 18, 2009).(англ.)
- Tom Macdonald, "C for Numerical Computing", Journal of Supercomputing vol. 5, pp. 31–48 (1991).(англ.)