# efficient replacing of certain diagonal parts of array/matrix

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!

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.

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]
]


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.

• Setting single entries of a SparseArray with Part is also very, very slow. – Henrik Schumacher Feb 2 '18 at 10:27
• Oh, I didn't realize that SparseArrays don't work well with Part assignment. Good to know. – Sjoerd Smit Feb 2 '18 at 10:31
• 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. – Henrik Schumacher Feb 2 '18 at 10:35
• 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? – Sjoerd Smit Feb 2 '18 at 10:40
• None that I know of... Good point to phrase it as the reverse of Extract! Wouldn't that be called Inject? – Henrik Schumacher Feb 2 '18 at 10:43