想象一下您在3个维度中赋予了n个点的s。任何2分之间的距离都是简单的欧几里得距离。您想从该集合中选择K点的子集Q,以使它们彼此之间最远。换句话说,没有其他k点的子集Q',因此Q中所有对距离的最小距离都小于Q'。
中的Q'。如果n大约是1600万,k大约是300,我们如何有效地执行此操作?
我的猜测是,这是NP-HARD,所以我们可能只想专注于近似值。我能想到的一个想法是使用多维缩放来在一行中对这些点进行排序,然后使用二进制搜索版本来获得该行中最远的点。
这被称为离散p分散(maxmin(问题。
最佳结合在White(1991(和Ravi等人中得到了证明。(1994(给出了问题2的近似,后者证明这种启发式是最好的(除非p = np(。
因子2近似
因子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.47,以获得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*-tree有效地找到候选点P和集合P之间的最小距离。但是,由于for
环仍涉及。<<<<<<<<<<<</p>
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)
我也很确定问题是NP-HARD,我发现最相似的问题是K-中心问题。如果运行时比正确性更重要,那么贪婪算法可能是您的最佳选择:
Q ={}
while |Q| < k
Q += p from S where mindist(p, Q) is maximal
旁注:在类似的问题中,例如,设定问题问题可以证明,贪婪算法的解决方案至少与最佳解决方案一样好。
为了加快事情的速度,我看到了3种可能性:
首先将数据集索引您的数据集,然后执行贪婪的搜索。R-Tree的构造是O(N log N(,但是尽管开发了最近的邻居搜索,但它也可以帮助您找到O(log n(中一组点的最远点。这可能比幼稚的O(k*n(算法更快。
从您的1600万点样本进行了一个子集,并在子集上执行贪婪算法。无论如何,您都是大约的,因此您可能能够更准确。您也可以将其与1.算法结合。
使用迭代方法,并在没有时间时停止。这里的想法是从s随机选择k点(让我们调用此set q'(。然后,在每个步骤中,您从q'切换点p_的最小距离,该距离的最小距离与q'中的另一个距离。与S的随机点相位。如果结果集q''更好地使用q'',则使用q'重复重复。。为了不卡住
如果您负担得起〜k*n距离计算,则可以
- 找到点分布的中心。
- 选择远离中心的点。(并将其从未选择点的集合中删除(。
- 从当前选择的所有点找到最远的点并选择它。
- 重复3.直到您以k点结尾。
找到所有点的最大程度。分成7x7x7素体。对于体素中的所有点,找到最接近其中心的点。返回这些7x7x7点。有些体素可能没有积分,希望不会太多。