PackedArrays有快速的产品操作吗？

Question

在Mathematica中,包含所有机器大小整数或浮点数的向量(或矩形数组)可以存储在打包数组中.这些对象占用的内存较少,而且某些操作的速度要快得多.

RandomReal尽可能生成一个打包数组.可以使用该Developer函数解压缩打包数组FromPackedArray

考虑这些时间

lst = RandomReal[1, 5000000];

Total[lst] // Timing
Plus @@ lst // Timing

lst = Developer`FromPackedArray[lst];

Total[lst] // Timing
Plus @@ lst // Timing

Out[1]= {0.016, 2.50056*10^6}

Out[2]= {0.859, 2.50056*10^6}

Out[3]= {0.625, 2.50056*10^6}

Out[4]= {0.64, 2.50056*10^6}

因此,在打包阵列的情况下,对于非打包阵列来说,Total快许多倍,Plus @@但大约相同.请注意,Plus @@在打包的数组上实际上要慢一点.

现在考虑

lst = RandomReal[100, 5000000];
Times @@ lst // Timing

lst = Developer`FromPackedArray[lst];
Times @@ lst // Timing

Out[1]= {0.875, 5.8324791357*10^7828854}

Out[1]= {0.625, 5.8324791357*10^7828854}

最后,我的问题是:Mathematica中有一个快速方法用于打包数组的列表产品,类似于Total？

我怀疑这可能是不可能的,因为数值误差与乘法复合的方式.此外,该函数需要能够返回非机器浮动才有用.

Answer 1

我也想知道是否有乘法相当于Total.

一个非常糟糕的解决方案是

In[1]:= lst=RandomReal[2,5000000];
        Times@@lst//Timing
        Exp[Total[Log[lst]]]//Timing
Out[2]= {2.54,4.370467929041*10^-666614}
Out[3]= {0.47,4.370467940*10^-666614}

只要数字是正数且不是太大或太小,那么舍入误差也不算太差.关于在评估期间可能发生的事情的猜测是:(1)如果数字是正浮点数,则Log操作可以快速应用于打包数组.(2)然后可以使用Total的打包数组方法快速添加数字.(3)然后,这只是需要出现非机器尺寸的浮子的最后一步.

请参阅此SO答案,了解适用于正浮动和负浮动的解决方案.

让我们快速检查这个解决方案是否适用于产生非机器大小答案的浮点数.与安德鲁相比(更快)compiledListProduct:

In[10]:= compiledListProduct = 
          Compile[{{l, _Real, 1}}, 
           Module[{tot = 1.}, Do[tot *= x, {x, l}]; tot], 
           CompilationTarget -> "C"]

In[11]:= lst=RandomReal[{0.05,.10},15000000];
         Times@@lst//Timing
         Exp[Total[Log[lst]]]//Timing
         compiledListProduct[lst]//Timing
Out[12]= {7.49,2.49105025389*10^-16998863}
Out[13]= {0.5,2.4910349*10^-16998863}
Out[14]= {0.07,0.}

如果你选择更大的(>1)实数,那么compiledListProduct会产生警告, CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation.并且需要一些时间来给出结果......

一个古玩是,无论Sum和Product可以采取任意列表.Sum工作良好

In[4]:= lst=RandomReal[2,5000000];
         Sum[i,{i,lst}]//Timing
         Total[lst]//Timing
Out[5]= {0.58,5.00039*10^6}
Out[6]= {0.02,5.00039*10^6}

但是对于long PackedArrays,例如此处的测试示例,Product由于自动编译的代码(在8.0版中)没有正确捕获下溢/溢出,因此失败:

In[7]:= lst=RandomReal[2,5000000];
         Product[i,{i,lst}]//Timing
         Times@@lst//Timing
Out[8]= {0.,Compile`AutoVar12!}
Out[9]= {2.52,1.781498881673*10^-666005}

有用的WRI技术支持提供的工作是使用关闭产品编译SetSystemOptions["CompileOptions" -> {"ProductCompileLength" -> Infinity}].另一种选择是使用lst=Developer`FromPackedArray[lst].