7
$\begingroup$

I have an array of sparse arrays (check update 3 below for a concrete example that has the same qualitative features as what I use in my real code) , say

sparse = Table[
  KroneckerProduct[RandomReal[{-10, 10}, {50, 50}], 
   IdentityMatrix[50, SparseArray]], {ii, 1, 5}, {jj, 1, 5}];

and I want to do 2 things:

1) Set the non-zero elements in every nth row (except for the last n+1 rows) in each sparse array to 0.

I currently do this by

n=3;
AbsoluteTiming[
 sparse[[All, All, 1 ;; -n-1 ;; n]] = 
   ConstantArray[0., Dimensions[sparse[[All, All, 1 ;; -n-1 ;; n]]]];]

{0.38174, Null}

and

2) update some of the rows of some of the sparse arrays (not all of them) with the rows from other arrays (that I generate separately)

for example

randMatrix1 = 
  KroneckerProduct[RandomReal[{-10, 10}, {50, 50}], 
   IdentityMatrix[50, SparseArray]];
randMatrix2 = 
  KroneckerProduct[RandomReal[{-10, 10}, {50, 50}], 
   IdentityMatrix[50, SparseArray]];

AbsoluteTiming[
 sparse[[1, 1]][[1 ;; -n - 1 ;; n + 1]] = 
  randMatrix1[[1 ;; -n - 1 ;; n + 1]];
 sparse[[2, 3]][[n + 1 ;; -1 ;; n + 1]] = 
  randMatrix2[[n + 1 ;; -1 ;; n + 1]];]

{0.001794, Null}

Ok, the last one is quite quick, but I have to do it many times, as I have different combinations of rows that should be replaced with rows from different matrices. Also in my real code the sizes are larger and randMatrix1 and randMatrix2 can be either sparse or dense.

How can these operations by sped up, or is there a completely different way of doing this?

UPDATE: Added the information that randMatrix1 and randMatrix2 can be either sparse or dense.

UPDATE 2: There is also the special case when randMatrix1 and randMatrix2 are the IdentitiyMatrix. Is that case any easier for some reason?

UPDATE 3: Per request in some of the answers, I provide code that produces an array of matrices similar to what I work with in my real problem

{dm1, dm2, dm3, dm4, dm5} = 
  Table[SparseArray[
    Flatten[Table[
       RandomReal[{0, 1}, 250][[ii]]*
        Block[{vec}, vec = RandomReal[{0, 10}, {10}]; 
         vec - Min[Abs[vec]]], {ii, 1, 250}]]*
     RandomReal[{0, 1}, {2500, 2500}]], 5];

sm = Table[
   SparseArray[
    Block[{vec}, vec = ConstantArray[0., 2500]; 
      vec[[3]] = RandomReal[{0., 10.}]; vec]*
     RandomReal[{0, 1}, {2500, 2500}]], 20];

sA = {{dm1, sm[[1]], sm[[2]], sm[[3]], sm[[4]]}, {sm[[5]], dm2, 
    sm[[6]], sm[[7]], sm[[8]]}, {sm[[9]], sm[[10]], dm3, sm[[11]], 
    sm[[12]]}, {sm[[13]], sm[[14]], sm[[15]], dm4, 
    sm[[16]]}, {sm[[17]], sm[[18]], sm[[19]], sm[[20]], dm5}};

Where sA is basically what I am working with in my real problem (apart from dimensions, of course, but that should be large enough). So the diagonal terms are very dense, whereas everything else is very sparse.

I am not sure I can ArrayFlatten the whole thing and then do the changes I want, as the changes to each block matrix are different (in the sense that the same rows get updated, but with different values) and implementing them in one go for all the blocks will be very hard, or at least my knowledge of Mathematica doesn't let me see how this can be done.

UPDATE 4: Partial solution So I think the problem is that the diagonal sparse arrays in my list of sparse arrays are very dense. Hence if I actually normalise them and then do the replacements things speed up. Here is how I do it

dim = 49; n = 49;
replace = Normal /@ Diagonal[sA];
(sA[[#, #]] = replace[[#]]) & /@ Range[5];
sA[[All, All, 1 ;; -n - 1 ;; n + 1]] *= 0.;(*and other replacements*)
replace = SparseArray /@ Diagonal[sA];
(sA[[#, #]] = replace[[#]]) & /@ Range[5];

and then I can ArrayFlatten the whole thing and continue with my calculation. Not very elegant but it helps a bit.

Improve this question
$\endgroup$
4
  • $\begingroup$ If someone edits this I’d like to upvote it! Sometimes I accidentally downvote questions that I actually like (like this one!) no clue how it happens. Thanks for asking a good question! $\endgroup$
    – CA Trevillian
    Nov 17, 2019 at 22:50
  • $\begingroup$ @CATrevillian what do you mean by 'someone edits it'? $\endgroup$
    – ThunderBiggi
    Nov 18, 2019 at 23:50
  • $\begingroup$ Just that I can’t change my vote until there is a useful/accepted edit provided. $\endgroup$
    – CA Trevillian
    Nov 18, 2019 at 23:52
  • 1
    $\begingroup$ @CATrevillian I made an edit $\endgroup$
    – ThunderBiggi
    Nov 19, 2019 at 1:19

2 Answers 2

Reset to default
8
$\begingroup$

The following did the first trick 40 times faster (on my machine):

sparse[[All, All, 1 ;; -n - 1 ;; n]] *= 0.

One of the major problems here is that 

ConstantArray[0., Dimensions[b[[All, All, 1 ;; -n - 1 ;; n]]]]

is a dense array although it conveys basically no information.

Improve this answer
$\endgroup$
5
  • $\begingroup$ I can't test it now, as I am traveling. However, I forgot to mention that in 2) randMatrix1 and randMatrix2 can be either dense or sparse. I updated the question with that information $\endgroup$
    – ThunderBiggi
    Oct 22, 2019 at 16:32
  • $\begingroup$ I tested the code and in this way it is way faster but applying it to my real code yields about the same times. The reason is that the structure of my array of sparse arrays is very specific (they all come after a lot of operations, starting with very sparse arrays). The diagonal ones are very dense, but all the rest are almost empty. In the end I 'ArrayFlatten' the whole thing and it is sparse as a whole. Shall I update my question with more details or I can provide a '{5,5,2601,2601}' examplary array if I am allowed to upload files here. $\endgroup$
    – ThunderBiggi
    Oct 22, 2019 at 23:24
  • $\begingroup$ @ThunderBiggi Maybe you can provide the code that would produce a typical example of your array. Moreover this: In contrast to in-place operations on dense matrices, "in-place operations" on sparse arrays have typically to rebuild the whole matrix structure if the sparsity pattern changes. So it might also make sense to ArrayFlatten the whole array and to perform the according row operations on the resulting sparse array. This might reduce overhead time (for rebuilding the sparsity pattern(s)). $\endgroup$
    – Henrik Schumacher
    Oct 22, 2019 at 23:31
  • $\begingroup$ If the sparsity pattern does not change, then one should also try to modify only the vector of "NonzeroValues"; that would be super fast. But from what you said, the sparsity pattern does change, right? $\endgroup$
    – Henrik Schumacher
    Oct 22, 2019 at 23:32
  • $\begingroup$ I updated my post as per your request. I was inactive as I have been travelling across continents for work. $\endgroup$
    – ThunderBiggi
    Oct 25, 2019 at 0:43
2
$\begingroup$

An alternative to Henrik solution:

Generate data:

dim = 50;
sparse = Table[
   KroneckerProduct[RandomReal[{-10, 10}, {dim, dim}], 
    IdentityMatrix[dim, SparseArray]], {ii, 1, 5}, {jj, 1, 5}];

My method:

n = 3;
AbsoluteTiming[
 sp = SparseArray[
   Drop[ArrayRules[IdentityMatrix[dim^2, SparseArray]], {1, -n, 
     n}], {dim^2, dim^2}];
 sparse1 = sp.# & /@ # & /@ sparse;
 ]

{0.170991, Null}

Henrik's method:

sparse2 = sparse;
n = 3;
AbsoluteTiming[
 sparse2[[All, All, 1 ;; -n - 1 ;; n]] *= 0.;
 ]

{0.092583, Null}

sparse1==sparse2

True

Unfortunately, my method takes twice the time of Henrik's one (tested with different dim values), and is much more verbose, but I think is an interesting comparison!

Improve this answer
$\endgroup$
1
  • $\begingroup$ I tested this exact implementation on my real code (which doesn't work that well with Henrik's solution) and it is slower than both my implementation and Henrik's, but this is due to the specifics of my real problem (see comment under Henrik's post) $\endgroup$
    – ThunderBiggi
    Oct 22, 2019 at 23:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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