3

While checking some matrix multiplication operations, I came across a strange behavior. I get different results when I perform the multiplication "by hand" (using the product and the sum) and when using the matrix multiplication operator %*%. 

c <- 1:10
a <- 100^(0:9)
p1 <- sum(a*c)
p2 <- a%*%c
p1==p2
      [,1]
[1,] FALSE
p1-p2
      [,1]
[1,] -2048

However, when I use any other value for a (e.g., a <- 101^0:9) , I do get the same results:

c <- 1:10
a <- 101^(0:9)
p1 <- sum(a*c)
p2 <- a%*%c
p1==p2
      [,1]
[1,] TRUE
p1-p2
      [,1]
[1,] 0

Any idea why this is happening?

Thank you, Pedro

Epsylon
  • 31
  • 2
  • 2
    See [this](http://stackoverflow.com/questions/9508518/why-are-these-numbers-not-equal). Use `all.equal` to test for equality. – Roland Dec 07 '12 at 11:48
  • While this is in the R-FAQ, seeing as it's really a "computer math for newbies" FAQ, is there somewhere in SO we could place the answer? – Carl Witthoft Dec 07 '12 at 12:36
  • "Any other value" is not necessarily OK. Try `f <- function(b) { a <- b^(0:9); c<- 1:10; sum(a*c) - a %*% c }; bvec <- 80:120; r <- sapply(bvec,f); plot(bvec,r,type="b")` – Ben Bolker Dec 07 '12 at 14:23

1 Answers1

2

%*% does compute its results in a slightly different way, which means that different rounding errors occur at different places, leading to a different overall result.

I'm just guessing, but I believe that this might be due to sum keeping its accumulator in a machine floating point register, which has 80 bit extended precision on Intel architectures. If you want to know for certain, you'd have to look at the assembly code of R.

MvG
  • 57,380
  • 22
  • 148
  • 276