Сложение большого множества чисел, существенно отличающихся по величине
Материал из MachineLearning.
Содержание |
Введение
Evan Manning
Задача вычисления суммы большого множества чисел с плавающей точкой часто возникает в практических приложениях. Самые простейшие примеры — численное интегрирование и вычисление среднего значения. Однако, какой бы тривиальной на первый взгляд эта задача не казалась, она не проста, если мы хотим надежно вычислить сумму.
Ниже будем предполагать, что числа, сумму которых нужно найти, записаны в массиве flt_arr. Длина массива ARR_SIZE достаточно большое число.
Тривиальный метод суммирования
Вначале рассмотрим простейший метод вычисления суммы чисел с плавающей точкой.
/* Simple Summation */
float f_add(float * flt_arr)
{
long i;
float sum = 0.0;
for (i = 0; i < ARR_SIZE; i++)
sum += flt_arr[i];
return sum;
}
Чаще всего результат применения данного метода оказывается неудовлетворительным. Проблема заключается в том, что складывание чисел с плавающей точкой сильно отличается от применения той же операции к целым числам. И основное отличие между ними — то, что при сложение вещественных чисел сначала их нужно преобразовать так, чтобы складывались цифры чисел одинакового порядка. При этом, если порядки чисел значительно различаются, тогда цифры младших разрядов отбрасываются. Если необходимо найти сумму сравнительно небольшого количества чисел, например, двух или трех, то данную ошибку можно и не заметить. Но при сложении большого множества чисел существенно отличающихся по величине, ошибка вычисления может оказаться неприемлемой.
Приведем простейший пример</br> Числа в компьютере представляются приближенно, что обусловлено конечностью разрядной сетки (см. Международный стандарт представления чисел с плавающей точкой в ЭВМ). Допустим, что порядок мантиссы не позволяет различить цифры порядка в 10 биллионов. Мы хотим сложить числа:
1.0 + 1.0e10 + −1.0e10
Выполняя суммирования чисел в различном порядке, мы получаем различный результат:
(1.0 + 1.0e10) + −1.0e10 = 1.0e10 + -1.0e10 = 0.0
но
1.0 + (1.0e10 + −1.0e10) = 1.0 + 0.0 = 1.0
В первом случае информация о первом слагаемом была полностью потеряна, поскольку два других слагаемых значительно больше по величине. При суммировании небольшого числа слагаемых ошибку можно и не заметить, но при суммировании большого множества чисел с плавающей точкой ошибка будет накапливаться за счет этих небольших погрешностей, и также за счет разницы между промежуточной суммой и текущем членом.
Можно посоветовать, складывать числа с большей точностью. Если, к примеру, необходимо сложить числа типа float, то результат лучше представлять в виде типа double, аналогично для изначального double — лучше использовать long double.
/* Extended Precision method */
float fx_add(float * flt_arr)
{
long i;
double sum = 0.0;
for (i = 0; i < ARR_SIZE; i++)
sum += flt_arr[i];
return sum;
}
Если есть такая возможность, то, конечно, лучше ей воспользоваться. Однако рассмотим другие методы повышения точности суммирования чисел.
Упорядоченное суммирование
Как было показано выше, ошибка вычисления сильно зависит от порядка выполнения операций. При сложении двух чисел, значительно различающихся по величине, точность вычисления сильно снижается. Представим, что мы суммируем большое множество чисел. Тогда одним из операндов всегда будет сумма, причем в большинстве случаев она по величине значительно превосходит числа, которые мы постепенно к ней прибавляем. Данный эффект потери точности вычислений можно минимизировать путем упорядоченного суммирования чисел от наименьшего к наибольшему по абсолютному значению.
/* Sorted Summation */
float fs_add(float * flt_arr)
{
long i;
float lcl_arr[ARR_SIZE], sum = 0.0;
for (i = 0; i < ARR_SIZE; i++)
lcl_arr[i] = flt_arr[i];
qsort(lcl_arr, ARR_SIZE, sizeof(float),
flt_abs_compare);
for (i = 0; i < ARR_SIZE; i++)
sum += lcl_arr[i];
return sum;
}
Сначала создаем локальную копию массива flt_arr — lcl_arr. Сортируем полученный массив с помощью qsort. Затем суммируем элементы отсортированного массива. Но результаты по-прежнему остаются неудовлетворительными. Это происходит из-за того, что сумма очень быстро становится значительно больше по порядку, чем прибавляемые числа. В результате этого, точность вычисления повышается незначительно.
На самом деле, особенности суммируемых чисел сильно влияют на сортировку. Лучший случай, когда сумма принимает значения около нуля, а распределение примерно симметричное. Это условие способствует тому, что сумма всегда имеет примерно тот же порядок, что и суммируемые числа. Этот случай рассмотрен на примерах. Однако на практике он встречается довольно редко. Иногда сортировку чисел невозможно предугадать. Упорядоченное суммирование всегда дает один и тот же результат на одних и тех же множеств чисел, даже если числа следуют в разном порядке. Но это не совсем верно в случае, например, если использовать функцию сравнения:
/* Simple comparison */
int flt_abs_compare(const void * a, const void * b)
{
float fa = fabs(*(float *)a);
float fb = fabs(*(float *)b);
if (fa == fb)
return 0;
else if (fa < fb)
return -1;
else
return 1;
}
В этом случае число и обратная (по знаку) ему величина считаются одинаковыми. А это означает, что они могут оказаться в любом порядке следования в упорядоченной последовательности. Когда же нам нужна уверенность, необходимо использовать функцию сравнения вида:
/* Tie-breaking comparison */
int flt_abs_compare(const void * a,
const void * b)
{
float va = *(float *)a;
float vb = *(float *)b;
float fa = fabs(va);
float fb = fabs(vb);
if (fa < fb)
return -1;
else if (fa > fb)
return 1;
else if (va < vb)
return -1;
else if (va > vb)
return 1;
else
return 0;
}
Данная функция добивается полной предсказуемости следования элементов.
Попарное суммирование
Ниже представлен алгоритм, в цикле обрабатывающий массив — на каждой итерации суммируются числа парами, при этом размер массива уменьшается вдвое.
/* Pairwise Summation */
float fp_add(float * flt_arr)
{
long i, j, limit;
float sum[ARR_SIZE / 2 + 1];
if (ARR_SIZE == 2)
return flt_arr[0] + flt_arr[1];
else if (ARR_SIZE == 1)
return flt_arr[0];
for (i = j = 0; i < ARR_SIZE / 2; i++, j += 2)
sum[i] = flt_arr[j] + flt_arr[j + 1];
if (ARR_SIZE & 1)
sum[ARR_SIZE / 2] = flt_arr[ARR_SIZE - 1];
limit = (ARR_SIZE + 1) / 2;
while (limit > 2) {
for (i = j = 0; i < limit / 2; i++, j += 2)
sum[i] = sum[j] + sum[j + 1];
if (limit & 1)
sum[limit / 2] = sum[limit - 1];
limit = (limit + 1) / 2;
}
return sum[0] + sum[1];
}
Данный алгоритм работает быстрее, чем упорядоченное суммирование, и при этом дает более аккуратный результат. Однако, нет предела совершенству…
Метод суммирования Кохена
D англоязычной литературе данный метод известен под названием Kahan summation algorithm или compensated summation.
Если перед нами стоит задача нахождения суммы большого множества чисел, то для достижения данной цели лучше всего использовать метод, известный под названием Kahan Summation. В данном методе коррекция промежуточной суммы производится на протяжении всей работы алгоритма. Ниже приведена реализация данного метода.
/* Kahan Summation */
float fk_add(float * flt_arr)
{
long i;
float sum, correction, corrected_next_term, new_sum;
sum = flt_arr[0];
correction = 0.0;
for (i = 1; i < ARR_SIZE; i++)
{
corrected_next_term = flt_arr[i] - correction;
/* sum большое, corrected_next_term маленькое, т.о. младшие биты corrected_next_term теряются */
new_sum = sum + corrected_next_term;
correction = (new_sum - sum) - corrected_next_term;
sum = new_sum;
}
return sum;
}
Каждый член вначале корректируется согласно ошибке, накопившаяся на предыдущих этапах. Затем новая сумма вычисляется путем суммирования скорректированного члена и промежуточной суммы. После этого вычисляется поправочный член как разница между изменением в сумме и прибавленного скорректированного члена.
Алгебраически, кажется, данная процедура не делает ничего особенного. Подставим выражение
new_sum = sum + corrected_next_term;
в строчку
correction = (new_sum — sum) — corrected_next_term;
Получим:
correction = ((sum + corrected_next_term) — sum) — corrected_next_term;
что можно упростить как:
correction = corrected_next_term — corrected_next_term = 0;
для точных чисел. Но для чисел с плавающей точкой поправочный член чаще всего ненулевой. Используя его, повышается точность вычисления.
Для иллюстрации сказанного рассмотрим пример десятичной арифметики чисел с плавающей точкой, в которой значащими являются 6 знаков после запятой.
Допустим, сумма sum достигла на каком-то шаге значения 100000 и нужно прибавить следующее число, равное 3.14159, correction равно нулю. На данной итерации будут произведены следующие высисления:
corrected_next_term = 3.14159 — 0
new_sum = 100000 + 3.14159
= 100003 /* Много цифр теряется */
correction = (100003 - 100000) - 3.14159 /* Должно вычисляться так, как написано! */
= 3.00000 - 3.14159
= -0.141590
sum = 100003
Значение суммы настолько велико, что в новое значение суммы вносят вклад лишь цифры числа большого порядка. Однако, на следующем этапе работы алгоритма, предположим, что слагаемое flt_arr[i] принимает значение 2.71828, но теперь correction ненулевое.
corrected_next_term = 2.71828 — −0.141590
= 2.85987
new_sum = 100003 + 2.85987
= 100006
correction = (100006 - 100003) - 2.85987
= 3.00000 - 2.85987
= .140130
sum = 100006
Сумма вычисляется в процессе выполнения двух операций — sum хранит текущую сумму, а correction накапливает части чисел, не вошедших в sum .
Список литературы
