wolfram-mathematica - 使用Compile高效处理小向量列表

Question

通常，我们需要处理由一系列坐标组成的数据:data = {{x1,y1}, {x2,y2}, ..., {xn,yn}}。它可以是2D或3D坐标，也可以是固定长度的小向量的任何其他任意长度列表。

让我通过汇总二维向量列表的简单示例来说明如何使用Compile解决此类问题:

data = RandomReal[1, {1000000, 2}];

首先，明显的版本:

fun1 = Compile[{{vec, _Real, 2}},
  Module[{sum = vec[[1]]},
   Do[sum += vec[[i]], {i, 2, Length[vec]}];
   sum
   ]
  ]

有多快？

In[13]:= Do[fun1[data], {10}] // Timing
Out[13]= {4.812, Null}

其次，不太明显的版本:

fun2 = Compile[{{vec, _Real, 1}},
  Module[{sum = vec[[1]]},
   Do[sum += vec[[i]], {i, 2, Length[vec]}];
   sum
   ]
  ]

In[18]:= Do[
  fun2 /@ Transpose[data],
  {10}
  ] // Timing

Out[18]= {1.078, Null}

如您所见，第二个版本要快得多。为什么？因为至关重要的操作，所以 sum += ...是 fun2中数字的加法，而它是 fun1中任意长度向量的加法。

您可以看到相同的“optimization” in this asnwer of mine的实际应用，但是在与此相关的情况下，可以给出许多其他示例。

现在，在这个简单的示例中，使用 fun2的代码不会比 fun1更长或更复杂，但在一般情况下可能会更好。

我如何告诉Compile它的参数之一不是一个任意的n*m矩阵，而是一个特殊的n*2或n*3矩阵，因此它可以自动进行这些优化，而不是使用通用矢量加法函数添加微小的length-2或length- 3个向量？

附录

为了更清楚地说明发生了什么，我们可以使用 CompilePrint:
CompilePrint[fun1]给

        1 argument
        5 Integer registers
        5 Tensor registers
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        T(R2)0 = A1
        I1 = 2
        I0 = 1
        Result = T(R1)3

1   T(R1)3 = Part[ T(R2)0, I0]
2   I3 = Length[ T(R2)0]
3   I4 = I0
4   goto 8
5   T(R1)2 = Part[ T(R2)0, I4]
6   T(R1)4 = T(R1)3 + T(R1)2
7   T(R1)3 = CopyTensor[ T(R1)4]]
8   if[ ++ I4 < I3] goto 5
9   Return

CompilePrint[fun2]给

        1 argument
        5 Integer registers
        4 Real registers
        1 Tensor register
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        T(R1)0 = A1
        I1 = 2
        I0 = 1
        Result = R2

1   R2 = Part[ T(R1)0, I0]
2   I3 = Length[ T(R1)0]
3   I4 = I0
4   goto 8
5   R1 = Part[ T(R1)0, I4]
6   R3 = R2 + R1
7   R2 = R3
8   if[ ++ I4 < I3] goto 5
9   Return

我选择包含此版本，而不是包含相当长的C版本，在C版本中，时间差异甚至更加明显。

Answer 1

实际上，您的附录几乎足以了解问题所在。对于第一个版本，您在内部循环中调用CopyTensor，这是效率低下的主要原因，因为必须在堆上分配许多小缓冲区，然后再释放它们。为了说明，这是一个不复制的版本:

fun3 =
 Compile[{{vec, _Real, 2}},
   Module[{sum = vec[[1]], len = Length[vec[[1]]]},   
     Do[sum[[j]] += vec[[i, j]], {j, 1, len}, {i, 2, Length[vec]}];
     sum], CompilationTarget -> "C"]

(顺便说一下，我认为速度比较在编译为C时更公平，因为例如Mathematica虚拟机确实更不鼓励嵌套循环)。对于这么小的向量，此函数仍然比您的函数慢，但比 fun1快3倍。

我认为，其余的效率低下是这种方法固有的。您可以将问题分解为单个组件的总和，这使您的第二个函数高效，这是因为您使用了 Transpose之类的结构操作，最重要的是，这使您可以从内部循环中挤出更多指令。因为这是最重要的-内部循环中的指令必须尽可能少。从 CompilePrint可以看到 fun1与 fun3的确如此。在某种程度上，您(针对此问题)找到了一种有效的高级方法来手动展开外循环(位于坐标索引上方的外循环)。您建议的替代方法是，根据矢量维的额外信息，要求编译器自动展开外循环。这听起来似乎是合理的优化，但可能尚未为Mathematica虚拟机实现。

还要注意的是，对于更长的向量(例如20)， fun1和 fun2之间的差异消失了，因为张量复制中的内存分配/解除分配的成本与大规模分配的成本相比微不足道(当您进行大规模分配时，实现效率更高)将向量分配给向量-可能是因为在这种情况下，您可以使用 memcpy之类的东西)。

总而言之，我认为虽然可以自动进行此优化会很好，但至少在此特定情况下，这是一种低级优化，很难指望它是全自动的-即使优化C编译器也并非总是如此执行它。您可以尝试做的一件事是将向量长度硬编码为已编译的函数，然后使用 SymbolicCGenerate(来自 CCodeGenerator`包)生成符号C，然后使用 ToCCodeString生成C代码(或通过其他任何方式获取C代码)代码，然后手动尝试创建和加载该库，从而通过 CreateLibrary的选项启用C编译器的所有优化。我不知道这是否行得通。编辑我实际上怀疑这是否有帮助，因为在生成C代码时，已经使用 goto -s实现了循环，以提高速度，并且这很可能会阻止编译器尝试展开循环。