algorithm - 如何向该示例添加并行计算？

Question

我有一个算法可以同步计算给定段上的某个积分。我想使用 Control.Parallel 库，或者更确切地说 par::a -> b -> b 来为这个算法添加并行计算。我该怎么做？

integrate :: (Double -> Double) -> Double -> Double -> Double
integrate f a b =
  let
    step     = (b - a) / 1000
    segments = [a + x * step | x <- [0..999]]
    area x   = step * (f x + f (x + step)) / 2
  in sum $ map area segments

Answer 1

从它的外观来看，您正在尝试使用梯形规则在从 b 到 a 的区域上逼近函数 f 的积​​分.您尝试并行化代码是正确的，但是尝试存在几个问题:

• 首先，您需要一个工作窃取调度程序才能获得任何好处，因为 par 不太可能给您提速
• 其次，它的实现方式是每个中间点 f(x) 计算两次，除了边界点 f(a) 和 f( b)

不久前我需要这个功能，所以我将它添加到 massiv 库中:trapezoidRule ，它方便地解决了上述两个问题并避免了列表的使用。

这是一个开箱即用的解决方案，但它不会自动并行计算，因为正在计算数组的一个元素(它旨在估计多个区域的积分)

integrate' :: (Double -> Double) -> Double -> Double -> Double
integrate' f a b = trapezoidRule Seq P (\scale x -> f (scale x)) a d (Sz1 1) n ! 0
  where
    n = 1000
    d = b - a

作为完整性检查:

λ> integrate (\x -> x * x) 10 20 -- implementation from the question
2333.3335
λ> integrate' (\x -> x * x) 10 20
2333.3335

这是一个将进行自动并行化并避免冗余评估的解决方案:

integrateA :: Int -> (Double -> Double) -> Double -> Double -> Double
integrateA n f a b =
  let step = (b - a) / fromIntegral n
      sz = size segments - 1
      segments = computeAs P $ A.map f (enumFromStepN Par a step (Sz (n + 1)))
      area y0 y1 = step * (y0 + y1) / 2
      areas = A.zipWith area (extract' 0 sz segments) (extract' 1 sz segments)
   in A.sum areas

由于列表融合，如果您的解决方案使用列表，则不会分配，因此对于简单的情况会非常快。在上面的解决方案中，将分配一个大小为 n+1 的数组，以促进共享并避免双重函数评估。由于调度，还会遇到额外的成本，因为 fork 线程不是免费的。但最后，对于非常昂贵的函数和非常大的 n，有可能在四核处理器上获得 ~x3 倍的加速。

以下是将高斯函数与 n = 100000 集成的一些基准:

benchmarking Gaussian1D/list
time                 3.657 ms   (3.623 ms .. 3.687 ms)
                     0.999 R²   (0.998 R² .. 1.000 R²)
mean                 3.627 ms   (3.604 ms .. 3.658 ms)
std dev              80.50 μs   (63.62 μs .. 115.4 μs)

benchmarking Gaussian1D/array Seq
time                 3.408 ms   (3.304 ms .. 3.523 ms)
                     0.987 R²   (0.979 R² .. 0.994 R²)
mean                 3.670 ms   (3.578 ms .. 3.839 ms)
std dev              408.0 μs   (293.8 μs .. 627.6 μs)
variance introduced by outliers: 69% (severely inflated)

benchmarking Gaussian1D/array Par
time                 1.340 ms   (1.286 ms .. 1.393 ms)
                     0.980 R²   (0.967 R² .. 0.989 R²)
mean                 1.393 ms   (1.328 ms .. 1.485 ms)
std dev              263.3 μs   (160.1 μs .. 385.6 μs)
variance introduced by outliers: 90% (severely inflated)

旁注建议。切换到辛普森规则将为您提供更好的近似值。 massiv 中提供了实现;)

编辑

这是一个非常有趣的问题，我决定看看在没有任何数组分配的情况下如何实现它。这是我想出的:

integrateS :: Int -> (Double -> Double) -> Double -> Double -> Double
integrateS n f a b =
  let step = (b - a) / fromIntegral n
      segments = A.map f (enumFromStepN Seq (a + step) step (Sz n))
      area y0 y1 = step * (y0 + y1) / 2
      sumWith (acc, y0) y1 =
        let acc' = acc + area y0 y1
         in acc' `seq` (acc', y1)
   in fst $ A.foldlS sumWith (0, f a) segments

上述方法在常量内存中运行，因为确实创建的少数数组不是由内存支持的真实数组，而是延迟数组。通过围绕折叠累加器的一些技巧，我们可以共享结果，从而避免双重函数评估。这导致惊人的加速:

benchmarking Gaussian1D/array Seq no-alloc
time                 1.788 ms   (1.777 ms .. 1.799 ms)
                     1.000 R²   (0.999 R² .. 1.000 R²)
mean                 1.787 ms   (1.781 ms .. 1.795 ms)
std dev              23.85 μs   (17.19 μs .. 31.96 μs)

上述方法的缺点是它不容易并行化，但并非不可能。拥抱你自己，这是一个可以在 8 种功能上运行的怪物(硬编码，在我的例子中是 4 个具有超线程的内核):

-- | Will not produce correct results if `n` is not divisible by 8
integrateN8 :: Int -> (Double -> Double) -> Double -> Double -> Double
integrateN8 n f a b =
  let k = 8
      n' = n `div` k
      step = (b - a) / fromIntegral n
      segments =
        makeArrayR D (ParN (fromIntegral k)) (Sz1 k) $ \i ->
          let start = a + step * fromIntegral n' * fromIntegral i + step
           in (f start, A.map f (enumFromStepN Seq (start + step) step (Sz (n' - 1))))
      area y0 y1 = step * (y0 + y1) / 2
      sumWith (acc, y0) y1 =
        let acc' = acc + area y0 y1
         in acc' `seq` (acc', y1)
      partialResults =
        computeAs U $ A.map (\(y0, arr) -> (y0, A.foldlS sumWith (0, y0) arr)) segments
      combine (acc, y0) (y1, (acci, yn)) =
        let acc' = acc + acci + area y0 y1
         in acc' `seq` (acc', yn)
   in fst $ foldlS combine (0, f a) partialResults

唯一真正分配的数组是用来保存 partialResults 的，它总共有 16 个 Double 元素。速度提升没有那么剧烈，但仍然存在:

benchmarking Gaussian1D/array Par no-alloc
time                 960.1 μs   (914.3 μs .. 1.020 ms)
                     0.968 R²   (0.944 R² .. 0.990 R²)
mean                 931.8 μs   (900.8 μs .. 976.3 μs)
std dev              129.2 μs   (84.20 μs .. 198.8 μs)
variance introduced by outliers: 84% (severely inflated)