如何根据坐标获取K个最远的点？

Question

我们有一个无聊的 CSV 文件，其中有10000行数据，包括年龄（浮点数）、标题（枚举/整数）、分数（浮点数）等。

• 我们在一张表中有N列，每列都含有int/float值，你可以将其想象为ND空间中的点。
• 我们想要选择K个点，使它们之间的距离最大。

如果我们有100个紧密聚集的点和远处一个点，对于三个点，我们会得到类似以下两种情况之一的结果： or this

当有4个点时，选择中间的某些点会更加有趣。

如何从N个点（包括任意复杂度）中选择K个最远的行（点）呢？这看起来像是具有给定分辨率的ND点云“三角剖分”，但不适用于3d点。

我正在寻找一种相对快速的方法（近似的方法-不需要精确解），用于K=200，N=100000和ND=6（可能是基于多重网格或KDTree、SOM或三角剖分的ANN）... 有人知道吗？

- DuckQueen

---

不，那并不是这样的。但这是一个非常棘手的问题。 - hrokr

你的第一个草图看起来不正确。它应该更像第二个，其中聚类已合并。 - Walter Tross

一个近似解是否可接受？ - Walter Tross

也许这个链接 https://flothesof.github.io/farthest-neighbors.html 可以帮助你。 - Ryabchenko Alexander

我很好奇你最终是如何解决你的问题的。 - Walter Tross

Answer 1

根据以往经验，对于类似问题，一种简单的解决方案是计算每组K个点之间所有配对的欧几里得距离的均值，然后取最大均值。正如上面有人指出的那样，很可能难以避免在所有组合（而不是所有配对）上进行循环。因此，这一切的一个可能的实现方式如下：

import itertools
import numpy as np
from scipy.spatial.distance import pdist

Npoints = 3 # or 4 or 5...
# making up some data:
data = np.matrix([[3,2,4,3,4],[23,25,30,21,27],[6,7,8,7,9],[5,5,6,6,7],[0,1,2,0,2],[3,9,1,6,5],[0,0,12,2,7]])
# finding row indices of all combinations:
c = [list(x) for x in itertools.combinations(range(len(data)), Npoints )]

distances = []
for i in c:    
    distances.append(np.mean(pdist(data[i,:]))) # pdist: a method of computing all pairwise Euclidean distances in a condensed way.

ind = distances.index(max(distances)) # finding the index of the max mean distance
rows = c[ind] # these are the points in question

- BossaNova

---

它能工作，但在c = [list(x) for x in itertools.combinations(range(len(data)), Npoints )]行上消耗了太多的内存。有没有一种更迭代的方法？ - DuckQueen

我知道...当行数N非常大+多维+K大于2时，贪婪搜索需要很长时间。您可以直接在itertool对象上循环，避免使用列表推导式，但这并不一定会减少运行时间。您所问的是一个非平凡的开放性问题...它有点像旅行商问题（虽然不完全相同！）。在这里可以看到一个不同解决方案的好方法（用R）：https://larssonjohan.com/post/2016-10-30-farthest-points/。 - BossaNova

Answer 2

我提出了一个近似解决方案。思路是从一组K个点开始，这些点是以我将在下面解释的方式选择的，并且反复循环这些点，用不属于该集合但包括当前点在内的N-K+1个点替换当前点，使得距离集合中点的距离之和最大化。这个过程会得出一组K个点，其中任何单个点的替换都会导致集合中点之间的距离之和减少。

为了启动这个过程，我们取到所有点的平均值最接近的K个点。这样，第一次循环时，K个点的集合将会分散在其最优解附近。随后的迭代会调整K个点集合，朝着距离之和的最大值方向调整，对于当前的N、K和ND的值，这个过程看起来只需要几秒钟就能达到。为了防止边界情况下的过度循环，尽管我们进行了限制，但我们仍然会在一定程度上迭代。

当迭代不再改善K个点之间的总距离时，我们停止迭代。当然，这是一个局部极大值。其他局部极大值将根据不同的初始条件或允许多个点的替换而达到，但我认为这是不值得的。

数据必须进行调整，以便每个维度上的单位位移具有相同的意义，即为了使欧几里得距离有意义。例如，如果您的维度是工资和孩子数量，则未经调整的算法可能会产生集中在极端工资区域的结果，而忽略那些有10个孩子的人。为了获得更现实的输出，您可以将工资和孩子数量除以它们的标准偏差或其他可比较工资差异和孩子数量差异的估计值。

为了能够为随机高斯分布绘制输出，我在代码中设置了ND = 2，但根据您的要求设置ND = 6没有问题（除非无法绘制）。

import matplotlib.pyplot as plt
import numpy as np
import scipy.spatial as spatial

N, K, ND = 100000, 200, 2
MAX_LOOPS = 20

SIGMA, SEED = 40, 1234
rng = np.random.default_rng(seed=SEED)
means, variances = [0] * ND, [SIGMA**2] * ND
data = rng.multivariate_normal(means, np.diag(variances), N)

def distances(ndarray_0, ndarray_1):
    if (ndarray_0.ndim, ndarray_1.ndim) not in ((1, 2), (2, 1)):
        raise ValueError("bad ndarray dimensions combination")
    return np.linalg.norm(ndarray_0 - ndarray_1, axis=1)

# start with the K points closest to the mean
# (the copy() is only to avoid a view into an otherwise unused array)
indices = np.argsort(distances(data, data.mean(0)))[:K].copy()
# distsums is, for all N points, the sum of the distances from the K points
distsums = spatial.distance.cdist(data, data[indices]).sum(1)
# but the K points themselves should not be considered
# (the trick is that -np.inf ± a finite quantity always yields -np.inf)
distsums[indices] = -np.inf
prev_sum = 0.0
for loop in range(MAX_LOOPS):
    for i in range(K):
        # remove this point from the K points
        old_index = indices[i]
        # calculate its sum of distances from the K points
        distsums[old_index] = distances(data[indices], data[old_index]).sum()
        # update the sums of distances of all points from the K-1 points
        distsums -= distances(data, data[old_index])
        # choose the point with the greatest sum of distances from the K-1 points
        new_index = np.argmax(distsums)
        # add it to the K points replacing the old_index
        indices[i] = new_index
        # don't consider it any more in distsums
        distsums[new_index] = -np.inf
        # update the sums of distances of all points from the K points
        distsums += distances(data, data[new_index])
    # sum all mutual distances of the K points
    curr_sum = spatial.distance.pdist(data[indices]).sum()
    # break if the sum hasn't changed
    if curr_sum == prev_sum:
        break
    prev_sum = curr_sum

if ND == 2:
    X, Y = data.T
    marker_size = 4
    plt.scatter(X, Y, s=marker_size)
    plt.scatter(X[indices], Y[indices], s=marker_size)
    plt.grid(True)
    plt.gca().set_aspect('equal', adjustable='box')
    plt.show()

输出： 将数据划分为 3 个等距的高斯分布，输出结果如下：

- Walter Tross

Answer 3

概括： 处理多个等距点和维度诅咒将是比仅仅找到这些点更大的问题。剧透：有一个惊喜结局。

我认为这是一个有趣的问题，但我对某些答案感到困惑。我认为这部分原因是由于提供的草图。您无疑已经注意到答案看起来相似--2D，带有聚类--即使您指出需要更广泛的范围。因为其他人最终也会看到这个问题，所以我会慢慢地解释我的思路，请在前面耐心等待。

从一个简化的例子开始是有意义的，以查看是否可以用易于理解的数据和线性2D模型来推广解决方案。

但我们不需要计算所有的距离。我们只需要极端值。因此，我们可以选择一些最高和最低值：

right = lin_2_D.nlargest(8, ['x'])
left = lin_2_D.nsmallest(8, ['x'])

graph = sns.scatterplot(x="x", y="y", data=lin_2_D, color = 'gray', marker = '+', alpha = .4)
sns.scatterplot(x = right['x'], y = right['y'],  color = 'red')
sns.scatterplot(x = left['x'], y = left['y'],  color = 'green')

fig = graph.figure
fig.set_size_inches(8,3)

我们目前所拥有的是：在100个点中，我们已经消除了需要计算84个点之间距离的需求。对于剩下的点，我们可以通过将结果排序并检查其与其他点之间的距离来进一步减少这一需求。

你可以想象一个情况，即您有几个数据点远离趋势线，可以通过取最大或最小的y值来捕捉，所有这些看起来都像Walter Tross的顶部图表。再添加几个额外的聚类，就会得到他底部的图表，看起来我们正在阐述同样的观点。

停留在这里的问题在于，您提到的要求是需要解决任意维度的问题。

不幸的是，我们遇到了四个挑战：

挑战1：随着维度的增加，您可能会遇到许多情况，在寻找中点时存在多个解决方案。因此，您正在寻找k个最远的点，但存在大量同样有效的可能解决方案，并且没有优先考虑它们的方法。以下是两个超级简单的例子说明这一点：

在这里，我们只有四个点和两个维度。你真的找不到比这更简单的了，对吗？从红色到绿色的距离微不足道。但是试着找到下一个最远的点，你会发现两个黑点都与红点和绿点等距离。想象一下，如果你使用第一个图表来找到最远的六个点，你可能会有20个或更多个点都等距离。

编辑：我刚注意到红色和绿色的点是在它们的圆的边缘而不是中心，我稍后会更新，但重点是相同的。

B）这很容易想象：想象一下D&D四面体骰子。在三维空间中有四个数据点，所有点距离相等，因此被称为三角形金字塔。如果你正在寻找最接近的两个点，那么哪两个？你有4个选择2（也就是6）种可能的组合。摆脱有效解决方案可能会有些问题，因为你不可避免地会面临问题，比如“为什么我们要摆脱这些而不是这一个？”

挑战2：维度诅咒。说的够清楚了。

挑战3 维度诅咒的复仇 因为你正在寻找最远的点，所以你必须对每个点进行x、y、z...n坐标或者你必须插值它们。现在，你的数据集更大，速度更慢。

挑战4 因为你正在寻找最远的点，因此降维技术如岭回归和套索将没有用处。

那么，对此该怎么办呢？

什么也不做。

等待。什么？！？

并非真正地、确切地和字面上的“什么也不做”。而是采用一种简单的启发式方法，这种方法易于理解和计算。保罗·凯宁表达得很好：

直观地说，当情况足够复杂或不确定时，只有最简单的方法才有效。然而，基于这些稳健适用技术的常识启发式方法可以产生几乎肯定是最优的结果。

在这种情况下，你没有遭受维数灾难，而是拥有维数祝福。的确，你有很多点，当你寻找其他等距点（k）时，它们会按线性比例增加，但空间的总维度体积将增加到维度的幂次方。你所关心的k个最远点与总点数相比微不足道。甚至k^2随着维数的增加也变得微不足道。

现在，如果你的维数较低，我会选择它们作为解决方案（除了使用嵌套for循环的解决方案……在NumPy或Pandas中）。

如果我处在你的位置，我会考虑我已经在其他答案中编写了可以用作基础的代码，并可能会想知道为什么我应该相信这个比其他答案更好，除了它提供了一个思考这个主题的框架。当然，应该有一些数学知识，也许是一些重要人物说同样的话。

让我引用《计算机控制和信号处理的计算密集方法》第18章和通过类比进行扩展的论证，其中包含一些繁重的数学知识。从上面的图表（边缘带有彩色点的图表）可以看出，中心被移除了，特别是如果你遵循了去除极端y值的想法。就像你把一个气球放在一个盒子里一样。你也可以把球放在一个立方体里。将其升级到多个维度，你就有了一个超球体在一个超立方体中。你可以在这里阅读更多关于这种关系的信息。

最后，让我们来谈谈一种启发式方法：

选择每个维度中具有最大或最小值的点。当/如果你用完了这些点，就选择靠近这些值的点，如果在最小/最大值处没有一个点。本质上，您正在选择一个框的角落。对于二维图形，您有四个点，对于三维图形，您有一个框的8个角（2^3）。

更准确地说，这将是一个4D或5D（取决于您如何分配标记形状和颜色）投影到3D。但您可以轻松地看到此数据云给出了完整的尺寸范围。

这里是一个快速的学习检查; 为了简单起见，忽略颜色/形状方面：你可以直观地判断，在决定可能稍微靠近的点之前，你最多可以处理 k 个点而没有问题。如果你有一个k < 2D，你可以看到你可能需要随机选择。如果你添加了另一个点，你可以看到它（k +1）会在一个质心中。所以这是检查：如果你有更多的点，它们会在哪里？我想我必须把这放在底部 - markdown的限制。

因此，对于一个6D数据云，小于64（实际上是65，我们马上就会看到）点的值相当容易。但是......

如果您没有数据云，而是有具有线性关系的数据，则将获得2^(D-1)个点。因此，在线性2D空间中，您有一条直线，在线性3D空间中，您将有一个平面。然后是菱形等等，即使您的形状是弯曲的，这也是正确的。我不想自己制作这个图表，而是使用来自Inversion Labs的优秀帖子Best-fit Surfaces for 3D Data中的图表。

• 如果点数k小于2^D，你需要一个过程来决定哪些点不使用。 线性判别分析 应该在你的备选名单上。尽管如此，你可以通过随机选择一个点来满足解决方案。

• 对于一个额外的点（k = 1 + 2^D），你要寻找一个距离边界空间中心最近的点。

• 当k > 2^D时，可能的解决方案将呈阶乘而非几何级数增长。这可能看起来不直观，所以让我们回到两个圆的例子。在二维空间中，只有两个点可能成为等距离候选点。但如果是三维空间并围绕着一条线旋转这些点，则现在任何在环状区域内的点都可以作为k的解决方案。对于三维例子，它们将是球体。从那里开始就是超球体（n-球体）。再次，2^D的比例。

最后一件事：如果您还不熟悉xarray，那么您应该认真考虑一下。

希望所有这些都有所帮助，也希望您会阅读这些链接。它们值得花时间去看。

*它将具有相同的形状，位于中心位置，顶点位于1/3处。就像有27个六面骰子形状的巨大立方体一样。每个顶点（或最近的点）将确定解决方案。您的原始k+1也必须重新定位。因此，您将选择8个顶点中的2个。最后一个问题：是否值得计算这些点之间的距离（请记住对角线略长于边缘），然后将它们与原始的2^D点进行比较？坦率地说，不值得。牺牲这个解决方案。

- hrokr

Answer 4

假设您将包含 N（10000）行和 D 维度（或特征）的 CSV 文件读入 N*D 矩阵 X 中。您可以按如下方式计算每个点之间的距离并将其存储在距离矩阵中：

import numpy as np
X = np.asarray(X) ### convert to numpy array
distance_matrix = np.zeros((X.shape[0],X.shape[0]))
for i in range(X.shape[0]):
    for j in range(i+1,X.shape[0]): 
    ## We compute triangle matrix and copy the rest. Distance from point A to point B and distance from point B to point A are the same. 
        distance_matrix[i][j]= np.linalg.norm(X[i]-X[j]) ## Here I am calculating Eucledian distance. Other distance measures can also be used.

        #distance_matrix = distance_matrix + distance_matrix.T - np.diag(np.diag(distance_matrix)) ## This syntax can be used to get the lower triangle of distance matrix, which is not really required in your case.
        K = 5 ## Number of points that you want to pick

        indexes = np.unravel_index(np.argsort(distance_matrix.ravel())[-1*K:], distance_matrix.shape)

        print(indexes)

- Ashwin Geet D'Sa

Answer 5

如果您想获取最远的点，可以利用为最近邻居开发的所有方法，只需提供不同的“度量标准”。例如，使用 scikit-learn 的最近邻居和距离度量工具，您可以像这样做。

import numpy as np
from sklearn.neighbors import BallTree
from sklearn.neighbors.dist_metrics import PyFuncDistance
from sklearn.datasets import make_blobs
from matplotlib import pyplot as plt


def inverted_euclidean(x1, x2):
    # You can speed this up using cython like scikit-learn does or numba
    dist = np.sum((x1 - x2) ** 2)
    # We invert the euclidean distance and set nearby points to the biggest possible
    # positive float that isn't inf
    inverted_dist = np.where(dist == 0, np.nextafter(np.inf, 0), 1 / dist)
    return inverted_dist

# Make up some fake data
n_samples = 100000
n_features = 200
X, _ = make_blobs(n_samples=n_samples, centers=3, n_features=n_features, random_state=0)

# We exploit the BallTree algorithm to get the most distant points
ball_tree = BallTree(X, leaf_size=50, metric=PyFuncDistance(inverted_euclidean))

# Some made up query, you can also provide a stack of points to query against
test_point = np.zeros((1, n_features))
distance, distant_points_inds = ball_tree.query(X=test_point, k=10, return_distance=True)
distant_points = X[distant_points_inds[0]]

# We can try to visualize the query results
plt.plot(X[:, 0], X[:, 1], ".b", alpha=0.1)
plt.plot(test_point[:, 0], test_point[:, 1], "*r", markersize=9)
plt.plot(distant_points[:, 0], distant_points[:, 1], "sg", markersize=5, alpha=0.8)
plt.show()

这将绘制类似以下内容的内容：

你可以改进许多方面：

1. 我使用numpy实现了inverted_euclidean距离函数，但你可以尝试像scikit-learn中的距离函数那样使用cython来实现它们。你也可以尝试使用numba进行jit编译。
2. 也许欧几里得距离不是你想用来找到最远点的度量，所以你可以自己实现或者简单地使用scikit-learn提供的。

使用Ball Tree算法（或KdTree算法）的好处在于，对于每个查询点，您只需要进行log(N)次比较即可找到训练集中最远的点。构建Ball Tree本身也需要log(N)次比较，因此最终如果您想为ball tree训练集 (X) 中的每个点找到k个最远的点，则它将具有几乎O(D N log(N))的复杂度（其中D是特征数量），随着k的增加而增加到O(D N^2)。

- lucianopaz