# Replace the diagonal of a list of matrices

I want to replace the diagonal of a square list of matrices with another list of matrices. What is the fastest and most efficient way to do so. (I will need to do this for very large arrays)

Small example below

m1 = {{{{1, 2}, {3, 4}}, {{5, 6}, {7, 8}}}, {{{9, 10}, {11,
12}}, {{13, 14}, {15, 16}}}};

v1 = {{{11, 3}, {7, 5}}, {{13, 97}, {3, 16}}};

m1[[1, 1]] = v1[[1]]; m1[[2, 2]] = v1[[2]];

m1

{{{{11, 3}, {7, 5}}, {{5, 6}, {7, 8}}}, {{{9, 10}, {11,
12}}, {{13, 97}, {3, 16}}}}

• One way is to use ReplacePart, e.g. ReplacePart[m1, {i_,i_} :> v1[[i]]] - not sure about efficiency tough – Lukas Lang Oct 26 '19 at 19:12

Let's generate sample data:

array = ConstantArray[
RandomInteger[{0, 10}, {100, 100}], {100, 100}];
replace = ConstantArray[RandomInteger[{0, 10}, {100, 100}], 100];

array1 = array2 = array3 = array4 = array;


and apply different methods:

RepeatedTiming[
Table[array1[[i, i]] = replace[[i]], {i, Length@replace}];]


{0.0064, Null}

RepeatedTiming[(array2[[#, #]] = replace[[#]]) & /@
Range[Length@replace];]


{0.0071, Null}

RepeatedTiming[array3=ReplacePart[array3, {i_, i_} :> replace[[i]]];]


{12.9, Null}

RepeatedTiming[
ParallelTable[array4[[i, i]] = replace[[i]], {i, Length@replace}];]


{0.016, Null}

array1 == array2 == array3 == array4


True

Suggested by Alx (the fastest so far):

RepeatedTiming[
Do[array1[[i, i]] = replace[[i]], {i, Length@replace}];]


{0.0018, Null}

• I love the second option! Thank you! – ThunderBiggi Oct 26 '19 at 19:49
• On my system Do instead of Table in the first method is 2--3 times faster. – Alx Oct 27 '19 at 0:53
• @Alx same on my machine, thx! I'll add it on my list. – Fraccalo Oct 27 '19 at 8:04

We can also use SparseArrays.

m1 = {{{{1, 2}, {3, 4}}, {{5, 6}, {7, 8}}}, {{{9, 10}, {11, 12}}, {{13, 14}, {15, 16}}}};
v1 = {{{11, 3}, {7, 5}}, {{13, 97}, {3, 16}}};
res = {{{{11, 3}, {7, 5}}, {{5, 6}, {7, 8}}}, {{{9, 10}, {11, 12}}, {{13, 97}, {3, 16}}}};

blockmatrixreplace[mat_, rep_] := Block[
{
diagonalmatrixindices = Flatten[Table[{i, i, j, k}, ##] & @@ ({#, Dimensions[mat][[1]]} & /@ {i, j, k}), 2],
emptydiag,
repten
},
emptydiag = SparseArray[mat] - SparseArray[diagonalmatrixindices -> Flatten@Diagonal@mat, Dimensions[mat]];
repten = SparseArray[diagonalmatrixindices -> Flatten@rep, Dimensions[mat]];
emptydiag + repten // Normal
]

{#1, #2 == res} & @@ RepeatedTiming@blockmatrixreplace[m1, v1]


{0.000055, True}

Or we could use Band in the SparseArray since

With an array a of the same rank as the whole sparse array, Band[start]->a by default inserts a at the position specified by start.

bandblockmatrixreplace[mat_, rep_] := Block[
{diagreps, offdiagreps, reps},
diagreps = MapIndexed[Band[{#2[[1]], #2[[1]], 1, 1}] -> {{#1}} &, rep];
offdiagreps = Table[Band[{j, k, 1, 1}] -> {{mat[[j, k]]}}, {j, #1}, {k, Range[#2]~Complement~{j}}] & @@ Dimensions[mat][[;; 2]] // Flatten;
reps = Join[diagreps, offdiagreps];
Normal@SparseArray[reps, Dimensions@mat]
]
{#1, #2 == res} & @@ RepeatedTiming@bandblockmatrixreplace[m1, v1]


{0.000070, True}

• This indeed answers my question exactly. However, in my real code mat is actually an array of sparse arrays, where the diagonal ones are very dense, but the off-diagonal ones are very sparse (so once array-flattened the whole thing is sparse, like here). So I can't really make them work with my real setup. I think for the 1st suggestion the list of rules doesn't work and in general taking SparseArray[mat] takes a lot of time – ThunderBiggi Nov 2 '19 at 17:25