使用Compile有效处理小向量列表

Question

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

让我Compile通过总结2D向量列表的简单示例来说明如何使用这些问题:

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.

你可以在我的asnwer中看到相同"优化"的实际应用,但是在相关的情况下可以给出许多其他例子.

现在在这个简单的例子中,使用的代码fun2不会比使用代码更复杂或更复杂fun1,但在一般情况下它可能很好.

我怎样才能知道Compile是它的一个参数是不是一个武断的n*m矩阵,而是一种特殊n*2或n*3一个,所以它可以做这些自动优化,而不是使用一个通用的向量加法函数中添加微小的长度,长度为二或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版本,其中时序差异更加明显.

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看出这确实是fun1vs 的情况fun3.在某种程度上,您发现(对于此问题)一种有效的高级方式来手动展开外部循环(坐标索引上的循环).您建议的另一种方法是要求编译器根据向量维度的额外信息自动展开外部循环.这听起来似乎是一种合理的优化,但尚未实现Mathematica虚拟机的实现.

还要注意,对于较大长度的向量(比如20),它们之间的差异fun1和fun2消失,因为与大规模分配的成本相比,张量复制中的内存分配/解除分配的成本变得无关紧要(当您分配向量时,仍然可以更有效地实现向量 - 也许是因为你可以使用像memcpy这种情况那样的东西.

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

Answer 2

寻找一个完全符合您要求的功能始终是一个很好的选择.

In[50]:= fun3=Compile[{{vec,_Real,2}},Total[vec]]

Out[50]= CompiledFunction[{vec},Total[vec],-CompiledCode-]

In[51]:= Do[fun3[data],{10}]//Timing

Out[51]= {0.121982,Null}

In[52]:= fun3[data]===fun1[data]

Out[52]= True

另一个选项,效率较低(*由于转置*)是使用Listable

fun4 = Compile[{{vec, _Real, 1}}, Total[vec], 
  RuntimeAttributes -> {Listable}]

In[63]:= Do[fun4[Transpose[data]],{10}]//Timing

Out[63]= {0.235964,Null}

In[64]:= Do[Transpose[data],{10}]//Timing

Out[64]= {0.133979,Null}

In[65]:= fun4[Transpose[data]]===fun1[data]

Out[65]= True