MATLAB 单精度浮点运算一个 bug （使用时要慎重）

最近用单精度运算处理一段大数据，无意间发现了一个奇葩的 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 :)

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

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

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

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

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

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

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)) 就都没问题了。

xp x86 2010b
结果都好着。。。
无论是command window还是m文件。。

winner245 发表于 2015-8-19 06:20
另外，即使输入带下标，只要没有输出赋值（或者说使用默认的ans做输出），就不会有问题，比如，直接写成 ...

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

winner245 发表于 2015-8-19 06:20
另外，即使输入带下标，只要没有输出赋值（或者说使用默认的ans做输出），就不会有问题，比如，直接写成 ...

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

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又会再现，再点运行结果又恢复正常

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))  
就没问题了

这个问题前几天楼主告诉我的时候讨论了一下，已经基本确定是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

honglei.chen 发表于 2015-8-19 05:40
我可以在15a里复现，但是15b的prerelease似乎已经没问题了

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

15b 的 prerelease 已经不存在这个bug了吗？
如果是这样的话，也就没有必要report这个bug了

winner245 发表于 2015-8-20 01:56
15b 的 prerelease 已经不存在这个bug了吗？
如果是这样的话，也就没有必要report这个bug了 ...

15b没有问题

2010b:

x =

    0.8147


x =

  3.6893e+019


y =

    0.8147


z =

    0.8147


ans =

     0     1