How can I effectively initialize this sparse array in Mathematica?

Question

I am trying to solve a rather large linear programming problem in Mathematica, but for some reason the bottleneck is setting up an array of linear constraints.

My code to initialize the matrix is ​​as follows:

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

Here n is 4455, iaa is a list of reals, and aa are lists of integers. The output I get for this line is

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

In other words, it takes 45 minutes to initialize this matrix, although it has only 111,742 non-zero entries. For comparison, in fact, solving a linear program takes only 17 seconds. What gives?

Edit: Also, can someone explain why this takes up as much memory as it works? Since in user time this calculation takes less than ten minutes ... most of the calculation time is spent on paging through memory.

Is Mathematica for some reason storing this matrix as a non-sparse matrix when it is created? Because that would be really stupid.

Answer 1

You can do much better. Here is the code based on the low-level sparse array of APIs placed here , which I will reproduce to make the code self-sufficient:

 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] 

Now, here is the timing:

 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} 

So, we have 100x acceleration for this array size. And, of course, the results are the same:

 In[254]:= Normal@res1 == Normal@res Out[254]= True 

The main idea of ​​the solution is to vectorize the problem ( Mod ) and stage-by-stage construction of the resulting sparse array in pieces using the above lower-level API.

EDIT

The code assumes that the lists are of the correct length - in particular, a must be of length n , and aa and iaa be 2n . So, in order to compare with other answers, the test code should be slightly modified (only for 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 

EDIT 2

The desired matrix size takes about 3 seconds to complete on my not-so-fast laptop (M8) and with pretty decent memory usage:

 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}]} 

Answer 2

Your design calls Divisible 2*n^2 ~= 40 000 000 times, so it doesn't scale too well! Depending on the size of the integers, this will be from one third to half the time. To test this, you just need to run

 AbsoluteTiming[Table[Divisible[a[[j]], aa[[i]]], {i, 2*n}, {j, n}]] 

Combine this with a pattern matching construct, which should call an evaluator for each element of a sparse array, and you can understand why this might be a little slow.

---

Edit: Here is some test code that compares timings (on my machine) for various constructs:

 In[1]:= n = 500; iaa = RandomReal[{0, 1}, 2 n]; a = Range[2 n]; aa = RandomInteger[{1, 4 n}, 2 n]; In[4]:= AbsoluteTiming[S = SparseArray[{{i_, 1} :> iaa[[i]], {i_, j_} /; Divisible[a[[j]], aa[[i]]] -> 1.}, {2*n, n}]] Out[4]= {3.423123, SparseArray[<2499>, {1000, 500}]} In[5]:= AbsoluteTiming[T = SparseArray[ReplacePart[Table[ Boole[Divisible[a[[i]], aa[[j]]]], {i, 1, n}, {j, 1, 2 n}], 1 -> iaa]]\[Transpose]] Out[5]= {1.524575, SparseArray[<2499>, {1000, 500}]} In[6]:= AbsoluteTiming[U = SparseArray[ReplacePart[Outer[ Boole[Divisible[#1, #2]]&, a[[1 ;; n]], aa], 1 -> iaa]]\[Transpose]] Out[6]= {0.916801, SparseArray[<2499>, {1000, 500}]} In[7]:= S == T == U Out[7]= True 

Unless you have other information about integers in a and aa , 2*n^2 tests using Divisible will effectively limit code execution speed.