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

如何解决提高在 2-D numpy 数组中查找最小元素的速度，该数组的许多条目设置为 np.inf

我有一个 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

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

解决方法

如果我们备份一个步骤，假设距离矩阵是对称的，并且基于一个 (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)

没有完全测试，但我认为应该可以。

,

我能想到的一件事可能是使用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]

第一次运行时一定不要计时。编译很慢。

另一种方法可能是以某种方式使用稀疏矩阵。

,

可以通过屏蔽选择数组的上三角，简单例子：

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

因此，不是将下三角形制作为 inf，您必须生成一个遮罩，其中下三角形为 False，上三角形为 True 并应用遮罩 arr[mask] 返回上三角形的一维数组，然后应用 min