Вы, безусловно, можете сделать намного лучше. Вот код, основанный на низкоуровневом 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] := Flatten@spart[s, 4][[2, 2]];
getSparseData[s_SparseArray] := spart[s, 4][[3]];
getDefaultElement[s_SparseArray] := spart[s, 3];
makeSparseArray[dims : {_, _}, jc : {__Integer}, ir : {__Integer}, data_List, defElem_: 0] :=
SparseArray @@ {Automatic, dims, defElem, {1, {jc, List /@ ir}, data}};
Clear[formSparseDivisible];
formSparseDivisible[a_, aa_, iaa_, chunkSize_: 100] :=
Module[{getDataChunkCode, i, start, ic, jr, sparseData, dims, dataChunk, res},
getDataChunkCode :=
If[# === {}, {}, SparseArray[1 - Unitize@(Mod[#, aa] & /@ #)]] &@
If[i*chunkSize >= Length[a],
{},
Take[a, {i*chunkSize + 1, Min[(i + 1)*chunkSize, Length[a]]}]];
i = 0;
start = getDataChunkCode;
i++;
ic = getIC[start];
jr = getJR[start];
sparseData = getSparseData[start];
dims = Dimensions[start];
While[True,
dataChunk = getDataChunkCode;
i++;
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]];
];
res = Transpose[makeSparseArray[dims, ic, jr, sparseData]];
res[[All, 1]] = N@iaa;
res]
Теперь вот время:
In[249]:=
n = 1500;
iaa = aa = Range[2 n];
a = Range[n];
AbsoluteTiming[res = formSparseDivisible[a, aa, iaa, 100];]
Out[252]= {0.2656250, Null}
In[253]:= AbsoluteTiming[
res1 = SparseArray[{{i_, 1} :>
iaa[[i]], {i_, j_} /; Divisible[a[[j]], aa[[i]]] -> 1.}, {2*n, n}];]
Out[253]= {29.1562500, Null}
Итак, у нас есть 100-кратное ускорение для этого размера массива. И, конечно же, результаты одинаковы:
In[254]:= Normal@res1 == Normal@res
Out[254]= True
Основная идея решения состоит в том, чтобы векторизовать проблему (Mod
) и постепенно создавать результирующий разреженный массив порциями, используя низкоуровневый API, описанный выше.
EDIT
Код предполагает, что списки имеют правильную длину - в частности, a
должен иметь длину n
, тогда как aa
и iaa
- 2n
. Таким образом, чтобы сравнить с другими ответами, тестовый код должен быть немного изменен (только для a
):
n = 500;
iaa = RandomReal[{0, 1}, 2 n];
a = Range[ n]; aa = RandomInteger[{1, 4 n}, 2 n];
In[300]:=
AbsoluteTiming[U=SparseArray[ReplacePart[Outer[Boole[Divisible[#1,#2]]&,
a[[1;;n]],aa],1->iaa]]\[Transpose]]
AbsoluteTiming[res = formSparseDivisible[a,aa,iaa,100]]
Out[300]= {0.8281250,SparseArray[<2838>,{1000,500}]}
Out[301]= {0.0156250,SparseArray[<2838>,{1000,500}]}
In[302]:= Normal@U==Normal@res
Out[302]= True
РЕДАКТИРОВАТЬ 2
Ваш желаемый размер матрицы устанавливается на моем не очень быстром ноутбуке (M8) примерно за 3 секунды, а также при достаточно приличном использовании памяти:
In[323]:=
n=5000;
iaa=RandomReal[{0,1},2 n];
a=Range[ n];aa=RandomInteger[{1,4 n},2 n];
AbsoluteTiming[res = formSparseDivisible[a,aa,iaa,200]]
Out[326]= {3.0781250,SparseArray[<36484>,{10000,5000}]}