Анализ числовой ошибки в функции C ++

Предположим, что у меня есть функция, которая принимает в качестве входных данных несколько значений с плавающей точкой (одинарные или двойные), выполняет некоторые вычисления и выдает выходные значения с плавающей точкой (также одинарные или двойные). Я работаю в основном с MSVC 2008, но также планирую работать с MinGW / GCC. Я программирую на C ++.

Каков типичный способ программного измерения количества ошибок в результатах? Предполагая, что мне нужно использовать библиотеку произвольной точности: какая библиотека лучше, если мне не важна скорость?

Если вы ищете надежную границу для ошибки округления, вам не обязательно нужна библиотека с точностью до aribtrary. Вместо этого вы можете использовать анализ ошибок.

Мне не удалось найти хорошую онлайн-справку, но все это описано в разделе 3.3 книги Ника Хайама «Точность и стабильность численных алгоритмов». Идея довольно проста:

1. Перефакторинг вашего кода, чтобы у вас было одно назначение одной арифметической операции в каждой строке.
2. Например, для каждой переменной xсоздайте переменную, x_errкоторая инициализируется нулем, когда ей xприсваивается константа.
3. Например, для каждой операции z = x * yобновите переменную, z_errиспользуя стандартную модель арифметики с плавающей запятой, а также возникающие zи текущие ошибки x_errи y_err.
4. Возвращаемое значение вашей функции должно также иметь соответствующее _errзначение. Это зависит от данных на вашей общей ошибки округления.

Сложная часть - шаг 3. Для самых простых арифметических операций вы можете использовать следующие правила:

• z = x + y -> z_err = u*abs(z) + x_err + y_err
• z = x - y -> z_err = u*abs(z) + x_err + y_err
• z = x * y -> z_err = u*abs(z) + x_err*abs(y) + y_err*abs(x)
• z = x / y -> z_err = u*abs(z) + (x_err*abs(y) + y_err*abs(x))/y^2
• z = sqrt(x) -> z_err = u*abs(z) + x_err/(2*abs(z))

где u = eps/2блок округления. Да, правила для +и -одинаковы. Правила для любой другой операции op(x)могут быть легко извлечены с использованием разложения в ряд Тейлора для примененного результата op(x + x_err). Или вы можете попробовать поискать в Google. Или используя книгу Ника Хайама.

В качестве примера рассмотрим следующий код Matlab / Octave, который оценивает полиномы в коэффициентах aв точке, xиспользуя схему Хорнера:

function s = horner ( a , x )
    s = a(end);
    for k=length(a)-1:-1:1
        s = a(k) + x*s;
    end

Для первого шага мы разделили две операции на s = a(k) + x*s:

function s = horner ( a , x )
    s = a(end);
    for k=length(a)-1:-1:1
        z = x*s;
        s = a(k) + z;
    end

Затем мы вводим _errпеременные. Обратите внимание, что входные данные aи xпредполагаются точными, но мы также можем также потребовать от пользователя передать соответствующие значения для a_errи x_err:

function [ s , s_err ] = horner ( a , x )
    s = a(end);
    s_err = 0;
    for k=length(a)-1:-1:1
        z = x*s;
        z_err = ...;
        s = a(k) + z;
        s_err = ...;
    end

Наконец, мы применяем правила, описанные выше, чтобы получить условия ошибки:

function [ s , s_err ] = horner ( a , x )
    u = eps/2;
    s = a(end);
    s_err = 0;
    for k=length(a)-1:-1:1
        z = x*s;
        z_err = u*abs(z) + s_err*abs(x);
        s = a(k) + z;
        s_err = u*abs(s) + z_err;
    end

Обратите внимание, что поскольку у нас нет a_errили x_err, например, они предполагаются равными нулю, соответствующие выражения просто игнорируются в выражениях ошибок.

И вуаля! Теперь у нас есть схема Хорнера, которая возвращает оценку ошибки в зависимости от данных (примечание: это верхняя граница ошибки) вместе с результатом.

В качестве примечания, поскольку вы используете C ++, вы можете рассмотреть возможность создания собственного класса для значений с плавающей запятой, который переносит _errтермин и перегружает все арифметические операции для обновления этих значений, как описано выше. Для больших кодов это может быть более простой, хотя и менее эффективный в вычислительном отношении маршрут. Сказав это, вы можете найти такой класс в Интернете. Быстрый поиск в Google дал мне эту ссылку .

PS Обратите внимание, что все это работает только на машинах, строго придерживающихся IEEE-754, то есть все арифметические операции с точностью до . Этот анализ также дает более жесткую, более реалистичную границу, чем использование интервальной арифметики, поскольку по определению вы не можете представлять число в плавающей точке, т.е. ваш интервал будет просто округляться до самого числа.x ( 1 ± u $\pm u$$x(1 \pm u)$

Хорошей переносимой библиотекой с открытым исходным кодом для арифметики с плавающей запятой произвольной точности (и многого другого) является NTL Виктора Шоупа , которая доступна в исходной форме C ++.

На более низком уровне находится библиотека GNU Multiple Precision (GMP) Bignum , также пакет с открытым исходным кодом.

NTL можно использовать с GMP, если требуется более высокая производительность, но NTL предоставляет свои собственные базовые подпрограммы, которые, безусловно, можно использовать, если вам «не безразлична скорость». GMP претендует на звание «самой быстрой библиотеки bignum». GMP в основном написан на C, но имеет интерфейс C ++.

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

Типичный способ определения размера ошибки, как для ошибок округления, так и для ошибок дискретизации и т. Д., Состоит в том, чтобы вычислить дополнительное значение точности и сравнить его со «стандартным» значением точности. Только умеренная дополнительная точность необходима для определения самого размера ошибки с разумной точностью, поскольку одни только ошибки округления существенно больше в «стандартной» точности, чем в вычислениях дополнительной точности.

Это можно проиллюстрировать, сравнив вычисления одинарной и двойной точности. Обратите внимание, что в C ++ промежуточные выражения всегда вычисляются с (по крайней мере) двойной точностью, поэтому, если мы хотим показать, на что похожи вычисления с «чистой» одинарной точностью, нам нужно хранить промежуточные значения с одинарной точностью.

Фрагмент кода C

    float fa,fb;
    double da,db,err;
    fa = 4.0;
    fb = 3.0;
    fa = fa/fb;
    fa -= 1.0;

    da = 4.0;
    db = 3.0;
    da = da/db;
    da -= 1.0;

    err = fa - da;
    printf("Single precision error wrt double precision value\n");
    printf("Error in getting 1/3rd is %e\n",err);
    return 0;

Выход сверху (цепочка инструментов Cygwin / MinGW32 GCC):

Single precision error wrt double precision value
Error in getting 1/3rd is 3.973643e-08

Таким образом, ошибка связана с тем, что можно ожидать при округлении 1/3 до одинарной точности. Нельзя было бы (я подозреваю) заботиться о том, чтобы получить правильную ошибку более, чем на пару знаков после запятой , поскольку измерение ошибки производится по величине, а не по точности.

GMP (то есть библиотека GNU Multiple Precision) - лучшая библиотека произвольной точности, которую я знаю.

Я не знаю ни одного программного способа измерения ошибки в результатах произвольной функции с плавающей запятой. Одна вещь, которую вы можете попробовать - это вычислить интервальное расширение функции, используя интервальную арифметику . В C ++ вам придется использовать какую-то библиотеку для вычисления расширений интервалов; одной из таких библиотек является арифметическая библиотека с интервалом увеличения, По сути, чтобы измерить ошибку, вы должны указать в качестве аргументов свои функциональные интервалы, которые имеют ширину округления в 2 раза (приблизительно), центрированную по интересующим значениям, и тогда ваш вывод будет представлять собой набор интервалов, ширину что даст вам некоторую консервативную оценку ошибки. Сложность этого подхода заключается в том, что интервальная арифметика, используемая таким образом, может значительно переоценить ошибку, но этот подход является наиболее «программным» из всех, которые я могу себе представить.

Точная и автоматическая оценка ошибок может быть достигнута с помощью интервального анализа . Вы работаете с интервалами вместо чисел. Например дополнение:

[a,b] + [c,d] = [min(a+c, a+d, b+c, b+d), max (a+c, a+d, b+c, b+d)] = [a+c, b+d]

С округлением также можно обращаться строго, см. Арифметику с закругленными интервалами .

Пока ваш вклад состоит из узких интервалов, оценки в порядке и очень дешевы для вычисления. К сожалению, ошибка часто переоценивается, см. Проблему зависимости .

Я не знаю ни одной арифметической библиотеки с произвольным интервалом точности.

От вашей проблемы зависит, может ли интервальная арифметика удовлетворить ваши потребности или нет.

Библиотека GNU MPFR - это библиотека с плавающей точкой произвольной точности, которая имеет высокую точность (в частности, правильное округление для всех операций, что не так просто, как кажется) в качестве одной из их основных точек фокусировки. Он использует GNU MP под капотом. Он имеет расширение под названием MPFI, которое выполняет интервальную арифметику, которая, как предполагает ответ Джеффа, может пригодиться для целей проверки: продолжайте увеличивать рабочую точность до тех пор, пока результирующий интервал не окажется в небольших пределах.

Это не всегда работает, хотя; в частности, это не обязательно эффективно, если вы делаете что-то вроде численного интегрирования, где каждый шаг несет «ошибку» независимо от проблем округления. В этом случае попробуйте специальный пакет, такой как COSY infinity, который делает это очень хорошо, используя специальные алгоритмы для ограничения ошибки интеграции (и используя так называемые модели Тейлора вместо интервалов).

Мне сказали, что MPIR - хорошая библиотека для использования, если вы работаете с Visual Studio:

http://mpir.org/