在给定的点集中选择最远点的子集
Choosing subset of farthest points in given set of points
假设给定了一个包含 3 个维度的 n 个点的集合 S。任意两点之间的距离是简单的欧氏距离。你想从这个集合中选择 k 个点的子集 Q，使它们彼此相距最远。换句话说，不存在 k 个点的其他子集 Q'，使得 Q 中所有成对距离的最小值小于 Q' 中的最小值。
如果 n 约为 1600 万，k 约为 300，我们如何有效地做到这一点？
我的猜测是，这个 NP-hard 问题可能是因为我们只想关注近似值。我能想到的一个想法是使用多维缩放将这些点排成一条直线，然后使用二进制搜索的版本来获取这条线上相距最远的点。
如果你有能力进行 ~ k*n 距离计算，那么你就可以
- 找到点分布的中心。
- Select离中心最远的点。 （并将其从 un-selected 点的集合中删除）。
- 找到距离所有当前 select 点最远的点并 select 它。
- 重复 3. 直到你得到 k 分。
我也很确定问题是 NP-Hard，我发现最相似的问题是 k-Center Problem。如果运行时比正确性更重要，贪心算法可能是您的最佳选择：
Q ={}
while |Q| < k
Q += p from S where mindist(p, Q) is maximal
旁注：在类似的问题中，例如 set-cover problem 可以证明贪心算法的解至少是最优解的 63%。
为了加快速度，我看到了 3 种可能性：
首先在 R-Tree 中索引您的数据集，然后执行贪婪搜索。 R-Tree 的构造是 O(n log n)，但尽管它是为最近邻搜索而开发的，但它也可以帮助您在 O(log n) 中找到一组点的最远点。这可能比朴素的 O(k*n) 算法更快。
从你的 1600 万个点中抽取一个子集并对子集执行贪心算法。无论如何你都是近似的，所以你可以节省更多的准确性。您也可以将其与 1. 算法结合使用。
使用迭代方法，时间到了就停止。这里的想法是从 S 中随机 select k 个点（我们称这个集合为 Q'）。然后在每个步骤中，将点 p_ 从 Q' 切换到 Q' 中具有最小距离的点 p_ 以及来自 S 的随机点。如果结果集 Q'' 更好地继续 Q''，否则重复 Q' .为了不被卡住，如果您在几次迭代中找不到足够的替代点，您可能想从 Q' 中选择另一个点而不是 p_。
找出所有点的最大范围。分成 7x7x7 体素。对于体素中的所有点，找到最接近其中心的点。 Return 这些 7x7x7 点。一些体素可能不包含点，希望不要太多。
这称为离散 p 色散 (maxmin) 问题。
最优性界限在 White (1991) and Ravi et al. (1994) 中得到证明，给出问题的因子 2 近似值，后者证明这种启发式是最好的（除非 P=NP）。
因子 2 近似值
factor-2 近似值如下：
Let V be the set of nodes/objects
Let i and j be two nodes at maximum distance
Let p be the number of objects to choose
p = set([i,j])
while size(P)<p:
Find a node v in V-P such that min_{v' in P} dist(v,v') is maximum
\That is: find the node with the greatest minimum distance to the set P
P = P.union(v)
Output P
您可以像这样在 Python 中实现它：
#!/usr/bin/env python3
import numpy as np
p = 50
N = 400
print("Building distance matrix...")
d = np.random.rand(N,N) #Random matrix
d = (d + d.T)/2 #Make the matrix symmetric
print("Finding initial edge...")
maxdist = 0
bestpair = ()
for i in range(N):
for j in range(i+1,N):
if d[i,j]>maxdist:
maxdist = d[i,j]
bestpair = (i,j)
P = set()
P.add(bestpair[0])
P.add(bestpair[1])
print("Finding optimal set...")
while len(P)<p:
print("P size = {0}".format(len(P)))
maxdist = 0
vbest = None
for v in range(N):
if v in P:
continue
for vprime in P:
if d[v,vprime]>maxdist:
maxdist = d[v,vprime]
vbest = v
P.add(vbest)
print(P)
精确解
您也可以将其建模为 MIP。对于 p=50，n=400 6000 秒后，最优性差距仍为 568%。近似算法用了 0.47s 来获得 100%（或更少）的最优性差距。朴素的 Gurobi Python 表示可能如下所示：
#!/usr/bin/env python
import numpy as np
import gurobipy as grb
p = 50
N = 400
print("Building distance matrix...")
d = np.random.rand(N,N) #Random matrix
d = (d + d.T)/2 #Make the matrix symmetric
m = grb.Model(name="MIP Model")
used = [m.addVar(vtype=grb.GRB.BINARY) for i in range(N)]
objective = grb.quicksum( d[i,j]*used[i]*used[j] for i in range(0,N) for j in range(i+1,N) )
m.addConstr(
lhs=grb.quicksum(used),
sense=grb.GRB.EQUAL,
rhs=p
)
# for maximization
m.ModelSense = grb.GRB.MAXIMIZE
m.setObjective(objective)
# m.Params.TimeLimit = 3*60
# solving with Glpk
ret = m.optimize()
缩放
显然，初始点的 O(N^2) 缩放很糟糕。通过认识到这对必须位于数据集的凸包上，我们可以更有效地找到它们。这给了我们一个 O(N log N) 的方法来找到对。一旦我们找到它，我们就会像以前一样进行（使用 SciPy 加速）。
缩放的最佳方法是使用 R* 树有效地找到候选点 p 和集合 P 之间的最小距离。但这在 Python 中无法有效地完成，因为 for 循环仍然存在。
import numpy as np
from scipy.spatial import ConvexHull
from scipy.spatial.distance import cdist
p = 300
N = 16000000
# Find a convex hull in O(N log N)
points = np.random.rand(N, 3) # N random points in 3-D
# Returned 420 points in testing
hull = ConvexHull(points)
# Extract the points forming the hull
hullpoints = points[hull.vertices,:]
# Naive way of finding the best pair in O(H^2) time if H is number of points on
# hull
hdist = cdist(hullpoints, hullpoints, metric='euclidean')
# Get the farthest apart points
bestpair = np.unravel_index(hdist.argmax(), hdist.shape)
P = np.array([hullpoints[bestpair[0]],hullpoints[bestpair[1]]])
# Now we have a problem
print("Finding optimal set...")
while len(P)<p:
print("P size = {0}".format(len(P)))
distance_to_P = cdist(points, P)
minimum_to_each_of_P = np.min(distance_to_P, axis=1)
best_new_point_idx = np.argmax(minimum_to_each_of_P)
best_new_point = np.expand_dims(points[best_new_point_idx,:],0)
P = np.append(P,best_new_point,axis=0)
print(P)
假设给定了一个包含 3 个维度的 n 个点的集合 S。任意两点之间的距离是简单的欧氏距离。你想从这个集合中选择 k 个点的子集 Q，使它们彼此相距最远。换句话说，不存在 k 个点的其他子集 Q'，使得 Q 中所有成对距离的最小值小于 Q' 中的最小值。
如果 n 约为 1600 万，k 约为 300，我们如何有效地做到这一点？
我的猜测是，这个 NP-hard 问题可能是因为我们只想关注近似值。我能想到的一个想法是使用多维缩放将这些点排成一条直线，然后使用二进制搜索的版本来获取这条线上相距最远的点。
如果你有能力进行 ~ k*n 距离计算，那么你就可以
- 找到点分布的中心。
- Select离中心最远的点。 （并将其从 un-selected 点的集合中删除）。
- 找到距离所有当前 select 点最远的点并 select 它。
- 重复 3. 直到你得到 k 分。
我也很确定问题是 NP-Hard，我发现最相似的问题是 k-Center Problem。如果运行时比正确性更重要，贪心算法可能是您的最佳选择：
Q ={}
while |Q| < k
Q += p from S where mindist(p, Q) is maximal
旁注：在类似的问题中，例如 set-cover problem 可以证明贪心算法的解至少是最优解的 63%。
为了加快速度，我看到了 3 种可能性：
首先在 R-Tree 中索引您的数据集，然后执行贪婪搜索。 R-Tree 的构造是 O(n log n)，但尽管它是为最近邻搜索而开发的，但它也可以帮助您在 O(log n) 中找到一组点的最远点。这可能比朴素的 O(k*n) 算法更快。
从你的 1600 万个点中抽取一个子集并对子集执行贪心算法。无论如何你都是近似的，所以你可以节省更多的准确性。您也可以将其与 1. 算法结合使用。
使用迭代方法，时间到了就停止。这里的想法是从 S 中随机 select k 个点（我们称这个集合为 Q'）。然后在每个步骤中，将点 p_ 从 Q' 切换到 Q' 中具有最小距离的点 p_ 以及来自 S 的随机点。如果结果集 Q'' 更好地继续 Q''，否则重复 Q' .为了不被卡住，如果您在几次迭代中找不到足够的替代点，您可能想从 Q' 中选择另一个点而不是 p_。
找出所有点的最大范围。分成 7x7x7 体素。对于体素中的所有点，找到最接近其中心的点。 Return 这些 7x7x7 点。一些体素可能不包含点，希望不要太多。
这称为离散 p 色散 (maxmin) 问题。
最优性界限在 White (1991) and Ravi et al. (1994) 中得到证明，给出问题的因子 2 近似值，后者证明这种启发式是最好的（除非 P=NP）。
因子 2 近似值
factor-2 近似值如下：
Let V be the set of nodes/objects
Let i and j be two nodes at maximum distance
Let p be the number of objects to choose
p = set([i,j])
while size(P)<p:
Find a node v in V-P such that min_{v' in P} dist(v,v') is maximum
\That is: find the node with the greatest minimum distance to the set P
P = P.union(v)
Output P
您可以像这样在 Python 中实现它：
#!/usr/bin/env python3
import numpy as np
p = 50
N = 400
print("Building distance matrix...")
d = np.random.rand(N,N) #Random matrix
d = (d + d.T)/2 #Make the matrix symmetric
print("Finding initial edge...")
maxdist = 0
bestpair = ()
for i in range(N):
for j in range(i+1,N):
if d[i,j]>maxdist:
maxdist = d[i,j]
bestpair = (i,j)
P = set()
P.add(bestpair[0])
P.add(bestpair[1])
print("Finding optimal set...")
while len(P)<p:
print("P size = {0}".format(len(P)))
maxdist = 0
vbest = None
for v in range(N):
if v in P:
continue
for vprime in P:
if d[v,vprime]>maxdist:
maxdist = d[v,vprime]
vbest = v
P.add(vbest)
print(P)
精确解
您也可以将其建模为 MIP。对于 p=50，n=400 6000 秒后，最优性差距仍为 568%。近似算法用了 0.47s 来获得 100%（或更少）的最优性差距。朴素的 Gurobi Python 表示可能如下所示：
#!/usr/bin/env python
import numpy as np
import gurobipy as grb
p = 50
N = 400
print("Building distance matrix...")
d = np.random.rand(N,N) #Random matrix
d = (d + d.T)/2 #Make the matrix symmetric
m = grb.Model(name="MIP Model")
used = [m.addVar(vtype=grb.GRB.BINARY) for i in range(N)]
objective = grb.quicksum( d[i,j]*used[i]*used[j] for i in range(0,N) for j in range(i+1,N) )
m.addConstr(
lhs=grb.quicksum(used),
sense=grb.GRB.EQUAL,
rhs=p
)
# for maximization
m.ModelSense = grb.GRB.MAXIMIZE
m.setObjective(objective)
# m.Params.TimeLimit = 3*60
# solving with Glpk
ret = m.optimize()
缩放
显然，初始点的 O(N^2) 缩放很糟糕。通过认识到这对必须位于数据集的凸包上，我们可以更有效地找到它们。这给了我们一个 O(N log N) 的方法来找到对。一旦我们找到它，我们就会像以前一样进行（使用 SciPy 加速）。
缩放的最佳方法是使用 R* 树有效地找到候选点 p 和集合 P 之间的最小距离。但这在 Python 中无法有效地完成，因为 for 循环仍然存在。
import numpy as np
from scipy.spatial import ConvexHull
from scipy.spatial.distance import cdist
p = 300
N = 16000000
# Find a convex hull in O(N log N)
points = np.random.rand(N, 3) # N random points in 3-D
# Returned 420 points in testing
hull = ConvexHull(points)
# Extract the points forming the hull
hullpoints = points[hull.vertices,:]
# Naive way of finding the best pair in O(H^2) time if H is number of points on
# hull
hdist = cdist(hullpoints, hullpoints, metric='euclidean')
# Get the farthest apart points
bestpair = np.unravel_index(hdist.argmax(), hdist.shape)
P = np.array([hullpoints[bestpair[0]],hullpoints[bestpair[1]]])
# Now we have a problem
print("Finding optimal set...")
while len(P)<p:
print("P size = {0}".format(len(P)))
distance_to_P = cdist(points, P)
minimum_to_each_of_P = np.min(distance_to_P, axis=1)
best_new_point_idx = np.argmax(minimum_to_each_of_P)
best_new_point = np.expand_dims(points[best_new_point_idx,:],0)
P = np.append(P,best_new_point,axis=0)
print(P)