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.