In this paper, we study parallel algorithms for the problem of maximizing a non-monotone DR-submodular function subject to a single cardinality constraint. The problem is a generalization of submodular maximization subject to a cardinality constraint. Many recent works have shown that DR-submodular maximization has a wide-range of applications beyond submodular maximization. These applications include maximum a-posteriori (MAP) inference for determinantal point processes (DPP), mean-field inference in log-submodular models, quadratic programming, and revenue maximization in social networks (Kulesza et al., 2012;Gillenwater et al., 2012;Bian et al., 2016;Ito & Fujimaki, 2016;Soma & Yoshida, 2017;Bian et al., 2017;2018).
The problem of maximizing a DR-submodular function subject to a convex constraint is a notable example of a non-convex optimization problem that can be solved with provable approximation guarantees. The continuous Greedy algorithm (Vondrák, 2008) developed in the context of the multilinear relaxation framework applies more generally to maximizing DR-submodular functions that are monotone increasing (if x ⃗ ≤ y ⃗ coordinate-wise then f (x ⃗ ) ≤ f (y ⃗)). Feldman et al. (Feldman et al., 2011) developed a variant of continuous Greedy for non-monotone objectives. Chekuri et al. (Chekuri et al., 2015) developed algorithms for both monotone and non-monotone DR-submodular maximization subject to packing constraints that are based on the continuous Greedy and multiplicative weights update framework. The work (Bian et al., 2017) generalized continuous Greedy for submodular functions to the DR-submodular case and developed Frank-Wolfe-style algorithms for maximimizing non-monotone DR-submodular function subject to general convex constraints.
A significant drawback of these algorithms is that they are inherently sequential and adaptive. In fact the highly adaptive nature of these algorithms goes back to the classical greedy algorithm for submodular functions: the algorithm sequentially selects the next element based on the marginal gain on top of previous elements. In certain settings such as feature selection (Khanna et al., 2017) evaluating the objective function is a time-consuming procedure and the main bottleneck of the optimization algorithm and therefore, parallelization is a must. Recent lines of work have focused on addressing these shortcomings and understanding the trade-offs between approximation guarantee, parallelization, and adaptivity. Starting with the work of Balkanski and Singer (Balkanski & Singer, 2018), there have been very recent efforts to understand the tradeoff between approximation guarantee and adaptivity for submodular maximization (Balkanski & Singer, 2018;Ene & Nguyen, 2019;Balkanski et al., 2019;Fahrbach et al., 2019;Chekuri & Quanrud, 2019;Balkanski et al., 2018). The adaptivity of an algorithm is the number of sequential rounds of queries it makes to the evaluation oracle of the function, where in every round the algorithm is allowed to make polynomially-many parallel queries. Recently, the work (Fahrbach et al., 2018) gave an algorithm for maximizing a submodular function subject to a cardinality constraint in O(log n/ϵ) rounds and 0.031ϵ approximation. For the general setting of DR-submodular functions with m packing constraints, the work (Ene et al., 2019) gave an algorithm with O(log(n/ϵ) log(1/ϵ) log(m + n)/ϵ 2 ) rounds and 1/e-ϵ approximation. In the special case of m = 1 constraint, this algorithm uses O(log 2 (n) log(1/ϵ)/ϵ 2 ) rounds.
In this work, we develop a new algorithm for DRsubmodular maximization subject to a single cardinality constraint using O(log n log(1/ϵ)/ϵ 3 ) rounds of adaptivity and obtaining 1/eϵ approximation. For small n compared to 1/ϵ, this is not an improvement over previous work. However, for large n, the number of rounds is almost a quadratic improvement from O(log 2 n) in the previous work to the nearly optimal O(log n) rounds. Our experimental evaluation shows that our algorithm outperforms the state of the art even in the small n regime that is advantageous for the existing algorithms.
Theorem 1. Let f : [0, 1] n → R + be a DR-submodular function and k ∈ R + . For every ϵ > 0, there is an algorithm for the problem max x ⃗ ∈[0,1] n : ∥x ⃗ ∥1≤k f (x ⃗ ) with the following guarantees:
• The algorithm is deterministic if provided oracle access for evaluating f and its gradient ∇f ;
• The algorithm achieves an approximation guarantee of 1 eϵ; • The number of rounds of adaptivity is
)︂ .
Let f : [0, 1] n → R + be a non-negative function. The function is diminishing returns submodular
where 1 ⃗ i is the i-th basis vector, i.e., the vector whose i-th entry is 1 and all other entries are 0.
Throughout the paper, we assume that f is differentiable. We assume that we are given black-box access to an oracle for evaluating f and its gradient ∇f . We extend the function
An example of a DR-submodular function is the multilinear extension of a submodular function g. The multilinear extension is defined as
where R(x ⃗ ) is a random subset of V where each i ∈ V is included independently at random with probability x ⃗ i .
Basic notation. We use e.g. x ⃗ = (x ⃗ 1 , . . . , x ⃗ n ) to denote a vector in R n . We use the following vector operations:
⃗ ) be the n-dimensional all-zeros (resp. all-ones) vector. Let 1 ⃗ S ∈ {0, 1} V denote the indicator vector of S ⊆ V , i.e., the vector that has a 1 in entry i if and only if i ∈ S.
We will use the following result that was shown in previous work (Chekuri et al., 2015).
Lemma 2 ( (Chekuri et al., 2015)
In this section, we present an idealized version of our algorithm where we assume that we can compute exactly the step size on line 16. The idealized algorithm is given in Algorithm 1. In the supplement (Section B), we show how to implement that step efficiently and incur only O(ϵ) additive error in the approximation.
The algorithm takes as input a target value M and it achieves the desired 1/e -O(ϵ) approximation if M is an (1 + ϵ) approximation of the optimal function value f (x ⃗ * ), i.e., we have f
As noted in previous work (Ene et al., 2019), it is straightforward to approximately guess such a value M using a single parallel round.
The algorithm maintains two solutions: the solution x ⃗ that will be returned, and a second solution z ⃗ ≥ x ⃗ . We use z ⃗ to handle the non-monotonicity of the objective. Throughout the algorithm, we have z ⃗ ≥ x ⃗ and thus ∇f (z ⃗) ≤ ∇f (x ⃗ ) by DR-submodularity. Instead of using the gradient at x ⃗ to decide which coordinates to update, we use the gradient at z ⃗. By doing so, we may under-estimate the potential gain of some of the coordinates, but it allows us to deal with the non-monotonicity of the function. This strategy is reminiscent of the gradient lookahead technique of (Ene et al., 2019), although we use it in a very different way.
The algorithm iteratively increases the two solutions over 1/ϵ phases (iterations of the outer for loop). During a single phase, the algorithm uses a descending thresholds approach.
Algorithm 1 Algorithm for max x ⃗ ∈[0,1] n : ∥x ⃗ ∥1≤k f (x ⃗ ), where f is a non-negative DR-submodular function. The algorithm takes as input a target value M such that f
while v > ϵv start and ∥z ⃗∥ 1 < ϵjk do 10:
else 15:
For a given η ∈ [0, ϵ 2 ], we define:
end if
end if
end while 27: end for 28: return x ⃗
Starting with an initial threshold v start , the algorithm only updates coordinates whose gradient is above the threshold.
The algorithm only updates a carefully chosen subset of these high-gain coordinates: the set S on line 11. The definition of S imposes two further restrictions on the fractional value of each coordinate in S. These restrictions allow us to control the ℓ ∞ norm of the solution and cope with the non-monotonicity of the objective function (cf. Lemma 2).
If the set S is empty, the algorithm lowers the threshold by an 1ϵ factor on line 13. Otherwise, the algorithm simultaneously increases a large fraction of the coordinates in S using an appropriately chosen step size η. The step size is carefully chosen to ensure several competing desiderata.
On one hand, we want to take large steps for fast convergence/high parallelism. On the other, we want to take small steps to ensure high function value gain, since the gradient approximates the function value only locally and we need to be conservative on how much we increase due to the non-monotonicity of the objective. We choose our step size so that, after the update, the size of S decreases by an 1ϵ factor or we reach the budget of ϵjk allocated to the phase or we reach our cap of ϵ 2 on the step size. Thus we converge fast by either shrinking S or filling up the budget, while simultaneously achieving high gains in function value.
In this section, we show that Algorithm 1 achieves a 1 e -O(ϵ) approximation. Recall that we assume that η 1 is computed exactly on line 16. In Section B of the supplement, we show how to extend the algorithm and the analysis so that the algorithm efficiently computes a suitable approximation to η 1 that suffices for obtaining a 1 e -O(ϵ) approximation. The omitted proofs can be found in the supplement (Section A).
In the following, we refer to each iteration of the outer for loop as a phase. We refer to each iteration of the inner while loop as an iteration. Note that the update vectors are non-negative in each iteration of the algorithm, and thus the vectors x ⃗ , z ⃗ remain non-negative throughout the algorithm and they can only increase. Additionally, since S(η) ⊆ T (η) ⊆ S, we have x ⃗ ≤ z ⃗ throughout the algorithm. We will also use the following observations repeatedly, whose straightforward proofs are deferred to Section A of the supplement. By DR-submodularity, since the relevant vectors can only increase in each coordinate, the relevant gradients can only decrease in each coordinate. This implies that, for every η ≤ η ′ , we have S(η) ⊇ S(η ′ ). Additionally, for every i ∈ T (η), we have
We will need an upper bound on the ℓ 1 and ℓ ∞ norms of x ⃗ and z ⃗. Since x ⃗ ≤ z ⃗, it suffices to upper bound the norms of z ⃗ (the ℓ 1 norm bound will be used to show that the final solution is feasible, and the ℓ ∞ norm bound will be used to derive the approximation guarantee). The omitted proofs can be found in the supplement.
Lemma 3. Consider phase j of the algorithm (the j-th iteration of the outer for loop). Throughout the phase, the algorithm maintains the invariant that ∥z ⃗∥ ∞ ≤ 1 -(1ϵ) j + ϵ 2 and ∥z ⃗∥ 1 ≤ ϵjk.
The following theorem is the heart of the analysis and it shows that each phase makes sufficient gain in function value. This is reminiscent of the gain made by the measured continuous greedy algorithm in the sequential setting and it implies the claimed approximation guarantee via induction. As discussed in Section 3, in order to achieve high parallelism, our algorithm proceeds very differently from the sequential continuous greedy algorithm. As a result, we need a very careful analysis to show the claimed gain.
Theorem 4. Consider phase j of the algorithm (the j-th iteration of the outer for loop). Let x ⃗ start and x ⃗ end be the vector x ⃗ at the beginning and end of the phase. We have
Proof. Recall that, within a single phase, the algorithm uses a descending thresholds approach and updates a carefully selected subset of the coordinates whose gradient/marginal gain is above the current threshold (line 11). The subset is chosen so that no coordinate increases too much. This allows us to control the ℓ ∞ norm of the solution throughout the phase, which in turn allows us to deal with the nonmonotonicity of the function. The main technical challenge is to show that we gain enough throughout the phase despite that:
(1) we are restricting the set of high-gain coordinates that are updated to ensure that no coordinate increases too much, and (2) we terminate phase j when we fill up ϵjk budget to ensure that we get a feasible solution after 1/ϵ phases. We divide the analysis into two cases, depending on whether the threshold is ever updated on line 13.
Case 1: we have v end = v start . Note that the phase terminates with ∥z ⃗ end ∥ 1 = ϵjk in this case.
We start with an overview of the analysis. All of the coordinates that were updated throughout the phase had gradient at least v start , and v start is chosen so that the gain is high.
The step size η is chosen so that we update an 1ϵ fraction of the coordinates in S. Thus, ignoring 1ϵ factors, an iteration of the phase increases the function value by v start times the increase in the ℓ 1 -norm of the solution. Since the phase fills up at least ϵk of the budget, the increase in value over all iterations is at least ϵkv start . Finally, the choice of v start ensures that we reached the desired gain.
We now give the formal argument. We fix an iteration of the phase that updates x ⃗ and z ⃗ on lines 20-23, and analyze the gain in function value in the current iteration. We let x ⃗ , z ⃗ denote the vectors right before the update on lines 20-23. Let x ⃗ ′ be the vector x ⃗ right after the update on line 20, and let z ⃗ ′ be the vector z ⃗ right after the update on line 21.
We have:
In (a), we used the fact that x ⃗ ′x ⃗ ≥ 0 and f is concave in non-negative directions.
We can show (b) as follows. We have x ⃗ ′ ≤ z ⃗ ′ = z ⃗(η) and thus ∇f (x ⃗ ′ ) ≥ ∇f (z ⃗(η)) by DR-submodularity. Additionally, for every coordinate i ∈ T (η), we have
In (c), we have used that S(η) ⊆ T (η), g ⃗(η) i > 0 for all i ∈ T (η), and g ⃗(η
We can show (d) as follows. Since η ≤ η 1 , we have
where the second inequality is by the choice of η 1 .
Let η t and S t denote η and S in iteration t of the phase (note that we are momentarily overloading η 1 and η 2 here, and they temporarily stand for the step size η in iterations 1 and 2, and not for the step sizes on lines 16 and 18). By summing up the above inequality over all iterations, we obtain:
We can show (a) as follows. Recall that we have ∥z
In (b), we used the definition of v start on line 7. In (c), we used that f (x ⃗ * ) ≤ (1 + ϵ)M .
Case 2: we have v end ̸ = v start . Note that this implies that v end ≤ (1ϵ)v start , since line 13 was executed at least once during the phase.
The analysis for this case is more involved. The difficulty stems from the fact that the algorithm updates only a subset of the high-gain coordinates to ensure that no coordinate increases too much. This means that, at the end of the phase, there could still be coordinates whose gain is above the threshold but we cannot increase them further during the current phase. We now formally define this problematic set of coordinates:
The coordinates in A are allowed to increase, and thus we can account for their gain similarly to the previous case: their contribution to the overall gain is at least v end times the increase in ℓ 1 -norm of the coordinates in A. For the coordinates in A, we lower bound their contribution to the gain by the gradient times the increase in those coordinates.
The following lemma precisely states the lower bound that we obtain. Note that, for the coordinates in A, we use the gradient value and not v end , to get a tighter bound.
Lemma 5. We have
Proof. Fix an iteration of the phase that updates x ⃗ and z ⃗ on lines 20-23. Let x ⃗ , z ⃗ denote the vectors right before the update on lines 20-23. Let x ⃗ ′ be the vector x ⃗ right after the update on line 20, and let z ⃗ ′ be the vector z ⃗ right after the update on line 21.
We have:
In (a), we used the fact that x ⃗ ′x ⃗ ≥ 0 and f is concave in non-negative directions.
We can show (b) as follows. We have x ⃗ ′ ≤ z ⃗ ′ = z ⃗(η) and thus ∇f (x ⃗ ′ ) ≥ ∇f (z ⃗(η)) by DR-submodularity. Additionally, for every coordinate i ∈ T (η), we have ∇ i f (z ⃗(η)) > 0. Therefore, for every i ∈ T (η), we have
We can show (c) as follows. Since η ≤ η 1 , we have
where the second inequality is by the choice of η 1 . By the definition of S(η), we have g ⃗(η) i ≥ v ≥ v end for every i ∈ S(η). By the definition of T (η), we have g ⃗(η) i > 0 for every i ∈ T (η).
Equality (d) follows from the fact that S ∩ A = T (η) ∩ A, which we can show as follows. We have T (η) ⊆ S, and S \ T (η) is the set of all coordinates with negative gradient g ⃗(η). For every i ∈ A, we have
where the first inequality is by DR-submodularity (since z ⃗(η) ≤ z ⃗ end ) and the second inequality is by the definition of A and the fact that (z ⃗ end ) i < 1 for all i ∈ [n] (Lemma 3). Moreover, we have z ⃗(η) i < 1 for all i ∈ [n] (Lemma 3). Thus g ⃗(η) i > 0 for all i ∈ A, and hence S ∩ A = T (η) ∩ A.
In (e), we have used that S(η) ⊆ T (η) and g ⃗(η) is nonnegative on the coordinates of T (η). In (f), we have used that g ⃗(η) i ≥ v ≥ v end for all i ∈ S(η) . In (g), we have used that z ⃗(η) ≤ z ⃗ end and thus ∇f (z ⃗(η)) ≥ ∇f (z ⃗ end ) by DR-submodularity. Let η t , S t , T t (η), z ⃗ t (η), g ⃗ t (η) denote η, S, T (η), z ⃗(η), g ⃗(η) in iteration t of the phase (note that we are overloading η 1 and η 2 in this proof only, and they temporarily stand for the step size η in iterations 1 and 2, and not for the step sizes on lines 16 and 18). By summing up the above inequality over all iterations, we obtain:
With the lower bound on the gain given by Lemma 5 in hand, the remaining work is to relate it to the optimal solution. Before doing so, we show that each coordinate in A increases sufficiently. Recall that the coordinates in A were prevented from increasing further by the two conditions on z ⃗ i in the definition of S: z
In the latter case, the lemma follows immediately. In the former case, we use that each iteration increases the coordinates by at most ϵ 2 . Lemma 6. For every i ∈ A, we have
Proof of Lemma 6. Since S was empty at the previous threshold v end /(1ϵ), we have
Therefore we may assume it is the former. By Lemma 3,
where in the second inequality we used that (1ϵ) j-1 ≥ 1/3 for sufficiently small ϵ (since
Recall that the goal of the remainder of the analysis is to relate the lower bound on the gain given in Lemma 5 to the optimal solution x ⃗ * . To this end, we use the inner product between the gradient at the end of the phase and the optimal solution as our potential. To complete the analysis, we sandwich this potential between the gain guaranteed by Lemma 5 and our target gain. We show that the potential is at least our target gain in Lemma 7. In Lemmas 8 and 9, we relate the gain from Lemma 5 to the potential.
As noted above, the following lemma shows that the potential is lower bounded by our target gain. The argument is similar to the analysis of the measured continuous greedy. In this part of the analysis, we use the fact that we carefully controlled the ℓ ∞ norm of the solution, which is crucial when working with non-monotone objectives.
Proof. We have
In (a), we used that (1a)b ≥ max{a, b}a for all a, b ∈ [0, 1].
In (b), we used the fact that f is concave in non-negative directions.
In (c), we used the fact that the algorithm maintains the invariant that f (x ⃗ ) ≥ f (z ⃗) via the update on line 23.
In (d), we used Lemma 2.
In (e), we used Lemma 3.
The next two lemmas relate the gain guaranteed by Lemma 5 to our potential and use Lemma 7 to complete the analysis.
Recall that the phase terminates with either v end ≤ ϵv start or ∥z ⃗ end ∥ 1 = ϵjk. We consider each of these cases in turn.
Lemma 8. Suppose that v end ≤ ϵv start . We have:
Proof. By Lemma 5, we have:
In (a), we used Lemma 6 and the fact that 0
In (b), we used that ∥x ⃗ * ∥ 1 ≤ k and the definition of A.
In (c), we used Lemma 7. Inequality (d) follows from the definition of v start and the fact that f (x ⃗ * ) ≥ M .
The analysis of the second case uses a similar but more involved argument.
Lemma 9. Suppose that ∥z ⃗ end ∥ 1 = ϵjk. We have:
Proof. By Lemma 5, we have:
In (a), we used Lemma 6 and the fact that 0 ⃗ ≤ x ⃗ * ≤ 1 ⃗ and the fact that ∇ i f (z ⃗ end ) ≥ v end for all i ∈ A.
We can show (b) as follows. Recall that we are in the case ∥z
In (c), we used that, for all i / ∈ A, we have
In (d), we used Lemma 7.
Thus in both cases we achieve the desired gain, and the proof is complete.
Using induction and Theorem 4, we can show that the final solution returned by the algorithm is a 1/e -O(ϵ) approximation. By construction, the final solution satisfies ∥x ⃗ ∥ 1 ≤ k, and thus it also satisfies the constraint.
Theorem 10. Let x ⃗ be the final solution returned by Algorithm 2. We have ∥x ⃗
Recall that we refer to each iteration of the outer for loop as a phase. We refer to each iteration of the inner while loop as an iteration.
Theorem 11. The total number of iterations of the algorithm is O(log(n) log(1/ϵ)/ϵ 3 ).
Proof. There are O(1/ϵ) phases. In each phase, there are O(log(1/ϵ)/ϵ) different thresholds v: the initial threshold is v start , the threshold right before the final one is at least ϵv start , and each update on line 13 decreases the threshold by a (1-ϵ) factor. Thus it only remains to bound the number of iterations with the same threshold.
In the following, we fix a single threshold and we consider only the iterations of the phase at that threshold. Over these iterations, z ⃗ is non-decreasing in every coordinate, g ⃗ is nonincreasing in every coordinate by DR-submodularity and 1 ⃗z ⃗ ≥ 0 ⃗ , and the set S can only lose coordinates and thus |S| is non-increasing. Additionally, for each coordinate i ∈ [n], the increase (z ⃗ end ) i -(z ⃗ start ) i over the entire phase is at most ϵ + ϵ 2 : the increase in each iteration is η if i ∈ S and 0 otherwise; since η ≤ ϵ 2 and z
We say that an iteration is a large-step iteration if η = ϵ 2 and it is a smaller-step iteration if η < ϵ 2 .
We first consider the large-step iterations. Let t be the last large-step iteration, and let S t be the set S in that iteration. Let i ∈ S t . Note that i ∈ S t ′ for all iterations t ′ ≤ t, since S t ⊆ S t ′ . Thus every large-step iteration increases z
Since z ⃗ i increases by at most ϵ + ϵ 2 over the entire phase, it follows that the number of large-step iterations is O(1/ϵ).
Next, we consider the smaller-step iterations. Note that, in every smaller-step iteration except possibly the last one, we have |S(η)| ≤ (1ϵ)|S| (if η = η 2 < ϵ 2 , at the end of the iteration we have ∥z ⃗∥ 1 = ϵjk and thus the phase ends; if η = η 1 , our choice of η 1 ensures that |S(η 1 )| ≤ (1ϵ)|S|). Thus every smaller-step iteration decreases |S| by at least an (1ϵ) factor. Now note that |S| ≤ n in the first iteration, |S| ≥ 1 in the last iteration, and |S| can only decrease with each iteration. Thus the number of smaller-step iterations is O(log n/ϵ).
In summary, there are O(1/ϵ) phases, O(log(1/ϵ)/ϵ) different thresholds per phase, and O(log(n)/ϵ) iterations per threshold. Thus the total number of iterations of the algorithm is O(log(n) log(1/ϵ)/ϵ 3 ).
We experimentally evaluate our parallel algorithm on instances of non-concave quadratic programming (NQP) and softmax extension of determinantal point processes (DPP) that were randomly generated similarly to previous work (Bian et al., 2017).
NQP instances are functions of the form f (
, where H ∈ R n×n is a matrix with non-positive entries, h ⃗ ∈ R n . We randomly generated such instances as follows: we sampled each entry of H uniformly at random from [-10, 0], and we set h ⃗ = -0.2H ⊤ 1 ⃗ .
DPP instances are functions of the form f (x ⃗ ) = log det(diag(x ⃗ )(L -I) + I), where L ∈ R n×n is a psd matrix and I is the identity matrix. We randomly generated such instances as follows. We sampled the eigenvalues of L as follows: the i-th eigenvalue is ℓ i = e ri , where r i was sampled uniformly from [-0.5, 1]. We sampled a random orthogonal matrix V . We set L = V diag(ℓ 1 , . . . , ℓ n )V ⊤ .
Algorithms, implementation details, and parameter choices.
We empirically compared our parallel algorithm with the sequential continuous greedy algorithm (Chekuri et al., 2015;Bian et al., 2017) and the parallel multiplicative weights update algorithm (Ene et al., 2019) (see Section C in the supplement for pseudocode descriptions). We implemented the sequential continuous greedy algorithm using a step size of ϵ/n, leading to O(n/ϵ) iterations and adaptive evaluations. We implemented our algorithm with a more aggressive update of the thresholds on line 13: instead of the update v ← (1ϵ)v, we performed the update v ← 0.75 • v.
We used error ϵ = 0.05 and budget k = 10 in all of the experiments.
Computing infrastructure. We implemented the algorithms in C++ and ran the experiments on an iMac with a 3.3 GHz Intel Core i5 processor and 8 GB of memory. The code used for generating the instances and evaluating the algorithms can be found in the supplement.
Results. The experimental results are shown in Figure 1. Each value is the average value for 5 independently sampled instances and the error bar is ±1 standard deviation. The sequential continuous greedy algorithm achieved the highest solution value in all of the runs, and we report the value obtained by the parallel algorithm as the fraction of the continuous greedy solution value.
In all of the runs, our parallel algorithm achieves higher function value than the parallel multiplicative weights update algorithm, while the number of evaluations is significantly lower. The running time of the MWU algorithm is prohibitive on larger instances, and thus we were only able to compare the algorithms in the small n regime, which is advantageous for the MWU algorithm.