Approximating Pareto Sum via Bounded Monotone Min-Plus Convolution

(Sketch.) We first shift all $x$-coordinates in $P$ and $Q$ so that the minimum $x$-value becomes 0; this preserves all dominance relations in the Pareto sum.

Next, we construct two sequences $A$ and $B$ of length $W$, initialized with value $W$. For each point $(x,y)\in P$ (resp. $(x,y)\in Q$), we set $A[x]=y$ (resp. $B[x]=y$) and then enforce monotonicity by a single left-to-right sweep: each entry is replaced by the minimum seen so far. After this step, $A[i]$ (resp. $B[i]$) equals the smallest $y$-value among all points of $P$ (resp. $Q$) with $x^{\prime }\le i$. This construction runs in $𝒪(W)$ time. Computing the min-plus convolution $C[k]={\min }_{i+j=k}A[i]+B[j]$ then yields, for each $k$, exactly the smallest achievable $y$-coordinate among all pairs $p\in P$, $q\in Q$ with $p.x+q.x=k$. Thus $C$ contains all candidate points of the Pareto sum. A final linear-time sweep removes dominated entries from $C$, leaving precisely the non-dominated pairs $(k,C[k])$. Hence the Pareto sum of $P$ and $Q$ can be obtained via a bounded monotone min-plus convolution instance computable in $𝒪(W)$ time. A more detailed proof can be found in the full version of the paper. $◀$

According to our theorem, if bounded monotone min-plus convolution can be solved in time $T(n,W)$ for arrays of size $n$ with entries in $[W]$, then we can solve bounded Pareto Sum in time $O(W+T(W,W))$. This makes algorithms for bounded monotone min-plus convolution applicable to Pareto Sum with little overhead to carry out the reduction and the back-transformation of the solution. By applying the algorithm of [11], we get:

We first argue that the set $Z$ computed in Algorithm 1 satisfies the relaxed domination condition with $\mathrm{\Delta }=2t$, specifically, that for all $p+q$ with $p\in P,q\in Q$, there exists some $z=\mathrm{rep}(p^{\prime })+\mathrm{rep}(q^{\prime })\in Z$ with $z\le p+q+(2t,2t)$. Indeed, for any $p\in P,q\in Q$, consider their rounded points $p^{\prime }=(\lfloor \frac{p.x}{t}\rfloor ,\lfloor \frac{p.y}{t}\rfloor )\in P^{\prime }$, $q^{\prime }=(\lfloor \frac{q.x}{t}\rfloor ,\lfloor \frac{q.y}{t}\rfloor )\in Q^{\prime }$. By definition of the Pareto sum $S^{\prime }$ of $P^{\prime },Q^{\prime }$, there must exist some $(\tilde{p}+\tilde{q})\in S^{\prime }$ with $\tilde{p}+\tilde{q}\le p^{\prime }+q^{\prime }$. We claim that $z:=\mathrm{rep}(\tilde{p})+\mathrm{rep}(\tilde{q})$ approximately dominates $p+q$, since

where we used $\mathrm{rep}(x^{\prime })\le t{x}^{\prime }+(t,t)$ for any rounded point $x^{\prime }$ in the first line, and $t{x}^{\prime }\le x$ for any $x$ and rounded point $x^{\prime }$ in the third line. We claim that the lemma follows: Filtering out dominated points in $Z$ ensures that $\tilde{S}$ is a Pareto set (establishing Condition 3.), without affecting the domination condition (if $z\in Z$ dominates $z^{\prime }\in Z$, then whenever $z^{\prime }\le p+q+(t,t)$, then $z\le z^{\prime }\le p+q+(t,t)$, so $z^{\prime }$ is never needed to ensure the condition). Finally, clearly $\tilde{S}\subseteq Z\subseteq P+Q$ as any element of $Z$ is explicitly given as $\mathrm{rep}(p^{\prime })+\mathrm{rep}(q^{\prime })$ for corresponding $\mathrm{rep}(p^{\prime })\in P$ and $\mathrm{rep}(q^{\prime })\in Q$. $◀$

Given Lemma 6, it remains to show that ApproximateParetoSum($P,Q,t)$ with $t:=\mathrm{\Delta }/2$ can be computed in the desired running time. After $O(n)$-time preprocessing, the main step in Algorithm 1 is to compute the Pareto sum $S^{\prime }$ of $P^{\prime },Q^{\prime }$. By our rounding, $P^{\prime },Q^{\prime }$ have integer coordinates in $\{0,\mathrm{\dots },\lfloor W/\mathrm{\Delta }\rfloor \}$, and Theorem 3 is applicable to obtain $S^{\prime }$ and thus $Z$ in time $T(\lfloor W/\mathrm{\Delta }\rfloor ,\lfloor W/\mathrm{\Delta }\rfloor )$. Observe that $Z$ is a set of size at most $O(W/\mathrm{\Delta })$. Thus, we can compute its Pareto front in time $O((W/\mathrm{\Delta })\log (W/\mathrm{\Delta }))$, concluding the claim. $◀$

Follows from the $\tilde{O}(n^{1.6})$ time algorithm of [11], and our witness reporting extension, described in the full version of the paper⁴⁴4We remark that it is very plausible that our witness reporting algorithm extends to the $\tilde{O}(n^{1.5})$ algorithm.. $◀$ Let us also now show a lower bound, which will nicely complement Algorithm 1.

The proof is an adaptation of the reduction from min-plus convolution to Pareto sum given in [13]. Consider an arbitrary min-plus convolution instance $A[0,\mathrm{\dots },n-1],B[0,\mathrm{\dots },n-1]$ with monotonicity (specifically, $A[i]>A[i+1],B[i]>B[i+1]$ for all $i$) and entries in $\{0,\mathrm{\dots },O(n)\}$. We construct the sets of points

Consider a $\mathrm{\Delta }$-approximation $\tilde{S}$ of the Pareto sum of $P,Q$ with $\mathrm{\Delta }=n^{\frac{\alpha }{1-\alpha }}$ . By construction of the point sets, if a point $\tilde{s}\in \tilde{S}\subseteq P+Q$ $\mathrm{\Delta }$-dominates a point $p+q$, i.e., $\tilde{s}\le p+q+(\mathrm{\Delta },\mathrm{\Delta })$, then it already dominates $p+q$, i.e., $\tilde{s}\le p+q$. Thus, $\tilde{S}$ is the true Pareto sum of $P,Q$, and we can extract the min-plus convolution $C$ of $A,B$: By construction, for every $k\in [2n-2]$, there must be a point in $PS(P,Q)$ with $x$-coordinate $k\cdot (2n^{\frac{\alpha }{1-\alpha }})$ and $y$-coordinate $({\min }_{i\in [n-1]}A[i]+B[k-i])\cdot (2n^{\frac{\alpha }{1-\alpha }})$, revealing $C[k]={\min }_{i\in [n-1]}(A[i]+B[k-i])$. The existence of a point with $x$-coordinate $k\cdot (2n^{\frac{\alpha }{1-\alpha }})$ for every $k\in [2n-2]$ follows by the decreasing values of $C[k]$. Also, witnesses for such points can be transferred.

Note that the constructed instance of Pareto sum Approximation consists of $n$ points in $P,Q$ with coordinates bounded by $W\le O(n^{1+\frac{\alpha }{1-\alpha }})=O(n^{\frac{1}{1-\alpha }})$ and $\mathrm{\Delta }=n^{\frac{\alpha }{1-\alpha }}$. Thus, if $T_{apx}(n,W,\mathrm{\Delta })\le O(n+{(W/\mathrm{\Delta })}^c)=O(n+{(n^{\frac{1}{1-\alpha }-\frac{\alpha }{1-\alpha }})}^c)$, then we can compute the min-plus convolution of $A,B$ in time $O(n^c)$, with witnesses. $◀$

We use three types of benchmarks in our evaluation. Two are synthetic to allow us to control instance properties, while the third benchmark consists of real data.

Benchmark 1: Instances with bounded range. First, we use an instance generator that allows us to fix $n$ and $W$. Here, it is important to ensure $W\ge n-1$ to guarantee that there are enough distinct natural numbers in the range to construct a Pareto set. We use $W=c\cdot n$ for $c\ge 1$. To efficiently create a strictly monotonically increasing sequence, we initialize an array with a counter ${\gamma }_i$ for each $i\in [W]$. Then, we increment a randomly selected counter in each of $n$ rounds. Afterwards, we sweep over this array, starting at position 0, and perform the following operations: If ${\gamma }_i>1$, we set ${\gamma }_i=1$ and ${\gamma }_{i+1}={\gamma }_{i+1}+{\gamma }_i-1$. Here, entry 0 is assumed to be the successor of ${\gamma }_W$. We stop when all entries are 0 or 1 which happens after each entry has been considered at most twice. Then, for all ${\gamma }_i=1$, element $i$ is added to the sequence. This produces a valid sequence in time $𝒪(n+W)$. Pareto sets can then be constructed by using two such sequences, one sorted increasingly for the $x$-coordinates and the other one decreasingly for the corresponding $y$-coordinates.

Benchmark 2: Function-based instance generator. For any strictly monotonically decreasing function $f$ defined on ${ℝ}^{+}$, sampling $n$ points $(x,f(x))$ provides us with a Pareto set. We use two base functions: the first uses $f(x)=-x$ and the second $f(x)=\frac{c}{x}$ with a constant $c\in {ℝ}^{+}$. We slightly perturb the points to have more input variability, while guaranteeing that pairwise non-dominance is upheld. We refer to the resulting instances as near-linear or near-curved.

See Figure 5 for an illustration. We note that by means of the perturbation, the resulting point sets are no longer convex.

Benchmark 3: Route planning instances. Finally, we also use real-world data sets for evaluation. The instances stem from bi-criteria route planning in line with [13]. Here, given a graph in which the edges are augmented with two cost values (in our case distance in meters and positive height difference in centimeters), the goal is to compute all Pareto-optimal paths from a source node $s$ to a target node $t$. Using bi-directional search to answer such queries, or constructing a data structure which accelerates query answering [23], demands to combine and filter the set of Pareto-optimal paths from $s$ to some intermediate node $v$, and from $v$ to the target $t$. The cost pairs of these two path sets form our inputs $P$ and $Q$. The benchmark contains roughly 100000 instances based on randomly chosen $s,v$ and $t$ triples in a directed road network with one million nodes and two million edges. Note that the sizes of $P$ and $Q$ might differ significantly in each instance.

For a general comparison of the different algorithms for Bounded Pareto sum discussed in this paper, we first applied a scalability study on the range-bound instances with all point coordinates being integers, see Figure 6.

We observe huge performance differences. SSS is by far the slowest approach. Its theoretical performance is only subquadratic in case the output size is in $o(n)$. However, for the considered instances there are roughly $4n$ output points. The BSC algorithm performs best and also significantly outperforms the plain SC algorithm. Even more remarkably, also the min-plus convolution based algorithms, naive and CP, outperform SC. While CDXZ beats SSS, it is a bit slower than SC on this instance type. As described in Section 4, our implementation of the CDXZ algorithm allows us to choose its parameters. We conducted a series of tests on synthetic instances of different sizes to determine a good parameter choice, see Figure 7 for an example.

We observe that the parameter values heavily impact the performance. Across all tested instances, scaling factors of $25$ and $p=2$ were among the top configurations. We therefore use those values in all subsequent experiments. We remark, though, that additional fine-tuning or choices adaptive to $n$ could potentially further improve the performance.

However, our main interest is in computing approximate solutions. In the following, we investigate how the algorithms perform when used as a plug-in in our approximation framework. In the remainder of the experiments, we focus on the BSC algorithm as well as the convolution based algorithms, as they outperform the baseline algorithms.

Next, we present comparative running time results of our approximation algorithm with different inner methods of Pareto sum computation applied to the scaled down inputs. For fair comparison, we use the non-witness reporting variants for all of them. Figures 8, 9, and 10 depict the results for range-bounded, near-linear and near-curve instances, respectively, using a value range of $2n$ and a scaling factor of $t=10$. Interestingly, for each input type, another algorithm performs best. For range-bounded instances, BSC is the fastest with CDXZ being second, while on near-linear instances CDXZ outperforms the other two methods. On near-curved instances, however, CP is clearly the best, beating the others by at least an order of magnitude.

The performance of CP relies on its ability to prune pairs, which on near-curved instances applies to up to 95% of total point pairs, see Figure 11.

Finally, we consider our third benchmark set, which consists of bi-criteria route planning instances of varying size. Figure 13 shows the running time and approximation quality of our algorithms for different scaling factors. We observe that BSC performs best with CP being a close second, especially for larger values of $t$. They both improve on exact computation by over an order of magnitude and also stay well below the theoretical additive approximation bound of $2t$.

Examples of approximate solutions can be found in Figures 14.

One can see that even when choosing a larger $t$, the Pareto sum is still well represented by the computed solution. Indeed, larger values of $t$ are quite reasonable in our application. For example, a route that is up to 100 meters longer than another one is almost indistinguishable to a user if the total route length is anyway in the orders of kilometers.

Furthermore, the user usually does not want to be presented with hundreds of routing options, but a concise and diverse set of reasonable options. This is achievable with our approach, as we can not only reduce the running time with larger $t$ but also the output size, see Figures 15 and 12. Clearly, our algorithm provides good trade-offs between running time and solution quality, which allows one to choose the best $t$ depending on the application.