在 Mathematica 中的稀疏数组上有效替代 Outer？

Posted 2023-02-24

【中文标题】在 Mathematica 中的稀疏数组上有效替代 Outer？

【英文标题】：Efficient alternative to Outer on sparse arrays in Mathematica?

【发布时间】：2012-01-25 15:18:11

【问题描述】：

假设我有两个非常大的列表 a1, a2, ... 和 b1, b2, ...，其中所有 ai 和 bj 都是大型稀疏数组。为了内存效率，我将每个列表存储为一个全面的稀疏数组。

现在我想在所有可能的 ai 和 bj 对上计算一些函数 f，其中每个结果 f[ai, bj] 又是一个稀疏数组。顺便说一下，所有这些稀疏数组都具有相同的维度。

虽然

Flatten[Outer[f, a1, a2, ..., b1, b2, ..., 1], 1]

返回所需的结果（原则上）它似乎消耗了过多的内存。尤其是因为返回值是一个稀疏数组列表，而在我感兴趣的案例中，一个全面的稀疏数组效率更高。

对于上述使用Outer，有没有一种有效的替代方法？

更具体的例子：

SparseArray[1, 1, 1, 1 -> 1, 2, 2, 2, 2 -> 1],
 SparseArray[1, 1, 1, 2 -> 1, 2, 2, 2, 1 -> 1],
 SparseArray[1, 1, 2, 1 -> 1, 2, 2, 1, 2 -> 1],
 SparseArray[1, 1, 2, 2 -> -1, 2, 2, 1, 1 -> 1],
 SparseArray[1, 2, 1, 1 -> 1, 2, 1, 2, 2 -> 1],
 SparseArray[1, 2, 1, 2 -> 1, 2, 1, 2, 1 -> 1],
 SparseArray[1, 2, 2, 1 -> -1, 2, 1, 1, 2 -> 1],
 SparseArray[1, 2, 2, 2 -> 1, 2, 1, 1, 1 -> 1];
ByteCount[%]

list = SparseArray[%%]
ByteCount[%]

Flatten[Outer[Dot, list, list, 1], 1];
ByteCount[%]
list1x2 = SparseArray[%%]
ByteCount[%]

Flatten[Outer[Dot, list1x2, list, 1], 1];
ByteCount[%]
list1x3 = SparseArray[%%]
ByteCount[%]

等等。不仅Outer（稀疏数组列表）的原始中间结果效率极低，Outer 在计算过程中似乎也消耗了太多内存。

【问题讨论】：

您可以制作一个独立的小示例来说明您的问题。因为任何愿意回答你的人都必须这样做。见sscce.org。他是你的例子吗？长度 = 1000; f[x_] := SparseArray[Table[x^2, len]]; list1 = SparseArray[Table[RandomReal[], len]]; list2 = SparseArray[Table[RandomReal[], len]];计时[外部[f, list1, list2]][[1]]

我不确定我是否理解你所说的“综合稀疏数组”是什么意思。你能详细说明一下吗？

@Nasser 你的例子中的函数 f 不应该是一个有两个参数的函数吗？

@SjoerdC.deVries，我相信您可能是正确的，但这就是我要问的重点。我不确定这个问题。如果用户发布了一个小的具体示例，那么回答的人不必弥补一些东西，并且可能会在这个过程中错过一些东西。

这里只是一些想法：Outer 似乎对Times 有特殊的支持：它保持了稀疏结构并且效率很高。一般来说，如果可以保证f 函数没有副作用，Outer[f, ...] 对于稀疏数组的效率会更高（因此f[a,b] always 返回相同的结果相同的a 和b）。那么稀疏数组的“未指定元素”（通常为0）只能传递给f一次，结果可以作为新稀疏数组中的新“未指定元素”。虽然无法自动检测到某些 ...

【参考方案1】：

如果 lst1 和 lst2 是您的列表，

Reap[
   Do[Sow[f[#1[[i]], #2[[j]]]],
       i, 1, Length@#1,
       j, 1, Length@#2
       ] &[lst1, lst2];
   ] // Last // Last

完成这项工作并且可能更节省内存。另一方面，也许不是。纳赛尔是对的，一个明确的例子会很有用。

编辑：使用 Nasser 随机生成的数组，对于 len=200，MaxMemoryUsed[] 表示此表单需要 170MB，而问题中的 Outer 表单需要 435MB。

【讨论】：

怎么样：Reap[Do[Sow[Dot[i, j]], i, r, j, r];] // Last // Last with r=SparseArray [],.. 或 Flatten[Table[Dot[i, j], i, r, j, r], 1]

【参考方案2】：

使用您的示例list 数据，我相信您会发现Append 到SparseArray 的能力非常有帮助。

acc = SparseArray[, 1, 2, 2, 2, 2, 2, 2]

Do[AppendTo[acc, i.j], i, list, j, list]

Rest[acc]

我需要Rest 删除结果中的第一个零填充张量。种子SparseArray 的第二个参数必须是带有前缀1 的每个元素的尺寸。您可能需要为种子 SparseArray 显式指定背景以优化性能。

【讨论】：

你可以使用 Internal`Bag 的东西，这取决于你有多少这样的数组。

@ruebenko，这与直接附加到 SparseArray 相比如何？

实际上没有；带有 SparseArray 的 Suffing Bag 使 SA 成为法线。我不知道。这至少解释了为什么 Internal`Bag 只是内部的。对于任意表达式，它没有完全实现。

使用acc = SparseArray[, 1, 2, 2, 2, 2, 2, 2]; Do[AppendTo[acc, Dot[i, j]], i, list, j, list]; list1x2 = Rest[acc]; Clear[acc] 等。实际上我的性能要差得多。一方面，稀疏数组list1x2、list1x3 等比通过上述Outer 构造的对应物占用更多的内存，尽管测试表明它们是相等的。此外，在计算 list1x6 时，MathKernel 内存不足，而 Outer 方法不会发生这种情况。知道为什么会这样吗？

@Mr.Wizard 试过这个，它在内存上的效率似乎高得离谱……但代价是运行时间非常慢。好吧，我想总会有一种权衡。感谢您的提示！

【参考方案3】：

我将提出一个相当复杂的解决方案，但只允许在计算期间使用大约两倍的内存 将最终结果存储为SparseArray。为此付出的代价将是执行速度要慢得多。

代码

稀疏数组构造/解构API

这里是代码。一、稍作修改（用于处理高维稀疏数组）稀疏数组构造——解构API，取自this answer：

ClearAll[spart, getIC, getJR, getSparseData, getDefaultElement, 
  makeSparseArray];
HoldPattern[spart[SparseArray[s___], p_]] := s[[p]];
getIC[s_SparseArray] := spart[s, 4][[2, 1]];
getJR[s_SparseArray] := spart[s, 4][[2, 2]];
getSparseData[s_SparseArray] := spart[s, 4][[3]];
getDefaultElement[s_SparseArray] := spart[s, 3];
makeSparseArray[dims_List, jc_List, ir_List, data_List, defElem_: 0] :=
   SparseArray @@ Automatic, dims, defElem, 1, jc, ir, data;

迭代器

以下函数产生迭代器。迭代器是封装迭代过程的好方法。

ClearAll[makeTwoListIterator];
makeTwoListIterator[fname_Symbol, a_List, b_List] :=
  With[indices = Flatten[Outer[List, a, b, 1], 1],
   With[len  = Length[indices],
    Module[i = 0,
      ClearAll[fname];
      fname[] := With[ind = ++i, indices[[ind]] /; ind <= len];
      fname[] := Null;
      fname[n_] := 
        With[ind = i + 1, i += n; 
           indices[[ind ;; Min[len, ind + n - 1]]] /; ind <= len];
      fname[n_] := Null;
 ]]];

请注意，我本可以更有效地实现上述功能，而不是在其中使用Outer，但就我们的目的而言，这不是主要问题。

这是一个更专业的版本，它为二维索引对生成交互器。

ClearAll[make2DIndexInterator];
make2DIndexInterator[fname_Symbol, i : iStart_, iEnd_, j : jStart_, jEnd_] :=
   makeTwoListIterator[fname, Range @@ i, Range @@ j];
 make2DIndexInterator[fname_Symbol, ilen_Integer, jlen_Integer] :=
   make2DIndexInterator[fname, 1, ilen, 1, jlen];

这是如何工作的：

In[14]:= 
makeTwoListIterator[next,a,b,c,d,e];
next[]
next[]
next[]

Out[15]= a,d
Out[16]= a,e
Out[17]= b,d

我们也可以使用它来获取批处理结果：

In[18]:= 
makeTwoListIterator[next,a,b,c,d,e];
next[2]
next[2]

Out[19]= a,d,a,e
Out[20]= b,d,b,e

，我们将使用第二种形式。

SparseArray - 构建函数

此函数将通过获取数据块（也是SparseArray 形式）并将它们粘合在一起来迭代地构建SparseArray 对象。它基本上是this答案中使用的代码，打包成一个函数。它接受用于生成下一个数据块的代码片段，包裹在 Hold 中（我也可以将其设为 HoldAll）

Clear[accumulateSparseArray];
accumulateSparseArray[Hold[getDataChunkCode_]] :=
  Module[start, ic, jr, sparseData, dims, dataChunk,
   start = getDataChunkCode;
   ic = getIC[start];
   jr = getJR[start];
   sparseData = getSparseData[start];
   dims = Dimensions[start];
   While[True, dataChunk = getDataChunkCode;
     If[dataChunk === , Break[]];
     ic = Join[ic, Rest@getIC[dataChunk] + Last@ic];
     jr = Join[jr, getJR[dataChunk]];
     sparseData = Join[sparseData, getSparseData[dataChunk]];
     dims[[1]] += First[Dimensions[dataChunk]];
   ];
   makeSparseArray[dims, ic, jr, sparseData]];

把它们放在一起

这个功能是主要的，放在一起：

ClearAll[sparseArrayOuter];
sparseArrayOuter[f_, a_SparseArray, b_SparseArray, chunkSize_: 100] := 
  Module[next, wrapperF, getDataChunkCode,
    make2DIndexInterator[next, Length@a, Length@b];
    wrapperF[x_List, y_List] := SparseArray[f @@@ Transpose[x, y]];
    getDataChunkCode :=
      With[inds = next[chunkSize],
         If[inds === Null, Return[]];
         wrapperF[a[[#]] & /@ inds[[All, 1]], b[[#]] & /@ inds[[All, -1]]]
      ];
    accumulateSparseArray[Hold[getDataChunkCode]]
  ];

在这里，我们首先生成迭代器，它将按需为我们提供索引对列表的部分，用于提取元素（也是SparseArrays）。请注意，我们通常会一次从两个大输入SparseArray-s 中提取一对以上的元素，以加快代码速度。我们一次处理多少对由可选的chunkSize 参数控制，默认为100。然后我们构造代码来处理这些元素并将结果放回SparseArray，我们使用辅助函数wrapperF。迭代器的使用不是绝对必要的（可以使用Reap-Sow，与其他答案一样），但允许我将迭代逻辑与稀疏数组的通用累积逻辑分离。

基准

首先我们准备大型稀疏数组并测试我们的功能：

In[49]:= 
arr = SparseArray[1,1,1,1->1,2,2,2,2->1],SparseArray[1,1,1,2->1,2,2,2,1->1],
SparseArray[1,1,2,1->1,2,2,1,2->1],SparseArray[1,1,2,2->-1,2,2,1,1->1],
SparseArray[1,2,1,1->1,2,1,2,2->1],SparseArray[1,2,1,2->1,2,1,2,1->1];

In[50]:= list=SparseArray[arr]
Out[50]= SparseArray[<12>,6,2,2,2,2]

In[51]:= larger = sparseArrayOuter[Dot,list,list]
Out[51]= SparseArray[<72>,36,2,2,2,2,2,2]

In[52]:= (large= sparseArrayOuter[Dot,larger,larger])//Timing
Out[52]= 0.047,SparseArray[<2592>,1296,2,2,2,2,2,2,2,2,2,2]

In[53]:= SparseArray[Flatten[Outer[Dot,larger,larger,1],1]]==large
Out[53]= True

In[54]:= MaxMemoryUsed[]
Out[54]= 21347336

现在我们进行功率测试

In[55]:= (huge= sparseArrayOuter[Dot,large,large,2000])//Timing
Out[55]= 114.344,SparseArray[<3359232>,1679616,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2]

In[56]:= MaxMemoryUsed[]
Out[56]= 536941120

In[57]:= ByteCount[huge]
Out[57]= 262021120

In[58]:= (huge1 = Flatten[Outer[Dot,large,large,1],1]);//Timing
Out[58]= 8.687,Null

In[59]:= MaxMemoryUsed[]
Out[59]= 2527281392

对于这个特定的示例，建议的方法比直接使用Outer 的内存效率高 5 倍，但速度慢了大约 15 倍。我不得不调整chunksize 参数（默认为100，但对于上面我使用2000，以获得最佳速度/内存使用组合）。我的方法仅用作存储最终结果所需内存的两倍的峰值。与基于Outer 的方法相比，内存节省的程度取决于所讨论的稀疏数组。

【讨论】：

在我有时间测试并理解这一点之前，我不会投票，但我希望你能回复。