Analyzing Numerical Error in C++ Function

Question

Suppose that I have a function that takes as input several floating-point values (single or double), does some computation, and produces output floating-point values (also single or double). I am working primarily with MSVC 2008, but also plan to work with MinGW/GCC. I'm programming in C++.

What is the typical way of programmatically measuring how much error is in the results? Assuming that I need to use an arbitrary precision library: what is the best such library if I don't care about speed?

Answer 1

If you're looking for a good bound on your rounding error, you don't necessarily need an aribtrary-precision library. You can use running error analysis instead.

I wasn't able to find a good online reference, but it's all described in Section 3.3 of Nick Higham's book "Accuracy and Stability of Numerical Algorithms". The idea is quite simple:

1. Re-factor your code so that you have a single assignment of a single arithmetic operation on each line.
2. For each variable, e.g. x, create a variable x_err which is initialized to zero when x is assigned a constant.
3. For each operation, e.g. z = x * y, update the variable z_err using the standard model of floating-point arithmetic and the resulting z and the running errors x_err and y_err.
4. The return value of your function should then also have a respective _err value attached to it. This is a data-dependent bound on your total roundoff error.

The tricky part is step 3. For the most simple arithmetic operations, you can use the following rules:

• 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))

where u = eps/2 is the unit roundoff. Yes, the rules for + and - are the same. Rules for any other operation op(x) can be easily extracted using the Taylor series expansion of the result applied to op(x + x_err). Or you can try googling. Or using Nick Higham's book.

As an example, consider the following Matlab/Octave code which evaluates a polynomials in the coefficients a at a point x using the Horner scheme:

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

For the first step, we split-up the two operations in 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

We then introduce the _err variables. Note that the inputs a and x are assumed to be exact, but we could just as well also require the user to pass corresponding values for a_err and 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

Finally, we apply the rules described above to get the error terms:

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

Note that since we have no a_err or x_err, e.g. they are assumed to be zero, the respective terms are simply ignored in the error expressions.

Et voilà! We now have a Horner scheme which returns a data-dependent error estimate (note: this is an upper bound on the error) alongside the result.

As a side note, since you're using C++, you might consider making your own class for floating-point values which carries around the _err term and overloading all the arithmetic operations to update these values as described above. For large codes, this may be the easier, albeit computationally less efficient, route. Having said that, you may be able to find such a class online. A quick Google search gave me this link.

P.S. Note that this all works only on machines adhering strictly to IEEE-754, i.e. all arithmetic operations are precise to $\pm u$. This analysis also gives a tighter, more realistic bound than using interval arithmetic since, by definition, you can not represent a number $x(1 \pm u)$ in floating point, i.e. your interval would just round to the number itself.

Answer 2

GMP (i.e., the GNU Multiple Precision library) is the best arbitrary precision library that I know of.

I don't know of any programmatic way to measure the error in the results of an arbitrary floating point function. One thing you could try is to calculate an interval extension of a function using interval arithmetic. In C++, you would have to use some sort of library to calculate the interval extensions; one such library is the Boost Interval Arithmetic Library. Basically, to measure the error, you would supply as arguments to your function intervals that have a width of 2 times unit roundoff (roughly), centered at the values of interest, and then your output would be a collection of intervals, the width of which would give you some conservative estimate of the error. A difficulty with this approach is that interval arithmetic used in this fashion can overestimate error by significant amounts, but this approach is the most "programmatic" one I can think of.

Answer 3

A nice portable and open source library for arbitrary precision floating point arithmetic (and much else) is Victor Shoup's NTL, which is available in C++ source form.

At a lower level is the GNU Multiple Precision (GMP) Bignum Library, also an open source package.

NTL can be used with GMP is faster performance is required, but NTL provides its own base routines that are certainly usable if you "don't care about speed". GMP claims to be the "fastest bignum library". GMP is largely written in C, but has a C++ interface.

Added: While interval arithmetic can give upper and lower bounds on the exact answer in an automated way, this doesn't accurately measure the error in a "standard" precision computation because the interval size typically grows with each operation (either in a relative or absolute error sense).

The typical way of finding error size, either for rounding errors or for errors of discretization, etc., is to compute an extra precision value and compare that to the "standard" precision value. Only modest extra precision is needed to determine the error size itself to reasonable accuracy, since the rounding errors alone are substantially larger in "standard" precision than they are in the extra precision computation.

The point can be illustrated by comparing single and double precision computations. Note that in C++ intermediate expressions are always computed in (at least) double precision, so if we want to illustrate what a computation in "pure" single precision would be like, we need to store out the intermediate values in single precision.

C code snippet

    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;

The output from above (Cygwin/MinGW32 GCC tool chain):

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

Thus the error is about what one expects in rounding 1/3rd to single precision. One would not (I'd suspect) care about getting more than a couple of decimal places in the error correct, since measurement of the error is for magnitude and not for exactitude.

Answer 4

Rigorous and automatic error estimation can be achieved by interval analysis. You work with intervals instead of numbers. For example addition:

[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]

Round-off can also be handled rigorously, see Rounded interval arithmetic.

As long as your input consists of narrow intervals, the estimates are OK and are faily cheap to compute. Unfortunately, the error is often overestimated, see the dependence problem.

I do not know of any arbitrary precision interval arithmetic library.

It depends on your problem at hand whether interval arithmetic can serve your needs or not.

Answer 5

The GNU MPFR library is an arbitrary-precision float library that has high accuracy (in particular, correct rounding for all operations, which is not as easy as it sounds) as one of their main focus points. It uses GNU MP under the hood. It has an extension called MPFI that does interval arithmetic, which - as Geoff's answer suggests - might come in handy for verification purposes: keep increasing the working precision until the resulting interval falls within small bounds.

This won't always work, though; in particular it's not necessarily effective if you're doing something like numerical integration, where every step carries an "error" independent of rounding problems. In that case, try a specialized package such as COSY infinity which does this very well using specific algorithms to bound the integration error (and using so-called Taylor models instead of intervals).

Answer 6

@Pedro mentioned running error analysis and introduced the math for propagating error while accounting for unit roundoff. They also gave a Matlab/Octave implementation. @plasmacel noticed that in this implementation, the error terms themselves were being computed using floating point values, which also have roundoff error which are not being accounted for.

This might be overkill for most purposes, but this C++ implementation strictly accounts for all round-offs in every operation, even when computing the error terms themselves:

#include <limits>
#include <cmath>

static_assert(std::numeric_limits<float>::is_iec559);

float unit_round_off(float result) {
    return std::numeric_limits<float>::epsilon() * result / 2;
}

void mul_w_err(float* r, float a, float b, float* r_err_output, float a_err, float b_err) {
    // r_err = a_err * b + b_err * a + a_err * b_err
    *r = a * b;
    float r_err = unit_round_off(*r);

    float term1 = a_err * b;
    term1 += unit_round_off(term1);

    float term2 = b_err * a;
    term2 += unit_round_off(term2);

    float term3 = a_err * b_err;
    term3 += unit_round_off(term3);

    float term12 = term1 + term2;
    term12 += unit_round_off(term12);

    float term3r = term3 + r_err;
    term3r += unit_round_off(term3r);

    float complete_err = term12 + term3r;
    complete_err += unit_round_off(complete_err);

    *r_err_output = complete_err;
}

void __add_sub_w_err_funct(float* r, float a, float b, float* r_err_output, float a_err, float b_err) {
    // being careful to round the floating point error up after every binary op
    float r_err = unit_round_off(*r);
    
    float ab_err_sum = a_err + b_err;
    ab_err_sum += unit_round_off(ab_err_sum);

    float abr_err_sum = r_err + ab_err_sum;
    abr_err_sum += unit_round_off(abr_err_sum);

    *r_err_output = abr_err_sum;
}

void add_w_err(float* r, float a, float b, float* r_err_output, float a_err, float b_err) {
    // r_err = a_err + b_err
    *r = a + b;
    __add_sub_w_err_funct(r, a, b, r_err_output, a_err, b_err);
}

void sub_w_err(float* r, float a, float b, float* r_err_output, float a_err, float b_err) {
    // r_err = a_err + b_err
    *r = a - b;
    __add_sub_w_err_funct(r, a, b, r_err_output, a_err, b_err);
}

Answer 7

I am told that MPIR is a good library to use if you are working with Visual Studio:

http://mpir.org/