Best of Both Worlds: Practical and Theoretically Optimal Submodular Maximization in Parallel
Year: 2022

Abstract:
For the problem of maximizing a monotone, submodular function with respect to a cardinality constraint k on a ground set of size n, we provide an algorithm that achieves the state-of-the-art in both its empirical performance and its theoretical properties, in terms of adaptive complexity, query complexity, and approximation ratio; that is, it obtains, with high probability, query complexity of O (n) in expectation, adaptivity of O (log(n)), and approximation ratio of nearly 1− 1/e. The main algorithm is assembled from two components which may be of independent interest. The first component of our algorithm, LINEARSEQ, is useful as a preprocessing algorithm to improve the query complexity of many algorithms. Moreover, a variant of LINEARSEQ is shown to have adaptive complexity of O(log(n/k)) which is smaller than that of any previous algorithm in the literature. The second component is a parallelizable thresholding procedure THRESHOLDSEQ for adding elements with gain above a constant threshold. Finally, we demonstrate that our main algorithm empirically outperforms, in terms of runtime, adaptive rounds, total queries, and objective values, the previous state-of-the-art algorithm FAST in a comprehensive evaluation with six submodular objective functions.

heading: 1 Introduction
text: The cardinality-constrained optimization of a monotone, submodular function f : 2N → R+, defined on subsets of a ground set N of size n, is a general problem formulation that is ubiquitous in wideranging applications, e.g. video or image summarization [30], network monitoring [26], information gathering [23], and MAP Inference for Determinantal Point Processes [20], among many others. The function f : 2N → R+ is submodular iff for all S ⊆ T ⊆ N , x 6∈ T , ∆ (x |T ) ≤ ∆ (x |S)1; and the function f is monotone if f(S) ≤ f(T ) for all S ⊆ T . In this paper, we study the following submodular maximization problem (SM)
maximizef(S), subject to |S| ≤ k, (SM) where f is a monotone, submodular function; SM is an NP-hard problem. There has been extensive effort into the design of approximation algorithms for SM over the course of more than 45 years, e.g. [32, 12, 10, 21, 24]. For SM, the optimal ratio has been shown to be 1− 1/e ≈ 0.63 [32].
1∆ (x |S) denotes the marginal gain of x to S: f(S ∪ {x})− f(S).
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
As instance sizes have grown very large, there has been much effort into the design of efficient, parallelizable algorithms for SM. Since queries to the objective function can be very expensive, the overall efficiency of an algorithm for SM is typically measured by the query complexity, or number of calls made to the objective function f [2, 9]. The degree of parallelizability can be measured by the adaptive complexity of an algorithm, which is the minimum number of rounds into which the queries to f may be organized, such that within each round, the queries are independent and hence may be arbitrariliy parallelized. Observe that the lower the adaptive complexity, the more parallelizable an algorithm is. To obtain a constant approximation factor, a lower bound of Ω(n) has been shown on the query complexity [24] and a lower bound of Ω(log(n)/ log log(n)) has been shown on the adaptive complexity [3].
Several algorithms have been developed recently that are nearly optimal in terms of query and adaptive complexities [14, 11, 17, 5]; that is, these algorithms achieve O (log n) adaptivity and O (npolylog(n)) query complexity (see Table 1). However, these algorithms use sampling techniques that result in very large constant factors that make these algorithms impractical. This fact is discussed in detail in Breuer et al. [9]; as an illustration, to obtain ratio 1− 1/e− 0.1 with 95% confidence, all of these algorithms require more than 106 queries of sets of size k/ log(n) in every adaptive round [9]; moreover, even if these algorithms are run as heuristics using a single sample, other inefficiencies preclude these algorithms of running even on moderately sized instances [9]. For this reason, the FAST algorithm of Breuer et al. [9] has been recently proposed, which uses an entirely different sampling technique called adaptive sequencing. Adaptive sequencing was originally introduced in Balkanski et al. [6], but the original version has quadratic query complexity in the size of the ground set and hence is still impractical on large instances. To speed it up, the FAST algorithm sacrifices theoretical guarantees to yield an algorithm that parallelizes well and is faster than all previous algorithms for SM in an extensive experimental evaluation. The theoretical sacrifices of FAST include: the adaptivity of FAST is Ω(log(n) log2(log n)), which is higher than the state-of-the-art, and more significantly, the algorithm obtains no approximation ratio for k < 8502; since many applications require small choices for k, this limits the practical utility of FAST. A natural question is thus: is it possible to design an algorithm that is both practical and theoretically optimal in terms of adaptivity, ratio, and total queries?

heading: 1.1 Contributions
text: In this paper, we provide three main contributions. The first contribution is the algorithm LINEARSEQ (LS, Section 2) that achieves with probability 1− 1/n a constant factor (4 +O(ε))−1 in expected linear query complexity and with O(log n) adaptive rounds (Theorem 1). Although the ratio of ≈ 0.25 is smaller than the optimal 1− 1/e ≈ 0.63, this algorithm can be used to improve the query
2The approximation ratio 1 − 1/e − 4ε of FAST holds with probability 1 − δ for k ≥ θ(ε, δ, k) = 2 log(2δ−1 log( 1
ε log(k)))/ε2(1− 5ε).
complexity of many extant algorithms, as we decribe in the related work section below. Interestingly, LINEARSEQ can be modified to have adaptivity O (log(n/k)) at a small cost to its ratio as discussed in Appendix F. This version of LINEARSEQ is a constant-factor algorithm for SM with smaller adaptivity than any previous algorithm in the literature, especially for values of k that are large relative to n.
Our second contribution is an improved parallelizable thresholding procedure THRESHOLDSEQ (TS, Section 3) for a commonly recurring task in submodular optimization: namely, add all elements that have a gain of a specified threshold τ to the solution. This subproblem arises not only in SM, but also e.g. in submodular cover [17] and non-monotone submodular maximization [4, 18, 15, 25]. Our TS accomplishes this task with probability 1− 1/n in O(log n) adaptive rounds and expected O(n) query complexity (Theorem 2), while previous procedures for this task only add elements with an expected gain of τ and use expensive sampling techniques [17]; have Ω(log2 n) adaptivity [22]; or have Ω(kn) query complexity [6].
Finally, we present in Section 3 the parallelized greedy algorithm PARALLELGREEDYBOOST (PGB), which is used in conjunction with LINEARSEQ and THRESHOLDSEQ to yield the final algorithm LS+PGB, which answers the above question affirmatively: LS+PGB obtains nearly the optimal 1− 1/e ratio with probability 1− 2/n in O(log n) adaptive rounds and O(n) queries in expectation; moreover, LS+PGB is faster than FAST in an extensive empirical evaluation (see Table 2). In addition, LS+PGB improves theoretically on the previous algorithms in query complexity while obtaining nearly optimal adaptivity (see Table 1).

heading: 1.2 Additional Related Work
text: Adaptive Sequencing. The main inefficiency of the adaptive sequencing method of Balkanski et al. [6] (which causes the quadratic query complexity) is an explicit check that a constant fraction of elements will be filtered from the ground set. In this work, we adopt a similar sampling technique to adaptive sequencing, except that we design the algorithm to filter a constant fraction of elements with only constant probability. This method allows us to reduce the quadratic query complexity of adaptive sequencing to linear query complexity while only increasing the adaptive complexity by a small constant factor. In contrast, FAST of Breuer et al. [9] speeds up adaptive sequencing by increasing the adaptive complexity of the algorithm through adaptive binary search procedures, which, in addition to the increasing the adaptivity by logarithmic factors, place restrictions on the k values for which the ratio can hold. This improved adaptive sequencing technique is the core of our THRESHOLDSEQ procedure, which has the additional benefit of being relatively simple to analyze.
Algorithms with Linear Query Complexity. Our LINEARSEQ algorithm also uses the improved adaptive sequencing technique, but in addition this algorithm integrates ideas from the Ω(n)-adaptive linear-time streaming algorithm of Kuhnle [24] to achieve a constant-factor algorithm with low adaptivity in expected linear time. Integration of the improved adaptive sequencing with the ideas of Kuhnle [24] is non-trivial, and ultimately this integration enables the theoretical improvement in query complexity over previous algorithms with sublinear adaptivity that obtain a constant ratio with
high probability (see Table 1). In Fahrbach et al. [17], a linear-time procedure SUBSAMPLEPREPROCESSING is described; this procedure is to the best of our knowledge the only algorithm in the literature that obtains a constant ratio with sublinear adaptive rounds and linear query complexity and hence is comparable to LINEARSEQ. However, SUBSAMPLEPREPROCESSING uses entirely different ideas from our LINEARSEQ and has much weaker theoretical guarantees: for input 0 < δ < 1, it obtains ratio δ 2
2×106 with probability 1 − δ in O(log(n)/δ) adaptive rounds and O(n) queries in expectation – the small ratio renders SUBSAMPLEPREPROCESSING impractical; also, its ratio holds only with constant probability. By contrast, with ε = 0.1, our LINEARSEQ obtains ratio ≈ 0.196 with probability 1− 1/n in O(log(n)) adaptive rounds and O(n) queries in expectation. Using LS for Preprocessing: Guesses of OPT. Many algorithms for SM, including FAST and all of the algorithms listed in Table 1 except for SM and our algorithm, use a strategy of guessing logarithmically many values of OPT. Our LINEARSEQ algorithm reduces the interval containing OPT from size k to a small constant size in expected linear time. Thus, LINEARSEQ could be used for preprocessing prior to running FAST or one of the other algorithms in Table 1, which would improve their query complexity without compromising their adaptive complexity or ratio; this illustrates the general utility of LINEARSEQ. For example, with this change, the theoretical adaptivity of FAST improves, although it remains worse than LS+PGB: the adaptive complexity of FAST becomes O ( 1 ε2 log(n) log ( 1 ε log(k) )) in contrast to the O ( 1 ε2 log (n/ε) ) of LS+PGB. Although SUBSAMPLEPREPROCESSING may be used for the same purpose, its ratio only holds with constant probability which would then limit the probability of success of any following algorithm.
Relationship of THRESHOLDSEQ to Existing Methods. The first procedure in the literature to perform the same task is the THRESHOLDSAMPLING procedure of Fahrbach et al. [17]; however, THRESHOLDSAMPLING only ensures that the expected marginal gain of each element added is at least τ and has large constants in its runtime that make it impractical [9]. In contrast, THRESHOLDSEQ ensures that added elements contribute a gain of at least τ with high probability and is highly efficient empirically. A second procedure in the literature to perform the same task is the ADAPTIVESEQUENCING method of Balkanski et al. [6], which similarly to THRESHOLDSEQ uses random permutations of the ground set; however, ADAPTIVE-SEQUENCING focuses on explicitly ensuring a constant fraction of elements will be filtered in the next round, which is expensive to check: the query complexity of ADAPTIVE-SEQUENCING is O(kn). In contrast, our THRESHOLDSEQ algorithm ensures this property with a constant probability, which is sufficient to ensure the adaptivity with the high probability of 1− 1/n in O(n) expected queries. Finally, a third related procedure in the literature is THRESHOLDSAMPLING of Kazemi et al. [22], which also uses random permutations to sample elements. However, this algorithm has the higher adaptivity of O (log(n) log(k)), in contrast to the O (log(n)) of THRESHOLDSEQ.
MapReduce Framework. Another line of work studying parallelizable algorithms for SM has focused on the MapReduce framework [13] in a distributed setting, e.g. [7, 8, 16, 28]. These algorithms divide the dataset over a large number of machines and are intended for a setting in which the data does not fit on a single machine. None of these algorithms has sublinear adaptivity and hence all have potentially large numbers of sequential function queries on each machine. In this work, our empirical evaluation is on a single machine with a large number of CPU cores; we do not evaluate our algorithms in a distributed setting.
Organization. The constant-factor algorithm LINEARSEQ is described and analyzed in Section 2; the details of the analysis are presented in Appendix C. The variant of LINEARSEQ with lower adaptivity is described in Appendix F. The algorithms THRESHOLDSEQ and PARALLELGREEDYBOOST are discussed at a high level in Section 3, with detailed descriptions of these algorithms and theoretical analysis presented in Appendices D and E. Our empirical evaluation is summarized in Section 4 with more results and discussion in Appendix H.

heading: 2 A Parallelizable Algorithm with Linear Query Complexity: LINEARSEQ
text: In this section, we describe the algorithm LINEARSEQ for SM (Alg. 1) that obtains ratio (4+O (ε))−1 in O (
1 ε3 log(n)
) adaptive rounds and expected O ( n ε3 ) queries. If ε ≤ 0.21, the ratio of LINEARSEQ
is lower-bounded by (4 + 16ε)−1 ≥ 0.135, which shows that a relatively large constant ratio is obtained even at large values of ε. An initial run of this algorithm is required for our main algorithm LS+PGB.
Algorithm 1 The algorithm that obtains ratio (4 +O (ε))−1 in O ( log(n)/ε3 ) adaptive rounds and
expected O ( n/ε3 ) queries.
1: procedure LINEARSEQ(f,N , k, ε) 2: Input: evaluation oracle f : 2N → R+, constraint k, error ε 3: a = arg maxu∈N f({u}) 4: Initialize A← {a} , V ← N , ` = d4(1 + 1/(βε)) log(n)e, β = ε/(16 log(8/(1− e−ε/2))) 5: for j ← 1 to ` do 6: Update V ← {x ∈ V : ∆ (x |A) ≥ f(A)/k} and filter out the rest 7: if |V | = 0 then break 8: V = {v1, v2, . . . , v|V |} ←random-permutation(V ) 9: Λ← {b(1 + ε)uc : 1 ≤ b(1 + ε)uc ≤ k, u ∈ N}
∪{bk + uεkc : bk + uεkc ≤ |V |, u ∈ N} ∪ {|V |} 10: B[λi] = false, for λi ∈ Λ 11: for λi ∈ Λ in parallel do 12: Tλi−1 ← {v1, v2, . . . , vλi−1} ; Tλi ← {v1, v2, . . . , vλi} ; T ′λi ← Tλi\Tλi−1 13: if ∆ ( T ′λi |A ∪ Tλi−1 ) /|T ′λi | ≥ (1− ε)f(A ∪ Tλi−1)/k then B[λi]← true 14: λ∗ ← max{λi ∈ Λ : B[λi] = false and ((λi ≤ k and B[1] to B[λi−1] are all true) or (λi > k and ∃m ≥ 1 s.t. | ⋃i−1 u=m T ′ λu | ≥ k and B[λm] to B[λi−1] are all true))}

heading: 15: A← A ∪ Tλ∗
text: 16: if |V | > 0 then return failure 17: return A′ ← last k elements added to A
Description of LS. The work of LS is done within iterations of a sequential outer for loop (Line 5); this loop iterates at most O (log(n)) times, and each iteration requires two adaptive rounds; thus, the adaptive complexity of the algorithm is O (log(n)). Each iteration adds more elements to the set A, which is initially empty. Within each iteration, there are four high-level steps: 1) filter elements from V that have gain less than f(A)/k (Line 6); 2) randomly permute V (Line 8); 3) compute in parallel the marginal gain of adding blocks of the sequence of remaining elements in V to A (for loop on Line 11); 4) select a prefix of the sequence V = ( v1, v2, . . . , v|V | ) to add to A (Line 14). The selection of the prefix to add is carefully chosen to approximately satisfy, on average, Condition 1 for elements added; and also to ensure that, with constant probability, a constant fraction of elements of V are filtered on the next iteration.
The following theorem states the theoretical results for LINEARSEQ. The remainder of this section proves this theorem, with intuition and discussion of the proof. The omitted proofs for all lemmata are provided in Appendix C.
Theorem 1. Let (f, k) be an instance of SM. For any constant 0 < ε < 1/2, the algorithm LINEARSEQ has adaptive complexity O ( log(n)/ε3 ) and outputs A′ ⊆ N with |A′| ≤ k such that the following properties hold: 1) The algorithm succeeds with probability at least 1 − 1/n. 2) There are O ( (1/(εk) + 1)n/ε3 ) oracle queries in expectation. 3) If the algorithm succeeds,[
4 + 4(2−ε)(1−ε)(1−2ε)ε ] f(A′) ≥ f(O), where O is an optimal solution to the instance (f, k).
Overview. The goal of this section is to produce a constant factor, parallelizable algorithm with linear query complexity. As a starting point, consider an algorithm1 that takes one pass through the ground set, adding each element e to candidate set A iff
∆ (e |A) ≥ f(A)/k. (1)
Condition 1 ensures two properties: 1) the last k elements in A contain a constant fraction of the value f(A); and 2) f(A) is within a constant fraction of OPT. By these two properties, the last k elements of A are a constant factor approximation to SM with exactly one query of the objective function per element of the ground set. For completeness, we give a pseudocode (Alg. 3) and proof in Appendix B. However, each query depends on all of the previous ones and thus there are n adaptive rounds. Therefore, the challenge is to approximately simulate Alg. 3 in a lowly adaptive (highly parallelizable) manner, which is what LINEARSEQ accomplishes.
1This algorithm is a simplified version of the streaming algorithm of Kuhnle [24].

heading: 2.1 Approximately Satisfying Condition 1
text: Discarding Elements. In one adaptive round during each iteration j of the outer for loop, all elements with gain to A of less than f (A) /k are discarded from V (Line 6). Since the size of A increases as the algorithm runs, by submodularity, the gain of these elements can only decrease and hence these elements cannot satisfy Condition 1 and can be safely discarded from consideration. The process of filtering thereby ensures the following lemma at termination.
Lemma 1. At successful termination of LINEARSEQ, f(O) ≤ 2f(A), where O ⊆ N is an optimal solution of size k.
Addition of Elements. Next, we describe the details of how elements are added to the set A. The random permutation of remaining elements on Line 8 constructs a sequence ( v1, v2, . . . , v|V | ) such that each element is uniformly randomly sampled from the remaining elements. By testing the marginal gains along the sequence in parallel, it is possible to determine a good prefix of the sequence (vi) to add to A to ensure the following: 1) Condition 1 is approximately satisfied; and 2) We will discard a constant fraction of V in the next iteration with constant probability. Condition 1 is important for the approximation ratio and discarding a constant fraction of V is important for the adaptivity and query complexity. Below, we discuss how to choose the prefix such that both are achieved. To speed up the algorithm, we do not test the marginal gain at each point in the sequence (vi), but rather test blocks of elements at once as determined by the index set Λ defined in the pseudocode.
Prefix Selection. Say a block is bad if this block does not satisfy the condition checked on Line 13 (which is an approximate, average form of Condition 1); otherwise, the block is good. At the end of an iteration, we select the largest block index λ∗, where this block is bad and the previous consecutive blocks which together have at least k elements are all good; or this block is bad and all the previous blocks are good blocks. Then, we add the prefix Tλ∗ = (v1, v2, . . . , vλ∗) into A. Now, the relevance of Condition 1 for the approximation ratio is that it implies f(A) ≥ 2f(A \A′), where A′ are the last k elements added to A. Lemma 2 shows that the conditions required on the marginal gains of blocks added imply an approximate form of this fact is satisfied by LINEARSEQ. Indeed, the proof of Lemma 2 informs the choice Λ of blocks evaluated and the computation of λ∗.
Lemma 2. Suppose LINEARSEQ terminates successfully. Then f(A) ≥ 2(1−ε+ε 2) 1+ε f(A\A ′).
Proof. If |A| ≤ k, the lemma is immediate, so assume |A| > k. For iteration j, let Tj,λ∗j denote the set added to A during iteration j; and let Tj,λ∗j = ∅ if the algorithm terminates before iteration j. Let Aj denote the value of A after iteration j. Define c = max{c ∈ N : A′ ⊆ (∪`j=cTj,λ∗j )}. Then, |Tc,λ∗c | > 0; and for any j > c, |Tj,λ∗j | < k. It holds that (∪ ` j=c+1Tj,λ∗j ) ⊂ A
′ ⊆ (∪`j=cTj,λ∗j ). Figure 1 shows how A is composed of these sets Tj,λ∗j and how each set is composed of blocks. The following claim is proven in Appendix C.3. Claim 1. It holds that ∆ ( Tc,λ∗c |A\A
′) ≥ (1 − ε) max{0, |Tc,λ∗c ∩ A′| − 2εk} · f(A\A′)/k. For j > c, it holds that ∆ ( Tj,λ∗j |Aj−1 ) ≥ 1−ε1+ε |Tj,λ∗j | · f(A\A ′)/k.
From Claim 1,
f(A)−f(A\A′) = ∆ ( Tc,λ∗c |A\A ′)+ ∑̀ j=c+1 ∆ ( Tj,λ∗j |Aj−1 ) ≥ (1− ε)(1− 2ε) 1 + ε ·f(A\A′), (2)
where Inequality 2 is proven in Appendix C.4.
In the remainder of this subsection, we will show that a (βε)-fraction of V is discarded at each iteration j of the outer for loop with probability at least 1/2, where β is a constant in terms of ε as defined on Line 4 in the pseudocode. The remainder of the proofs in this section are implicitly conditioned on the behavior of the algorithm prior to iteration j. The next lemma describes the behavior of the number of elements that will be filtered at iteration j + 1. Observe that the set Si defined in the next lemma is the set of elements that would be filtered at the next iteration if prefix Ti is added to A. Lemma 3. Let Si = {x ∈ V : ∆ (x |A ∪ Ti) < f(A ∪ Ti)/k} . It holds that |S0| = 0, |S|V || = |V |, and |Si| ≤ |Si+1|.
By Lemma 3, we know the number of elements in Si increases from 0 to |V | with i. Therefore, there exists a t such that t = min{i ∈ N : |Si| ≥ βε|V |}. If λ∗ ≥ t, |Sλ∗ | ≥ βε|V |, and we will successfully filter out more than (βε)-fraction of V at the next iteration. In this case, we say that the iteration j succeeds. Otherwise, if λ∗ < t, the iteration may fail. The remainder of the proof bounds the probability that λ∗ < t, which is an upper bound on the probability that iteration j fails. Let λt = max{λ ∈ Λ : λ < t}, and let λ′t = max({λ′ ∈ Λ : ∑ λ∈Λ,λ′≤λ≤λt |Tλ| ≥ k} ∪ {1}). If λ∗ < t, there must be at least one index λ between λ′t and λt such that the block T ′ λ is bad. The next lemma bounds the probability that any block T ′λ, with λ < λt, is bad. Lemma 4. Let t = min{i ∈ N : |Si| ≥ βε|V |}; λt = max{λ ∈ Λ : λ < t}; (Yi) be a sequence of independent and identically distributed Bernoulli trials, where the success probability is βε. Then for any λ < λt, Pr (B[λ] = false) ≤ Pr (∑|T ′λ| i=1 Yi > ε|T ′λ| ) .
Finally, we bound the probability that an iteration j of the outer for loop fails. Let B1 = {λ ∈ Λ : λ ≤ k and λ < λt}, B2 = {λ ∈ Λ : |Λ ∩ [λ, λt]| ≤ d1/εe}. Then
Pr (iteration j fails) ≤ Pr (∃λ ∈ B1 ∪B2 with B[λ] = false) ≤ 1/2, (3) where the proof of Inequality 3 is in Appendix C.7.

heading: 2.2 Proof of Theorem 1
text: From Section 2.1, the probability at any iteration of the outer for loop of successful filtering of an (βε)-fraction of V is at least 1/2. We can model the success of the iterations as a sequence of dependent Bernoulli random variables, with success probability that depends on the results of previous trials but is always at least 1/2.
Success Probability of LINEARSEQ. If there are at leastm = dlog1−βε(1/n)e successful iterations, the algorithm LINEARSEQ will succeed. The number of successful iterations X` up to and including the `-th iteration is a sum of dependent Bernoulli random variables. With some work (Lemma 6 in Appendix A), the Chernoff bounds can be applied to ensure the algorithm succeeds with probability at least 1− 1/n, as shown in Appendix C.2. Adaptivity and Query Complexity. Oracle queries are made on Lines 6 and 13 of LINEARSEQ. The filtering on Line 6 is in one adaptive round, and the inner for loop is also in one adaptive round. Thus, the adaptivity is proportional to the number of iterations of the outer for loop, O (`) = O ( log(n)/ε3 ) . For the query complexity, let Yi be the number of iterations between the (i− 1)-th and i-th successful iterations of the outer for loop. By Lemma 6 in Appendix A, E [Yi] ≤ 2. From here, we show in Appendix C.8 that there are at most O ( n/ε3 ) queries in expectation.
Approximation Ratio. Suppose LINEARSEQ terminates successfully. We have the approximation ratio as follows:
f(A′) (a) ≥ f(A)− f(A\A′) (b) ≥ f(A)− 1 + ε 2(1− ε+ ε2) f(A) (c) ≥ 1 4 + 4(2−ε)(1−ε)(1−2ε) · ε f(O),
where Inequality (a) is from submodularity of f , Inequality (b) is from Lemma 2, and Inequality (c) is from Lemma 1.

heading: 3 Improving to Nearly the Optimal Ratio
text: In this section, we describe how to obtain the nearly optimal ratio in nearly optimal query and adaptive complexities (Section 3.2). First, in Section 3.1, we describe THRESHOLDSEQ, a parallelizable procedure to add all elements with gain above a constant threshold to the solution. In Section 3.2, we describe PARALLELGREEDYBOOST and finally the main algorithm LS+PGB. Because of space constraints, the algorithms are described in the main text at a high level only, with detailed descriptions and proofs deferred to Appendices D and E.

heading: 3.1 The THRESHOLDSEQ Procedure
text: In this section, we discuss the algorithm THRESHOLDSEQ, which adds all elements with gain above an input threshold τ up to accuracy ε in O(log n) adaptive rounds and O(n) queries in expectation. Pseudocode is given in Alg. 4 in Appendix D.
Overview. The goal of this algorithm is, given an input threshold τ and size constraint k, to produce a set of size at most k such that the average gain of elements added is at least τ . As discussed in Section 1, this task is an important subroutine of many algorithms for submodular optimization (including our final algorithm), although by itself it does not produce any approximation ratio for SM. The overall strategy of our parallelizable algorithm THRESHOLDSEQ is analagous to that of LINEARSEQ, although THRESHOLDSEQ is considerably simpler to analyze. The following theorem summarizes the theoretical guarantees of THRESHOLDSEQ and the proofs are in Appendix D. Theorem 2. Suppose THRESHOLDSEQ is run with input (f, k, ε, δ, τ). Then, the algorithm has adaptive complexity O(log(n/δ)/ε) and outputs A ⊆ N with |A| ≤ k such that the following properties hold: 1) The algorithm succeeds with probability at least 1− δ/n. 2) There are O(n/ε) oracle queries in expectation. 3) It holds that f(A)/|A| ≥ (1 − ε)τ/(1 + ε). 4) If |A| < k, then ∆ (x |A) < τ for all x ∈ N .

heading: 3.2 The PARALLELGREEDYBOOST Procedure and the Main Algorithm
text: Algorithm 2 The PARALLELGREEDYBOOST procedure. 1: Input: evaluation oracle f : 2N → R+, constraint k,
constant α, value Γ such that Γ ≤ f(O) ≤ Γ/α, accuracy parameter ε 2: Initialize τ ← Γ/(αk), δ← 1/(log1−ε(α/3) + 1), A← ∅ 3: while τ ≥ Γ/(3k) do 4: τ ← τ(1− ε) 5: S ← THRESHOLDSEQ(fA,N , k − |A|, δ, ε/3, τ) 6: A← A ∪ S 7: if |A| = k then 8: return A 9: return A
In this section, we describe the greedy algorithm PARALLELGREEDYBOOST (PGB, Alg. 2) that uses multiple calls to THRESHOLDSEQ with descending thresholds. Next, our state-of-the-art algorithm LS+PGB is specified.
Description of PARALLELGREEDYBOOST. This procedure takes as input the results from running an α-approximation algorithm on the instance (f, k) of SM; thus, PARALLELGREEDYBOOST is not meant to be used as a standalone algorithm. Namely, PARALLELGREEDYBOOST takes as input Γ, the solution value of an α-approximation algorithm for SM; this solution value Γ is then boosted to ensure the ratio 1− 1/e− ε on the instance. The values of Γ and α are used to produce an initial threshold value τ for THRESHOLDSEQ. Then, the threshold value is iteratively decreased by a factor of (1− ε) and the call to THRESHOLDSEQ is iteratively repeated to build up a solution, until a minimum value for the threshold of Γ/(3k) is reached. Therefore, THRESHOLDSEQ is called at most O (log(1/α)/ε) times. We remark that α is not required to be a constant approximation ratio. Theorem 3. Let (f, k) be an instance of SM. Suppose an α- approximation algorithm for SM is used to obtain Γ, where the approximation ratio α holds with probability 1− pα. For any constant ε > 0,
the algorithm PARALLELGREEDYBOOST has adaptive complexity O ( logα−1
ε2 log
( n log(α−1)
ε
))
and outputs A ∈ N with |A| ≤ k such that the following properties hold: 1) The algorithm succeeds with probability at least 1− 1/n− pα. 2) If the algorithm succeeds, there are O ( n log ( α−1 ) /ε2 ) oracle queries in expectation. 3) If the algorithm succeeds, f(A) ≥ (1− 1/e− ε)f(O), where O is an optimal solution to the instance (f, k).
Proof. Success Probability. For the while loop in Line 3-8, there are no more than dlog1−ε(α/3)e iterations. If THRESHOLDSEQ completes successfully at every iteration, Algorithm 2 also succeeds. The probability that this occurs is lower bounded in Appendix E.1.1. For the remainder of the proof of Theorem 3, we assume that every call to THRESHOLDSEQ succeeds.
Adaptive and Query Complexity. There are at most dlog1−ε(α/3)e iterations of the while loop. Since log(x) ≤ x − 1, dlog1−ε(α/3)e = d log(α/3) log(1−ε)e, and ε < 1 − 1/e, it holds that dlog1−ε(α/3)e ≤ d log(3/α)
ε e. And for each iteration, queries to the oracle happen only on Line 5, the call to THRESHOLDSEQ. Since the adaptive and query complexity of THRESHOLDSEQ is O (log(n/δ)/ε) and O (n/ε), the adaptive and query complexities for Algorithm 2 are
O
( logα−1
ε2 log
( n log(α−1)
ε
)) , O ( logα−1 ε2 n ) , respectively.
Approximation Ratio. Let Aj be the set A we get after Line 6, and let Sj be the set returned by THRESHOLDSEQ in iteration j of the while loop. Let ` be the number of iterations of the while loop.
First, in the case that |A| < k at termination, THRESHOLDSEQ returns 0 ≤ |S`| < k − |A`−1| at the last iteration. From Theorem 2, for any o ∈ O, ∆ (o |A) < τ < Γ/(3k). By submodularity and monotonicity, f(O)− f(A) ≤ f(O ∪A)− f(A) ≤ ∑ o∈O\A ∆ (o |A) ≤ ∑ o∈O\A Γ/(3k) ≤ f(O)/3, and the ratio holds.
Second, consider the case that |A| = k. Suppose in iteration j + 1, THRESHOLDSEQ returns a nonempty set Sj+1. Then, in the previous iteration j, THRESHOLDSEQ returns a set Sj that 0 ≤ |Sj | < k − |Aj−1|. From Theorem 2,
f(O)− f(Aj+1) ≤ (
1− (1− ε/3)(1− ε) (1 + ε/3)k
|Aj+1\Aj | ) (f(O)− f(Aj)). (4)
The above inequality also holds when Aj+1 = Aj . Therefore, it holds that
f(O)− f(A) ≤ e− (1−ε/3)(1−ε) 1+ε/3 f(O) ≤ (1/e+ ε)f(O). (5)
The detailed proof of Inequality 4 and 5 can be found in Appendix E.1.2.
Main Algorithm: LS+PGB. To obtain the main algorithm of this paper (and its nearly optimal theoretical guarantees), we use PARALLELGREEDYBOOST with the solution value Γ and ratio α given by LINEARSEQ. Because this choice requires an initial run of LINEARSEQ, we denote this algorithm by LS+PGB. Thus, LS+PGB integrates LINEARSEQ and PARALLELGREEDYBOOST to get nearly the optimal 1− 1/e ratio with query complexity of O (n) and adaptivity of O (log(n)).

heading: 4 Empirical Evaluation
text: In this section, we demonstrate that the empirical performance of LS+PGB outperforms that of FAST for the metrics of total time, total queries, adaptive rounds, and objective value across six applications of SM: maximum cover on random graphs (MaxCover), twitter feed summarization (TweetSumm), image summarization (ImageSumm), influence maximization (Influence), revenue maximization (RevMax), and Traffic Speeding Sensor Placement (Traffic). See Appendix H.2 for the definition of the objectives. The sizes n of the ground sets range from n = 1885 to 100000.
Implementation and Environment. We evaluate the same implementation of FAST used in Breuer et al. [9]. Our implementation of LS+PGB is parallelized using the Message Passing Interface (MPI) within the same Python codebase as FAST (see the Supplementary Material for source code). Practical optimizations to LINEARSEQ are made, which do not compromise the theoretical guarantees, which are discussed in Appendix G. The hardware of the system consists of 40 Intel(R) Xeon(R) Gold 5218R CPU @ 2.10GHz cores (with 80 threads available), of which up to 75 threads are
made available to the algorithms for the experiments. On each instance, the algorithms are repeated independently for five repetitions, and the mean and standard deviation of the objective value, total queries, adaptive rounds and parallel (wall clock) runtime to the submodular function is plotted.
Parameters. The parameters ε, δ of FAST are set to enforce the nominal ratio of 1−1/e−0.1 ≈ 0.53 with probability 0.95; these are the same parameter settings for FAST as in the Breuer et al. [9] evaluation. The ε parameter of LS+PGB is set to enforce the same ratio with probability 1−2/n. With these parameters, FAST ensures its ratio only if k ≥ θ(ε, δ, k) = 2 log(2δ−1 log( 1ε log(k)))/ε
2(1− 5ε) ≥ 7103. Since k < 7103 on many of our instances, FAST is evaluated in these instances as a theoretically motivated heuristic. In contrast, the ratio of LS+PGB holds on all instances evaluated. We use exponentially increasing k values from n/1000 to n/10 for each application to explore the behavior of each algorithm across a broad range of instance sizes.
Overview of Results. Figure 2 illustrates the comparison with FAST across the ImageSumm and RevenueMax application; results on other applications are shown in Appendix H. Runtime: LS+PGB is faster than FAST by more than 1% on 80% of instances evaluated; and is faster by an order of magnitude on 14% of instances. Objective value: LS+PGB achieves higher objective by more than 1% on 50% of instances, whereas FAST achieves higher objective by more than 1% on 8% of instances. Adaptive rounds: LS+PGB achieves more than 1% fewer adaptive rounds on 75% of instances, while FAST achieves more than 1% fewer adaptive rounds on 22% of instances. Total queries: LS+PGB uses more than 1% fewer queries on 84% of scenarios with FAST using more than 1% fewer queries on 9% of scenarios. In summary, LS+PGB frequently gives substantial improvement in objective value, queries, adaptive rounds, and parallel runtime. Comparison of the arithmetic means of the metrics over all instances is given in Table 2. Finally, FAST and LS+PGB show very similar linear speedup with the number of processors employed: as shown in Fig. 7.

heading: 5 Concluding Remarks
text: In this work, we have introduced the algorithm LS+PGB, which is highly parallelizable and achieves state-of-the-art empirical performance over any previous algorithm for SM; also, LS+PGB is nearly optimal theoretically in terms of query complexity, adaptivity, and approximation ratio. An integral component of LS+PGB is our preprocessing algorithm LINEARSEQ, which reduces the interval containing OPT to a small constant size in expected linear time and low adaptivity, which may be independently useful. Another component of LS+PGB is the THRESHOLDSEQ procedure, which adds all elements with gain above a threshold in a parallelizable manner and improves existing algorithms in the literature for the same task.

heading: Acknowledgements
text: The work of Yixin Chen, Tonmoy Dey, and Alan Kuhnle was partially supported by Florida State University. The authors have received no third-party funding in direct support of this work. The authors have no additional revenues from other sources related to this work.

Review summaries:
In recent years there has been a theoretical breakthrough in low adaptivity O(log(n)) parallel algorithms for submodular maximization, but these algorithms are very slow in practice. Very recently, Breuer et al. proposed an algorithm, FAST, that trades off a small factor in asymptotic adaptive complexity to obtain the practically fastest algorithm. Building on this line of research, the authors propose a new algorithm with slightly better practical performance absent the tradeoff in asymptotic adaptive complexity. The authors demonstrate the practical performance advantage by benchmarking their approach against the FAST algorithm on numerous experiments.  

Submodular maximization is an important problem with many relevant applications in machine learning. With the huge amount of existing data and increase in the available computational power (increase in the number of machines), it is important to design scalable and parallelizable algorithms with provable guarantees from both practical and theoretical points of view. Many recent works have studied the submodular maximization in the context of adaptivity (“The degree of parallelizability can be measured by the adaptive complexity of an algorithm, which is the minimum number of rounds into which the queries to f may be organized”). In this paper, the authors propose an algorithm that achieves the state-of-the-art in terms of adaptive complexity, query complexity, and approximation ratio.

This work introduces the LinearSeq subroutine, which achieves a 1/4 approximation ratio for monotone submodular maximization subject to a cardinality constraint w.h.p. while using $O(n)$ queries in expectation and $O(\log n)$ adaptive rounds. The authors then show how to use LinearSeq together with a thresholding subroutine to achieve an algorithm with a $(1-1/e-\varepsilon)$ approximation for monotone submodular maximization in $O(1/\varepsilon^2 \log(n / \varepsilon))$ adaptive rounds and an expected $O(n / \varepsilon^2)$ queries. This is an improvement over the previous works of [Balkanski-Rubinstein-Singer, SODA 2019], [Ene-Nguyen, SODA 2019], [Fahrbach-Mirrokni-Zadimoghaddam, SODA 2019] in terms of $\varepsilon$ factors (for the query complexity) and practicability.  The authors compare their work extensively to the recent FAST algorithm of [Breuer-Balkanski-Singer, ICML 2019], which sacrifices some theoretical guarantees in favor of being extremely practical. Compared to all previous works, the main result here improves the query complexity without sacrificing the approximation quality or (much of) the adaptivity. The experiments in this paper compare against FAST and demonstrate that the algorithm in this paper is substantially faster than FAST without compromising on the solution quality (see Figure 3 in the full version of the paper).