Suppose I have a tensor A (maybe it’s a free-variable covariance matrix). I would like to construct a block diagonal matrix with k repetitions of the matrix along the diagonal.

If I were working with numpy objects, I could achieve this with the code: np.kron(np.eye(k), A)

Kroneker products are available in aesara, inside aesara.tensor.linalg:

import aesara.tensor as at
x = at.linalg.kron(at.eye(k), A)

There is also a helper function in PyMC for block-concatenation of a list of matrices, pm.math.block_diagonal, that functions similarly to scipy.linalg.block_diag