Numpy: how to find the unique local minimum of sub matrixes in matrix A?

Question

Given a matrix A of dimension MxN (4x4), how would one find the next-best minimum of each 2x2 submatrix?

A = array([[ 32673.    ,  15108.2   ,  26767.2   ,   9420.   ],
           [ 32944.2   ,  14604.01  ,  26757.01  ,   9127.2  ],
           [ 26551.2   ,  9257.01   ,  26595.01  ,   9309.2  ],
           [ 26624.    ,   8935.2   ,  26673.2   ,   8982.   ]])

The next-best minimum of a set of submatrixes, is the minimum of that submatrix that does not conflict with the local position of other minima:

Example Algorithm:

1. Find the minimum in A: 8935.2 global coords[3,1], local coords [1,1]
2. No other matrix has been evaluated so no conflict yet.
3. Find the next submatrix min: 8982. gc [3,3], lc [1,1]
4. Conflict exists, find next min in same submatrix: 9309.2 gc [2,3], lc [0,1]
5. Find next submatrix min: 9420 gc [0,3] lc[0,1]
6. Conflict exists, find next min: 26757.01 gc [1,2] lc [1,0]
7. Find next submatrix min: 14604 -- conflict with lc[1,1]
8. Find next submatrix min: 15108.2 -- conflict with lc [0,1]
9. Find next submatrix min: 32673. gc [0,0], lc [0,0]

one approach I have thought of trying is to follow the algorithm above, but instead of exhaustively searching each submatrix again, I globally update each submatrix local position with a 'high' value (>> max(A)), which is incremented on each successful find of a minima.

The expected output would be a list:

[((0, 0), (0, 0), 32673), ((0, 1), (1, 0), 26757.01), ((1, 0), (1, 1), 8935.2), ((1, 1), (0, 1), 9309.2)]

of the form [((t1), (t2), value) ... ], where t1 is the coordinates of the submatrix in A, and t2 is the coordinates of the selected minimum in the submatrix.

Edit: the submatrices are defined as ZxZ, where MxN modulo ZxZ == 0, and are non-overlapping starting at (0,0), and tiled to match the dimensions of MxN.

Edit: Below is a solution I've constructed, but it is slow. I suspect that that if I delete submatrices from the matrix on each iteration, then the performance might be improved, but I am not sure how to do this.

    def get_mins(self, result):
    # result is the 2d array
    dim = 2  # 2x2 submatrix
    mins = []
    count = 0
    while count < dim**2:
        a, b = result.shape
        M4D = result.reshape(a//dim, dim, b//dim, dim)
        lidx = M4D.transpose(0, 2, 1, 3).reshape(-1, b//dim, dim**2).argmin(-1)
        r, c = numpy.unravel_index(lidx, [dim, dim])

        yy = M4D.min(axis=(1, 3))
        ww = numpy.dstack((r, c))

        super_min = numpy.unravel_index(numpy.argmin(yy), (dim, dim))

        rows = super_min[0]
        cols = super_min[1]

        # ww[rows,cols] g_ves us 2x2 position
        offset_r, offset_c = ww[rows, cols]
        # super_min gives us submatrix position

        mins.append((tuple(super_min), (offset_r, offset_c), yy.min()))

        if dim > 1:
            # update all other positions with inf >> max(result)
            result[numpy.ix_([offset_r + (d * dim) for d in range(dim)], [offset_c + (d * dim) for d in range(dim)])] = numpy.inf
            # update the submatrix to all == numpy.inf
            result[rows*dim:((rows*dim)+dim), cols*dim:((cols*dim)+dim)] = numpy.inf
        count += 1
    return mins

Answer 1

You still weren't very clear about the submatrix definition, but from your expected output I've deduced that you want to divide it into 4 non-overlapping arrays - which I can create with reshape and transpose:

In [113]: A1=A.reshape(4,2,2).transpose(0,2,1)
In [114]: A1
Out[114]: 
array([[[ 32673.  ,  26767.2 ],
        [ 15108.2 ,   9420.  ]],

       [[ 32944.2 ,  26757.01],
        [ 14604.01,   9127.2 ]],

       [[ 26551.2 ,  26595.01],
        [  9257.01,   9309.2 ]],

       [[ 26624.  ,  26673.2 ],
        [  8935.2 ,   8982.  ]]])

argmin gives the position in each (in flattened coor)

In [115]: np.argmin(A1[1])
Out[115]: 3
In [116]: [np.argmin(a) for a in A1]
Out[116]: [3, 3, 2, 2]

So no real advantage to working with 2x2 subarrays - let's just ravel them and stick with simpler 1d ones - and single argmin values

In [117]: A2=A1.reshape(4,4)
In [118]: A2
Out[118]: 
array([[ 32673.  ,  26767.2 ,  15108.2 ,   9420.  ],
       [ 32944.2 ,  26757.01,  14604.01,   9127.2 ],
       [ 26551.2 ,  26595.01,   9257.01,   9309.2 ],
       [ 26624.  ,  26673.2 ,   8935.2 ,   8982.  ]])
In [119]: [np.argmin(a) for a in A2]
Out[119]: [3, 3, 2, 2]

At the end, I could convert those indices back to 2d ones:

In [123]: [np.unravel_index(np.argmin(a),(2,2)) for a in A2]
Out[123]: [(1, 1), (1, 1), (1, 0), (1, 0)]

I think the rest is just an iterative search on this A2 structure.

In [124]: A2[1:,3]=np.inf
In [125]: [np.argmin(a) for a in A2]
Out[125]: [3, 2, 2, 2]
In [126]: A2[2:,2]=np.inf
In [127]: [np.argmin(a) for a in A2]
Out[127]: [3, 2, 0, 0]
In [128]: A2[3:,0]=np.inf
In [129]: [np.argmin(a) for a in A2]
Out[129]: [3, 2, 0, 1]

In [139]: A2
Out[139]: 
array([[ 32673.  ,  26767.2 ,  15108.2 ,   9420.  ],
       [ 32944.2 ,  26757.01,  14604.01,       inf],
       [ 26551.2 ,  26595.01,       inf,       inf],
       [      inf,  26673.2 ,       inf,       inf]])

Oops, I thought I'd figured out how you defined the submatrices, but it doesn't look right. But I'll leave this answer. It may help you clear up your question.

Answer 2

Hmm, was a bit more work than expected ^^. The algorithm is roughly:

• sort the values of the matrix a in a 1d array aSrt
• loop through aSrt and identify the 2x2 submatrix according to this value
• if we do not have this submatrix already in lstSubMat, add it by global coordinate of submatrix(0,0)
• now we have a list lstSubMat, which contains the submatrices ordered by there minimum value
• for each submatrix we find now the minimum value under the precondition that the local coordinate is still available (i.e. is not listet in msk). This is stored in result (by global coordinate)

What you can nicely see in the code below is:

• How to find elements in ndarrays
• How to order an adarray rowwise first by second column index, then by first column index
• how to transform lists in tupels and the other way round

Code:

#lc: local coordinates
#gc: global coordinates
#sc: submatrix coordinates


import numpy as np
a = np.array(
    [[ 32673.    ,  15108.2   ,  26767.2   ,   9420.   ],
    [ 32944.2   ,  14604.01  ,  26757.01  ,   9127.2  ],
    [ 26551.2   ,  9257.01   ,  26595.01  ,   9309.2  ],
    [ 26624.    ,   8935.2   ,  26673.2   ,   8982.   ]]
    )
#print(a)



#sort values of a in 1d array
aSrt=np.sort(a.flatten())
#print(aSrt)

#list of submatrix coordinates ordered by their minimum
lstSubMat=[]
for ii in range(0,len(aSrt)):
    #print('just to make things clear:',np.where(a==aSrt[ii]))
    gc=[elem[0] for elem in list(np.where(a==aSrt[ii]))]
    lc = [elem%2 for elem in gc]
    sc = [gc[jj]-lc[jj] for jj in range(0,2)]
    #print('gc:',gc,'sc',sc,'lc:',lc, 'value:',aSrt[0])
    if not sc in lstSubMat:
        lstSubMat.append(sc)
        #lstSubMat[1].append(lc)
        #lstSubMat[2].append(value)

# result is list of gc
result=np.empty((4,2),dtype=int)
#result=np.empty([4,2])
nmbFound=0

#check list with lc
msk=[]

while nmbFound<4:
    sc=lstSubMat[0]
    subMat=a[sc[0]:sc[0]+2,sc[1]:sc[1]+2]
    #print('subMat:',subMat)
    valSubMatSrt=np.sort(subMat.flatten())
    for ii in range(0,4):
        lc=[elem[0] for elem in list(np.where(subMat==valSubMatSrt[ii]))]
        if not lc in msk:
            msk.append(lc)
            #result.append([sc[jj]+lc[jj] for jj in range(0,2)])
            #result[nmbFound]=[sc[jj]+lc[jj] for jj in range(0,2)]
            result[nmbFound,0]=sc[0]+lc[0]
            result[nmbFound,1]=sc[1]+lc[1]
            nmbFound+=1
            #print('gc:',result[-1],'sc',sc,'lc:',lc, 'value:',aSrt[0])
            lstSubMat=lstSubMat[1:]
            break

#print(result)

#sort first by row then by col index of submatrix -> //2
result=result[(result[:,1]//2).argsort()] 
result=result[(result[:,0]//2).argsort()] 
#print(result)

print('\n\nresult:')
for ii in range(0,len(result)):
    sc=tuple([elem//2 for elem in result[ii,:]])
    lc=tuple([result[ii,jj]%2 for jj in range(0,2)])
    print(sc,lc,a[tuple(result[ii,:])])

Output:

result:
(0, 0) (0, 0) 32673.0
(0, 1) (1, 0) 26757.01
(1, 0) (1, 1) 8935.2
(1, 1) (0, 1) 9309.2

Answer 3

Given the dependency between iterations in choosing the global minimum, here's an approach with one-loop -

def unq_localmin(A, dim):
    m, n = A.shape
    M4D = A.reshape(m//dim, dim, n//dim, dim)
    M2Dr = M4D.swapaxes(1,2).reshape(-1,dim**2)
    a = M2Dr.copy()

    N = M2Dr.shape[0]
    R = np.empty(N,dtype=int)
    C = np.empty(N,dtype=int)
    shp = M2Dr.shape
    for i in range(N):
        r,c = np.unravel_index(np.argmin(a),shp)
        a[r] = np.inf
        a[:,c] = np.inf
        R[i], C[i] = r, c
    out = M2Dr[R,C]
    idr = np.column_stack(np.unravel_index(R,(dim,dim)))
    idc = np.column_stack(np.unravel_index(C,(dim,dim)))
    return zip(map(tuple,idr),map(tuple,idc),out)

Let's verify results with a random bigger 9x9 array and 3x3 submatrix/subarray to test out variety against OP's implementation get_mins -

In [66]: A   # Input data array
Out[66]: 
array([[ 927.,  852.,   18.,  949.,  933.,  558.,  519.,  118.,   82.],
       [ 939.,  782.,  178.,  987.,  534.,  981.,  879.,  895.,  407.],
       [ 968.,  187.,  539.,  986.,  506.,  499.,  529.,  978.,  567.],
       [ 767.,  272.,  881.,  858.,  621.,  301.,  675.,  151.,  670.],
       [ 874.,  221.,   72.,  210.,  273.,  823.,  784.,  289.,  425.],
       [ 621.,  510.,  303.,  935.,   88.,  970.,  278.,  125.,  669.],
       [ 702.,  722.,  620.,   51.,  845.,  414.,  154.,  154.,  635.],
       [ 600.,  928.,  540.,  462.,  772.,  487.,  196.,  499.,  208.],
       [ 654.,  335.,  258.,  297.,  649.,  712.,  292.,  767.,  819.]])

In [67]: unq_localmin(A, dim = 3) # Using proposed approach
Out[67]: 
[((0, 0), (0, 2), 18.0),
 ((2, 1), (0, 0), 51.0),
 ((1, 0), (1, 2), 72.0),
 ((1, 1), (2, 1), 88.0),
 ((0, 2), (0, 1), 118.0),
 ((2, 2), (1, 0), 196.0),
 ((2, 0), (2, 2), 258.0),
 ((1, 2), (2, 0), 278.0),
 ((0, 1), (1, 1), 534.0)]

In [68]: out = np.empty((9,9))

In [69]: get_mins(out,A) # Using OP's soln with dim = 3 edited
Out[69]: 
[((0, 0), (0, 2), 18.0),
 ((2, 1), (0, 0), 51.0),
 ((1, 0), (1, 2), 72.0),
 ((1, 1), (2, 1), 88.0),
 ((0, 2), (0, 1), 118.0),
 ((2, 2), (1, 0), 196.0),
 ((2, 0), (2, 2), 258.0),
 ((1, 2), (2, 0), 278.0),
 ((0, 1), (1, 1), 534.0)]

Simplification

The above solution gets us the row and col indices that could be used to construct the indices tuples as printed with get_mins. If you don't need those, we could simplify the proposed approach a bit, like so -

def unq_localmin_v2(A, dim):
    m, n = A.shape
    M4D = A.reshape(m//dim, dim, n//dim, dim)
    M2Dr = M4D.swapaxes(1,2).reshape(-1,dim**2)    
    N = M2Dr.shape[0]
    out = np.empty(N)
    shp = M2Dr.shape
    for i in range(N):
        r,c = np.unravel_index(np.argmin(M2Dr),shp)
        out[i] = M2Dr[r,c]
        M2Dr[r] = np.inf
        M2Dr[:,c] = np.inf        
    return out

Runtime test -

In [52]: A = np.random.randint(11,999,(9,9)).astype(float)

In [53]: %timeit unq_localmin_v2(A, dim=3)
10000 loops, best of 3: 93.1 µs per loop

In [54]: out = np.empty((9,9))

In [55]: %timeit get_mins(out,A)
1000 loops, best of 3: 907 µs per loop