5
$\begingroup$

I make a SparseArray and then modify it by adding rows together. This seems to take a long time. Can it be speeded up?

na = 1000;
nb = 60000;
nc = 50;
c = 0.8;
data = Flatten[RandomReal[{-1, 1}, {na, na}]];
loc = RandomInteger[{1, 5386/2}, {na^2, 2}];
mat = SparseArray[loc -> data, {nb, nb}, 0.0];
ri = RandomInteger[{1, nb/2}, nc];

Now to do the adding of rows. This is what takes the time.

Timing[Do[mat[[n]] = mat[[n]] + c mat[[nb/2 + n]],
   {n, ri}];]

This takes 1.09 seconds on my machine Version; 11.3. Could it be done faster? Is there a workaround?

Thanks

$\endgroup$
8
$\begingroup$

Put as many operations as possible into one call:

mat[[ri]] += c mat[[nb/2 + ri]]; // AbsoluteTiming // First

0.010061

Edit

I also thought about a more general situation where rows enumerated by a list readidx are read from an input matrix, multiplied by a factor, and then added to rows enumerated by a writeidx. If readidx or writeidx contains any duplicates, one can create a lower triangular matrix to achieve the goal much quicker than with the Do loop:

SeedRandom[666];
nnz = 100000;
nb = 60000;
nc = 5000;

data = Flatten[RandomReal[{-1, 1}, nnz]];
loc = RandomInteger[{1, nb}, {nnz, 2}];
A = mat0 = SparseArray[loc -> data, {nb, nb}, 0.0];
readidx = RandomInteger[{1, nb}, nc];
writeidx = RandomInteger[{1, nb}, nc];
factors = RandomReal[{-1, 1}, nc];
DuplicateFreeQ[readidx]
DuplicateFreeQ[writeidx]

False

False

B1 = A;
Do[B1[[writeidx[[i]]]] += factors[[i]] A[[readidx[[i]]]], {i, 1, Length[readidx]}]; // AbsoluteTiming // First

B2 = With[{L = With[{spopt = SystemOptions["SparseArrayOptions"]},
        Internal`WithLocalSettings[
         SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}],
         SparseArray[
          Transpose[{writeidx, readidx}] -> factors, 
          {Length[A], Length[A]}
         ],
         SetSystemOptions[spopt]]
        ]},
     A + L.A
     ]; // AbsoluteTiming // First
B1 == B2

18.1406

0.008159

True

Note that we have to SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}] so that the matrix L gets assempled correctly, even if Transpose[{writeidx, readidx}] contains duplicates.

$\endgroup$
  • $\begingroup$ How do I get the times c into it? $\endgroup$ – Hugh Apr 25 '18 at 21:55
  • $\begingroup$ This is wonderful. I have implemented in my code. Now a calculation that took 30 seconds takes 0.03 seconds. A factor of 1000. Your further problem is also enlightening. This is the second time you have helped me in 24 hours. Many thanks. $\endgroup$ – Hugh Apr 26 '18 at 9:45
  • $\begingroup$ Glad to hear that! You're always welcome! $\endgroup$ – Henrik Schumacher Apr 26 '18 at 11:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.