On Constrained Boolean Pareto Optimization

1 Introduction

Optimization tasks usually come with constraints [Bertsekas, 1999], which must be satisfied by the final solutions. For general constrained optimization, i.e., we don't make assumptions on its objective function and constraints, the most widely used method is the penalty function method [Ben Hadj-Alouane and Bean, 1997; Coello Coello, 2002]. This method forms an unconstrained optimization problem by adding a penalty term to the objective function of the original problem, where the penalty term usually measures the degree of the constraints violation. When the penalty term is well balanced with the objective function, the unconstrained optimization problem will yield the same optimal solutions as the original constrained optimization problem.

Pareto optimization [Coello Coello, 2002; Cai and Wang, 2006] solves a constrained optimization problem in a new way: it treats the constraints violation degree as an additional objective, and reformulates the original constrained optimization problem as a bi-objective optimization problem, which is to optimize the original objective function and to minimize the violation degree simultaneously; the bi-objective optimization problem is then solved by a multi-objective evolutionary algorithm; and finally the produced solutions are transformed back to the original problem. Previously, the performance of this approach was only tested in experiments on some benchmark problems [Venkatraman and Yen, 2005; Cai and Wang, 2006]. Recently it has been applied to the ensemble pruning problem [Qian et al., 2015], where Pareto optimization is shown to be superior to state-of-the-art methods both theoretically and empirically [Qian et al., 2015]. However, the theoretical understanding of Pareto optimization is still quite insufficient. It is also important to theoretically understand how widely Pareto optimization is applicable.

∗This research was supported by NSFC (61375061, 61333014) and JiangsuSF (BK2012303).

This paper tries to shed some light on the power of Pareto optimization. We theoretically compare the efficiency of Pareto optimization with a penalty function method on two large Boolean optimization problem classes. The efficiency is measured by the expected running time for achieving an optimal or approximate solution. Given a problem class F, we compare the worst case running time [Cormen et al., 2001]: denoting POM(f) and PFM(f) as the expected running time of the Pareto optimization and the penalty function method for solving f ∈ F, respectively, we compare POM(F) = max{POM(f)|f ∈ F} with PFM(F) = max{PFM(f)|f ∈ F}. Our main results are:

• On the P-solvable minimum matroid optimization problem (denoted as M) for finding an optimal solution, POM(M) ∈ O(rn(log n + log wmax + r)) (Theorem 1) and P FM(M) ∈ Ω(r2n(log n+log wmax)) (Theorem 2), where n, wmax and r are parameters of problem size, maximum weight and matroid ranking, respectively. Thus, Pareto optimization is faster than the penalty function method by at least a factor of min{log n + log wmax, r}.
• On the NP-hard minimum cost coverage problem (denoted as C) for finding an approximate solution, POM(C) ∈ O(Nn(log n+log wmax+N)) (Theorem 3) and PFM(C) is at least exponential w.r.t. n, N and logwmax (Theorem 4), where N is a submodular function parameter. Thus, Pareto optimization is exponentially faster than the penalty function method.
• On a specific minimum cost coverage instance C ∈ C, the greedy algorithm finds a local optimal solution (Proposition 1), while POM(C) ∈ O(n2 logn) (Proposition 2) and PFM(C) ∈ O(n3 logn) (Proposition 3) for finding

the global optimal solution. Thus, Pareto optimization can be better than the greedy algorithm.

The rest of the paper starts with a section of preliminaries. Section 3 and 4 present the studies on minimum matroid optimization and minimum cost coverage, respectively. We further compare Pareto optimization with a greedy algorithm in section 5. Section 6 concludes the paper.

2 Preliminaries

2.1 Constrained Optimization

Constrained optimization [Bertsekas, 1999] is to optimize an objective function with constraints that must be satisfied. It can be formally stated as follows.

Definition 1 (Constrained Optimization)

where f(x) is the objective function, gi(x) and hi(x) are the equality and inequality constraints, respectively.

Note that we consider binary spaces {0, 1}n in this paper. We also only consider minimization since maximizing f is equivalent to minimizing -f. A solution is (in)feasible if it does (not) satisfy the constraints. Thus, constrained optimization is to find a feasible solution minimizing the objective f .

2.2 The Penalty FunctionMethod

The penalty function method [Ben Hadj-Alouane and Bean, 1997; Coello Coello, 2002] is the most widely employed method to solve constrained optimization. For the optimization problem in Eq. (1), the penalty function method turns to solve an unconstrained optimization problem:

where f is the objective in Eq. (1), λ is the penalty coefficient, and fi is the penalty term expressing a violation degree of the i-th constraint. A straightforward way to set fi is

For the penalty coefficient λ, a static value has been shown to be simple and robust [Homaifar et al., 1994; Michalewicz and Schoenauer, 1996]. Furthermore, Zhou and He [2007] proved that a large enough static λ value makes the penalty function method equivalent to the "superiority of feasible points" strategy [Deb, 2000] that first prefers a smaller constraints violation degree, and then compares the objective value if having the same violation degree. Such a strategy ensures that a feasible solution is always better than an infeasible one, and thus the optimal solutions must be feasible and are the same as that of the constrained optimization problem.

Then, an unconstrained optimization algorithm will be employed to solve this unconstrained problem. For tackling various constrained optimization problems, we consider the context of general purpose optimization algorithms, such as randomized local search and evolutionary algorithms [B¨ack, 1996]. The penalty function method with a general purpose optimization algorithm [He and Yao, 2001; Yu and Zhou, 2008] is presented in Algorithm 1.

Algorithm 1 (The Penalty Function Method) Given a constrained optimization problem as in Eq. (1), it contains:

1. Let h(x)=f(x)+λ∑mi=1 fi(x) according to Eq. (2).
2. x = selected from {0, 1}n uniformly at random.
3. repeat until the termination condition is met
4. x′ =flip each bit of x independently with prob. 1n .
5. if h(x′ ) ≤ h(x)
6. x = x′ .
7. return x

Sharing a common structure with other general purpose optimization algorithms, it starts from a random solution, generates a new solution from the current solution by some heuristic, decides if the new solution is good enough to be kept, and repeats the generation and decision until a good solution is found. It employs a typical heuristic for generating new solutions in evolutionary algorithms, that flips each bit (0 to 1 or inverse) of a solution with probability 1n . This algorithm does not make any assumption on the optimization problem, but only require that two solutions can be compared for the goodness (i.e., step 5). Thus it can be widely applied.

2.3 The Pareto OptimizationMethod

By Pareto optimization [Coello Coello, 2002; Cai and Wang, 2006], the constrained optimization problem is reformulated into a multi-objective (bi-objective) problem:

where fi is the constraint violation degree and can be set as Eq. (2). In this bi-objective formulation, a solution has an objective vector instead of a scalar objective value. Unlike single-objective optimization, the objective vector makes the comparison between two solutions not straightforward, because it is possible that one solution is better on the first dimension while the other is better on the second dimension. For this situation, the domination relationship between solutions is usually used, which is introduced in Definition 2 for the bi-objective (i.e., two objectives) minimization case.

Definition 2 (Domination) Let g = (g1, g2) : X → R2 be the objective vector. For two solutions x, x′ ∈ X: (1) x weakly dominates x′ if g1(x) ≤ g1(x′ ) and g2(x) ≤ g2(x′ ), denoted as ⪰g; (2) x dominates x′ if x ⪰g x′ and either g1(x) < g1(x′ ) or g2(x) < g2(x′ ), denoted as ≻g.

Thus, a bi-objective optimization problem may not have a single optimal solution, but instead have a set of Pareto optimal solutions. A solution x is Pareto optimal if there is no other solution in X that dominates x. The set of objective vectors of all the Pareto optimal solutions is called the Pareto front.

For the fairness of comparison, a general purpose biobjective optimization algorithm [Laumanns et al., 2004;

Qian et al., 2013] with the same initialization and random perturbation as Algorithm 1 is then employed to solve this bi-objective problem, presented in Algorithm 2.

Algorithm 2 (The Pareto OptimizationMethod) Given a constrained optimization problem as in Eq. (1), it contains:

1. Let g(x) = (f(x), ∑mi=1 fi(x)). 2. x = selected from {0, 1}n uniformly at random. 3. P = {x}. 4. repeat until the termination condition is met 5. x = selected from P uniformly at random. 6. x′ =flip each bit of x independently with prob. 1n . 7. if ∄z ∈ P such that z ≻g x′ 8. P = (P - {z ∈ P |x′ ⪰g z}) ∪ {x′ }. 9. return x ∈ P with ∑mi=1 fi(x) = 0

It firstly generates a random solution, and puts it into the candidate solution set P ; and then follows a cycle to improve the solutions in P iteratively. In each iteration, a solution x is randomly selected from P; and is then perturbed to generate a new solution x′ ; if x′ is not dominated by any solution in P, x′ is added into P and at the same time solutions in P that are weakly dominated by x′ get removed. Finally, the solution in P that violates no constraints is selected.

Note that Pareto optimization here is different from traditional multi-objective optimization [Deb et al., 2002]. The latter aims at finding a uniform approximation of the Pareto front, while Pareto optimization aims at finding the optimal solutions for the original constrained optimization problem.

2.4 Running Time Analysis Comparison

To measure the performance of an algorithm solving a problem, we employ the expected running time complexity, which counts the number of objective function evaluations (the most time-consuming step) until finding a desired solution, since it has been a leading theoretical aspect for randomized search heuristics [Neumann and Witt, 2010; Auger and Doerr, 2011]. Because both the two analyzed methods perform only one objective evaluation (i.e., h(x′ ) and g(x′ )) in each iteration, the expected running time equals the average number of iterations. For P-solvable problems, we use exact analysis [He and Yao, 2001], i.e., the running time until an optimal solution is found. For NP-hard problems, we use approximate analysis [Yu et al., 2012], i.e., the running time until finding an αapproximate solution having the objective value ≤ α · OPT , where α > 1 and OPT is the optimal function value.

3 OnMinimumMatroid Optimization

We first analyze the P-solvable minimum matroid optimization problem [Edmonds, 1971], which covers several combinatorial optimization problems, e.g., minimum spanning tree.

Let | · | denote the size (i.e., cardinality) of a set. A matroid is a pair (U, S), where U is a finite set and S ⊆ 2U , satisfying

1. ∅ ∈ S; (2) ∀A ⊆ B ∈ S,A ∈ S;
2. ∀A,B ∈ S, |A| > |B| : ∃e ∈ A-B,B ∪ {e} ∈ S.

The elements of S are called independent. For any A ⊆ U, a maximal independent subset of A is called a basis of A;

the rank of A is the maximal cardinality of a basis of A, i.e., r(A) = max{|B| | B ⊆ A, B ∈ S}. Note that for a matroid, all bases of A have the same cardinality.

Given a matroid (U, S) and a weight function w : U → N where N denotes the positive integer set, the minimum matroid optimization problem is to find a basis of U with the minimum weight. Let U = {e1, e2, . . . , en} and wi = w(ei). Let x ∈ {0, 1}n represent a subset of U, where xi = 1 means that ei is selected. For notational convenience, we will not distinguish x and its corresponding subset {ei ∈ U|xi = 1}. The problem then can be formally stated as follows.

Definition 3 (MinimumMatroid Optimization) Given a matroid (U, S) and a weight function w : U → N, it is to find a solution such that

Note that, a rank oracle that computes the rank of a subset is polynomially equivalent to an independence oracle that tells whether a subset is independent [Korte and Vygen, 2012].

For Pareto optimization, noting r(U) ≥ r(x), the objective vector g is implemented as

Theorem 1 proves the running time bound, where wmax = max{wi | 1 ≤ i ≤ n} is the maximum element weight in U.

Theorem 1 For minimum matroid optimization, the expected running time of the Pareto optimization method for finding an optimal solution is O(r(U)n(log n + log wmax + r(U))).

The proof idea is to divide the optimization process into two phases: (1) starts after initialization and finishes until finding the special solution {0}n ; (2) starts from {0}n and follows the greedy behavior [Edmonds, 1971] to find an optimal solution. The running time of phase (1) shown in Lemma 2 is derived by applying multiplicative drift analysis (i.e., Lemma 1), a recently proposed approach for analyzing the hitting time of a random process.

Lemma 1 (Multiplicative Drift Analysis) [Doerr et al., 2012] Let S ⊆ R+ be a finite set of positive numbers with minimum smin. Let {Xt}t∈N be a sequence of random variables over S∪{0}. Let T be the random variable that denotes the first point in time t ∈ N for which Xt = 0. Suppose that there exists a real number δ > 0 such that

holds for all s ∈ S. Then for all s0 ∈ S, we have

Lemma 2 For minimum matroid optimization, the expected running time of the Pareto optimization method for finding {0}n is O(r(U)n(log n + log wmax)).

Proof. We prove it by using Lemma 1. Let Xt = min{w(x) | x ∈ P} (i.e., the minimum weight of solutions in the candidate solution set P) after t iterations of Algorithm 2. Xt = 0 implies that the solution {0}n has been found, and thus T is the running time for finding {0}n .

Then, we are to analyze E[[Xt -Xt+1|Xt]]. Let x be the corresponding solution with w(x) = Xt. It is easy to see that Xt cannot increase (i.e., Xt+1 ≤ Xt) because a newly generated solution x′ with a larger weight cannot weakly dominate x. Let Pmax denote the largest size of P during the optimization. In the (t + 1)-th iteration, we consider that x is selected in step 5 of Algorithm 2, which happens with probability at least 1 Pmax . Then, in step 6, we denote xi as the solution generated by flipping the i-th bit (i.e., xi) of x and keeping other bits unchanged, which happens with probability 1n (1 - 1n )n-1 ≥ 1 en . If xi = 1, the generated solution xi will be included into P since it has the smallest weight now and Xt+1 will decrease to Xt - wi. Thus, we have

This implies that E[[Xt -Xt+1|Xt]] ≥ Xt/(enPmax).

The procedure of Algorithm 2 ensures that the solutions maintained in P must be incomparable, which implies that there exists at most one solution in P for each value of one objective. Since the objective r(U) - r(x) ∈ {0, 1, . . . , r(U)}, Pmax ≤ r(U) + 1. Note that X0 ≤ nwmax and smin ≥ 1 by wi ∈ N. Thus, by Lemma 1, for any X0,

Proof of Theorem 1. We analyze phase (2) after finding {0}n . Let wk = min{w(x) | r(x) = k} denote the smallest weight in all the solutions with rank k. It is easy to see that wk increases with k. Then the Pareto front of g is {(wk , r(U) - k) | 0 ≤ k ≤ r(U)}. We are to show that staring from {0}n with (w0 = 0, r(U)), the algorithm can find the solutions corresponding to the Pareto front from k = 0 to r(U) sequentially, where the last solution with (wr(U) , 0) is just the the optimal solution of Eq. (3) we are to find.

We first prove that for any x corresponding to (wk , r(U)- k), adding the lightest element that makes r(x) increase by 1 will produce a solution x′ corresponding to (wk+1 , r(U) - k - 1). Assume that x = {ei1 , . . . , eik }, x′ = x ∪ {eik+1 }, and there exists another solution ˆx = {ej1 , . . . , ejk , ejk+1 } with rank k+1 and w(ˆx) < w(x′ ). Assume that wj1 ≤ . . . ≤ wjk+1 . If wjk+1 ≥ wik+1 , w(ˆx)-wjk+1 < w(x′ )-wik+1 = w(x), which contradicts with that x has the smallest weight in all the solutions with rank k. If wjk+1 < wik+1 , then ∃e ∈ ˆx - x, x ∪ {e} has rank k + 1 and w(e) ≤ wjk+1 < wik+1 , which contradicts with that eik+1 is the lightest element that increases r(x) by 1. Thus, our claim holds.

Assume that {(wk , r(U) - k) | 0 ≤ k ≤ i} has been found. This is well defined because i ≥ 0 after finding {0}n .

Based on the above analysis, for finding (wi+1 , r(U)-i-1) in one iteration, it is sufficient to select the solution corresponding to (wi , r(U) - i) in step 5 of Algorithm 2 and flip only the lightest element that increases its rank by 1 in step 6, which happens with probability at least 1 Pmax · 1n (1- 1n )n-1 ≥ 1 en(r(U)+1) . Then, the expected running time is O(r(U)n) for finding (wi+1 , r(U) - i - 1). Since r(U) such processes are sufficient to find the optimal solution with (wr(U) , 0), the expected running time of phase (2) is O(r2 (U)n).

By combining the two phases, the total expected running time is O(r(U)n(log n + log wmax + r(U))). □

For the penalty function method, we first show that minimum spanning tree (MST) is an instance of minimum matroid optimization and then give a concrete MST example where the penalty function method needs more running time than Pareto optimization.

Given an undirected connected graph G = (V,E) on m vertices and n edges with weights w : E → N, the MST problem is to find a connected subgraph with the minimum weight. Let E = {e1, e2, . . . , en}. Let x ∈ {0, 1}n represent a subgraph, where xi = 1 means that the edge ei is selected; let c(x) denote its connected components. The MST problem can be formally stated as follows.

Definition 4 (Minimum Spanning Tree) Given an undirected connected graph G = (V,E) with a weight function w : E → N, it is to find a solution such that

Let U = E and S = {subgraphs without cycles}. Then, (U, S) is a matroid, c(x)+r(x) = m and r(U) = m-1. It is easy to verify that MST is a minimum matroid optimization instance (i.e., Definition 3) with such a configuration.

Neumann and Wegener [2007] have analyzed the running time of the penalty function method on a specific MST example, as shown in Lemma 3 (i.e., Theorem 6 in their paper). Note that the parameters m and n in [Neumann and Wegener, 2007] are opposite to ours, and we have exchanged them here for consistency. Thus, we have Theorem 2 by r(U) = m-1.

Lemma 3 [Neumann andWegener, 2007] The expected running time until the penalty function method finds a minimum spanning tree for an example graph (n ∈ Θ(m2 ), wmax = 3m2 ) equals Θ(n2 logm).

Theorem 2 There exists a minimum matroid optimization instance, where the expected running time of the penalty function method for finding an optimal solution is Θ(r2 (U )n(log n + log wmax)).

4 OnMinimumCost Coverage

We then analyze the minimum cost coverage problem, which is NP-hard. A representative instance is the submodular set cover problem [Wolsey, 1982], which arises in many applications, e.g., social networks [Kempe et al., 2003] and sensor placement [Krause et al., 2008].

Let U = {e1, e2, . . . , en} be a finite set. A set function f : 2U → R is monotone and submodular iff ∀A,B ⊆ U, f(A) ≤ f(B)+ ∑e∈A-B (f(B∪{e})-f(B)) [Nemhauser et al., 1978]. Let x ∈ {0, 1}n represent a subset of U. Minimum cost coverage can be formally stated as follows.

Definition 5 (Minimum Cost Coverage) Given a monotone submodular function f : 2U →R, some value q ≤ f(U) and a weight function w : U → N, it is to find a solution such that

We consider f to be non-negative. Since the rank of a matroid is a monotone submodular function, minimum matroid optimization is actually a minimum cost coverage instance with q = r(U). We analyze them separately because minimum matroid optimization is P-solvable and minimum cost coverage is NP-hard in general. We use exact and approximate analysis, respectively.

For minimum cost coverage by Pareto optimization, the objective vector g is implemented as

The running time bound is proved in Theorem 3. Let Hi = ∑ij=1 1/j be the i-th harmonic number, and let c denote the minimum real number making c · f(x) for all x and c · q to be integer. Let N denote the number of distinct f values in [0, q). The proof idea is similar to that for Theorem 1 except that it follows the greedy behavior [Wolsey, 1982] to find an approximate solution.

Theorem 3 For minimum cost coverage, the expected running time of the Pareto optimization method for finding a Hcqapproximate solution is O(Nn(log n + log wmax + N)).

Proof. Because the objective max{0, q - f(x)} has N + 1 number of distinct values and the solutions in P are incomparable, the largest size of P is not larger than N +1. By using the same proof as Lemma 2 except Pmax ≤ N + 1 here, we can derive that the expected running time for finding {0}n is O(Nn(log n + log wmax)).

Let Rk = Hcq - Hc(q-k). Let x∗ and OPT denote an optimal solution and the optimal function value of Eq. (4), respectively. Then, w(x∗ ) = OPT. Let Kmax denote the maximum value of k such that there exists one solution x in P with min{q, f(x)} = k and w(x) ≤ Rk · OPT . That is, Kmax = max{k | ∃x ∈ P,min{q, f(x)} = k ∧ w(x) ≤ Rk · OPT}. Then, we are to analyze the expected running time until Kmax = q, which implies that it finds a Rq (i.e., Hcq)-approximate solution.

After finding {0}n , Kmax ≥ 0. Assume that currently Kmax = k < q. Let x denote the corresponding solution with the value k, i.e., min{q, f(x)} = k and w(x) ≤ Rk · OPT . We are first to show that Kmax cannot decrease. If x is kept in P, Kmax obviously will not decrease. If x is deleted, by step 7 and 8 of Algorithm 2, a newly generated solution x′ weakly dominating x (i.e., x′ is not worse on both the two objectives) must be included. Note that max{0, q f(x)} = q -min{q, f(x)}. Thus, min{q, f(x′ )} = k′ ≥ min{q, f(x)} = k and w(x′ ) ≤ w(x). Because Rk increases with k, w(x′ ) ≤ Rk · OPT ≤ Rk′ · OPT. Thus, Kmax ≥ k′ , i.e., Kmax will not decrease.

Then, we are to show that Kmax can increase by flipping a specific 0 bit of x. Let f ′ (x) = min{q, f(x)}. Let xi denote the solution generated by flipping the i-th bit of x. Let I = {i ∈ [1, n] | f ′ (xi ) - f ′ (x) > 0} denote the 0 bit positions of x where the flipping can generate a solution with a positive increment on f ′ . Let δ =min{ wi f ′ (xi )-f ′ (x) | i ∈ I}. Then, δ ≤ OPT/(q - k). Otherwise, for any ei ∈ x∗ - x, wi > (f ′ (xi )-f ′ (x))·OPT/(q-k). Thus, ∑ei∈x∗-x wi > (∑ei∈x∗-x(f ′ (xi ) - f ′ (x))) · OPT/(q - k). Since f is monotone and submodular, f ′ is also monotone and submodular, then f ′ (x∗ ) - f ′ (x) ≤ ∑ei∈x∗-x(f ′ (xi ) - f ′ (x)). Thus, ∑ei∈x∗-x wi > (f ′ (x∗ ) - f ′ (x)) · OPT/(q - k) = OPT, which contradicts with ∑ei∈x∗-x wi ≤ w(x∗ ) = OPT. Thus, by selecting x in step 5 of Algorithm 2 and flipping only the 0 bit corresponding to δ in step 6, it can generate a new solution x′ with min{q, f(x′ )} = k′ > k and

Once generated, x′ will be included into P. Otherwise, there must exist a solution in P dominating x′ which has a larger f ′ and a smaller w; this implies that Kmax has already been larger than k, which contradicts with the assumption Kmax = k. After including x′ , Kmax increases from k to k′ .

The probability of flipping a specific 0 bit of x is at least 1 N+1 · 1n (1- 1n )n-1 ≥ 1 en(N+1) . Thus, the expected running time for such a step of increasing Kmax is at most en(N+1). Since N such steps are sufficient to make Kmax = q, the expected running time of this phase is O(N 2n).

By combining the two phases, the expected running time of the whole process is O(Nn(logn+ logwmax +N)). □

For the penalty function method, we first show that minimum set cover (MSC) is an instance of minimum cost coverage and then give a concrete MSC example where the penalty function method is less efficient than Pareto optimization.

Definition 6 (Minimum Set Cover) Given a set S = {e1, . . . , em}, a collection C = {C1, . . . , Cn} of subsets of S with corresponding costs w : C → N, it is to find a subset of C (represented by x ∈ {0, 1}n ) with the minimum cost such that all the elements of S are covered, that is

Let U = C and f(x) = | ⋃i:xi=1 Ci|. Then f is monotone and submodular. Since f(x) ≤ m and |S| = m, it is easy to verify that the MSC problem is an instance of minimum cost coverage (i.e., Definition 5) with q = m. Obviously, the parameters c ≤ 1 and N ≤ m for MSC.

Friedrich et al. [2010] have analyzed the running time of the penalty function method on a specific MSC example, as shown in Lemma 4 (i.e., Theorem 8 in their paper). Thus, we have Theorem 4 by letting ϵ being a constant and wmax = 2n .

Lemma 4 [Friedrich et al., 2010] Let δ > 0 be a constant and nδ-1 ≤ ϵ < 1/2. The expected running time of the penalty function method on a MSC example (m = ϵ(1-ϵ)n2 ) for finding an approximation better than ((1 - ϵ)wmax)/ϵ is exponential w.r.t. n.

Theorem 4 There exists a minimum cost coverage instance, where the expected running time of the penalty function method for finding a Hcq-approximate solution is exponential w.r.t. n, N and logwmax.

5 Pareto Optimization vs. Greedy Algorithm

For the above two problems, there exist greedy algorithms which efficiently find an optimal solution [Edmonds, 1971] and a Hcq-approximate solution [Wolsey, 1982], respectively. Thus, a natural question is whether Pareto optimization can be better than greedy algorithms. We give a positive answer by analyzing a concrete instance of minimum set cover (also minimum cost coverage).

Example 1 [Yu et al., 2012] As shown in Figure 1, the ground set S contains m = k(k - 1) elements, and the collection C contains n = (k + 1)(k - 1) subsets {S∗ 1 , . . . , S∗ k-1 , S1,1, . . . , S1,k, . . . , Sk-1,1, . . . , Sk-1,k} with costs ∀i : w(S∗ i ) = 1 + 1/k2 and ∀i, j : w(Si,j) = 1/j.

For Example 1, the global optimal solution is Cglobal = {S∗ 1 , . . . , S∗ k-1} with cost (k - 1)(1 + 1/k2 ), and Clocal = {S1,1, . . . , S1,k, . . . , Sk-1,1, . . . , Sk-1,k} is a local optimal solution with cost (k - 1)Hk. We show the performance of the greedy algorithm [Chvatal, 1979], the Pareto optimization method and the penalty function method on this example in Propositions 1-3, respectively. The greedy algorithm iteratively selects a subset from C with the smallest ratio of its cost and the number of newly covered elements, and its performance is directly from Proposition 1 in [Yu et al., 2012].

Proposition 1 [Yu et al., 2012] For Example 1, the greedy algorithm [Chvatal, 1979] finds Clocal.

Proposition 2 For Example 1, the Pareto optimization method finds Cglobal in O(n2 log n) expected running time.

Proof. The objective vector g is implemented as (w(x),m| ⋃i:xi=1 Ci|). By using the same proof as Lemma 2 except that wmax = 1 + 1/k2 and Pmax ≤ m + 1 in this example, we can derive an upper bound O(mn log n) = O(n2 log n) on the expected running time for finding {0}n .

Then, we are to show that starting from {0}n , it can find the subset of Cglobal = {S∗ 1 , . . . , S∗ k-1} with size from 0 to k - 1 sequentially. Note that any subset of Cglobal is Pareto optimal, since the cost of covering k elements using subsets only from {Si,j} is larger than 1 + 1 k2 . After finding x ⊆ Cglobal with |x| = i, it can generate x′ ⊆ Cglobal with |x′ | = i + 1 in one iteration with probability at least k-1-i (m+1)n (1- 1n )n-1 ≥ k-1-i e(m+1)n , since it is sufficient to select x and flip one 0 bit corresponding to any unselected S∗ i . Thus, the expected running time to find the optimal solution Cglobal is at most ∑k-2 i=0 e(m+1)n k-1-i ∈ O(n2 logn). □

Proposition 3 For Example 1, the penalty function method finds Cglobal in O(n3 log n) expected running time.

Proof. Let u denote the number of uncovered elements of the current solution. Obviously, u will not increase. Since each uncovered element has a corresponding Si,j, flipping the corresponding bit of Si,j can decrease u by 1. Thus, u can decrease by 1 in one iteration with probability at least un (1 - 1n )n-1 ≥ u en . Hence, for finding a set cover (i.e., u decreases to 0), the expected running time is at most ∑1u=m en u ∈ O(n logn). Note that after finding a set cover, the solution will always keep being a set cover.

Then, we apply Lemma 1 to derive the running time for finding Cglobal. Let Xt = w(xt ) - wopt, where xt is the solution after t iterations of Algorithm 2 and wopt is the optimal weight, i.e., w(Cglobal). Then, Xt = 0 implies that Cglobal is found, and T is just the running time for finding Cglobal.

We are to analyze E[[Xt -Xt+1|Xt]]. It is easy to see that w(xt+1 ) ≤ w(xt ). For xt , there are three covering cases for each column of Figure 1. For case (1) which contains S∗ i and several Si,j , by deleting one Si,j , w(xt+1 ) = w(xt ) - 1/j with probability 1n (1 - 1n )n-1 . For case (2) which contains all Si,j but not S∗ i , by including S∗ i and deleting Si,1 and another Si,j , w(xt+1 ) = w(xt ) - (1/j - 1/k2 ) with probability 1 n3 (1- 1n )n-3 ; by including S∗ i and deleting Si,1 and any two from {Si,2, . . . , Si,k}, w(xt+1 ) < w(xt ) - (k - 2)/k2 with probability (k-1 2 ) n4 (1 - 1n )n-4 . For case (3) which contains only S∗ i , it reaches the optimal case. By unifying these probabilities to a lower bound k-2 2e(k+1)n3 , the sum of the improvements is just w(xt ) -wopt. Thus,

Then, we have E[[Xt - Xt+1|Xt]]

6 Conclusion

Pareto optimization employs a bi-objective optimization as an intermediate step for solving constrained optimization problems. Recent application [Qian et al., 2015] has shown its theoretical and empirical advantages. This work further theoretically investigates its optimization ability. Pareto optimization is compared with the widely used penalty function method on two large problem classes, the P-solvable minimum matroid optimization and the NP-hard minimum cost coverage problem. On both problems, Pareto optimization is proved to be more efficient. Moreover, Pareto optimization is shown to be better than the greedy algorithm. With these theoretical results, we expect that Pareto optimization will attract more attentions for constrained optimization tasks.