Increase speed of finding minimum element in a 2-D numpy array which has many entries set to np.inf

user_15426295 0
Viewed 77

I have a 16000*16000 matrix and want to find the minimum entry. This matrix is a distance matrix, so it is symmetric about diagonal. In order to get exactly one minimum at each time, I set the lower triangle and the diagonal to np.inf. Below is an 5*5 matrix example:

inf a0  a1  a2  a3
inf inf a4  a5  a6
inf inf inf a7  a8
inf inf inf inf a9
inf inf inf inf inf

I want to find the index of the minimum entry only in the upper triangle. However, when I use np.argmin(), it will still go through the whole matrix. Is there any way to "ignore" the lower triangle and increase speed?

I have tried many methods, such as:

  1. Use masked array
  2. Use triu_indices() to extract the upper triangle and then find the minimum
  3. Set the entries in the lower triangle and diagonal to None instead of np.inf, then use np.nanargmin() to find the minimum

However, all of the methods I tried are slower the using np.argmin() directly.

Thank you for your time, I would appreciate it if you can help me.

UPDATE 1: Some background of my problem

In fact, I am implementing a modified version of agglomerative clustering from scratch. The original dataset is 16000*64 (I have 16000 points, each is 64-dimensional). At first, I build 16000 clusters and each contains exactly one point. In each iteration, I find the nearest 2 clusters and merge them, until meet the terminate condition.

To avoid repeated calculation of distances, I store the distances in a 16000*16000 distance matrix. I set the diagonal and lower triangle to np.inf. In each iteration, I will find the smallest entry in the distance matrix, and the index of this entry corresponds to the 2 nearest clusters, say c_i and c_j. Afterwards, in the distance matrix, I fill the 2 rows and 2 columns corresponding to c_i and c_j to np.inf, which means that these 2 clusters are merged and do not exist anymore. Then I will calculate an array of the distances between the new cluster and all other clusters, then put the array in the 1 row and 1 column corresponding to c_i.

Let me make it clear: in the whole process, the size of the distance matrix never change. In each iteration, for the 2 rows and 2 columns correspond to the 2 nearest clusters I found, I fill 1 row and 1 column with np.inf and put the distance array of the new cluster in the other 1 row and 1 column.

Now the bottleneck of the performance is finding the smallest entry in the distance matrix, which takes 0.008s. The run time of the whole algorithm is about 40 minutes.

UPDATE 2: How I compute distance matrix

Below is the code I used to generate distance matrix:

from sklearn.metrics import pairwise_distances

dis_matrix = pairwise_distances(dataset)

for i in range(num_dim):
    for j in range(num_dim):
        if i >= j or (cluster_list[i].contain_reference_point and cluster_list[j].contain_reference_point):
            dis_matrix[i][j] = np.inf

Nevertheless, I need to say that generating the distance matrix is not the bottleneck in the algorithm now, because I generate it only once, and then I just update the distance matrix (as mentioned above).

minimum matrix numpy python performance
3 Answers

If we back up a step, assuming the distance matrix is symmetric and based on an (i, n) shaped array with i points in n dimensions, and the distance metric is cartesian, this can be done very efficiently with a KDTree data structure:

i = 16000
n = 3
points = np.random.rand(i, n) * 100

from scipy.spatial import cKDTree
tree = cKDTree(points)
close = tree.sparse_distance_matrix(tree, 
                                    max_distance = 1, #can tune for your application
                                    output_type  = "coo_matrix") 
close.eliminate_zeros()
ix = close.data.argmin()
i, j = (close.row[ix], close.col[ix])

This is pretty blazing fast, but it depends on your application and distance metric if it's useful for you.

If you don't need the distance matrix at all (and only need indices), you can do:

d, ix = tree.query(points, 2)
j, i = ix[d[:, 1].argmin()]

EDIT: this doesn't work well for high-dimensionality data. Since you're up against the curse of dimensionality, you'll probably need to brute force. I recommend scipy.spatial.distance.pdist for this:

from scipy.spatial.distance import pdist
D = pdist(points, metric = 'seuclidean')  # this only returns the upper diagonal
ix = np.argmin(D)

def ix_to_ij(ix, n):
    sorter = np.arange(n-1)[::-1].cumsum()
    j = np.searchsorted(sorter, ix)
    i = ix - sorter[j]
    return i, j

ix_to_ij(ix, 16000)

Not completely tested but I think that should work.

user_44277770

One thing I can think of that might give you a boost is using numba.njit:

@njit
def upper_min(m):
    x = np.inf
    for r in range(0, m.shape[0] - 1):
        for c in range(r + 1, m.shape[1] + 1):
            if x < m[r, c]:
                x = m[r, c]

Be sure not to time it the first time you run it. The compilation is slow.

Another way could be to use sparse matrices somehow.

user_29887300

You can select upper triangle of array by masking, simple example:

import numpy as np
arr = np.array([[0, 1], [2, 3]])
# Mask of upper triangle
mask = np.array([[True, True],[False, True]])
# Masking returns only upper triangle as 1D array
min_val = np.min(arr[mask]) # Equal to np.min([0, 1, 3])

So instead of making lower triangle as inf, you have to generate a mask where lower triangle is False and upper triangle is True and apply masking arr[mask] which returns 1D array of upper triangle, then you apply min

user_82405100
Related