R Precision for Double - Why code returns negative why positive outcome expected?

Question

I am testing 2 ways of calculating Prod(b-a), where a and b are vectors of length n. Prod(b-a)=(b1-a1)(b2-a2)(b3-a3)*... (bn-an), where b_i>a_i>0 for all i=1,2,3, n. For some special cases, another way (Method 2) of calculation this prod(b-a) is more efficient. It uses the following formula, which is to expand the terms and sum them:

Here is my question is: When it happens that a_i very close to b_i, the true outcome could be very, very close 0, something like 10^(-16). Method 1 (substract and Multiply) always returns positive output. Method 2 of using the formula some times return negative output ( about 7~8% of time returning negative for my experiment). Mathematically, these 2 methods should return exactly the same output. But in computer language, it apparently produces different outputs.

Here are my codes to run the test. When I run the testing code for 10000 times, about 7~8% of my runs for method 2 returns negative output. According to the official document, the R double has the precision of "2.225074e-308" as indicated by R parameter: ".Machine$double.xmin". Why it's getting into the negative values when the differences are between 10^(-16) ~ 10^(-18)? Any help that sheds light on this will be apprecaited. I would also love some suggestions concerning how to practically increase the precision to higher level as indicated by R document.

##########  Testing code 1.  
ftest1case<-function(a,b)  {
  n<-length(a)
  if (length(b)!=n) stop("---------  length a and b are not right.")
  if ( any(b<a) ) stop("----------  b has to be greater than a all the time.")
  out1<-prod(b-a)
  
  out2<-0
  N<-2^n
  for ( i in 1:N ) {
    tidx<-rev(as.integer(intToBits(x=i-1))[1:n])
    tsign<-ifelse( (sum(tidx)%%2)==0,1.0,-1.0)
    out2<-out2+tsign*prod(b[tidx==0])*prod(a[tidx==1])
  }
  c(out1,out2)
}


##########  Testing code 2.  
ftestManyCases<-function(N,printFreq=1000,smallNum=10^(-20))
{
  tt<-matrix(0,nrow=N,ncol=2)
  n<-12
  for ( i in 1:N) {
    a<-runif(n,0,1)
    b<-a+runif(n,0,1)*0.1
    tt[i,]<-ftest1case(a=a,b=b)
  
    if ( (i%%printFreq)==0 ) cat("----- i = ",i,"\n")
    if ( tt[i,2]< smallNum ) cat("------ i = ",i, " ---- Negative summation found.\n")
  }
  
  tout<-apply(tt,2,FUN=function(x)  { round(sum(x<smallNum)/N,6) } )
  names(tout)<-c("PerLess0_Method1","PerLee0_Method2")
  list(summary=tout, data=tt)  
}  
 

########  Step 1.  Test for 1 case.
n<-12
a<-runif(n,0,1)
b<-a+runif(n,0,1)*0.1
ftest1case(a=a,b=b)   

########  Step 2   Test Code 2 for multiple cases.
N<-300
tt<-ftestManyCases(N=N,printFreq = 100) 
tt[[1]]

Answer 1

It's hard for me to imagine when an algorithm that consists of generating 2^n permutations and adding them up is going to be more efficient than a straightforward product of differences, but I'll take your word for it that there are some special cases where it is.

As suggested in comments, the root of your problem is the accumulation of floating-point errors when adding values of different magnitudes; see here for an R-specific question about floating point and here for the generic explanation.

First, a simplified example:

n <- 12
set.seed(1001)
a <- runif(a,0,1)
b <- a + 0.01
prod(a-b) ## 1e-24

out2 <- 0
N <- 2^n
out2v <- numeric(N)
for ( i in 1:N ) {
    tidx <- rev(as.integer(intToBits(x=i-1))[1:n])
    tsign <- ifelse( (sum(tidx)%%2)==0,1.0,-1.0)
    j <- as.logical(tidx)
    out2v[i] <- tsign*prod(b[!j])*prod(a[j])
}
sum(out2v)  ## -2.011703e-21

Using extended precision (with 1000 bits of precision) to check that the simple/brute force calculation is more reliable:

library(Rmpfr)
a_m <- mpfr(a, 1000)
b_m <- mpfr(b, 1000)
prod(a_m-b_m) 
##  1.00000000000000857647286522936696473705868726043995807429578968484409120647055193862325070279593735821154440625984047036486664599510856317884962563644275433171621778761377125514191564456600405460403870124263023336542598111475858881830547350667868450934867675523340703947491662460873009229537576817962228e-24

This proves the point in this case, but in general doing extended-precision arithmetic will probably kill any performance gains you would get.

Redoing the permutation-based calculation with mpfr values (using out2 <- mpfr(0, 1000), and going back to the out2 <- out2 + ... running summation rather than accumulating the values in a vector and calling sum()) gives an accurate answer (at least to the first 20 or so digits, I didn't check farther), but takes 6.5 seconds on my machine (instead of 0.03 seconds when using regular floating-point).

Why is this calculation problematic? First, note the difference between .Machine$double.xmin (approx 2e-308), which is the smallest floating-point value that the system can store, and .Machine$double.eps (approx 2e-16), which is the smallest value such that 1+x > x, i.e. the smallest relative value that can be added without catastrophic cancellation (values a little bit bigger than this magnitude will experience severe, but not catastrophic, cancellation).

Now look at the distribution of values in out2v, the series of values in out2v:

hist(out2v)

There are clusters of negative and positive numbers of similar magnitude. If our summation happens to add a bunch of values that almost cancel (so that the result is very close to 0), then add that to another value that is not nearly zero, we'll get bad cancellation.

It's entirely possible that there's a way to rearrange this calculation so that bad cancellation doesn't happen, but I couldn't think of one easily.