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) 那么这只是最后一步，需要出现非机器尺寸的浮子。

请参阅此答案了解解决方案适用于正浮点数和负浮点数。

让我们快速检查一下该解决方案是否适用于产生非机器大小答案的浮点数。与 Andrew 的（更快）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：遇到数字错误；继续进行未编译的评估。 并且需要一些时间才能给出结果...

---

一个好奇的是 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}

但对于长 PackedArray 来说，例如这里的测试示例，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”->无穷大}]。另一种选择是使用 lst=Developer`FromPackedArray[lst]。

Answer 2

首先，为了避免混淆，请看一个示例，其结果均可以表示为硬件机器精度数字，这些数字必须全部小于

In[1]:= $MaxMachineNumber

Out[1]= 1.79769*10^308

您的 Total 示例，该示例已具有此良好（且快速）的属性。这是使用机器编号的 Times 示例的变体：

In[2]:= lst = RandomReal[{0.99, 1.01}, 5000000];
Times @@ lst // Timing

Out[3]= {1.435, 1.38851*10^-38}

现在我们可以使用 编译使编译函数能够有效地执行此操作：

In[4]:= listproduct = 
 Compile[{{l, _Real, 1}}, 
  Module[{tot = 1.}, Do[tot *= x, {x, l}]; tot]]

Out[4]= CompiledFunction[{l},Module[{tot=1.},Do[tot*=x,{x,l}];tot],-CompiledCode-]

速度要快得多：

In[5]:= listproduct[lst] // Timing

Out[5]= {0.141, 1.38851*10^-38}

假设您有 C 编译器和 Mathematica 8，您还可以自动编译一路 C 代码。创建一个临时 DLL，并在运行时链接回 Mathematica。

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

Out[6]= CompiledFunction[{l},Module[{tot=1.},Do[tot*=x,{x,l}];tot],-CompiledCode-]

这提供的性能与内置 Mathematica 函数的性能没有太大区别：

In[7]:= compiledlistproduct[lst] // Timing

Out[7]= {0.015, 1.38851*10^-38}

请注意，如果您的产品确实会超出 $MaxMachineNumber （或 $MinMachineNumber)，那么你最好坚持使用Apply[Times, list]。如果您的结果可以达到那么大，同样的评论也适用于 Total：

In[11]:= lst = RandomReal[10^305, 5000000];
Plus @@ lst // Timing

Out[12]= {1.435, 2.499873364498981*10^311}

In[13]:= lst = RandomReal[10^305, 5000000];
Total[lst] // Timing

Out[14]= {1.576, 2.500061580905602*10^311}

Answer 3

西蒙的方法速度很快，但在负值时失败。结合他对我的其他问题的回答，这里是处理底片的快速解决方案。谢谢，西蒙。

功能

f = (-1)^(-1 /. Rule @@@ Tally@Sign@# /. -1 -> 0) * Exp@Total@Log@Abs@# &;

---

测试

lst = RandomReal[{-50, 50}, 5000000];

Times @@ lst // Timing
f@lst // Timing

lst = Developer`FromPackedArray[lst];

Times @@ lst // Timing
f@lst // Timing

---

{0.844, -4.42943661963*10^6323240}

{0.062, -4.4294366*10^6323240}

{0.64, -4.42943661963*10^6323240}

{0.203, -4.4294366*10^6323240}