Say I have a set of $n \times n$ matrices $A_1, ..., A_m$ as numpy arrays. I'd like to create the block matrix defined below.
I'm looking for a clean, elegant, and easy-to-interpret way of doing this in numpy. I tried this with
a1, a2 = np.full((2, 2), 1), np.full((2, 2), 2) out = np.block([[a1, (a1+a2)/2], [(a1+a2)/2, a2]])
but that approach doesn't generalize to an arbitrary number of $m$ matrices.
An approach that I found which is general is the following:
A = np.array([a1, a2]) out = (A[:, :, None, :] + A.transpose(1, 0, 2)[None, :, :, :]).reshape(n * m, -1)
but that one, while being efficient, is fairly hard-to-read (this code will be read much more often than written).
scipy.linalg.block_diag gets me halfway there, but I don't get the off-diagonals.
Can anyone think of a good alternative solution? I was thinking of looking into numpy's array-generation-from-function routines, but haven't found a good way of going about that yet.