2019年8月更新
这里有另一种简单的实现方式，适用于较小的维度并且速度相对较快。输入点被假定为唯一的。
def keep_efficient(pts):
'returns Pareto efficient row subset of pts'
pts = pts[pts.sum(1).argsort()[::-1]]
undominated = np.ones(pts.shape[0], dtype=bool)
for i in range(pts.shape[0]):
n = pts.shape[0]
if i >= n:
break
undominated[i+1:n] = (pts[i+1:] >= pts[i]).any(1)
pts = pts[undominated[:n]]
return pts
我们首先按照坐标之和对点进行排序。这很有用，因为：
- 对于许多数据分布，具有最大坐标和的点将支配大量点。
- 如果点x的坐标和比点y大，则y无法支配x。
以下是与Peter答案相对应的一些基准测试，使用np.random.randn。
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)
凸包启发式算法
最近我研究了这个问题，并找到了一种有用的启发式算法，如果点分布独立且维度很少，则效果很好。
这个想法是计算点的凸包。对于维度较少且独立分布的点，凸包的顶点数量将很小。直观地说，我们可以期望凸包的一些顶点支配许多原始点。此外，如果凸包中的一个点没有被凸包中的任何其他点支配，那么它也不会被原始集合中的任何点支配。
这给出了一个简单的迭代算法。我们重复执行以下步骤：
- 计算凸包。
- 保存凸包中未被支配的 Pareto 点。
- 过滤点以删除被凸包元素支配的点。
我为三维空间添加了一些基准测试。对于某些点分布情况，该方法似乎产生更好的渐近性。
import numpy as np
from scipy import spatial
from functools import reduce
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 = []
while pts.shape[0]:
if pts.shape[0] < 10:
pareto_groups.append(get_pareto_undominated_by(pts))
break
hull_vertices = spatial.ConvexHull(pts).vertices
hull_pts = pts[hull_vertices]
nonhull_mask = np.ones(pts.shape[0], dtype=bool)
nonhull_mask[hull_vertices] = False
pts = pts[nonhull_mask]
pareto = get_pareto_undominated_by(hull_pts)
pareto_groups.append(pareto)
pts = get_pareto_undominated_by(pts, pareto)
return np.vstack(pareto_groups)
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)
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
pts_list = [list(pt) for pt in pts]
import timeit
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())))
结果
原始帖子
我尝试用一些调整来重新编写相同的算法。我认为你遇到的大部分问题来自于inputPoints.remove(row)。这需要搜索点列表；通过索引删除将更加高效。
我还修改了dominates函数，以避免一些冗余比较。在更高的维度中，这可能会很方便。
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
numpy.array([[1,2], [3,4], [2,1], [1,1]])它返回以下结果：[ True False True True]但根据帕累托前沿的定义，它应该返回以下结果：[False False False True]- hyperionbis_efficient[is_efficient] = np.any(costs[is_efficient]<c, axis=1) # 删除支配点is_efficient[i] = True- hyperionb