In matrix of sub-matrices, select diagonal indices of off-diagonal submatrices

Question

Let's say I have a square matrix of size (n*m) x (n*m) that is composed of n x n sub-matrices, with each submatrix being a square of size m x m.

I want to select the diagonal indices of the off-diagonal submatrices. This is the code that I have so far, but I would like to do this without the calls to itertools.permutations if possible.

n_submatrix = 3
submatrix_size = 2
my_matrix = np.zeros([n_submatrix * submatrix_size,
                      n_submatrix * submatrix_size])

for submatrix_index in itertools.permutations(range(n_submatrix),2):
    my_slice0 = slice(submatrix_index[0]*submatrix_size, (submatrix_index[0]+1)*submatrix_size)
    my_slice1 = slice(submatrix_index[1]*submatrix_size, (submatrix_index[1]+1)*submatrix_size)
    np.fill_diagonal(my_matrix[my_slice0,my_slice1],1)

my_matrix

#array([[ 0.,  0.,  1.,  0.,  1.,  0.],
#       [ 0.,  0.,  0.,  1.,  0.,  1.],
#       [ 1.,  0.,  0.,  0.,  1.,  0.],
#       [ 0.,  1.,  0.,  0.,  0.,  1.],
#       [ 1.,  0.,  1.,  0.,  0.,  0.],
#       [ 0.,  1.,  0.,  1.,  0.,  0.]])

Answer 1

That blocky format looks like a job for Kronecker product and luckily we have a NumPy built-in for the same in np.kron. Additionally, we need to use np.eye to create such blocky arrays and feed to np.kron. Thus, to solve our case, we can do something like this -

np.kron(1-np.eye(n_submatrix),np.eye(submatrix_size))

Here's a sample run with the parameters same as listed in the question -

In [46]: n_submatrix = 3
    ...: submatrix_size = 2
    ...: 

In [47]: np.kron(1-np.eye(n_submatrix),np.eye(submatrix_size))
Out[47]: 
array([[ 0.,  0.,  1.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  1.],
       [ 1.,  0.,  0.,  0.,  1.,  0.],
       [ 0.,  1.,  0.,  0.,  0.,  1.],
       [ 1.,  0.,  1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  1.,  0.,  0.]])

Answer 2

I'm sure there are other ways, but in the meantime, how about this:

my_matrix = np.tile(np.eye(submatrix_size), (n_submatrix, n_submatrix))
np.fill_diagonal(my_matrix, 0)

Then again, the concatenation of matrixes might not be particularly efficient, so instead I'd rather compute the target indexes directly, fill them, and possibly clear out the diagonal again (all destructive on a previously allocated matrix):

cells = n_submatrix * submatrix_size
my_matrix = np.zeros((cells, cells))

for x in range(n_submatrix):
    for y in range(submatrix_size):
        r = range(y, n_submatrix * submatrix_size, submatrix_size)
        my_matrix[x * submatrix_size + y, r] = 1

np.fill_diagonal(my_matrix, 0)

That is still easily explained: we step first by number of submatrixes (because that's where the pattern repeats), then for each row we set all indexes in that row at once, offset by the current step in the submatrix.

Again, that seems less efficient than calculating indexes at once. I ended up with the following:

cells = n_submatrix * submatrix_size
my_matrix = np.zeros((cells, cells))

for y in range(submatrix_size):
    r = range(y, n_submatrix * submatrix_size, submatrix_size)
    my_matrix[np.ix_(r, r)] = 1

np.fill_diagonal(my_matrix, 0)

To be fair, the documentation for numpy.ix_ and indexing isn't particularly intuitive, but basically I was looking for a way to compress the notion of "for each of these rows, select all the indexes in that row" and then setting them all at once. Now there is only one loop for each row in the submatrix and I'm out of ideas.

There might be a clever formula to generate all of the indexes in one go (by calculating them modulo some number) for all the diagonals though.

---

Btw. for your code I'd suggest extracting expressions as variables (like cells above) and possibly a bit more succinct names, but that's just me.