# How to accelerate replacing of matrix rows in a loop as quickly as possible?

I am working with a very large sparse matrix (for example) given in what follows:

m = 50; n = 40; o = 30; size = m*n*o;
B = SparseArray[{
{i_, i_} -> RandomReal[], {size, size - 1} ->
2., {i_, j_} /; Abs[i - j] == 5 ->
1., {i_, j_} /; Abs[3 i - j] == 2 -> 2.
}, {size, size}, 0.]


This is a matrix of the size $60000\times 60000$ as a simple instance (in practice I have several of such matrices of higher sizes).

I wish to replace many rows of the matrix $B$ with special rows (here with zero rows) coming from boundary conditions. My list of rows' numbers are given by

index1 = Flatten[Table[{i, j, k}, {i, 1, m}, {j, 1, n}, {k, 1, o}], 2];
Table[ind = index1[[l]];
If[ind[] == 1 || ind[] == n, var[l] = 0, var[l] = l], {l,
Length@index1}];
bounindex = Table[var[l], {l, Length@index1}];


And my zero row is produced by

vector1 = SparseArray@ConstantArray[0, size];


Now for doing the replacement only for the 0 numbers obtained in bounindex, I write the following replacing rule in a loop as simply as possible:

Table[Which[bounindex[[i]] == 0, B[[i]] = vector1], {i,
size}]; // AbsoluteTiming


This works OK but it takes a tremendous time (in my real problem I have larger sizes)! So I am wondering if there is a speed-up technique to replace the special rows of the matrix $B$ with appropriate rows (that I want).

I would be thankful if some hints for accelerating such a replacement be given.

(updated)

For later comparisons, I will store the value of B:

B0 = B;


Original post

Set all the rows to zero at once. First, use Pick to find the rows to be zeroed:

delete = Pick[Range, bounindex, 0];
len = Length[delete]


2670

Here is the number of nonzero values before zeroing out rows:

Length @ B["NonzeroValues"]


219989

Next, use Part to set all of the rows to 0:

B[[delete]] = SparseArray[{}, {len, 60000}]; //AbsoluteTiming


{0.00192, Null}

Here is the number of nonzero values after zeroing out rows:

Length @ B["NonzeroValues"]


208689

Update

The OP asked about replacements with a more complicated RHS (in his example the RHS could be constructed from an appropriate Dot product of a sparse matrix with B). Here is an example with a random, but fairly sparse RHS:

sub = SparseArray @ RandomChoice[{.9999,.0001} -> {0, 1}, {2670, 60000}];


Using the same procedure as before:

B = B0;
Length @ B["NonzeroValues"]
B[[delete]] = sub; //AbsoluteTiming


219989

{0.009667, Null}

Check:

Length @ sub["NonzeroValues"]
Length @ B["NonzeroValues"]

%-%%


16024

224713

208689

The timing above depends heavily on the sparsity of the substitute matrix.

One final note. @george2079 tried multiplying 2 sparse matrices as an alternative. A similar idea is to multiply a sparse array with a spare vector. This is also quite fast, but not as flexible (one can only zero out rows, not replace rows with some other desired row):

B = B0;
Length @ B["NonzeroValues"]

B = B SparseArray @ Clip @ bounindex; //AbsoluteTiming
Length @ B["NonzeroValues"]


219989

{0.006694, Null}

208689

• tricky use of sparsearray! – Ulrich Neumann Dec 28 '17 at 20:58
• This is cool in my viewpoint. Now, I have an inquiry to see if it is possible to extend your response. Is there a way to generalize this way by a way at which several rows are replaced by a combination of some other rows? For instance, in the following piece of code: Table[B[[i]] = (2 B[[i - 2 (n*o)]] + 4 B[[i - (n*o)]]), {i, size, size - (n*o) + 1, -1}]; // AbsoluteTiming – Fazlollah Dec 28 '17 at 21:00
• @Fazlollah I'm afk for a bit, but the rhs looks like a dot product of a matrix with 2s and 4s in the right place with B. If this matrix is Sparse, then the Dot product will be sparse, and then the same kind of Part assignment should work. – Carl Woll Dec 28 '17 at 21:09
• This response is great and fruitful for anyone solves high dimensional PDEs and wishes to impose the boundaries directly into the coefficient matrix. Thanks a lot. – Fazlollah Dec 30 '17 at 16:18

Perhaps ReplacePart[] is the answer you are looking for:

index0 = Position[bounindex, 0];
Bnew=ReplacePart[B, index0 -> vector1];


It is as fast as the Table-version...

• Does this do the replacement properly? Please re-check it. – Fazlollah Dec 28 '17 at 20:53
• ReplacePart command doesn't change B, Bnew(see my edit/answer is the 'replaced' matrix! – Ulrich Neumann Dec 29 '17 at 11:03

Edit no this is not good. It seems the matrix multiplication expands the sparse arrays ( evidenced by memory usage).

--

I dont know if this is good performance wise or not, but here is a simple idea:

construct a 1-base sparse array with a zero row and multiply:

B = B SparseArray[ {3, i_} -> 0, Dimensions[base] , 1]