Я предложу решение, которое является довольно сложным, но позволяет использовать только примерно вдвое больше памяти во время вычислений, сколько необходимо для сохранения окончательного результата как SparseArray. Цена за это будет гораздо медленнее.
код
API для создания / деконструкции разреженных массивов
Вот код. Во-первых, немного измененная (для решения многомерных разреженных массивов) конструкция разреженного массива - API деконструкции, взятый из этого ответа :
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) и склеивая их вместе. Это в основном код, используемый в этом ответе, упакованный в функцию. Он принимает фрагмент кода, использованный для создания следующего фрагмента данных, обернутый в 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 -ов за один раз, чтобы ускорить код. Сколько пар мы обрабатываем одновременно, определяется необязательным параметром 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, будет зависеть от рассматриваемых разреженных массивов.