11
$\begingroup$

I have a numerical function f(x, y) returning a double floating point number that implements some formula and I want to check that it is correct against analytic expressions for all combination of the parameters x and y that I am interested in. What is the proper way to compare the computed and analytical floating point numbers?

Let's say the two numbers are a and b. So far I've been making sure that both absolute (abs(a-b) < eps) as well as relative (abs(a-b)/max(abs(a), abs(b)) < eps) errors are less than eps. That way it will catch numerical inaccuracies even if the numbers are let's say around 1e-20.

However, today I discovered a problem, the numerical value a and analytic value b were: 

In [47]: a                                                                     
Out[47]: 5.9781943146790832e-322

In [48]: b                                                                     
Out[48]: 6.0276008792632078e-322

In [50]: abs(a-b)                                                              
Out[50]: 4.9406564584124654e-324

In [52]: abs(a-b) / max(a, b)                                                  
Out[52]: 0.0081967213114754103

So the absolute error [50] is (obviously) small, but the relative error [52] is large. So I thought that I have a bug in my program. By debugging, I realized, that these numbers are denormal. As such, I wrote the following routine to do the proper relative comparison:

real(dp) elemental function rel_error(a, b) result(r)
real(dp), intent(in) :: a, b
real(dp) :: m, d
d = abs(a-b)
m = max(abs(a), abs(b))
if (d < tiny(1._dp)) then
    r = 0
else
    r = d / m
end if
end function

Where tiny(1._dp) returns 2.22507385850720138E-308 on my computer. Now everything works and I simply get 0 as the relative error and all is ok. In particular, the above relative error [52] is wrong, it's simply caused by insufficient accuracy of the denormal numbers. Is my implementation of the rel_error function correct? Should I just check that abs(a-b) is less than tiny (=denormal), and return 0? Or should I check some other combination, like max(abs(a), abs(b))?

I would just like to know what the "proper" way is.

Improve this question
$\endgroup$

3 Answers 3

Reset to default
11
$\begingroup$

You can directly check for denormals using isnormal() from math.h (C99 or later, POSIX.1 or later). In Fortran, if the module ieee_arithmetic is available, you can use ieee_is_normal(). To be more precise about fuzzy equality, you have to consider the floating point representation of denormals and decide what you mean for the results to be good enough.

More to the point, to believe that either result is correct, you have to be sure that you didn't lose too many digits at an intermediate step. Computing with denormals is generally unreliable and should be avoided by having your algorithm rescale internally. To ensure that your internal scaling was successful, I recommend activating floating point exceptions using feenableexcept() with C99 or the ieee_arithmetic module in Fortran.

Although you can have your application catch the signal that is raised on floating point exceptions, all kernels I have tried reset the hardware flag so fetestexcept() does not return a useful result. When run with -fp_trap, PETSc programs will (by default) print a stack trace when a floating point error is raised, but will not identify the offending line. If you run in a debugger, the debugger preserves the hardware flag and breaks on the offending expression. You can check the precise reason by calling fetestexcept from the debugger where the result is a bitwise OR of the following flags (values may vary by machine, see fenv.h; these values are for x86-64 with glibc).

  • FE_INVALID = 0x1
  • FE_DIVBYZERO = 0x4
  • FE_OVERFLOW = 0x8
  • FE_UNDERFLOW = 0x10
  • FE_INEXACT = 0x20
Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for the excellent answer. The analytic expression that I compare against in the asymptotic regime is exp(log_gamma(m+0.5_dp) - (m+0.5_dp)*log(t)) / 2 for m=234, t=2000. It goes to zero quickly as I increase m. All I want to make sure that my numerical routine returns "correct" numbers (to return zero is perfectly fine too) to at least 12 significant digits. So if the calculation returns a denormal number, then it's simply zero, and there should be no problem. So just the comparison routine needs to be robust against this. $\endgroup$
    – Ondřej Čertík
    Sep 24, 2012 at 4:55
5
$\begingroup$

Donald Knuth has a proposal for a floating point comparison algorithm in volume 2 "Seminumerical algorithms" of "The Art of Computer Programming". It was implemented in C by Th. Belding (see fcmp package) and is available in the GSL.

Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ Here is my Fortran implementation: gist.github.com/3776847, notice that I need to explicitly handle denormal numbers anyway in it. Otherwise I think it's pretty much equivalent to the relative error, the only difference is that instead of doing abs(a-b)/max(a, b) < eps, we do abs(a-b)/2**exponent(max(a, b)) < eps, which pretty much just drops the mantissa in the max(a, b), so in my opinion the difference is negligible. $\endgroup$
    – Ondřej Čertík
    Sep 24, 2012 at 16:31
5
$\begingroup$

Optimally rounded denormalized numbers may indeed have a high relative error. (Flushing that to zero while still calling it a relative error is misleading.)

But close to zero, computing relative erros is meaningless.

Therefore, even before you reach denormalized numbers, you should probably switch to absolute accuracy (namely the one you want to guarantee in this case).

I'd therefore suggest to test the computed $y$ against the true $x$ by checking the validity of a formula such as $|y-x|\le absacc+relacc*\max(|x|,|y|)$. say with relacc=1e-12 and absacc=1e-150.

Then users of your code know precisely how much accuracy they really have.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Are you sure that it is meaningless to compute relative errors close to zero? I think it's meaningless only if there is loss of accuracy (for whatever reason). If for example there is loss of accuracy for x < 1e-150 due to some numerical issues (like subtracting two big numbers), then you are right. In my case however the numbers seem to be accurate all the way down to zero, except when it hits the denormal numbers. So in my case absacc=1e-320 or so and I can just check abs(a-b) < tiny(1._dp) as I do above. $\endgroup$
    – Ondřej Čertík
    Sep 24, 2012 at 16:44
  • $\begingroup$ @OndřejČertík: In that case replace the 1e-150 by 1e-300 or whatever bound you can verify. In any case very close to zero you make an absolute error, and your error claim should reflect this rather than declare the relative error to be zero. $\endgroup$
    – Arnold Neumaier
    Sep 25, 2012 at 7:15
  • $\begingroup$ I see. I can verify, that all works for numbers higher than tiny(1._dp)=2.22507385850720138E-308 (I made a mistake in my previous comment, it's 2e-308, not 1e-320). So this is my absolute error. Then I need to compare the relative error. I see your point, I think you are right. Thanks! $\endgroup$
    – Ondřej Čertík
    Sep 25, 2012 at 17:59
  • 1
    $\begingroup$ @OndřejČertík: To find the additional relative error given absacc, monitor the maximum of $\frac{|y-x|-absacc}{\max(|x|,|y|)}$. $\endgroup$
    – Arnold Neumaier
    Sep 25, 2012 at 18:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.