nondominated sortを高速に計算する

はじめに

表題の通り，高速なnondominated sortの計算に必要な関数を作りました．
nondominated sortとは主に多目的最適化においてサンプルの優劣をつけるために利用されるSort方法です．

追記 (2022年4月21日)

Githubにコードを公開し，pip installで利用できるようにしました．

pip install fast-pareto

以下の方法でimport可能です．

from fast_pareto import is_pareto_front, nondominated_rank

実装

以下の nondominated_rank 関数では入力 costs (shape := (n_points, n_costs)) に対して，ranks (shape := (n_points, )) を返します．

import numpy as np

def is_pareto_front(costs):
    on_front_indices = np.arange(costs.shape[0])
    (n_points, _) = costs.shape
    next_index = 0

    while next_index < len(costs):
        nd_mask = np.any(costs < costs[next_index], axis=1)
        nd_mask[next_index] = True
        # Remove dominated points
        on_front_indices, costs = on_front_indices[nd_mask], costs[nd_mask]
        next_index = np.sum(nd_mask[:next_index]) + 1

    mask = np.zeros(n_points, dtype=np.bool8)
    mask[on_front_indices] = True
    return mask

    
def nondominated_rank(costs):
    ranks = np.zeros(len(costs), dtype=np.int32)
    rank = 0
    indices = np.arange(len(costs))
    while indices.size > 0:
        on_front = is_pareto_front(costs)
        ranks[indices[on_front]] = rank
        # Remove pareto front points
        indices, costs = indices[~on_front], costs[~on_front]
        rank += 1

    return ranks

ベンチマーク

Input

for n_points in [100, 1000, 10000]:
    for n_costs in [1, 5, 10, 50]:
        print(f"n_points={n_points}, n_costs={n_costs}")
        %time nondominated_rank(np.random.normal(size=(n_points, n_costs)))
        print("\n")

Output

n_points=100, n_costs=1
CPU times: user 10.7 ms, sys: 0 ns, total: 10.7 ms
Wall time: 10.2 ms


n_points=100, n_costs=5
CPU times: user 3.29 ms, sys: 0 ns, total: 3.29 ms
Wall time: 3.3 ms


n_points=100, n_costs=10
CPU times: user 3 ms, sys: 0 ns, total: 3 ms
Wall time: 3 ms


n_points=100, n_costs=50
CPU times: user 3.57 ms, sys: 0 ns, total: 3.57 ms
Wall time: 3.57 ms


n_points=1000, n_costs=1
CPU times: user 105 ms, sys: 0 ns, total: 105 ms
Wall time: 105 ms


n_points=1000, n_costs=5
CPU times: user 37.8 ms, sys: 0 ns, total: 37.8 ms
Wall time: 37.8 ms


n_points=1000, n_costs=10
CPU times: user 53.5 ms, sys: 0 ns, total: 53.5 ms
Wall time: 53.5 ms


n_points=1000, n_costs=50
CPU times: user 90.4 ms, sys: 0 ns, total: 90.4 ms
Wall time: 90.4 ms


n_points=10000, n_costs=1
CPU times: user 3.36 s, sys: 0 ns, total: 3.36 s
Wall time: 3.36 s


n_points=10000, n_costs=5
CPU times: user 1.22 s, sys: 0 ns, total: 1.22 s
Wall time: 1.22 s


n_points=10000, n_costs=10
CPU times: user 2.72 s, sys: 0 ns, total: 2.72 s
Wall time: 2.72 s


n_points=10000, n_costs=50
CPU times: user 14.4 s, sys: 0 ns, total: 14.4 s
Wall time: 14.4 s

(non-dominant set or 非劣解に関する) 補足

\newcommand{\xv}{\boldsymbol{x}}

多出力関数 $f: \mathbb{R}^D \rightarrow \mathbb{R}^M$ の各目的最小化問題を考える．
任意の $i \in [1, M]$ に関して，$f_i(\xv) \leq f_i(\xv^\prime)$ であり，ある $i \in [1, M]$ が存在し，$f_i(\xv) < f_i(\xv^\prime)$ であるとき $f(\xv)$ は$f(\xv^\prime)$を dominate (優越) するという．
仮に$f(\xv)$ を優越するような $f(\xv^\prime)$ が存在しない場合にそのような $f(\xv)$ はパレート最適と呼ばれて，non-domination rankが 1 であると定義される．
また，non-domination rankが $n$ もしくは $n$ 以下であるような $f(\xv^\prime)$ を全て除外した集合において$f(\xv)$がパレート最適を達成する場合にそのような $f(\xv)$ は non-domination rankが $n + 1$ であると定義される．