Re: st: Trapping Underflows and Overflows in Mata

[jlinhart@stata.com (Jean Marie Linhart, StataCorp LP)]

Joseph Coveney, <jcoveney at bigplanet dot com>, asks about underflow with
Mata:

> But I need some help in trapping a potential underflow.
> 
> : x = smallestdouble()
> 
> : x = x / 10
> 
> : x
>  4.4501e-309
> 
> : if (x) printf("Not underflowed");;
> Not underflowed
> :

Numeric underflow is when you have a tiny (possibly zero) result and a loss of
accuracy.  The result above *is* underflowed, i.e., it is not of full double
precision accuracy, though it is not zero.  In some underflows, you will indeed
get an answer of 0.

Numbers with less than full double precision accuracy (1 part in 2^52) are
called denormalized numbers.

> It seems as if Mata is handling underflows gradually (denormalized numbers)
> and not abruptly (set to zero).  So, instead of
>
>    if (!result) . . .
>
> it should be
> 
>    if ((result > 0 && result < smallestdouble()) ||
>      (result < 0 && result  > -smallestdouble())) . . .
>
> or
>
>    if (result && result > -smallestdouble()) &&
>      result < smallestdouble()) . . .
>
> or perhaps
> 
>    if (result && abs(result) < smallestdouble()) . . .
>
> Am I on track here?

Yes.  Joseph is testing for what are likely denormalized results, and that
seems to be what he wants to do.

In addition to this gradual underflow, there are cases where you lose precision
and get a zero result.

> I guess that I could also use some guidance on when I should expect
> underflows to become problematic in Mata--for example, the following seem
> hold up well:
> 
> : y = x = smallestdouble()
> 
> : (y / 10) / (x / 10)
>   1
> 
> : (y / 100) / (x / 10)
>  .1
>
> But holding the exponent at -308 (or rather its equivalent in base 2) and
> successively shifting the fraction rightward will eventually show
> unacceptably degraded precision for a typical application.

In order to answer this, Joseph needs to think about what is "unacceptably
degraded precision" for him?  And what is he trying to do?  The answer may
change with the question posed; problems with loss of precision must be thought
through carefully.

Double precision numbers generally have 52 bits of precision, so accuracy to 1
part in 2^52, or 15-16 base 10 digits of precision.  For most of us, 7-8 digits
of base 10 precision are plenty, this is a denormalized result of about the
order of magnitude of 5x10^(-316).  This is roughly midway in order of
magnitude between the smallest full precision double and the smallest
denormalized double.

Underflow is not the only potential culprit when it comes to loss of precision.
An ill-thought out subtraction call kill off quite a bit of precision.  E.g.,

: x = 2*10^20
: y = 1
: x - (x + y)
  0

Due to precision problems in this example, x+y is the same as x.  x and y each
have about 16 decimal digits of precision (52 binary digits of precision), but
after adding y and subtracting, there is nothing left.  How often do most
people think hard about subtraction and precision?  We don't normally need to
-- our results rarely need to have a full 16 digits of base 10 precision.
Double precision calculations give us a comfortable safety margin.

Below my signature, I'll add some suggestions for further reading.

Hope that helps,

--Jean Marie
jlinhart@stata.com

If you want to know more, try:

Wikipedia on the IEEE standard for binary floating point arithmetic:
http://en.wikipedia.org/wiki/IEEE_754

What every computer scientist should know about floating-point arithmetic:
http://docs.sun.com/source/806-3568/ncg_goldberg.html

An Interview with the Old Man of Floating-Point (interview with William Kahan)
(with a section on the controversial history of gradual underflow)
http://www.cs.berkeley.edu/~wkahan/ieee754status/754story.html
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/