Re: [R] Question about converting from square roots to decimals

From: Ted Harding <Ted.Harding_at_manchester.ac.uk>
Date: Sat 20 Jan 2007 - 19:53:07 GMT


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:

  1. 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.
  2. 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@manchester.ac.uk> Fax-to-email: +44 (0)870 094 0861
Date: 20-Jan-07                                       Time: 19:53:03
------------------------------ XFMail ------------------------------

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun Jan 21 06:59:04 2007

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Sun 21 Jan 2007 - 02:30:22 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.