R语言双精度计算异常：Prod(b-a)方法2为何返回负值？

问题背景 #

我正在测试两种计算$\text{Prod}(b-a)$的方法，其中$a$、$b$是长度为$n$的向量，$\text{Prod}(b-a)$定义为所有对应元素差的乘积：
$$\text{Prod}(b-a)=(b_1-a_1) \times (b_2-a_2) \times \dots \times (b_n-a_n)$$
已知对所有$i$满足$b_i>a_i>0$。部分场景下，展开为求和形式的方法2效率更高，其数学公式为：

$$\prod_{i=1}^n (b_i - a_i) = \sum_{S\subseteq{1,\dots,n}} (-1)^{|S|} \left( \prod_{i\notin S} b_i \right) \left( \prod_{i\in S} a_i \right)$$

异常现象 #

当$a_i$与$b_i$非常接近时，真实结果趋近于0（如$10^{-16}$量级）：

• 方法1（直接相减后相乘）始终返回正值；
• 方法2（展开求和）约7~8%的测试案例返回负值。

数学上两种方法结果完全一致，但实际计算存在差异，我疑惑为何在结果差异处于$10^{-16} \sim 10^{-18}$区间时会出现符号反转，同时希望得到提升计算精度的实用建议。

测试代码 #

单案例测试函数 #

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)
}

多批量测试函数 #

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. 单案例测试
n<-12
a<-runif(n,0,1)
b<-a+runif(n,0,1)*0.1
ftest1case(a=a,b=b)   

# Step 2. 多案例批量测试
N<-300
tt<-ftestManyCases(N=N,printFreq = 100) 
tt[[1]]

原因分析 #

1. 浮点数舍入误差累积：方法2需要计算$2^n$项的加减（当$n=12$时为4096项），每一步浮点数运算都会产生微小舍入误差。当最终结果趋近于0时，这些误差的累积会打破正负项的完美抵消，导致结果符号反转。
2. 大数相消问题：方法2中的项包含$\prod b_i$（大数）和$\prod a_i$（接近大数），当$a_i$与$b_i$接近时，这些项的量级远大于最终结果。相减时会丢失大量有效数字，微小的舍入误差就足以改变结果的符号。
3. 方法1的稳定性优势：直接计算$b_i-a_i$后相乘，每一步都是小数相乘，不会出现大数相消的场景，因此符号始终正确，精度表现更稳定。

精度提升建议 #

• 优先选择方法1：除非有明确的性能需求且能接受精度损失，否则直接计算元素差再相乘是更可靠的方案。
• 优化方法2的求和逻辑：
  • 对求和项按绝对值从小到大排序，先累加小项再累加大项，减少误差累积；
  • 使用Kahan补偿求和算法，降低加减运算中的舍入误差。
• 使用高精度数据类型：借助R的Rmpfr包，用多精度浮点数替代默认双精度类型，大幅提升计算精度。
• 符号修正机制：若必须使用方法2，可在计算后添加符号修正——由于数学上结果必为正，若计算结果为负且绝对值极小（如小于$10^{-15}$），可将其修正为正的极小值。

内容的提问来源于stack exchange，提问作者Max577