11

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.

Ondřej Čertík
  • 2,930
  • 18
  • 40

3 Answers3

11

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
Jed Brown
  • 25,650
  • 3
  • 72
  • 130
  • 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. – Ondřej Čertík Sep 24 '12 at 04:55
5

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.

GertVdE
  • 6,179
  • 1
  • 21
  • 36
  • 2
    Here is my Fortran implementation: https://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. – Ondřej Čertík Sep 24 '12 at 16:31
5

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.

Arnold Neumaier
  • 11,318
  • 20
  • 47
  • 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. – Ondřej Čertík Sep 24 '12 at 16:44
  • @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. – Arnold Neumaier Sep 25 '12 at 07:15
  • 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! – Ondřej Čertík Sep 25 '12 at 17:59
  • 1
    @OndřejČertík: To find the additional relative error given absacc, monitor the maximum of $\frac{|y-x|-absacc}{\max(|x|,|y|)}$. – Arnold Neumaier Sep 25 '12 at 18:06