为什么这个 Monte Carlo Haskell 程序这么慢?

Why is this Monte Carlo Haskell program so slow?

更新:

我的编译命令是ghc -O2 Montecarlo.hs。我的运行dom版本是random-1.1,ghc版本是8.6.4,我的系统是macOS Big Sur 11.1(Intel芯片)。我用来测试速度的命令是time ./Montecarlo 10000000,它returns的结果是real 0m17.463s, user 0m17.176s, sys 0m0.162s.


下面是一个Haskell程序,使用Monte Carlo计算pi。但是,当输入为1000万时,程序运行持续了20秒。同样逻辑编写的C程序只用了0.206秒。为什么会这样,我怎样才能加快速度?谢谢大家

这是 Haskell 版本:

import System.Random
import Data.List
import System.Environment

montecarloCircle :: Int -> Double
montecarloCircle x
   = 4*fromIntegral
        (foldl' (\x y -> if y <= 1 then x+1 else x) 0
           $ zipWith (\x y -> (x**2 + y**2))
                     (take x $ randomRs (-1.0,1) (mkStdGen 1) :: [Double])
                     (take x $ randomRs (-1.0,1) (mkStdGen 2) :: [Double]) )
        / fromIntegral x

main = do
    num <- getArgs
    let n = read (num !! 0)
    print $ montecarloCircle n

这是C版本:

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

#define N       10000000

#define int_t   long // the type of N and M

// Rand from 0.0 to 1.0
double rand01()
{
    return rand()*1.0/RAND_MAX;
}

int main()
{
        srand((unsigned)time(NULL));

        double x,y;
        int_t M = 0;

        for (int_t i = 0;i < N;i++)
        {
                x = rand01();
                y = rand01();
                if (x*x+y*y<1) M++;
        }
        double pi = (double)4*M/N;
        printf("%lf\n", pi);
}

在数值应用程序中,Haskell 代码往往比用经典命令式语言(例如 C/C++/Fortran)编写的等效代码慢 2 到 5 倍。然而,Haskell 代码慢了 100 倍是非常意外的!

完整性检查:重现结果

使用带有 Intel Celeron N2840 64 位 cpu、运行ning Linux 内核 5.3、GLIBC 2.28、GCC 8.3、GHC 8.2.2、GHC 随机的有点旧的笔记本-1.1,我们得到时间:

C code,        gcc -O3:    950 msec
Haskell code,  ghc -O3:    104 sec

确实如此,使用这些配置,Haskell 代码 运行 比 C 代码慢大约 100 倍。

这是为什么?

首先要注意的是随机数生成之上的 π 计算算法看起来非常简单,因此 运行 次可能主要由 2*1000 万个伪随机数的生成主导。此外,我们完全有理由期待 Haskell 和 C 使用完全不相关的随机数生成算法。

因此,与其将 Haskell 与 C 进行比较,不如比较这两种语言碰巧使用的随机数生成算法,以及它们各自的实现。

在 C 中,语言标准没有指定应该使用哪个算法函数 srand(),但它往往是旧的、简单且快速的算法。

另一方面,在 Haskell,我们通常会看到很多关于 StdGen 效率低下的抱怨,就像 2006 年的这里:

https://gitlab.haskell.org/ghc/ghc/-/issues/427

其中一位领先的 Haskell 名人提到 StdGen 可能会快 30 倍。

幸运的是,帮助已经在路上了。 recent blog post explains how the new Haskell random-1.2 package 解决了 StdGen 速度不足的问题,使用了一种完全不同的算法 splitmix

伪随机数生成 (PRNG) 领域是一个相当活跃的领域。算法通常会被更新更好的算法淘汰。为了透视起见,这是相对较新的 (2018) review paper on the subject.

转向更新、更好的 Haskell 组件:

使用另一台功能稍强的机器,配备 Intel Core i5-4440 3GHz 64 位 cpu、运行ning Linux 内核 5.10、GLIBC 2.32、GCC 10.2,关键是Haskell 随机包 1.2:

C code,        gcc -O3:    164 msec
Haskell code,  ghc -O3:    985 msec

所以 Haskell 现在“只是”慢了 6 倍,而不是 100 倍。

而且我们仍然必须解决 Haskell 不得不使用 x**2+y**2 而不是 C 得到 x*x+y*y 的不公平,这个细节在 之前并不重要随机 1.1。这给了我们 379 毫秒! 所以我们又回到了通常的 2X-5X 范围,以进行 Haskell 与 C 速度的比较。

请注意,如果我们向 Haskell 可执行文件请求 运行 统计信息,我们将得到以下输出:

$ time q66441802.x +RTS -s -RTS 10000000
3.1415616
     925,771,512 bytes allocated in the heap
     ...
  Alloc rate    2,488,684,937 bytes per MUT second

  Productivity  99.3% of total user, 99.3% of total elapsed

real   0m0,379s
user   0m0,373s
sys    0m0,004s

所以发现Haskell沿途分配了接近1GB的内存,这有助于理解速度差异。

代码修复:

我们注意到 C 代码使用单个随机序列,而 Haskell 代码使用两个随机序列,两次调用 mkStdGen 以数字 1 和 2 作为种子。这不仅不公平,而且不正确。设计 PRNG 算法的应用数学家非常注意确保任何单个序列都具有正确的统计属性,但基本上不保证不同序列之间可能存在的相关性。甚至使用种子作为单个全局序列的偏移量也不是闻所未闻。

这已在以下代码中修复,不会显着改变性能:

computeNorms :: [Double] -> [Double]
computeNorms []  = []
computeNorms [x] = []
computeNorms (x:y:xs2) = (x*x + y*y) : (computeNorms xs2)

monteCarloCircle :: Int -> Double
monteCarloCircle nn =
    let
         randomSeed   =  1
         gen0         =  mkStdGen randomSeed
                      -- use a single random serie:
         uxys         =  (randomRs  (-1.0, 1.0)  gen0) :: [Double]
         norms        =  take nn (computeNorms uxys)
         insiderCount =  length  $  filter (<= 1.0) norms
    in
         (4.0::Double) * ((fromIntegral insiderCount) / (fromIntegral nn))

旁注:

新的 Haskell random-1.2 包已在此 recent SO question 中提及,尽管是在新的 monadic 接口的上下文中。

一句话总结:

假设练习的目标是衡量 C 和 Haskell 语言的相对 运行 时间速度,基本上有两种可能性:

一个是完全避免使用 PRNG,因为使用的算法不同。

另一种是通过手动为两种语言提供一种相同的算法来控制随机数的生成。例如,可以使用公开可用的 MRG32k3a algorithm designed by Pierre L'Ecuyer. Candidate Haskell MRG32k3a implementation here(截至 2021 年 3 月)。留作 reader.

的练习