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

Posted 2023-03-16

【中文标题】提高在具有许多条目设置为 np.inf 的二维 numpy 数组中查找最小元素的速度

【英文标题】：Increase speed of finding minimum element in a 2-D numpy array which has many entries set to np.inf

【发布时间】：2021-06-16 15:59:49

【问题描述】：

我有一个 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() 时，它仍然会遍历整个矩阵。有什么办法可以“忽略”下三角并提高速度？

我试过很多方法，比如：

使用掩码数组 使用triu_indices()提取上三角，然后求最小值 将下方三角形和对角线中的条目设置为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.008 秒。整个算法的运行时间约为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

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

【问题讨论】：

你是如何创建距离矩阵的？如果它是对称的，它可能是自我参照的，对吧？你能用scipy.spatial.distance.pdist代替你现在做的吗？那只输出（并且只计算）上三角形。然后你可以使用argmin 的结果来对抗triu_indices，或者找到一些方法直接计算它（因为所有这些索引都会很大）。

展示你如何计算距离。我想我可以通过完全重写这一步来帮助你

【参考方案1】：

如果我们后退一步，假设距离矩阵是对称的，并且基于具有i 维度的i 点的(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)

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

【讨论】：

您好，首先非常感谢您抽出宝贵时间！整个算法仍然需要 40 分钟才能完成。我刚刚添加了一些我的问题的背景，你能抽出时间看看吗？非常感谢您的帮助！

啊，是的。要使 KDTree 高效，您需要 i > 2**n 。使用 64 个维度，您会遇到暴力破解。是否可以在聚类之前对您的数据进行 PCA 降维？

我会尝试 PCA 并检查准确性是否下降。我希望准确性不会发生太大变化......

【参考方案2】：

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

第一次运行时，请确保不要计时。编译很慢。

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

【讨论】：

首先非常感谢您的宝贵时间！我尝试了numba，它让我的速度提高了 28%。但是，整个算法仍然需要 40 分钟才能完成。我刚刚添加了一些我的问题的背景，你能抽出时间看看吗？非常感谢您的帮助！

你甚至可以使用 prange 来代替 range 来启用多线程。

我更新了计算距离矩阵的方法

【参考方案3】：

可以通过掩码选择数组的上三角，简单示例：

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]，它返回上三角的一维数组，然后您应用最小值

【讨论】：

triu_indices 这样做效率更高，但速度更慢