你如何在 Repa 中计算 a[i] = f(a[i-1]) ？答案

【问题标题】：How do you compute a[i] = f(a[i-1]) in Repa?你如何在 Repa 中计算 a[i] = f(a[i-1]) ？
【发布时间】：2012-06-27 14:36:53
【问题描述】：

是否可以在 Repa 中计算取决于过去值（即较小的索引）的数组？给出了数组的初始部分（例如，a[0]）。（请注意，我使用类似 C 的符号来表示数组的元素；请不要混淆。）

我阅读了tutorial 并快速检查了hackage，但我找不到执行此操作的函数。

（我想在一维数组中进行这种计算在 Repa 中没有意义，因为你不能并行化它。但我认为你可以在二维或更多维情况下并行化它。）

编辑：可能我应该更具体地说明我想使用什么样的f。由于在a[i] 是标量的情况下无法并行化，所以让我们关注a[i] 是N 暗向量的情况。我不需要a[i] 更高维（例如矩阵），因为您可以将其“展开”为向量。所以，f 是一个将 R^N 映射到 R^N 的函数。

大部分情况是这样的：

b = M a[i-1]
a[i][j] = g(b)[j]

其中b 是一个 N 暗向量，M 是一个 N × N 矩阵（不假设稀疏性），g 是一些非线性函数。我想在给定a[0]、g 和M 的情况下为i=1,..N-1 计算它。我希望有一些通用的方法来（1）并行化这种类型的计算和（2）使b等中间变量的分配有效（在类C语言中，你可以重用它，这会很好如果 Repa 或类似的库可以在不破坏纯度的情况下像魔术一样做到这一点）。

【问题讨论】：

对于关联f，可以并行化，称为“扫描”。 en.wikipedia.org/wiki/Prefix_sum 不过，我在 Repa 文档中找不到 scan。
您可以使用repa 模板来做到这一点，hackage.haskell.org/packages/archive/repa/2.0.2.1/doc/html/…。但另请参阅stackoverflow.com/questions/6170008/…
@Heatsink 扫描不需要一个系列作为输入吗？对我来说，这更像是unfold。
哦，这不是扫描。我将等式误读为a[i] = f(a[i], a[i-1])。其实更像take n $ iterate f z。
@Heatsink 是的，这正是我的意思。您认为是否可以在 Repa 中进行这种计算，尤其是当 f 是多维的时？另外，我担心使用无限列表并比较类 C 语言会产生某种开销，您可以在计算之前分配内存。与裸 Haskell 列表相比，我希望 Repa 减少这种开销。

标签： haskell repa

【解决方案1】：

我看不到 Repa 的做法。但是有 Vector：Data.Vector.iterateN 构建您想要的向量。然后Data.Array.Repa.fromUnboxed 将其从 Vector 转换为 Repa。

iterateN :: Int -> (a -> a) -> a -> Vector aSource

O(n) 对值应用 n 次函数。第零个元素是原始值。

【讨论】：

谢谢。这在进行一维计算时应该足够了。可能我不需要关心雷帕。我只想要在 Haskell 的速度方面与（未优化的）C 相媲美的东西。这可以是多维的吗？（a 可以是另一个向量吗？）如果可以，可以并行化吗？
嗯，这取决于您要询问的是哪种多维概括。在 1D 情况下，显然不可能并行化，因为函数本身本质上是顺序的。您能否更具体地说明您要询问的多维概括是什么样的？
当然。假设a[i] 是一个向量（所以a 是一个二维数组），我想计算a[i] = M a[i-1]，给定一个初始值a[0] 和一个矩阵M，直到i=N-1。使用我用于标题的符号，它的f(x) = M x 现在f 被矢量化了。你可以想到任何一种f，但基本上它只使用前一步的值。
好吧，如果你是幸运的f 是线性的，有几个优化矩阵乘法的选项，我会先尝试。但这取决于M 的形式。由于我不认为您的 M 是稀疏的（或者是吗？），您可以寻找适合您的矩阵类型的 this algorithm 之类的东西。不过，一般来说，每个a[i][j] 都可以依赖于整个a[i-1]，所以我猜你几乎无法并行化任意矩阵。
你是对的，我没有假设稀疏M。我应该澄清一下。我认为 Repa 能够并行化任何类型的向量到向量计算，而不仅仅是线性计算。我错了吗？我知道你不能跨步骤并行计算，但我认为可以在步骤内并行计算（例如，向量到向量计算）。

【解决方案2】：

编辑：实际上，我想我误解了这个问题。我会在这里留下我的答案，以防它对其他人有用......

你可以使用traversehttp://hackage.haskell.org/packages/archive/repa/3.2.1.1/doc/html/Data-Array-Repa.html#v:traverse:

Prelude Data.Array.Repa R> let x = fromListUnboxed (Z :. 10 :: DIM1) [1..10]
Prelude Data.Array.Repa R> R.computeUnboxedS $ R.traverse x (\ (Z :. i) -> (Z :. (i - 1))) (\f  (Z :. i) -> f (Z :. (i + 1)) - f (Z :. i))
AUnboxed (Z :. 9) (fromList [1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0])

剖析它：

    R.computeUnboxedS $                            -- force the vector to be "real"
    R.traverse x                                   -- traverse the vector
    (\ (Z :. i) -> (Z :. (i - 1)))                 -- function to get the shape of the result
    (\f (Z :. i) -> f (Z :. (i + 1)) - f (Z :. i)) -- actual "stencil"

将它扩展到多维数组应该是微不足道的。

【讨论】：

您仅根据源数组计算结果数组。我要计算的是一个数组，其元素取决于其前一步。你可能认为标题中的i 是数组的大小。我的意思是索引。
结合使用Data.Vector.Unboxed.scanl1' 和Repa.toUnboxed 怎么样？例如，cumsum 将是：R.fromUnboxed (R.extent y) $ U.scanl1' (+) $ (R.toUnboxed y)。 [丑得要命，而且可能无法扩展到更高的维度。]
cumsum 仍然是从源数组到目标数组的计算。它取决于 source 数组的前面部分，而不是目标。正如上述 cmets 中提到的 Heatsink 和 phg ，展开是我所需要的。我想知道如何在 Repa（或 Haskell 中的任何其他面向数组的库）中展开。