问题 浮点比较不可再现性


我和我的博士学生在物理数据分析环境中遇到了一个问题,我可以使用一些洞察力。我们的代码可以分析来自其中一个LHC实验的数据,这些实验可以产生不可重复的结果。特别是从中得到的计算结果 相同的二进制,跑上去 同一台机器 连续执行可能会有所不同。我们知道许多不同的不可再生性来源,但已排除了通常的嫌疑人。

当比较两个名义上具有相同值的数字时,我们已经将问题跟踪到(双精度)浮点比较操作的不可再现性。由于分析中的先前步骤,偶尔会发生这种情况。我们刚刚找到一个例子来测试一个数字是否小于0.3(注意我们从不测试浮点值之间的相等性)。事实证明,由于探测器的几何形状,计算有时可能产生精确到0.3(或其最接近的双精度表示)的结果。

我们非常清楚比较浮点数以及FPU中过度精度可能影响比较结果的缺陷。我想回答的问题是“为什么结果不可复制?”是因为FPU寄存器加载或其他FPU指令没有清除多余的位,因此先前计算中的“剩余”位会影响结果吗? (这似乎不太可能)我在另一个论坛上看到一个建议,即进程或线程之间的上下文切换也可能导致浮点比较结果的变化,因为FPU的内容存储在堆栈中,因此被截断。任何关于这些=或其他可能的解释的评论将不胜感激。


6984
2018-02-15 19:35


起源

你能否加入关于上下文切换的建议的参考?虽然我可以想象一个处理器移动累加器数据并丢弃位,但这种机制对我来说似乎不是一个好的解释,而且一些更详细的信息可能会很有趣。 - Coffee on Mars
也许使用不同的编译器优化标志可能会解决此问题 - tkerwin
@Coffee on Mars:这将是我的建议,所以我想我可以解释:)问题是FPU可以在其寄存器中使用更多位,在一些最近的处理器中多达80位用于双打。现在,在单线程环境中,FPU将能够以该精度执行所有操作,您将获得一个结果。如果将其他线程/进程添加到混合中,当OS执行上下文切换时,它必须以80位双精度存储80位寄存器的值,从而失去精度。 - David Rodríguez - dribeas
不是答案,而是'猜测'。如果它以float的整数形式进行比较,那么工作软件会更明智,比如将float转换为unsigned long的sign-expoent-mantissa形式。这也意味着以该形式存储实验的值,因此当数字彼此相邻时,您不会遇到问题。 - Giuliano
为了测试你的“多余位”理论,我可能会建议在比较之前声明并设置一个新变量。这会降低性能,并且可能无论如何都不是“解决方案”,但这将是拒绝您的假设的有趣方式。 - Matthew


答案:


猜测,正在发生的事情是,您的计算通常在FPU内部执行一些额外的精度,并且仅在特定点处舍入(例如,当您将结果赋值给某个值时)。

然而,当存在上下文切换时,必须保存和恢复FPU的状态 - 并且这些额外位是至少非常公平的。  在上下文切换中保存和恢复。当它发生时,这可能不会引起重大变化,但是如果(例如)你稍后从每个变量中减去固定数量并乘以剩下的数量,那么差异也会成倍增加。

需要明确的是:我怀疑“遗留下来”是罪魁祸首。相反,它会丢失额外的位,导致计算中稍微不同点的舍入。


6
2018-02-15 19:53



我敢打赌,操作系统正在使用FSAVE / FRSTOR指令,它实际上转储并恢复x87状态的所有位。这包括内部寄存器的所有80位。假设x86具有32位操作系统,Brian忘记告诉我们他是否正在使用。 :-) - Bo Persson
@Bo Persson:FP寄存器有80个 可见 比特,但大多数也具有至少一个或两个“保护位”,以在80中给出可见的正确圆形LSB。 - Jerry Coffin
所有这些转换和状态保存/恢复都是确定性的。保存的FP状态是完美的副本。 - zvrba


什么平台?

大多数FPU可以在内部存储比ieee双重表示更高的精度 - 以避免中间结果中的舍入误差。通常有一个编译器切换来交换速度/准确性 - 请参阅 http://msdn.microsoft.com/en-us/library/e7s85ffb(VS.80).aspx


4
2018-02-15 19:42



这将如何产生非确定性结果? - zvrba
上下文切换如Jerry所述 - Martin Beckett


我做的:

#include <stdio.h>
#include <stdlib.h>

typedef long double ldbl;

ldbl x[1<<20];

void hexdump( void* p, int N ) {
  for( int i=0; i<N; i++ ) printf( "%02X", ((unsigned char*)p)[i] );
}

int main( int argc, char** argv ) {

  printf( "sizeof(long double)=%i\n", sizeof(ldbl) );

  if( argc<2 ) return 1;

  int i;
  ldbl a = ldbl(1)/atoi(argv[1]);

  for( i=0; i<sizeof(x)/sizeof(x[0]); i++ ) x[i]=a;

  while(1) {
    for( i=0; i<sizeof(x)/sizeof(x[0]); i++ ) if( x[i]!=a ) {
      hexdump( &a, sizeof(a) );
      printf( " " );
      hexdump( &x[i], sizeof(x[i]) );
      printf( "\n" );
    }
  }

}

使用/ Qlong_double与IntelC编译,以便它产生:

;;;     for( i=0; i<sizeof(x)/sizeof(x[0]); i++ ) if( x[i]!=a ) {

        xor       ebx, ebx                                      ;25.10
                                ; LOE ebx f1
.B1.9:                          ; Preds .B1.19 .B1.8
        mov       esi, ebx                                      ;25.47
        shl       esi, 4                                        ;25.47
        fld       TBYTE PTR [?x@@3PA_TA+esi]                    ;25.51
        fucomp                                                  ;25.57
        fnstsw    ax                                            ;25.57
        sahf                                                    ;25.57
        jp        .B1.10        ; Prob 0%                       ;25.57
        je        .B1.19        ; Prob 79%                      ;25.57
[...]
.B1.19:                         ; Preds .B1.18 .B1.9
        inc       ebx                                           ;25.41
        cmp       ebx, 1048576                                  ;25.17
        jb        .B1.9         ; Prob 82%                      ;25.17

并开始了10个具有不同“种子”的实例。如你所见,它比较了 10字节长从内存加倍,FPU堆栈上有一个,所以在这种情况下 操作系统不能保持完全精确,我们肯定会看到错误。 好吧,他们仍在运行而没有检测到任何东西......这不是真的 令人惊讶的是,因为x86具有立即保存/恢复整个FPU状态的命令, 无论如何,一个不能保持完全精度的操作系统将完全被破坏。

所以要么是它的一些独特的OS / cpu /编译器,要么是不同的比较结果 在更改程序中的某些内容并重新编译它之后生成它 程序中的错误,例如。缓冲区溢出。


2
2018-02-15 20:48





该程序是多线程的吗?

如果是的话,我会怀疑是一场竞争。

如果不是,程序执行是确定性的。在给定相同输入的情况下获得不同结果的最可能的原因是未定义的行为,即程序中的错误。读取未初始化的变量,过时的指针,覆盖堆栈上某些FP编号的最低位等。可能性是无穷无尽的。如果你在linux上运行它,尝试在valgrind下运行它,看看它是否揭示了一些错误。

顺便说一下,你是如何将问题缩小到FP比较的?

(远射:硬件故障?例如,RAM芯片出现故障可能会导致数据在不同情况下以不同方式被读取。但是,这可能会很快崩溃操作系统。)

任何其他解释都是不可信的 - 操作系统或硬件中的错误不会被发现很长时间。


2
2018-02-15 22:01



虽然程序在多进程上下文中运行,但它不是多线程的。它已经通过Valgrind(包括Memcheck)运行,没有任何问题。由于没有确定不可再生性的来源而感到沮丧,我们采用了低技术调试 - 将价值与分析“削减”(0.3)相比较倾向于cout。两个单独的执行都打印0.3到15有效数字精度(... << std :: setw(15)<< dR)但是执行dR与0.3比较的后续行产生不同的结果。 - Brian Cole
为什么不记录变量值的十六进制转储? log2(10 ^ 15)是~49.82,并且在long double中有更多的尾数位。 - Shelwien
好吧,虽然valgrind并不总能检测到所有问题。我建议你在机器代码级别上执行调试。在print语句上设置断点,按指令执行指令,并在每条指令后检查所有寄存器的状态。一定要查看XMM寄存器。你可以在这里粘贴有问题的代码(print语句和比较本身是什么)? - zvrba
我同意这种分析。我真诚地怀疑上下文切换的事情,除非我当然被证明是错误的。我不相信寄存器状态被错误地压入堆栈。 - Alexandre C.
在2016年,如果这仍然是一个问题,你可以重新编译 -fsanitize=undefined 代替。这应该可以捕获更多的算术错误。 - Jean-Michaël Celerier


CPU的内部FPU可以以高于double或float的精度存储浮点。只要寄存器中的值必须存储在其他地方,包括当存储器被换出到缓存中时,这些值就必须被转换(我知道这是事实),并且该核心上的上下文切换或OS中断听起来像另一个简单的来源。当然,OS中断或上下文切换的时序或非热存储器的交换是完全不可预测的,并且应用程序无法控制。

当然,这取决于平台,但您的描述听起来像是在现代桌面或服务器上运行(因此x86)。


1
2018-02-15 19:59



很好的尝试,但所有这些转换都是 确定性。 FP状态通过专用指令保存,这种保存/恢复机制不会丢失信息。 - zvrba
@zvrba:确定性并不意味着无损。如果非确定性地调用,则确定性的指令不具有确定性。 - Puppy
咦?所有FP指令在给定相同输入的情况下产生完全相同(偶数按位!)的结果,无论它们在何种情况下被调用。 (可能的例外是其中一个参数是NaN的无穷大。) - zvrba


我将合并David Rodriguez和Bo Persson的一些评论并做出疯狂猜测。

在使用SSE3指令时是否可以进行任务切换?基于此 英特尔关于使用SSE3指令的文章 保存寄存器状态FSAVE和FRESTOR的命令已被FXSAVE和FXRESTOR取代,FXSAVE和FXRESTOR应该处理 累加器的全长。

在x64机器上,我认为“错误”指令可能包含在某些外部编译库中。


0
2018-02-15 22:33



这些指令通常不被用户模式程序使用,并且如果linux内核存在这样的错误(即,使用错误的指令来保存/恢复FP状态),那么很久以前就会发现它。 - zvrba
是的,我相信你是对的。另一方面,良好的任务切换不应该在寄存器中留下位,并且在问题的评论中已经讨论了这种可能性。我将在这里留下答案,作为可能在问题的装配方面查询的参考。 - Coffee on Mars
这里没有人证明任务切换会在寄存器中留下位。用于FP状态切换的硬件指令当然不会留下任何垃圾。 - zvrba


你当然会打 GCC Bug n°323,其他指出是由于FPU的精度过高。

方案:

  • 使用SSE(或AVX,它是2016 ......)来执行计算
  • 使用 -ffloat-store 编译开关。来自海湾合作委员会的文件。

不要将浮点变量存储在寄存器中,并禁止可能会改变从寄存器或内存中获取浮点值的其他选项。
   这个选项可以防止诸如68000之类的机器上出现不希望的过度精度,其中浮动寄存器(68881)的精度比double应该具有的精度高。   同样适用于x86架构。   对于大多数程序来说,多余的精度只会很好,但有些程序依赖于IEEE浮点的精确定义。   在修改它们以将所有相关的中间计算存储到变量中之后,对这些程序使用-ffloat-store。


0
2017-10-08 06:21