第一轮回 浮点陷阱

在穿过黄泉路之后,我们到了第一站:善良异教徒之家.这些人供奉于”对浮点数无知之神”,这些人认为:

                                               .1 == .3 / 3

是正确的.

他们同样确切地认为: 

                                          seq(0, 1, by=.1) == .3

只有一个值是真值.

但是你不应该认为:               

                              unique(c(.3, .4 - .1, .5 - .2, .6 - .3, .7 - .4))

的长度为1.

我的第一个程序写于石器时代,一段关于求二次方程的代码.当时代码的媒介还是纸带.而打孔机上是没有后退键的---一旦你打下了洞,洞就在那,无法修改.当你在纸带的最后犯了点小错误,你就得扔掉,然后重新打洞,打洞,打洞.我对这些都信手拈来的.

 

煞费苦心地请求额外的一些纸带充满乐趣,短暂的快乐.接下来就是把这些纸带放入一个由计算器系统监控的篮子里.几个小时后,程序的输出会放在归档文件中.当然,程序中会有错误的.经过又一段时间的煎熬(比之前的相对短一点啦),我的纸带又放入篮子里了.

 

没几次和计算机交互,我就会意识到计算机仅仅会告诉我程序中的第一个错误.终于,在第三天,结果中不再有错误提示信息.但是,这是一个答案---一个错误的答案.我的程序是一个简单的二次方程求解,答案很简单,就是2和3,可是程序跑出来的结果是:1.999997和3.000001.我经过了这么多的曲折,却得到了一个错误的答案.

 

我现在可以很快地写出一个二次方程式求解:

> quadratic.formula <- function (a, b, c)
{
  rad <- b^2 - 4 * a * c
  if(is.complex(rad) || all(rad >= 0)) 
  {
    rad <- sqrt(rad)
  }   
  else 
  {
   rad <- sqrt(as.complex(rad))
  }
  cbind(-b - rad, -b + rad) / (2 * a)
}

> quadratic.formula(1, -5, 6)
    [,1] [,2]
[1,] 2    3

> quadratic.formula(1, c(-5, 1), 6)
     [,1]            [,2]
[1,] 2.0+0.000000i   3.0+0.000000i
[2,] -0.5-2.397916i  -0.5+2.397916i</span>

这个程序比我之前写的那个更有普遍意义,更能说明得到正确答案2和3的问题.R仅仅是输出,所以大多数数值错误是被隐藏起来的.我们可以通过减去正确答案来看它是怎么错误的.

> quadratic.formula(1, -5, 6) - c(2, 3)
    [,1] [,2]
[1,] 0    0

好极了,这种情况下我们得到的答案是正确的.但是如果我们将问题稍微改变一点,错误就出现了:

> quadratic.formula(1/3, -5/3, 6/3)
    [,1] [,2]
[1,] 2    3

> print(quadratic.formula(1/3, -5/3, 6/3), digits=16) [,1] [,2]
[1,] 1.999999999999999 3.000000000000001

> quadratic.formula(1/3, -5/3, 6/3) - c(2, 3)
     [,1]          [,2]
[1,] -8.881784e-16 1.332268e-15

R输出答案有时是美好的祝福,同时也是一个诅咒.R在隐藏数值错误方面非常棒,以至于我们很容易忘记错误的存在.不要忘记.

无论何时,只要你在做浮点运算,即使是很简单的,你都应该假定这里将会有数值错误.如果恰好没有错误,将它视为快乐的意外---不是你原因.你可以用函数 all.equal()代替’==’ 去检查浮点数是否相等.

如果有时一个逻辑整数被计算过,你需要用函数round()来确保它真的是一个整数.

不要把数值错误和一个错误搞混了,一个错误是一次计算被错误地执行了.数值错误是由于有限的数字表达(输出)造成的明显的噪声.例如:1/3被表达成33%.

我们再看了善良异教徒信仰的另一面---输出的即是全部.

> 7/13 - 3/31
[1] 0.4416873

R的输出---默认情况下---是一个手工的截取,不是计算机运行结果的全部.

> print(7/13 - 3/31, digits=16)
[1] 0.4416873449131513

很多总结类函数甚至对它们的输出限制的更多.

> summary(7/13 - 3/31)
Min.    1st Qu.  Median  Mean    3rd Qu.  Max.
0.4417  0.4417   0.4417  0.4417  0.4417   0.4417</span>

因有限的算术表达造成的数值错误不仅使我们对问题感到困惑,有时也会让我们对问题感到困惑.在数学中,矩阵的秩是一些明确的整数.在计算中,矩阵的秩是一个模糊的概念.由于特征值不必是明确的零或者非零,秩不必是一个确定的数.

我们下降到第一站的边缘,这里由米诺斯守卫,他咬着自己的牙齿.他的尾巴在自己周围摇动着,次数标示着他面前的罪人的级别。

完.