Python Matrix Diagonal Of Inf Without Fill_diagonal

I need to set the diagonal elements of a matrix to Inf. An easy way to do it is to use np.fill_diagonal. np.fill_diagonal(my_matrix, float('inf') However fill_diagonal modifies th

Solution 1:

Approach #1

The magic you are looking for is in NumPy strides that gets us a view of the array without the diagonal elements and as such doesn't occupy anymore of the memory space. Here's the implementation to get such a view -

def nodiag_view(a):
    m = a.shape[0]
    p,q = a.strides
    return np.lib.stride_tricks.as_strided(a[:,1:], (m-1,m), (p+q,q))

Let's take a look at a sample run to verify that its a view -

In [124]: a  # Input array
Out[124]: 
array([[ 0, 61, 43, 26, 21],
       [20,  0, 78, 29, 64],
       [34, 49,  0, 64, 60],
       [36, 96, 67,  0, 75],
       [36, 85, 40, 74,  0]])

# Get the no-diag view
In [125]: a_nodiag = nodiag_view(a)

# Lets's verify by changing elements in the view and that should change# elements in the original array too
In [126]: a_nodiag[:] = 999

In [127]: a
Out[127]: 
array([[  0, 999, 999, 999, 999],
       [999,   0, 999, 999, 999],
       [999, 999,   0, 999, 999],
       [999, 999, 999,   0, 999],
       [999, 999, 999, 999,   0]])

Finally, let's see how we can setup it up to solve your entire problem -

def argmin_without_diag(a):
    a_nodiag = nodiag_view(a)
    idx_nodiag = np.argmin(a_nodiag)
    m = a.shape[0]
    idx = idx_nodiag + np.unravel_index(idx_nodiag, (m-1,m))[0]+1
    return  np.unravel_index(idx, a.shape)

Sample run -

In[142]: aOut[142]: 
array([[ 0, 60, 79, 55, 77],
       [62,  0, 86, 84, 25],
       [32, 96,  0, 74, 89],
       [24, 33, 64,  0, 93],
       [14, 74, 30, 44,  0]])

In[143]: argmin_without_diag(a)
Out[143]: (4, 0)

---

Approach #2

If you worry about both memory and performance, you could temporarily set the diagonal ones as infnite, then get the argmin index and then put back the original diagonal values. Thus, effectively the input array is unchanged. The implementation would look something like this -

defargmin_without_diag_replacement(a):
    # Store diagonal values
    vals = a.ravel()[::a.shape[1]+1].copy()

    # Set diag ones as infinites
    a.ravel()[::a.shape[1]+1] = np.inf

    # Get argmin index
    idx = np.argmin(a)

    # Put back the original diag values
    a.ravel()[::a.shape[1]+1] = vals
    return np.unravel_index(idx, a.shape)

Thus, for (n x n) shaped array, the temporary array would have just n elements.

Sample run -

In[237]: aOut[237]: 
array([[  0.,  95.,  57.,  75.,  92.],
       [ 37.,   0.,  69.,  71.,  62.],
       [ 42.,  72.,   0.,  30.,  57.],
       [ 41.,  80.,  94.,   0.,  26.],
       [ 36.,  45.,  71.,  76.,   0.]])

In[238]: argmin_without_diag_replacement(a)
Out[238]: (3, 4)

Runtime test

In [271]: a = np.random.randint(11,99,(1000,1000)).astype(float)

In [272]: np.fill_diagonal(a,0)

In [273]: %timeit argmin_without_diag(a)
1000 loops, best of 3: 1.76 ms per loop

In [274]: %timeit argmin_without_diag_replacement(a)
1000 loops, best of 3: 688 µs per loop

Solution 2:

orig_mat = np.array([[1.2,2,3],[4,5,6],[7,8,9]])

#set diagonal to inf without making a copy of the array.
orig_mat + np.where(np.eye(orig_mat.shape[0])>0,np.inf,0)
array([[ inf,   2.,   3.],
       [  4.,  inf,   6.],
       [  7.,   8.,  inf]])

#the original array remains untorched.
print(orig_mat)
[[ 1.2  2.   3. ]
 [ 4.   5.   6. ]
 [ 7.   8.   9. ]]

Solution 3:

From what I understand you want to find the index of the min of the matrix, but not include the diagonal elements:

orig_mat = np.array([[1.2,2,3],[4,5,6],[7,8,9]])
min_idx = np.argmin(orig_mat+np.multiply(np.eye(orig_mat.shape[0]),1e9))

Solution 4:

My nomination is the explicit copy and fill, wrapped in a function

def foo(mat):
    a = mat.copy()
    np.fill_diagonal(a,np.inf)
    return a

This does make a copy, but only one.

orig_mat + np.where(np.eye(orig_mat.shape[0])>0,np.inf,0))  

returns a new array as well,but in the process has to make two added temporary ones. So it ends up being slower:

for orig_mat=np.ones((1000,1000))

In [67]: timeit np.argmin(orig_mat + np.where(np.eye(orig_mat.shape[0])>0, np.inf, 0))
100 loops, best of 3: 18.7 ms per loop
In [68]: timeit np.argmin(foo(orig_mat))
100 loops, best of 3: 6.94 ms per loop

Surprisingly @Divakar's solution isn't much faster:

In [69]: timeit argmin_without_diag(orig_mat)
100 loops, best of 3: 5.97 ms per loop