# Efficient way to compose a SparseArray from system of linear equations [closed]

I have several relations between some variables, e.g. id[1] - id[2] = 0, id[2] = 0. I need to compose a SparseArray that contains the coefficients. In that example,{1,-1},{0,1}. But I have several thousand terms and my method takes almost 1h. I use a Do loop and I imagine Mathematica has a more efficient way.

Any suggestion or help is extremely welcome.

Minimal example:

nmax = 10000;
list = Table[{RandomInteger[{1, 10}]*id[RandomInteger[{1, nmax}]] +
RandomInteger[{1, 10}]*id[RandomInteger[{1, nmax}]]}, {j,
nmax}];
Do[list[[j]] = list[[j]] /. id[n_] :> id[n, j], {j, Length[list]}]; (*<- too slow*)
SparseArray[
list /. Plus -> List /. a_. id[n__] :> ({n} -> a) // Flatten]

(*Time:2.5s*)

• Generation of the initial list already takes half of final computation time. I think you have to look at SymmetrizedArray[ ] if you want to impose symmetries in an efficient way.
– Acus
Commented Apr 12, 2023 at 12:08
• Thanks for the comment. What takes time in my machine is the Do loop. Generating the first list takes only 0.04s, but anyway that I already have in my real case. Commented Apr 12, 2023 at 12:21
• Check document of CoefficientArrays Commented Apr 12, 2023 at 12:25
• Thank you for this. CoefficientArrays[list, id[#] & /@ Range[1, nmax]][[2]] // Transpose does the same as my function in 0.023 s. Commented Apr 12, 2023 at 13:01
• But testing I discovered that the bottleneck was in reassigning elements of list: list[[j]] = list[[j]] // .... . Simply doing newlist = list and nlist[[j]] = list[[j]] // ... speed it up to 0.048s Commented Apr 12, 2023 at 13:03

You can add indices to your matrix instead of doing iteration and then a simple Replace can do all the work.

Assuming:

nmax = 10000;

list = Table[{RandomInteger[{1, 10}]*id[RandomInteger[{1, nmax}]] +
RandomInteger[{1, 10}]*id[RandomInteger[{1, nmax}]]},
{j, nmax}];

(* list: {{id[5] + 4 id[6]} , ... } *)


Solution:

Block[{temp},

temp =
ArrayPad[list, {0, {0, 1}}, Transpose@{Range[Length[list]]}];

(* temp: {{id[5] + 4 id[6], 1} , ... } *)

r2 = SparseArray@
Replace[temp, {
{m1_.*id[v1_] + m2_.*id[v2_], i_} :> Sequence[{v1, i} -> m1, {v2, i} -> m2],
{m_.*id[v_], i_} :> ({v, i} -> m)
}, {1}]

]


If we store your solution's result in r1, you'll see r1 == r2 and the calculation is ~60 times faster (3.49 vs 0.057 on my computer).

• Thank you so much. I had tried to reduce everything to a simple replace but didn't succeed. This is very clever. However the bottleneck was simply in list[[j]] = list[[j]] //... (see the comments below the question) Commented Apr 12, 2023 at 13:08

You may use "MapIndexed" for this task. This command takes every sublist of list and applies a function to it. This function receives as first argument the sublist and as second element its position.

We first define the function that transforms expressions like: num1 id[num2] into {row number, num2} -> num 1, needed to create the sparse array:

fun[sub_, pos_] = sub /. { n1_?NumericQ : 1  id[n2_] :> {pos[[1]], n2} -> n1, Plus -> List };


With this we apply SparseArray to the list to create the definition of the array:

sa = MapIndexed[fun, list] //Flatten;
SparseArray[sa]


• Very cool. Mathematica has so many useful functions like this. The bottleneck of my function was actually somewhere else, but thank you so much for your solution Commented Apr 12, 2023 at 13:10