4
$\begingroup$

in writing some code for FEM i need an efficient way to replace parts of the diagonal of a sparse array. Certain positions in the diagonal of the array (given in a list) should be changed to "1".

An example would be:

mat = RandomInteger[3, {10, 10}];
replacelist = {{1, 1}, {3, 3}, {7, 7}, {10, 10}};
ReplacePart[mat, replacelist -> 999] // MatrixForm
(*
999 0   2   1   1   3   2   3   2   1
1   2   0   3   1   1   0   0   0   0
0   0   999 0   2   0   2   3   3   2
1   2   1   1   0   0   3   0   0   1
3   3   2   3   1   1   0   3   1   1
1   2   1   3   0   2   3   1   0   3
2   2   3   1   2   1   999 0   3   3
0   0   1   2   2   1   2   3   3   1
3   1   3   1   2   3   2   3   2   2
2   1   1   0   2   2   0   0   2   999
*)

This works, but it is way to slow, since I assume it is done elementwise. It has to be orders of magnitude faster. My skills in using pure functions is unfortunately limited, but I guess there is some clever way to do it with those. Also I have to construct a list first (replacelist) containing all positions twice as row and column for which I guess is a more elegant way possible.

Thank you very much in advance!

$\endgroup$
4
$\begingroup$

I gave a presentation on this subject a few years back. You can find efficient matrix assembly and modification codes in the talk Applications of Numerical Computation.

$\endgroup$
4
$\begingroup$

This should work with arbitrary lists replacelist and valuelist.

A = SparseArrya[mat];
A += SparseArray[
    replacelist -> Subtract[valuelist, Extract[A, replacelist]], 
    Dimensions[A]
   ]

Packaged into a function, this could look like this:

Inject[A_SparseArray, replacelist_ -> valuelist_]:= A + SparseArray[
    replacelist -> Subtract[valuelist, Extract[A, replacelist]], 
    Dimensions[A]
   ]
$\endgroup$
1
$\begingroup$

I just did a little check and was surprised to see how slow ReplacePart really is. I guess that you'll want to use Part assignment here since it's really efficient (e.g., mat[[1, 1]] = 999), but I don't know of a way to do Part assignment to the diagonal at once. However, you can easily do it by Scanning over the elements you want to assign or using Do.

mat = RandomInteger[3, {1000, 1000}];
replacelist = {{1, 1}, {3, 3}, {7, 7}, {10, 10}};
RepeatedTiming[ReplacePart[mat, replacelist -> 999];]
RepeatedTiming[
 Scan[Function[mat[[#, #]] = 999], replacelist[[All, 1]]];]
RepeatedTiming[Do[mat[[i, i]] = 999, {i, replacelist[[All, 1]]}];]

As you can see, it's much faster than ReplacePart, even if the elements are done one at a time. The difference is that in Part assignment the matrix mat does not have to be copied over while ReplacePart will create a new copy of the matrix internally.

$\endgroup$
  • $\begingroup$ Setting single entries of a SparseArray with Part is also very, very slow. $\endgroup$ – Henrik Schumacher Feb 2 '18 at 10:27
  • $\begingroup$ Oh, I didn't realize that SparseArrays don't work well with Part assignment. Good to know. $\endgroup$ – Sjoerd Smit Feb 2 '18 at 10:31
  • $\begingroup$ Yeah, the reason is that (in the worst case), the sparsity pattern (internally, that is A["ColumnIndices"] and A["RowPointers"]) have to be recomputed after each assignment. Best strategy is to manipulate as many entries as possible at once, like setting whole rowa, colums, or rectangular subarrays in one go. $\endgroup$ – Henrik Schumacher Feb 2 '18 at 10:35
  • $\begingroup$ So just to make sure I never missed anything: I know you can assign to whole rows, columns or blocks at a time, but I couldn't find a mechanism to assign to arbitrary sub-parts of an array (other than ReplacePart, which isn't doesn't do in-place assignment). Kind of like the reverse of Extract. I take it that WL doesn't have anything for that? $\endgroup$ – Sjoerd Smit Feb 2 '18 at 10:40
  • $\begingroup$ None that I know of... Good point to phrase it as the reverse of Extract! Wouldn't that be called Inject? $\endgroup$ – Henrik Schumacher Feb 2 '18 at 10:43

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.