我将提出一种相当复杂的解决方案,但只允许使用大约两倍的内存during将最终结果存储为所需的计算SparseArray
。为此付出的代价将是执行速度慢得多。
The code
稀疏数组构造/解构API
这是代码。首先,稍微修改一下(以解决高维稀疏数组)稀疏数组构造 - 解构 API,取自这个答案 https://stackoverflow.com/questions/7525782/import-big-files-arrays-with-mathematica/7527064#7527064:
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 https://stackoverflow.com/questions/7525782/import-big-files-arrays-with-mathematica/7527064#7527064答案,封装成一个函数。它接受用于生成下一个数据块的代码段,包装在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
对于这个特定的例子,建议的方法比直接使用的内存效率高 5 倍Outer
,但慢了大约 15 倍。我不得不调整chunksize
参数(默认为100
,但是对于上面我使用的2000
,以获得最佳速度/内存使用组合)。我的方法仅使用了存储最终结果所需内存两倍的峰值内存。与相比,节省内存的程度Outer
- 基于方法将取决于所讨论的稀疏数组。