Эффективная альтернатива Outer для разреженных массивов в Mathematica?

Я предложу решение, которое довольно сложное, но позволяет использовать только вдвое больше памяти во время вычислений, чем необходимо для хранения конечного результата в виде 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 будет зависеть от рассматриваемых разреженных массивов.

Если lst1 и lst2 ваши списки,

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

выполняет свою работу и может быть более эффективным с точки зрения использования памяти. С другой стороны, может быть, и нет. Насер прав, явный пример был бы полезен.

РЕДАКТИРОВАТЬ: Используя случайно сгенерированные массивы Нассера, и для len=200, MaxMemoryUsed[] указывает, что для этой формы требуется 170 МБ, а форма Outer в вопросе занимает 435 МБ.

Используя данные вашего примера 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, чтобы оптимизировать производительность.