MATLAB中文论坛

 找回密码
 注册

QQ登录

只需一步,快速开始

扫一扫,访问微社区

查看: 2622|回复: 12|关注: 0

[我分享] MATLAB 单精度浮点运算一个 bug (使用时要慎重)

[复制链接]

论坛优秀回答者

24

主题

1万

帖子

1633

最佳答案
  • 关注者: 647
发表于 2015-8-19 02:36:29 | 显示全部楼层 |阅读模式
最近用单精度运算处理一段大数据,无意间发现了一个奇葩的 bug:

x = rand('single')
y = x;
z = x;
x(1) = x(1) + eps('single')   
y = y + eps('single')
z(end) = z(end) + eps(single(1.0))  
[isequal(x,y), isequal(y,z)]

如果没有 bug 的话,将得到 x、y、z 完全相等(这可以在双精度运算下得到验证,只需将上面代码里所有 single 换成 double 即可),但最终的结果却是 x 不等于 y,y 和 z 相等。之所以说这个 bug 有点奇葩,主要是因为以下几个现象有些奇怪:

1:多次运行上述代码,不论每次产生的初值 x 如何,最终得到的 x 只能取几个特定的值,如 1.1921e-07、2、-2、3.6893e+19、 -3.6893e+19 (对!连负数都出来了!)

2:x 与 y 的唯一区别在于,x 使用了数组索引而 y 没有,采用了完全相同的运算,然而不带索引的 y 的结果无误,而带索引的 x 明显错误 (索引还能造成差别?)

3:上面看到带索引会出错,可是 z 明明也带了索引(索引换成1或者end其实无所谓),可是 z 的结果却偏偏无误,而这里唯一的不同的是将eps('single') 换成了 eps(single(1.0)) ,可是 help 文档里清清楚楚地记载: eps('single') is equivalent to eps(single(1.0)) ,可他们造成的结果却截然不同!

如果这还不算奇怪,那我们继续往下看:

4:随便在上述代码里加一个断点后运行,x、y、z 的结果都准确无误了 (断点也能make a difference了?)

5:如果在程序里关掉 jit 加速器后,x、y、z 的结果也都准确无误了

通过上面的测试(以及一些额外的测试),似乎可以断定,只要代码里使用的是 eps(single(x)) (如果 x 原本不是single类型)或者 eps(x) (如果 x 已经是single类型),则不会有问题;但一旦代码里写的是 eps('single')  的形式,则很可能会出现一些奇怪的现象。所以,建议在新版本修正这一 bug 之前,尽量不要使用  eps('single') 的写法 。

最后交代一下测试平台:

64 WIN 8 + 64 MATLAB R2015a
64 WIN 10 + 64 MATLAB R2015a

不清楚在其它平台下表现如何,有兴趣的朋友不妨也测试一下,让我们一起来找 bug :)

论坛优秀回答者

3

主题

1万

帖子

859

最佳答案
  • 关注者: 230
发表于 2015-8-19 05:40:25 | 显示全部楼层
本帖最后由 honglei.chen 于 2015-8-19 05:51 编辑

我可以在15a里复现,但是15b的prerelease似乎已经没问题了

从问题的描述来看,这个应该和编译的优化有关,所以就是JIT的问题。这也是为什么关掉JIT就可以的原因。至于debug模式,这个时候是JIT是关掉的,所以这个时候结果正确。

还有就是从我的测试来看,问题更可能和下标的优化有关,而不是eps('single'),不过在优化的时候,这些都很难说

x(1) = x + eps('single') 没问题
x = x(1) + eps('single') 有问题

论坛优秀回答者

24

主题

1万

帖子

1633

最佳答案
  • 关注者: 647
 楼主| 发表于 2015-8-19 06:05:39 | 显示全部楼层
honglei.chen 发表于 2015-8-19 05:40
我可以在15a里复现,但是15b的prerelease似乎已经没问题了

从问题的描述来看,这个应该和编译的优化有关, ...

我觉得下标只是一个可能的因素,实际的情况可能非常复杂。

之前还测试过很多种情况,比如,如果事先将 eps('single') 赋值了,则无论下标情形如何,都不会有问题:

a = eps('single');
x(1) = x + a  
x = x(1) + a  
x(1) = x(1) + a
x = x + a  

但如果直接使用 eps('single'), 则只要输入 x 带下标就有问题,比如以下两种
x = x(1) + eps('single')  
x(1) = x(1) + eps('single')  
只要输入不带下标就没问题,比如以下两种:
x(1) = x + eps('single')   
x = x + eps('single')  

但比上述更复杂的情况还有很多种,感觉规律远没有这么简单,很难简单的归纳全部出错的情形。

论坛优秀回答者

24

主题

1万

帖子

1633

最佳答案
  • 关注者: 647
 楼主| 发表于 2015-8-19 06:20:10 | 显示全部楼层
honglei.chen 发表于 2015-8-19 05:40
我可以在15a里复现,但是15b的prerelease似乎已经没问题了

从问题的描述来看,这个应该和编译的优化有关, ...

另外,即使输入带下标,只要没有输出赋值(或者说使用默认的ans做输出),就不会有问题,比如,直接写成:
x(1) + eps('single')
或者 ans = x(1) + eps('single')
就都没问题,但只要是自己给定了输出变量(不是ans),如
t = x(1) + eps('single')
x = x(1) + eps('single')
x(1) = x(1) + eps('single')
就都会有问题。

另外,在上述所有出错的情形里,只要将 eps('single') 换成 2^-23,或者换成预先赋值的a,a = eps('single') , 或者换成 eps(single(1)) 就都没问题了。



论坛优秀回答者

入门

423 麦片

财富积分


50500


1

主题

1022

帖子

90

最佳答案
  • 关注者: 12
发表于 2015-8-19 08:20:07 | 显示全部楼层
xp x86 2010b
结果都好着。。。
无论是command window还是m文件。。

论坛优秀回答者

3

主题

1万

帖子

859

最佳答案
  • 关注者: 230
发表于 2015-8-19 09:23:42 | 显示全部楼层
winner245 发表于 2015-8-19 06:20
另外,即使输入带下标,只要没有输出赋值(或者说使用默认的ans做输出),就不会有问题,比如,直接写成 ...

因为这是在优化时候出的错,所以一定是一些组合在一起的原因。我不知道具体的细节,但是我觉得你先定义一个变量等于eps('single')实际上就改变了优化的条件,因为整个程序的节点图就不一样了。不过不管怎么说这是一个bug。

论坛优秀回答者

3

主题

1万

帖子

859

最佳答案
  • 关注者: 230
发表于 2015-8-19 09:40:01 | 显示全部楼层
winner245 发表于 2015-8-19 06:20
另外,即使输入带下标,只要没有输出赋值(或者说使用默认的ans做输出),就不会有问题,比如,直接写成 ...

还有我比较好奇你是怎么碰到这个问题的,如果知道你真正要做什么的话可以看看有什么解决的方法。

论坛优秀回答者

权威

3488 麦片

财富积分



19

主题

3730

帖子

753

最佳答案
  • 关注者: 306
发表于 2015-8-19 09:55:56 | 显示全部楼层
64 WIN 7 + 64 MATLAB R2015a:
第一次运行代码的时候也是一样的bug:
x =

    0.1712


x =

  1.1921e-07


y =

    0.1712


z =

    0.1712


ans =

     0     1
但是第二次点运行的时候结果就正常了,如果这时候重启MATLAB或者是修改一下代码,比如加个分号,bug又会再现,再点运行结果又恢复正常

论坛优秀回答者

24

主题

1万

帖子

1633

最佳答案
  • 关注者: 647
 楼主| 发表于 2015-8-19 13:06:44 | 显示全部楼层
honglei.chen 发表于 2015-8-19 09:40
还有我比较好奇你是怎么碰到这个问题的,如果知道你真正要做什么的话可以看看有什么解决的方法。 ...

这个问题的背景大致是这样的:我要处理的数组很长,担心内存不够用,于是将double类型的数组用single来代替了,所以降低了一些精度,使得原本用cumsum对该数组求和后最后一个元素本该为1的,而在单精度计算下,最后一个元素与1的误差比双精度要大一些,于是通过归一化将最后一个元素修正为了1。但考虑到代码的努棒性,人为希望在这个1的基础上加上一个单精度eps间隔,本来可以直接写成
x(end) = 1 + eps(single(1))   % 没有问题
但考虑到eps(single(1))和eps('single') 等价,一开始写的形式是:
x(end) = x(end) + eps('single')
结果就发现出错了。最后就改成了
x(end) = x(end) + eps(x(end))  
就没问题了

论坛优秀回答者

1

主题

9782

帖子

1463

最佳答案
  • 关注者: 262
发表于 2015-8-19 15:27:27 | 显示全部楼层
这个问题前几天楼主告诉我的时候讨论了一下,已经基本确定是JIT的问题
除了1L说的几种情况外还有可能(但是不同机器上测试结果不同,即使matlab版本和系统相同)出现以下情况:
t = isequal(x(1)+eps('single'),x(1)+eps('single'))
t为0,如果放在m文件中执行,第一次(有JIT过程)t为0,重复执行t为1,对文件做出任何修改(例如增加一个空格或换行)后的第一次运行t=0

isequal(x(1)+eps('single'),x(1)+eps('single'))
ans为1

t = isequal(sqrt(x(1)+eps('single')),sqrt(x(1)+eps('single')))
在m文件中执行,文件每次被修改后的前两次运行t=0,从第三次开始t=1

eq1 = @isequal;
t = eq1(sqrt(x(1)+eps('single')),sqrt(x(1)+eps('single')))
t为1
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /2 下一条

快速回复 返回顶部 返回列表