python - 提高在二维 numpy 数组中查找最小元素的速度，该数组有许多条目设置为 np.inf

Question

我有一个 16000*16000 的矩阵，想要找到最小的条目。该矩阵是距离矩阵，因此关于对角线对称。为了每次都能精确地获得一个最小值，我将下三角形和对角线设置为 np.inf。下面是一个 5*5 矩阵示例:

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

我想仅在上三角形中找到最小条目的索引。但是，当我使用 np.argmin() 时，它仍然会遍历整个矩阵。有没有办法“忽略”下三角并提高速度？

我尝试过很多方法，例如:

1. 使用掩码数组
2. 使用triu_indices()提取上三角形，然后找到最小值
3. 将下三角形和对角线中的条目设置为 None 而不是 np.inf，然后使用 np.nanargmin() 查找最小值

但是，我尝试的所有方法都比直接使用 np.argmin() 慢。

感谢您的宝贵时间，如果您能帮助我，我将不胜感激。

更新 1:我的问题的一些背景

事实上，我正在从头开始实现凝聚聚类的修改版本。原始数据集是16000*64(我有16000个点，每个点都是64维)。首先，我构建了 16000 个簇，每个簇只包含一个点。在每次迭代中，我找到最近的 2 个簇并将它们合并，直到满足终止条件。

为了避免重复计算距离，我将距离存储在16000*16000的距离矩阵中。我将对角线和下三角形设置为 np.inf。在每次迭代中，我都会找到距离矩阵中的最小条目，并且该条目的索引对应于 2 个最近的簇，例如 c_i 和 c_j。然后，在距离矩阵中，我将c_i和c_j对应的2行2列填充到np.inf中，这意味着这2个簇被合并并且不合并不再存在了。然后我将计算新簇与所有其他簇之间的距离的数组，然后将该数组放入与 c_i 对应的 1 行 1 列中。

我要明确一点:在整个过程中，距离矩阵的大小始终没有改变。在每次迭代中，对于与我找到的最近的2个簇相对应的2行和2列，我用np.inf填充1行和1列，并将新簇的距离数组放在另一个中1 行 1 列。

现在性能的瓶颈是找到距离矩阵中的最小条目，这需要 0.008s。整个算法的运行时间约为40分钟。

更新 2:如何计算距离矩阵

下面是我用来生成距离矩阵的代码:

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

尽管如此，我需要说的是，生成距离矩阵现在并不是算法的瓶颈，因为我只生成一次，然后我只需更新距离矩阵(如上所述)。

Answer 1

如果我们备份一个步骤，假设距离矩阵是对称的，并且基于 (i, n) 形状的数组，其中 i 点位于 n 维度，距离度量是笛卡尔坐标，这可以使用 KDTree 数据结构非常有效地完成:

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])

这非常快，但它是否对您有用取决于您的应用程序和距离指标。

如果您根本不需要距离矩阵(只需要索引)，您可以这样做:

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

编辑:这对于高维数据来说效果不佳。由于你面临着维度的诅咒，你可能需要使用蛮力。我推荐 scipy.spatial.distance.pdist 为此:

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)

尚未完全测试，但我认为应该可行。