[R] Question about converting from square roots to decimals

[(Ted Harding)]

On 20-Jan-07 Robert Barber wrote:
> Hi,
> 
> I apologize if there is a simple answer to this question that I've
> missed.  I did search the mailing list but I might not have used the
> right keywords. Why does sum(A3^2) give the result of 1, but
> sum(A3^2)==1 give the result of FALSE? 
> 
>> A3<-matrix(nrow=3,c(1/(2^.5),1/(2^.5),0))
>> A3
>           [,1]
> [1,] 0.7071068
> [2,] 0.7071068
> [3,] 0.0000000
>> sum(A3^2)
> [1] 1
>> sum(A3^2)^.5
> [1] 1
>> sum(A3^2)==1         # here's the part I don't understand
> [1] FALSE
>> sum(A3^2)^.5==1      # here's the part I don't understand
> [1] FALSE
> 
> I realize that it has something to do with the conversion of the square
> roots into decimals.  But shouldn't it then give me some number other
> than 1 as the result for sum(A3^2)?  Are there other ways to do this
> than what I've tried?  I'm trying to confirm that A3 is a unit vector.

This is an instance of what must be a candidate for the MFAQAT
(most frequently asked qustion of all time).

The nub of the matter can be found in FAQ 7.31:

  http://www.r-project.org/  -->  FAQs

where, at 7.31, it says:

  7.31 Why doesn't R think these numbers are equal?

  The only numbers that can be represented exactly in R's
  numeric type are integers and fractions whose denominator
  is a power of 2. Other numbers have to be rounded to
  (typically) 53 binary digits accuracy. As a result, two
  floating point numbers will not reliably be equal unless
  they have been computed by the same algorithm, and not
  always even then. For example

     R> a <- sqrt(2)
     R> a * a == 2
     [1] FALSE
     R> a * a - 2
     [1] 4.440892e-16

  The function all.equal() compares two objects using a
  numeric tolerance of .Machine$double.eps ^ 0.5. If you
  want much greater accuracy than this you will need to
  consider error propagation carefully.

  For more information, see e.g. David Goldberg (1991),
  "What Every Computer Scientist Should Know About Floating-Point
  Arithmetic", ACM Computing Surveys, 23/1, 5-48, also
  available via
  http://docs.sun.com/source/806-3568/ncg_goldberg.html.

However, what this FAQ does not point out is that "==" tests
for exact equality, and even the 'help' page

  ?"=="

is not as explicit about this:

  For numerical values, remember '==' and '!=' do not allow
  for the finite representation of fractions, nor for
  rounding error. Using 'all.equal' with 'identical' is
  almost always preferable.  See the examples.

Nor does the help

  ?all.equal

give the explanation:

  'all.equal(x,y)' is a utility to compare R objects 'x'
  and 'y' testing "near equality". ...

and a user who is not aware of the imprecision problems inherent
in computer arithmetic may well not see the point of testing for
"near equality" when the user expects exact equality mathematically.

Statistically speaking, the frequency of occurrence of variants
of this question is a phenomenon!

I hypothesise that it arises from a combination of two main
circumstances:

a) Many users who are unfamiliar with the technical details of
   finite-length binary arithmetic will not be expecting that
   there could be a problem of this kind in the first place.
   So, when it occurs, they will simply be puzzled.

b) It's actually quite difficult to find your way to the above
   explanation in the FAQs. First, you need to anticipate that
   this is the sort of thing that will be a FAQ. If you're subject
   to (a) above, you may not be thinking on these lines.

   Secondly, even if you get as far as looking at the FAQs at
   the above URL, you have a lot of scrolling down to do before
   you find the question

     7.31 Why doesn't R think these numbers are equal?

   and even then, if you blink you'll miss it: trying it just
   now, I in fact didn't spot it in the list of questions
   immediately below the header

     7 R Miscellanea

   even though I already knew it was there somewhere and was
   actively looking for it. It was only when I landed on the
   full FAQ itself that I recognised it.

   It would have been quicker (I just tried that too) to start
   at the top of the FAQs page, and use the browser's Search
   tool to search for == : the sixth occurrence of "==" was it!

   But, even then, you still need to be thinking that the answer
   is to be found in connection with "==". If you're subject to
   (a), you'll be wondering instead why the answer was wrong.

This is such a frequently occurring issue that I feel there is
a case for a prominently displayed link, maybe in FAQs but maybe
better at the top level of r-project.org, to a brief but adequate
discourse with a title like

  Arithmetic [im]precision and apparent errors in R

What do others think?

Ted.

--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 20-Jan-07                                       Time: 19:53:03
------------------------------ XFMail ------------------------------