# Re: st: Another precision problem with Mata

 From wgould@stata.com (William Gould, Stata) To statalist@hsphsun2.harvard.edu Subject Re: st: Another precision problem with Mata Date Mon, 22 Aug 2005 09:55:42 -0500

```Ben Jann <ben.jann@soz.gess.ethz.ch> asks,

> in a Mata program I am using the expression
>
>                abs((x-y)/z) < 1
>
> where x, y and z are real scalars. Note that in some cases y has been
> set to x + z earlier in the program. Although the expression should
> evaluate to "false" in these cases, it sometimes does evaluate to "true"
> due to roundoff error.
>
> My question now is: Would it be reasonable to code
>
>             abs((x-y)/z) < 1-epsilon(1)
>
> Or how would one approach such a problem?

The first step in approaching this problem is to understand the source
of the error.  Ben reports that when y = x + z, he can have problems.
In the y=x+z case,

abs((x-y)/z)  = abs((x-(x+z))/z)
= abs(-z/z)            (inifinite precision)
= 1                    (inifinite precision)

Ben is hoping that x-(x+z) == z, and that is not necessarily true when
we perform finite-precision arithmetic.

So now let's consider (x-(x+z)), the source of the problem.
I am going to use the following notation

-         infinite-precision subtraction

SUB       finite-precision subtraction

In this case

x PLUS z = x+z + error

The error can be so large so that

x PLUS z = x               (1)
or
x PLUS z = z               (2)

All that is required is that (1) x be very large relative to z, or
(2) vice versa.

The importance of that point is that we cannot bound the error by
epsilon(1), as Ben was hoping.

I cannot go further with this problem, in fact, without knowing something
about the relative magnitudes of x and z.

However, Ben has just revealed that it is important to him to identify
abs((x-y)/z) == 1 when y = x+z.  In that case, I suggest Ben remove y from his
code and substitue (d+z) everywhere for it.  Then abs((x-y)/z) == 1 is
equivalent to abs(d)==abs(x).

Perhaps there are other possible code-specific solutions, or perhaps,
Ben does know something about magnitudes and we can bound the error, and
Ben will be willing to accept abs((x-y)/z)<1 will evaluate to false in
some cases where, in infinite precision, it would evaluate to true.

But one thing I can guarantee:  There is no reasons to think about how
one might code

abs((x-y)/z) < 1

so that it is guaranteed to be false when y = x PLUS z, for all values
of x and z.

-- Bill
wgould@stata.com
*
*   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/
```