我在3D空间中有一组点,我需要从中找到Pareto前沿.执行速度在这里非常重要,并且当我添加测试点时,时间会非常快.
点集看起来像这样:
[[0.3296170319979843, 0.0, 0.44472108843537406], [0.3296170319979843,0.0, 0.44472108843537406], [0.32920760896951373, 0.0, 0.4440408163265306], [0.32920760896951373, 0.0, 0.4440408163265306], [0.33815192743764166, 0.0, 0.44356462585034007]]
现在,我正在使用这个算法:
def dominates(row, candidateRow):
return sum([row[x] >= candidateRow[x] for x in range(len(row))]) == len(row)
def simple_cull(inputPoints, dominates):
paretoPoints = set()
candidateRowNr = 0
dominatedPoints = set()
while True:
candidateRow = inputPoints[candidateRowNr]
inputPoints.remove(candidateRow)
rowNr = 0
nonDominated = True
while len(inputPoints) != 0 and rowNr < len(inputPoints):
row = inputPoints[rowNr]
if dominates(candidateRow, row):
# If it is worse on all features remove the row from the array
inputPoints.remove(row)
dominatedPoints.add(tuple(row))
elif dominates(row, candidateRow):
nonDominated = False
dominatedPoints.add(tuple(candidateRow))
rowNr += 1
else:
rowNr += 1
if nonDominated:
# add the non-dominated point to the Pareto frontier
paretoPoints.add(tuple(candidateRow))
if len(inputPoints) == 0:
break
return paretoPoints, dominatedPoints
在此处找到:http://code.activestate.com/recipes/578287-multidimensional-pareto-front/
在一组解决方案中找到一组非支配解决方案的最快方法是什么？或者,简而言之,Python能比这个算法做得更好吗？
Pet*_*ter 21
如果你担心实际速度,你肯定想要使用numpy(因为聪明的算法调整可能比使用数组操作获得的效果要小).以下是三种解决方案,它们都计算相同的功能.该is_pareto_efficient_dumb解决方案是在大多数情况下慢,但在随着成本的增加快一些is_pareto_efficient_simple解决方案比许多点哑解决方案更加有效,并最终is_pareto_efficient功能是不易阅读,但最快的(所以都是帕累托有效的!).
import numpy as np
# Very slow for many datapoints. Fastest for many costs, most readable
def is_pareto_efficient_dumb(costs):
"""
Find the pareto-efficient points
:param costs: An (n_points, n_costs) array
:return: A (n_points, ) boolean array, indicating whether each point is Pareto efficient
"""
is_efficient = np.ones(costs.shape[0], dtype = bool)
for i, c in enumerate(costs):
is_efficient[i] = np.all(np.any(costs[:i]>c, axis=1)) and np.all(np.any(costs[i+1:]>c, axis=1))
return is_efficient
# Fairly fast for many datapoints, less fast for many costs, somewhat readable
def is_pareto_efficient_simple(costs):
"""
Find the pareto-efficient points
:param costs: An (n_points, n_costs) array
:return: A (n_points, ) boolean array, indicating whether each point is Pareto efficient
"""
is_efficient = np.ones(costs.shape[0], dtype = bool)
for i, c in enumerate(costs):
if is_efficient[i]:
is_efficient[is_efficient] = np.any(costs[is_efficient]<c, axis=1) # Keep any point with a lower cost
is_efficient[i] = True # And keep self
return is_efficient
# Faster than is_pareto_efficient_simple, but less readable.
def is_pareto_efficient(costs, return_mask = True):
"""
Find the pareto-efficient points
:param costs: An (n_points, n_costs) array
:param return_mask: True to return a mask
:return: An array of indices of pareto-efficient points.
If return_mask is True, this will be an (n_points, ) boolean array
Otherwise it will be a (n_efficient_points, ) integer array of indices.
"""
is_efficient = np.arange(costs.shape[0])
n_points = costs.shape[0]
next_point_index = 0 # Next index in the is_efficient array to search for
while next_point_index<len(costs):
nondominated_point_mask = np.any(costs<costs[next_point_index], axis=1)
nondominated_point_mask[next_point_index] = True
is_efficient = is_efficient[nondominated_point_mask] # Remove dominated points
costs = costs[nondominated_point_mask]
next_point_index = np.sum(nondominated_point_mask[:next_point_index])+1
if return_mask:
is_efficient_mask = np.zeros(n_points, dtype = bool)
is_efficient_mask[is_efficient] = True
return is_efficient_mask
else:
return is_efficient
分析测试(使用从正态分布中提取的点):
10,000个样品,2个成本:
is_pareto_efficient_dumb: Elapsed time is 1.586s
is_pareto_efficient_simple: Elapsed time is 0.009653s
is_pareto_efficient: Elapsed time is 0.005479s
有1,000,000个样品,2个成本:
is_pareto_efficient_dumb: Really, really, slow
is_pareto_efficient_simple: Elapsed time is 1.174s
is_pareto_efficient: Elapsed time is 0.4033s
10,000个样品,15个成本:
is_pareto_efficient_dumb: Elapsed time is 4.019s
is_pareto_efficient_simple: Elapsed time is 6.466s
is_pareto_efficient: Elapsed time is 6.41s
请注意,如果效率是一个问题,您可以通过预先重新排序数据获得2倍左右的加速,请参阅此处.
我最近最后查看了这个问题,发现了一个有用的启发式方法,如果有很多点独立分布且维度很少,那么它很有效.
想法是计算点的凸包.由于具有很少的尺寸和独立分布的点,凸包的顶点数量将很小.直观地,我们可以预期凸包的一些顶点支配许多原始点.此外,如果凸包中的点不受凸包中任何其他点的支配,那么它也不受原始集中任何点的支配.
这给出了一个简单的迭代算法.我们反复
我为维度3添加了一些基准.似乎对于某些点的分布,这种方法会产生更好的渐近性.
def keep_efficient(pts):
'returns Pareto efficient row subset of pts'
# sort points by decreasing sum of coordinates
pts = pts[pts.sum(1).argsort()[::-1]]
# initialize a boolean mask for undominated points
# to avoid creating copies each iteration
undominated = np.ones(pts.shape[0], dtype=bool)
for i in range(pts.shape[0]):
# process each point in turn
n = pts.shape[0]
if i >= n:
break
# find all points not dominated by i
# since points are sorted by coordinate sum
# i cannot dominate any points in 1,...,i-1
undominated[i+1:n] = (pts[i+1:] > pts[i]).any(1)
# keep points undominated so far
pts = pts[undominated[:n]]
return pts
N=10000 d=2
keep_efficient
1.31 ms ± 11.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
is_pareto_efficient
6.51 ms ± 23.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
N=10000 d=3
keep_efficient
2.3 ms ± 13.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
is_pareto_efficient
16.4 ms ± 156 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
N=10000 d=4
keep_efficient
4.37 ms ± 38.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
is_pareto_efficient
21.1 ms ± 115 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
N=10000 d=5
keep_efficient
15.1 ms ± 491 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
is_pareto_efficient
110 ms ± 1.01 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
N=10000 d=6
keep_efficient
40.1 ms ± 211 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
is_pareto_efficient
279 ms ± 2.54 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
N=10000 d=15
keep_efficient
3.92 s ± 125 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
is_pareto_efficient
5.88 s ± 74.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
我通过几次调整重新编写了相同的算法.我认为你的大部分问题都来自于此x.这需要搜索点列表; 按索引删除会更有效率.我还修改了y函数以避免一些冗余的比较.这可能在更高的维度上得心应手.
import numpy as np
from scipy import spatial
from functools import reduce
# test points
pts = np.random.rand(10_000_000, 3)
def filter_(pts, pt):
"""
Get all points in pts that are not Pareto dominated by the point pt
"""
weakly_worse = (pts <= pt).all(axis=-1)
strictly_worse = (pts < pt).any(axis=-1)
return pts[~(weakly_worse & strictly_worse)]
def get_pareto_undominated_by(pts1, pts2=None):
"""
Return all points in pts1 that are not Pareto dominated
by any points in pts2
"""
if pts2 is None:
pts2 = pts1
return reduce(filter_, pts2, pts1)
def get_pareto_frontier(pts):
"""
Iteratively filter points based on the convex hull heuristic
"""
pareto_groups = []
# loop while there are points remaining
while pts.shape[0]:
# brute force if there are few points:
if pts.shape[0] < 10:
pareto_groups.append(get_pareto_undominated_by(pts))
break
# compute vertices of the convex hull
hull_vertices = spatial.ConvexHull(pts).vertices
# get corresponding points
hull_pts = pts[hull_vertices]
# get points in pts that are not convex hull vertices
nonhull_mask = np.ones(pts.shape[0], dtype=bool)
nonhull_mask[hull_vertices] = False
pts = pts[nonhull_mask]
# get points in the convex hull that are on the Pareto frontier
pareto = get_pareto_undominated_by(hull_pts)
pareto_groups.append(pareto)
# filter remaining points to keep those not dominated by
# Pareto points of the convex hull
pts = get_pareto_undominated_by(pts, pareto)
return np.vstack(pareto_groups)
# --------------------------------------------------------------------------------
# previous solutions
# --------------------------------------------------------------------------------
def is_pareto_efficient_dumb(costs):
"""
:param costs: An (n_points, n_costs) array
:return: A (n_points, ) boolean array, indicating whether each point is Pareto efficient
"""
is_efficient = np.ones(costs.shape[0], dtype = bool)
for i, c in enumerate(costs):
is_efficient[i] = np.all(np.any(costs>=c, axis=1))
return is_efficient
def is_pareto_efficient(costs):
"""
:param costs: An (n_points, n_costs) array
:return: A (n_points, ) boolean array, indicating whether each point is Pareto efficient
"""
is_efficient = np.ones(costs.shape[0], dtype = bool)
for i, c in enumerate(costs):
if is_efficient[i]:
is_efficient[is_efficient] = np.any(costs[is_efficient]<=c, axis=1) # Remove dominated points
return is_efficient
def dominates(row, rowCandidate):
return all(r >= rc for r, rc in zip(row, rowCandidate))
def cull(pts, dominates):
dominated = []
cleared = []
remaining = pts
while remaining:
candidate = remaining[0]
new_remaining = []
for other in remaining[1:]:
[new_remaining, dominated][dominates(candidate, other)].append(other)
if not any(dominates(other, candidate) for other in new_remaining):
cleared.append(candidate)
else:
dominated.append(candidate)
remaining = new_remaining
return cleared, dominated
# --------------------------------------------------------------------------------
# benchmarking
# --------------------------------------------------------------------------------
# to accomodate the original non-numpy solution
pts_list = [list(pt) for pt in pts]
import timeit
# print('Old non-numpy solution:s\t{}'.format(
# timeit.timeit('cull(pts_list, dominates)', number=3, globals=globals())))
print('Numpy solution:\t{}'.format(
timeit.timeit('is_pareto_efficient(pts)', number=3, globals=globals())))
print('Convex hull heuristic:\t{}'.format(
timeit.timeit('get_pareto_frontier(pts)', number=3, globals=globals())))
| 归档时间： |
|
| 查看次数： |
11528 次 |
| 最近记录： |