user1658296 user1658296 - 1 month ago 11
Python Question

NumPy: Finding minimum in submatrices/blocks and corresponding row, column indices

Given a square matrix M (mxm), composed of submatrices of dimension (nxn), where mxm modulo nxn == 0, return a matrix of all minima from each submatrix.

e.g. M: (4x4)

M = array([[19834, 29333, 29333, 33120],
[ 5148, 25560, 29393, 8083],
[ 5148, 29393, 25560, 8083],
[19958, 25748, 25748, 24506]], dtype=uint64)


Each submatrix m is 2x2, with minima returned as follows:

M* = array([[ 5148, 8083],
[ 5148, 8083]], dtype=uint64)


I've tried to reshape the arrays and invoke min(), but it only allows you to do so across one axis.

Additionally I'd like to get the original index of the minima as well, but I've assumed I'd have to search for it after generating M*.

Answer

You could split the two axes into two more resulting in a 4D array with a length of 2 each as the second one of those two split axes, which would be the second and fourth axes in the resulting 4D array. Then, simply find minimum along those two axes for the desired output.

Thus, the implementation would look something like this -

m,n = M.shape
out = M.reshape(m//2,2,n//2,2).min(axis=(1,3))

Sample run -

In [41]: M
Out[41]: 
array([[33, 26, 15, 53, 72, 53],
       [12, 64, 28, 27, 58, 51],
       [61, 42, 70, 92, 61, 95],
       [35, 62, 48, 27, 53, 33]])

In [42]: m,n = M.shape

In [43]: M.reshape(m//2,2,n//2,2).min(axis=(1,3))
Out[43]: 
array([[12, 15, 51],
       [35, 27, 33]])

Getting argmins in each submatrix corresponding to original array

For getting those argmins, we need to do some additional work as listed below -

B = 2 # Blocksize 
m,n = M.shape

# Reshape into 4D array as discussed for the previous problem
M4D = M.reshape(m//B,B,n//B,B)

# Bring the second and fourth axes together and then merge.
# Then, get linear indices within each submatrix
lidx = M4D.transpose(0,2,1,3).reshape(-1,n//B,B**2).argmin(-1)

# Convert those linear indices into row, col indices corresponding
# to submatrix and finally corresponding to the original array
r,c = np.unravel_index(lidx,[B,B])
row = r + (B*np.arange(m//B)[:,None])
col = c + B*np.arange(n//B)

Sample input, output -

In [170]: M
Out[170]: 
array([[40, 91, 90, 72, 86, 44],
       [63, 56, 20, 95, 60, 41],
       [28, 50, 32, 89, 69, 46],
       [41, 41, 33, 81, 30, 63]])

In [171]: np.column_stack((row.ravel(),col.ravel()))
Out[171]: 
array([[0, 0],
       [1, 2],
       [1, 5],
       [2, 0],
       [2, 2],
       [3, 4]])
Comments