如何在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}]}

换句话说,初始化此矩阵需要45分钟,即使它只有111,742个非零项.为了比较,实际上解决线性程序只需要17秒.是什么赋予了？

编辑:此外,任何人都可以解释为什么这会占用如此多的内存,因为它运行？因为在用户时间,这个计算花费的时间不到十分钟......大部分计算时间都花费在内存中进行分页.

Mathematica出于某种原因将这个矩阵存储为非稀疏矩阵,而它正在构建它吗？因为那真的很蠢.

Answer 1

你肯定可以做得更好.这是一个基于此处发布的低级稀疏数组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}]}