Быстрый расчет фронта Парето в Python

Ответ 1

Если вас беспокоит реальная скорость, вам определенно нужно использовать 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 раза, предварительно упорядочив данные, см. здесь.

Ответ 2

Обновлено августа 2019 года

Вот еще одна простая реализация, которая довольно быстрая для скромных размеров. Точки ввода предполагаются уникальными.

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

Мы начинаем с сортировки точек по сумме координат. Это полезно, потому что

• Для многих распределений данных точка с наибольшей суммой координат будет доминировать над большим количеством точек.
• Если точка x имеет большую сумму координат, чем точка y, то y не может доминировать над x.

Вот несколько тестов относительно ответа Питера с использованием 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)

Эвристическая выпуклая оболочка

Я недавно остановился на этой проблеме и обнаружил полезную эвристику, которая хорошо работает, если много точек распределены независимо, а измерений мало.

Идея состоит в том, чтобы вычислить выпуклую оболочку точек. С небольшими размерами и независимо распределенными точками количество вершин выпуклой оболочки будет небольшим. Интуитивно можно ожидать, что некоторые вершины выпуклой оболочки будут доминировать во многих исходных точках. Более того, если точка в выпуклой оболочке не доминирует ни в одной другой точке в выпуклой оболочке, то в ней также не доминирует ни одна точка в исходном наборе.

Это дает простой итеративный алгоритм. Мы неоднократно

1. Вычислить выпуклый корпус.
2. Сохраните неоплаченные точки Парето из выпуклой оболочки.
3. Отфильтруйте точки, чтобы удалить те, в которых преобладают элементы выпуклой оболочки.

Я добавляю несколько тестов для измерения 3. Кажется, что для некоторого распределения точек этот подход дает лучшую асимптотику.

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())))

Результаты

# >>= python temp.py # 1,000 points
# Old non-numpy solution:      0.0316428339574486
# Numpy solution:              0.005961259012110531
# Convex hull heuristic:       0.012369581032544374
# >>= python temp.py # 1,000,000 points
# Old non-numpy solution:      70.67529802105855
# Numpy solution:              5.398462114972062
# Convex hull heuristic:       1.5286884519737214
# >>= python temp.py # 10,000,000 points
# Numpy solution:              98.03680767398328
# Convex hull heuristic:       10.203076395904645

Исходное сообщение

Я сделал попытку переписать тот же алгоритм с помощью нескольких настроек. Я думаю, что большинство ваших проблем связано с 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

Ответ 3

Определение dominates неверно. A доминирует B тогда и только тогда, когда он лучше или равен B для всех измерений, и строго лучше, по крайней мере, для одного измерения.

Ответ 4

Питер, хороший ответ.

Я просто хотел обобщить для тех, кто хочет выбрать максимизацию по умолчанию для минимизации. Это тривиальное исправление, но приятно документировать здесь:

def is_pareto(costs, maximise=False):
    """
    :param costs: An (n_points, n_costs) array
    :maximise: boolean. True for maximising, False for minimising
    :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]:
            if maximise:
                is_efficient[is_efficient] = np.any(costs[is_efficient]>=c, axis=1)  # Remove dominated points
            else:
                is_efficient[is_efficient] = np.any(costs[is_efficient]<=c, axis=1)  # Remove dominated points
    return is_efficient