双精度浮点数在接近零的值上计算速度变慢

3

我收到了一位朋友的请求,要分享我在过去某个时候偶然发现的东西。原始帖子可从此处获取。问题陈述可以在这里找到。基本上是一个算法竞赛的网站。

我被放在一个算法问题前,我使用以下代码解决了这个问题:

double dp[80002][50];
class FoxListeningToMusic {
public:
    vector <double> getProbabilities(vector <int> length, int T)  {    
        memset(dp, 0, sizeof(dp));
        int n = length.size();
        for(int i = 0; i < n; i++)
            dp[0][i] = 1.0 / (double)n;

        double mul = 1.0 / (double)n;
        int idx ;
        for(int i = 1; i <= T; i++) {
            for(int j = 0; j < n; j++)  {
                idx = i - length[j];
                if(idx >= 0)  {
                    for(int k = 0; k < n; k++)
                        dp[i][k] += mul * dp[idx][k];
                }
                else
                    dp[i][j] += mul;
                }
            }
        }

        vector<double> v(n);
        for(int i = 0; i < n; i++)
            v[i] = dp[T][i];
        return v;
    }

};

不管代码是否能正确解决问题,至少对于我要讨论的内容来说并不重要。事实是,这段代码在一些测试用例上执行超过了2秒的时间限制。这种情况有些意料之中,因为这里的复杂度是O(T * length.size() ^ 2),如果考虑到问题的约束条件,就会变成2 * 108。然而,有趣的是,我特别针对时间限制测试了我的解决方案。我使用的测试用例似乎是我的解决方案的“最坏情况”:给定长度为50个1和T = 80000。这段代码运行了0.75秒。这远低于2秒的时间限制。
我说我使用的测试用例是最坏情况,因为将执行的指令数取决于内部for循环中的分支条件idx >= 0。如果这是真的,就需要再执行一次循环(复杂度为O(n))。否则,只会执行单个操作O(1)。正如你所看到的,长度越短,这种情况就越多。
尽管有这样的推理,我的问题在以下情况下测试失败:
length = {1, 1, 1, 1, 3, 3, 3, 3, 1, 3, 3, 2, 3, 2, 3, 3, 1, 2, 3, 1, 2, 3, 2,
          1, 3, 1, 1, 1, 2, 3, 2, 3, 2, 2, 1, 3, 1, 1, 3, 1, 3, 1, 3, 2, 3, 1,
          1, 3, 2, 76393} T= 77297.
For this case my program runs for 5.204000 seconds.

我的第一个假设是这个意外的运行时比率(在第一种情况下我们应该期望执行更少的处理器指令)的原因是处理器以某种方式缓存了相似的计算:在我的例子中,计算与长度的所有元素关于对称都是相同的,并且非常聪明的处理器可以利用这一点节省重复相同的指令序列。因此,我尝试组成另一个例子:这次在长度数组中使用不同的值:

length = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
          21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
          39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 77943}
T=80000 runs for  0.813000 seconds. 

在这个例子之后,我不再能够解释为什么这些时间度量如此 - 我的第二个例子似乎需要比测试失败的更多的处理器指令,并且不允许我认为在第一个示例中发生的缓存。实际上,我无法定义这种行为的原因,但我非常确定它应该与处理器缓存或传送带有关。我非常好奇那些实验在不同芯片组上的行为,所以请随意在这里发表评论。
此外,如果有任何比我更了解硬件的人可以解释这种行为,将不胜感激。
在此之前,我应该给自己做个注释 - 在估计算法复杂性时,不要低估处理器优化。有时,它们似乎会显着减少/增加特定示例的分摊速度。

可能是重复的问题:为什么将0.1f更改为0会使性能降低10倍? - Eric Postpischil
2个回答

7
这种奇怪的行为原因是非规格化数。将代码处理这样的数字作为纯零可以极大地加快在这些边缘情况下的代码运行速度。

提示:在这种情况下,非规格化数是指非常接近于零的数字(例如浮点数的10-38;由于@PascalCuoq的更正)。对于这样的数字,处理器在处理时会变得非常慢,因为这个原因(摘自维基百科):

有些系统在硬件上处理非规格化值,就像处理规格化值一样。其他系统将非规格化值的处理留给系统软件,在硬件上仅处理规格化值和零。在软件中处理非规格化值总是导致性能显著降低。

编辑我还在SO上找到了此建议,可以检查数字是否变成了非规格化数。


你使用的操作系统和编译器是什么? - Violet Giraffe
那你做了什么?启用 Denormals Are Zero 标志吗?还是 Flush To Zero?两个都有吗?(另外,使用了什么微架构?) - harold
@harold 我在最内层循环之后放置了一个检查 dp[i][k] = (dp[i][k] < VERY_SMALL_EPSYLON) ? 0 : dp[i][k];。这是比赛的一部分,我没有权利更改标志。 - Boris Strandjev
1
“例如 10^-20”?!单精度浮点数中最大的非规格化数约为10^-38。你对于单精度非规格化数的估计偏差了近10^18倍(双精度则是10^288倍)。 - Pascal Cuoq
1
@BorisStrandjev,我没有提到这个,因为Harold已经提过了,但通常的解决方案是让FPU将非规格化值映射为零。你的解决方案有其优点(可移植性,例如如果集成了几个数值库,它们对FPU状态有自己的要求才能正常工作),但配置FPU刷新为零意味着不必担心极限是什么,还有其他一些事情。 - Pascal Cuoq
显示剩余3条评论

1
另一种解决这种情况的选项是使用定点运算,完全避免使用浮点数。问题说明要求答案精确到1e-9,由于2^64约为10^19,并且最多只进行80000次迭代,因此具有足够的精度。其工作方式是定义一个大常量,例如:
const uint64_t ONE = pow(10,17);

你需要将 uint64_t 数组初始化为 ONE/n 而不是 1.0/double(n),主循环应该像这样:
  for(int i = 1; i <= T; i++) {
    for(int j = 0; j < n; j++)  {
      idx = i - length[j];

      if(idx >= 0)  {
        for(int k = 0; k < n; k++){
          dpi[i][k] += dpi[idx][k];
        }
      }    
      else
        dpi[i][j] += ONE;

    }
    for(int k = 0; k < n; k++){
      dpi[i][k] = dpi[i][k]/n;
    }
  }

理论上,这应该更快,因为您避免了主循环中的浮点运算,内部循环仅由整数加法组成。在我的机器上,性能提升仅约为10%,这表明真正的瓶颈可能是内存访问。但在其他情况下,您可能会看到更大的性能提升。


嗯,我对“/n”的计算有些担忧,它们在你的方法中可能不够精确。你尝试用这种方法解决问题了吗? - Boris Strandjev
是的,我尝试了几个输入,包括您提供的输入。这种方法与“double”解决方案的差异大约在1e-15的数量级上。 - mrip
非常抱歉回复晚了,但我现在已经测试了您的方法,并确认它有效,问题也得到了解决。感谢您分享另一种解决这个特定问题的方法。 - Boris Strandjev

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接