wolfram-mathematica - 如何在Mathematica中有效地初始化此稀疏数组？

Question

我正在尝试解决Mathematica中的一个很大的线性编程问题，但是由于某种原因，瓶颈正在建立线性约束数组。

我用于初始化矩阵的代码如下所示:

AbsoluteTiming[S = SparseArray[{{i_, 1} -> iaa[[i]],
    {i_, j_} /; Divisible[a[[j]], aa[[i]]] -> 1.}, {2*n, n}]]

这里n是4455，iaa是实数列表，而a，aa是整数列表。我在这一行得到的输出是

{2652.014773,SparseArray[<111742>,{8910,4455}]}

换句话说，即使仅具有111,742个非零条目，也需要45分钟来初始化此矩阵。为了进行比较，实际求解线性程序仅需17秒。是什么赋予了？

编辑:另外，谁能解释为什么它占用这么多的内存，因为它正在运行？因为在用户时间内，此计算只需不到十分钟的时间……大部分计算时间都花在了通过内存分页上。

Mathematica是否出于某种原因在构建该矩阵时将其存储为非稀疏矩阵？因为那真的很愚蠢。

Answer 1

您肯定可以做得更好。这是基于发布于here的低级稀疏数组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以块的形式逐步构建所得的稀疏数组。

编辑

该代码假定列表的长度正确-特别地，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}]}