# distributed_submodular_maximization__aada57a6.pdf Journal of Machine Learning Research 17 (2016) 1-44 Submitted 10/14; Published 12/16 Distributed Submodular Maximization Baharan Mirzasoleiman baharanm@inf.ethz.ch Department of Computer Science ETH Zurich Universitaetstrasse 6, 8092 Zurich, Switzerland Amin Karbasi amin.karbasi@yale.edu School of Engineering and Applied Science Yale University New Haven, USA Rik Sarkar rsarkar@inf.ed.ac.uk Department of Informatics University of Edinburgh 10 Crichton St, Edinburgh EH8 9AB, United Kingdom Andreas Krause krausea@ethz.ch Department of Computer Science ETH Zurich Universitaetstrasse 6, 8092 Zurich, Switzerland Editor: JeffBilmes Many large-scale machine learning problems clustering, non-parametric learning, kernel machines, etc. require selecting a small yet representative subset from a large dataset. Such problems can often be reduced to maximizing a submodular set function subject to various constraints. Classical approaches to submodular optimization require centralized access to the full dataset, which is impractical for truly large-scale problems. In this paper, we consider the problem of submodular function maximization in a distributed fashion. We develop a simple, two-stage protocol Gree Di, that is easily implemented using Map Reduce style computations. We theoretically analyze our approach, and show that under certain natural conditions, performance close to the centralized approach can be achieved. We begin with monotone submodular maximization subject to a cardinality constraint, and then extend this approach to obtain approximation guarantees for (not necessarily monotone) submodular maximization subject to more general constraints including matroid or knapsack constraints. In our extensive experiments, we demonstrate the effectiveness of our approach on several applications, including sparse Gaussian process inference and exemplar based clustering on tens of millions of examples using Hadoop. Keywords: distributed computing, submodular functions, approximation algorithms, greedy algorithms, map-reduce 1. Introduction Numerous machine learning tasks require selecting representative subsets of manageable size out of large datasets. Examples range from exemplar based clustering (Dueck and c 2016 Baharan Mirzasoleiman, Amin Karbasi, Rik Sarkar and Andreas Krause. Mirzasoleiman, Karbasi, Sarkar and Krause Frey, 2007) to active set selection for non-parametric learning (Rasmussen, 2004), to viral marketing (Kempe et al., 2003), and data subset selection for the purpose of training complex models (Lin and Bilmes, 2011). Many such problems can be reduced to the problem of maximizing a submodular set function subject to cardinality or other feasibility constraints such as matroid, or knapsack constraints (Krause and Gomes, 2010; Krause and Golovin, 2012; Lee et al., 2009a). Submodular functions exhibit a natural diminishing returns property common in many well known objectives: the marginal benefit of any given element decreases as we select more and more elements. Functions such as entropy or maximum weighted coverage are typical examples of functions with diminishing returns. As a result, submodular function optimization has numerous applications in machine learning and social networks: viral marketing (Kempe et al., 2003; Babaei et al., 2013; Mirzasoleiman et al., 2012), information gathering (Krause and Guestrin, 2011), document summarization (Lin and Bilmes, 2011), and active learning (Golovin and Krause, 2011; Guillory and Bilmes, 2011). Although maximizing a submodular function is NP-hard in general, a seminal result of Nemhauser et al. (1978) states that a simple greedy algorithm produces solutions competitive with the optimal (intractable) solution (Nemhauser and Wolsey, 1978; Feige, 1998). However, such greedy algorithms or their accelerated variants (Minoux, 1978; Badanidiyuru and Vondr ak, 2014; Mirzasoleiman et al., 2015a) do not scale well when the dataset is massive. As data volumes in modern applications increase faster than the ability of individual computers to process them, we need to look at ways to adapt our computations using parallelism. Map Reduce (Dean and Ghemawat, 2008) is arguably one of the most successful programming models for reliable and efficient parallel computing. It works by distributing the data to independent machines: map tasks redistribute the data for appropriate parallel processing and the output then gets sorted and processed in parallel by reduce tasks. To perform submodular optimization in Map Reduce, we need to design suitable parallel algorithms. The greedy algorithms that work well for centralized submodular optimization do not translate easily to parallel environments. The algorithms are inherently sequential in nature, since the marginal gain from adding each element is dependent on the elements picked in previous iterations. This mismatch makes it inefficient to apply classical algorithms directly to parallel setups. In this paper, we develop a distributed procedure for maximizing submodular functions, that can be easily implemented in Map Reduce. Our strategy is to partition the data (e.g., randomly) and process it in parallel. In particular: We present a simple, parallel protocol, called Gree Di for distributed submodular maximization subject to cardinality constraints. It requires minimal communication, and can be easily implemented in Map Reduce style parallel computation models. We show that under some natural conditions, for large datasets the quality of the obtained solution is provably competitive with the best centralized solution. We discuss extensions of our approach to obtain approximation algorithms for (notnecessarily monotone) submodular maximization subject to more general types of constraints, including matroid and knapsack constraints. Distributed Submodular Maximization We implement our approach for exemplar based clustering and active set selection in Hadoop, and show how our approach allows to scale exemplar based clustering and sparse Gaussian process inference to datasets containing tens of millions of points. We extensively evaluate our algorithm on several machine learning problems, including exemplar based clustering, active set selection and finding cuts in graphs, and show that our approach leads to parallel solutions that are very competitive with those obtained via centralized methods (98% in exemplar based clustering, 97% in active set selection, 90% in finding cuts). This paper is organized as follows. We begin in Section 2 by discussing background and related work. In Section 3, we formalize the distributed submodular maximization problem under cardinality constraints, and introduce example applications as well as naive approaches toward solving the problem. We subsequently present our Gree Di algorithm in Section 4, and prove its approximation guarantees. We then consider maximizing a submodular function subject to more general constraints in Section 5. We also present computational experiments on very large datasets in Section 6, showing that in addition to its provable approximation guarantees, our algorithm provides results close to the centralized greedy algorithm. We conclude in Section 7. 2. Background and Related Work 2.1 Distributed Data Analysis and Map Reduce Due to the rapid increase in dataset sizes, and the relatively slow advances in sequential processing capabilities of modern CPUs, parallel computing paradigms have received much interest. Inhabiting a sweet spot of resiliency, expressivity and programming ease, the Map Reduce style computing model (Dean and Ghemawat, 2008) has emerged as prominent foundation for large scale machine learning and data mining algorithms (Chu et al., 2007; Ekanayake et al., 2008). A Map Reduce job takes the input data as a set of < key; value > pairs. Each job consists of three stages: the map stage, the shuffle stage, and the reduce stage. The map stage, partitions the data randomly across a number of machines by associating each element with a key and produce a set of < key; value > pairs. Then, in the shuffle stage, the value associated with all of the elements with the same key gets merged and send to the same machine. Each reducer then processes the values associated with the same key and outputs a set of new < key; value > pairs with the same key. The reducers output could be an input to another Map Reduce job, and a program in Map Reduce paradigm can consist of multiple rounds of map and reduce stages (Karloffet al., 2010). 2.2 Centralized and Streaming Submodular Maximization The problem of centralized maximization of submodular functions has received much interest, starting with the seminal work of Nemhauser et al. (1978). Recent work has focused on providing approximation guarantees for more complex constraints (for a more detailed account, see the recent survey by Krause and Golovin, 2012). Golovin et al. (2010) consider an algorithm for online distributed submodular maximization with an application to sensor selection. However, their approach requires k stages of communication, which is unrealistic Mirzasoleiman, Karbasi, Sarkar and Krause for large k in a Map Reduce style model. Krause and Gomes (2010) consider the problem of submodular maximization in a streaming model; however, their approach makes strong assumptions about the way the data stream is generated and is not applicable to the general distributed setting. Recently, Badanidiyuru et al. (2014) provide a single pass streaming algorithm for cardinality-constrained submodular maximization with 1/2 ε approximation guarantee to the optimum solution that makes no assumptions on the data stream. There has also been new improvements in the running time of the standard greedy solution for solving SET-COVER (a special case of submodular maximization) when the data is large and disk resident (Cormode et al., 2010). More generally, Badanidiyuru and Vondr ak (2014) and Mirzasoleiman et al. (2015a) improve the running time of the greedy algorithm for maximizing a monotone submodular function by reducing the number of oracle calls to the objective function. Very recently, Mirzasoleiman et al. (2016) provided a fast algorithm for maximizing non-monotone submodular functions under general constraints. In a similar spirit, Wei et al. (2014) propose a multi-stage framework for submodular maximization. In order to reduce the memory and computation cost, they apply an approximate greedy procedure to maximize surrogate (proxy) submodular functions instead of optimizing the target function at each stage. The above approaches are sequential in nature and it is not clear how to parallelize them. However, they can be naturally integrated into our distributed framework to achieve further acceleration. 2.3 Scaling Up: Distributed Algorithms Recent work has focused on specific instances of submodular optimization in distributed settings. Such scenarios often occur in large-scale graph mining problems where the data itself is too large to be stored on one machine. In particular, Chierichetti et al. (2010) address the MAX-COVER problem and provide a (1 1/e ϵ) approximation to the centralized algorithm at the cost of passing over the dataset many times. Their result is further improved by Blelloch et al. (2011). Lattanzi et al. (2011) address more general graph problems by introducing the idea of filtering, namely, reducing the size of the input in a distributed fashion so that the resulting, much smaller, problem instance can be solved on a single machine. This idea is, in spirit, similar to our distributed method Gree Di. In contrast, we provide a more general framework, and characterize settings where performance competitive with the centralized setting can be obtained. The present version is a significant extension of our previous conference paper (Mirzasoleiman et al., 2013), providing theoretical guarantees for both monotone and non-monotone submodular maximization problems subject to more general types of constraints, including matroid and knapsack constraints (described in Section 5), and additional empirical results (Section 6). Parallel to our efforts (Mirzasoleiman et al., 2013), Kumar et al. (2013) has taken the approach of adapting the sequential greedy algorithm to distributed settings. However, their method requires knowledge of the ratio between the largest and smallest marginal gains of the elements, and generally requires a non-constant (logarithmic) number of rounds. We provide empirical comparisons in Section 6.4. Distributed Submodular Maximization Figure 1: Cluster exemplars (left column) discovered by our distributed algorithm Gree Di described in Section 4 applied to the Tiny Images dataset (Torralba et al., 2008), and a set of representatives from each cluster. 3. Submodular Maximization In this section, we first review submodular functions and how to greedily maximize them. We then describe the distributed submodular maximization problem, the focus of this paper. Finally, we discuss two naive approaches towards solving this problem. 3.1 Greedy Submodular Maximization Suppose that we have a large dataset of images, e.g., the set of all images on the Web or an online image hosting website such as Flickr, and we wish to retrieve a subset of images that best represents the visual appearance of the dataset. Collectively, these images can be considered as exemplars that summarize the visual categories of the dataset as shown in Fig. 1. One way to approach this problem is to formalize it as the k-medoid problem. Given a set V = {e1, e2, . . . , en} of images (called ground set) associated with a (not necessarily symmetric) dissimilarity function, we seek to select a subset S V of at most k exemplars or cluster centers, and then assign each image in the dataset to its least dissimilar exemplar. If an element e V is assigned to exemplar v S, then the cost associated with e is the dissimilarity between e and v. The goal of the k-medoid problem is to choose exemplars that minimize the sum of dissimilarities between every data point e V and its assigned cluster center. Solving the k-medoid problem optimally is NP-hard, however, as we discuss in Section 3.4, we can transform this problem, and many other summarization tasks, to the problem of maximizing a monotone submodular function subject to a cardinality constraint max S V f(S) s.t. |S| k. (1) Submodular functions are set functions which satisfy the following natural diminishing returns property. Definition 1 (c.f., Nemhauser et al. (1978)) A set function f : 2V R is submodular, if for every A B V and e V \ B f(A {e}) f(A) f(B {e}) f(B). Furthermore, f is called monotone ifffor all A B V it holds that f(A) f(B). Mirzasoleiman, Karbasi, Sarkar and Krause We will generally additionally require that f is nonnegative, i.e., f(A) 0 for all sets A. Problem (1) is NP-hard for many classes of submodular functions (Feige, 1998). A fundamental result by Nemhauser et al. (1978) establishes that a simple greedy algorithm that starts with the empty set and iteratively augments the current solution with an element of maximum incremental value v = arg max v V \A f(A {v}), (2) continuing until k elements have been selected, is guaranteed to provide a constant factor approximation. Theorem 2 (Nemhauser et al., 1978) For any non-negative and monotone submodular function f, the greedy heuristic always produces a solution Agc[k] of size k that achieves at least a constant factor (1 1/e) of the optimal solution. f(Agc[k]) (1 1/e) max |A| k f(A). This result can be easily extended to f(Agc[l]) (1 e l/k) max|A| k f(A), where l and k are two positive integers (see, Krause and Golovin, 2012). 3.2 Distributed Submodular Maximization In many today s applications where the size of the ground set |V | = n is very large and cannot be stored on a single computer, running the standard greedy algorithm or its variants (e.g., lazy evaluations, Minoux, 1978; Leskovec et al., 2007; Mirzasoleiman et al., 2015a) in a centralized manner is infeasible. Hence, we seek a solution that is suitable for largescale parallel computation. The greedy method described above is in general difficult to parallelize, since it is inherently sequential: at each step, only the object with the highest marginal gain is chosen and every subsequent step depends on the preceding ones. Concretely, we consider the setting where the ground set V is very large and cannot be handled on a single machine, thus must be distributed among a set of m machines. While there are several approaches towards parallel computation, in this paper we consider the following model that can be naturally implemented in Map Reduce. The computation proceeds in a sequence of rounds. In each round, the dataset is distributed to m machines. Each machine i carries out computations independently in parallel on its local data. After all machines finish, they synchronize by exchanging a limited amount of data (of size polynomial in k and m, but independent of n). Hence, any distributed algorithm in this model must specify: 1) how to distribute V among the m machines, 2) which algorithm should run on each machine, and 3) how to communicate and merge the resulting solutions. In particular, the distributed submodular maximization problem requires the specification of the above steps in order to implement an approach for submodular maximization. More precisely, given a monotone submodular function f, a cardinality constraint k, and a number of machines m, we wish to produce a solution Ad[m, k] of size k such that f(Ad[m, k]) is competitive with the optimal centralized solution max|A| k,A V f(A). Distributed Submodular Maximization 3.3 Naive Approaches Towards Distributed Submodular Maximization One way to solve problem (1) in a distributed fashion is as follows. The dataset is first partitioned (randomly, or using some other strategy) onto the m machines, with Vi representing the data allocated to machine i. We then proceed in k rounds. In each round, all machines in parallel compute the marginal gains of all elements in their sets Vi. Next, they communicate their candidate to a central processor, who identifies the globally best element, which is in turn communicated to the m machines. This element is then taken into account for computing the marginal gains and selecting the next elements. This algorithm (up to decisions on how break ties) implements exactly the centralized greedy algorithm, and hence provides the same approximation guarantees on the quality of the solution. Unfortunately, this approach requires synchronization after each of the k rounds. In many applications, k is quite large (e.g., tens of thousands), rendering this approach impractical for Map Reduce style computations. An alternative approach for large k would be to greedily select k/m elements independently on each machine (without synchronization), and then merge them to obtain a solution of size k. This approach that requires only two rounds (as opposed to k), is much more communication efficient, and can be easily implemented using a single Map Reduce stage. Unfortunately, many machines may select redundant elements, and thus the merged solution may suffer from diminishing returns. It is not hard to construct examples for which this approach produces solutions that are a factor Ω(m) worse than the centralized solution. In Section 4, we introduce an alternative protocol Gree Di, which requires little communication, while at the same time yielding a solution competitive with the centralized one, under certain natural additional assumptions. 3.4 Applications of Distributed Submodular Maximization In this part, we discuss two concrete problem instances, with their corresponding submodular objective functions f, where the size of the datasets often requires a distributed solution for the underlying submodular maximization. 3.4.1 Large-scale Nonparametric Learning Nonparametric learning (i.e., learning of models whose complexity may depend on the dataset size n) are notoriously hard to scale to large datasets. A concrete instance of this problem arises from training Gaussian processes or performing MAP inference in Determinantal Point Processes, as considered below. Similar challenges arise in many related learning methods, such as training kernel machines, when attempting to scale them to large data sets. Active Set Selection in Sparse Gaussian Processes (GPs). Formally a GP is a joint probability distribution over a (possibly infinite) set of random variables XV , indexed by the ground set V , such that every (finite) subset XS for S = {e1, . . . , es} is distributed according to a multivariate normal distribution. More precisely, we have P(XS = x S) = N(XS; µS, ΣS,S), where µ = (µe1, . . . , µes) and ΣS,S = [Kei,ej] are prior mean and covariance matrix, respectively. The covariance matrix is parametrized via a positive definite kernel K( , ). As a Mirzasoleiman, Karbasi, Sarkar and Krause concrete example, when elements of the ground set V are embedded in a Euclidean space, a commonly used kernel in practice is the squared exponential kernel defined as follows: K(ei, ej) = exp( ||ei ej||2 2/h2). Gaussian processes are commonly used as priors for nonparametric regression. In GP regression, each data point e V is considered a random variable. Upon observations y A = x A + n A (where n A is a vector of independent Gaussian noise with variance σ2), the predictive distribution of a new data point e V is a normal distribution P(Xe | y A) = N(µe|A, Σ2 e|A), where mean µe|A and variance σ2 e|A are given by µe|A = µe + Σe,A(ΣA,A + σ2I) 1(x A µA), (3) σ2 e|A = σ2 e Σe,A(ΣA,A + σ2I) 1ΣA,e. (4) Evaluating (3) and (4) is computationally expensive as it requires solving a linear system of |A| variables. Instead, most efficient approaches for making predictions in GPs rely on choosing a small so called active set of data points. For instance, in the Informative Vector Machine (IVM) one seeks a set S such that the information gain, defined as f(S) = I(YS; XV ) = H(XV ) H(XV |YS) = 1 2 log det(I + σ 2ΣS,S) is maximized. It can be shown that this choice of f is monotone submodular (Krause and Guestrin, 2005a). For medium-scale problems, the standard greedy algorithms provide good solutions. For massive data however, we need to resort to distributed algorithms. In Section 6, we will show how Gree Di can choose near-optimal subsets out of a dataset of 45 million vectors. Inference for Determinantal Point Processes. A very similar problem arises when performing inference in Determinantal Point Processes (DPPs). DPPs (Macchi, 1975) are distributions over subsets with a preference for diversity, i.e., there is a higher probability associated with sets containing dissimilar elements. Formally, a point process P on a set of items V = {1, 2, ..., N} is a probability measure on 2V (the set of all subsets of V ). P is called determinantal point process if for every S V we have: P(S) det(KS), where K is a positive semidefinite kernel matrix, and KS [Kij]i,j S, is the restriction of K to the entries indexed by elements of S (we adopt that det(K ) = 1). The normalization constant can be computed explicitly from the following equation X S det(KS) = det(I + K), where I is the N N identity matrix. Intuitively, the kernel matrix determines which items are similar and therefore less likely to appear together. In order to find the most diverse and informative subset of size k, we need to find arg max|S| k det(KS) which is NP-hard, as the total number of possible subsets is exponential (Ko et al., 1995). However, the objective function is log-submodular, i.e. f(S) = log det(KS) is a submodular function (Kulesza, 2012). Hence, MAP inference in large DPPs is another potential application of distributed submodular maximization. Distributed Submodular Maximization 3.4.2 Large-scale Exemplar Based Clustering Suppose we wish to select a set of exemplars, that best represent a massive dataset. One approach for finding such exemplars is solving the k-medoid problem (Kaufman and Rousseeuw, 2009), which aims to minimize the sum of pairwise dissimilarities between exemplars and elements of the dataset. More precisely, let us assume that for the dataset V we are given a nonnegative function l : V V R (not necessarily assumed symmetric, nor obeying the triangle inequality) such that l( , ) encodes dissimilarity between elements of the underlying set V . Then, the cost function for the k-medoid problem is: L(S) = 1 |V | v V min e S l(e, υ). (5) Finding the subset S = arg min |S| k L(S) of cardinality at most k that minimizes the cost function (5) is NP-hard. However, by introducing an auxiliary element e0, a.k.a. phantom exemplar, we can turn L into a monotone submodular function (Krause and Gomes, 2010) f(S) = L({e0}) L(S {e0}). (6) In words, f measures the decrease in the loss associated with the set S versus the loss associated with just the auxiliary element. We begin with a phantom exemplar and try to find the active set that together with the phantom exemplar reduces the value of our loss function more than any other set. Technically, any point e0 that satisfies the following condition can be used as a phantom exemplar: max v V l(v, v ) l(v, e0), v V \ S. This condition ensures that once the distance between any v V \ S and e0 is greater than the maximum distance between elements in the dataset, then L(S {e0}) = L(S). As a result, maximizing f (a monotone submodular function) is equivalent to minimizing the cost function L. This problem becomes especially computationally challenging when we have a large dataset and we wish to extract a manageable-size set of exemplars, further motivating our distributed approach. 3.4.3 Other Examples Numerous other real world problems in machine learning can be modeled as maximizing a monotone submodular function subject to appropriate constraints (e.g., cardinality, matroid, knapsack). To name a few, specific applications that have been considered range from efficient content discovery for web crawlers and multi topic blog-watch (Chierichetti et al., 2010), over document summarization (Lin and Bilmes, 2011) and speech data subset selection (Wei et al., 2013), to outbreak detection in social networks (Leskovec et al., 2007), online advertising and network routing (De Vries and Vohra, 2003), revenue maximization in social networks (Hartline et al., 2008), and inferring network of influence (Gomez Rodriguez et al., 2010). In all such examples, the size of the dataset (e.g., number of webpages, Mirzasoleiman, Karbasi, Sarkar and Krause size of the corpus, number of blogs in the blogosphere, number of nodes in social networks) is massive, thus Gree Di offers a scalable approach, in contrast to the standard greedy algorithm, for such problems. 4. The Gree Di Approach for Distributed Submodular Maximization In this section we present our main results. We first provide our distributed solution Gree Di for maximizing submodular functions under cardinality constraints. We then show how we can make use of the geometry of data inherent in many practical settings in order to obtain strong data-dependent bounds on the performance of our distributed algorithm. 4.1 An Intractable, yet Communication Efficient Approach Before we introduce Gree Di, we first consider an intractable, but communication efficient two-round parallel protocol to illustrate the ideas. This approach, shown in Algorithm 1, first distributes the ground set V to m machines. Each machine then finds the optimal solution, i.e., a set of cardinality at most k, that maximizes the value of f in each partition. These solutions are then merged, and the optimal subset of cardinality k is found in the combined set. We denote this distributed solution by f(Ad[m, k]). As the optimum centralized solution Ac[k] achieves the maximum value of the submodular function, it is clear that f(Ac[k]) f(Ad[m, k]). For the special case of selecting a single element k = 1, we have f(Ac[1]) = f(Ad[m, 1]). Furthermore, for modular functions f (i.e., those for which f and f are both submodular), it is easy to see that the distributed scheme in fact returns the optimal centralized solution as well. In general, however, there can be a gap between the distributed and the centralized solution. Nonetheless, as the following theorem shows, this gap cannot be more than 1/ min(m, k). Furthermore, this result is tight. Theorem 3 Let f be a monotone submodular function and let k > 0. Then, f(Ad[m, k])) 1 min(m,k)f(Ac[k]). In contrast, for any value of m and k, there is a monotone submodular function f such that f(Ac[k]) = min(m, k) f(Ad[m, k]). The proof of all the theorems can be found in the appendix. The above theorem fully characterizes the performance of Algorithm 1 in terms of the best centralized solution. In practice, we cannot run Algorithm 1, since there is no efficient way to identify the optimum subset Ac i[k] in set Vi, unless P=NP. In the following, we introduce an efficient distributed approximation Gree Di. We will further show, that under some additional assumptions, much stronger guarantees can be obtained. 4.2 Our Gree Di Approximation Our efficient distributed method Gree Di is shown in Algorithm 2. It parallels the intractable Algorithm 1, but replaces the selection of optimal subsets, i.e., Ac i[k], by greedy solutions Agc i [k]. Due to the approximate nature of the greedy algorithm, we allow it to pick sets slightly larger than k. More precisely, Gree Di is a two-round algorithm that takes the ground set V , the number of partitions m, and the cardinality constraint κ. It first distributes the ground set over m machines. Then each machine separately runs the Distributed Submodular Maximization Algorithm 1 Inefficient Distributed Submodular Maximization Input: Set V , #of partitions m, constraints k. Output: Set Ad[m, k]. 1: Partition V into m sets V1, V2, . . . , Vm. 2: In each partition Vi find the optimum set Ac i[k] of cardinality k. 3: Merge the resulting sets: B = m i=1Ac i[k]. 4: Find the optimum set of cardinality k in B. Output this solution Ad[m, k]. Algorithm 2 Greedy Distributed Submodular Maximization (Gree Di) Input: Set V , #of partitions m, constraints κ. Output: Set Agd[m, κ]. 1: Partition V into m sets V1, V2, . . . , Vm (arbitrarily or at random). 2: Run the standard greedy algorithm on each set Vi to find a solution Agc i [κ]. 3: Find Agc max[κ] = arg max A{F(A) : A {Agc 1 [κ], . . . , Agc m[κ]}} 4: Merge the resulting sets: B = m i=1Agc i [κ]. 5: Run the standard greedy algorithm on B to find a solution Agc B [κ]. 6: Return Agd[m, κ] = arg max A{F(A) : A {Agc max[κ], Agc B [κ]}}. standard greedy algorithm by sequentially finding an element e Vi that maximizes the discrete derivative (2). Each machine i in parallel continues adding elements to the set Agc i [ ] until it reaches κ elements. We define Agc max[κ] to be the set with the maximum value among {Agc 1 [κ], Agc 2 [κ], . . . , Agc m[κ]}. Then the solutions are merged, i.e., B = m i=1Agc i [κ], and another round of greedy selection is performed over B until κ elements are selected. We denote this solution by Agc B [κ]. The final distributed solution with parameters m and κ, denoted by Agd[m, κ], is the set with a higher value between Agc max[κ] and Agc B [κ] (c.f., Figure 2 shows Gree Di schematically). The following result parallels Theorem 3. Theorem 4 Let f be a monotone submodular function and κ k. Then f(Agd[m, κ]) (1 e κ/k) min(m, k) f(Ac[k]). For the special case of κ = k the result of 4 simplifies to f(Agd[m, κ]) (1 1/e) min(m,k)f(Ac[k]). Moreover, it is straightforward to generalize Gree Di to multiple rounds (i.e., more than two) for very large datasets. In light of Theorem 3, one can expect that in general it is impossible to eliminate the dependency of the distributed solution on min(k, m)1. However, as we show in the sequel, in many practical settings, the ground set V exhibits rich geometrical structure that can be used to obtain stronger guarantees. 1. It has been very recently shown by Mirzasoleiman et al. (2015b) that the tightest dependency is Θ( p min(m, k)). Mirzasoleiman, Karbasi, Sarkar and Krause A1gc[k] = greedy(V1, k) Vm Vi V1 V2 A2gc[k] = greedy(V2, k) Aigc[k] = greedy(Vi, k) Amgc[k] = greedy(Vm, k) B = A1gc[k] U A2gc[k] U ... U Amgc[k] ABgc[k] = greedy(B, k) Amaxgc[k] = argmax{f(A1gc[k]), ... , f(Amgc[k])} Figure 2: Illustration of our two-round algorithm Gree Di 4.3 Performance on Datasets with Geometric Structure In practice, we can hope to do much better than the worst case bounds shown previously by exploiting underlying structure often present in real data and important set functions. In this part, we assume that a metric d : V V R exists on the data elements, and analyze performance of the algorithm on functions that vary slowly with changes in the input. We refer to these as Lipschitz functions: Definition 5 Let λ > 0. A set function f : 2V R is λ-Lipschitz w.r.t. metric d on V , if for any integer k, any equal sized sets S = {e1, e2, . . . , ek} V and S = {e 1, e 2, . . . , e k} V and any matching of elements: M = {(e1, e 1), (e2, e 2) . . . , (ek, e k)}, the difference between f(S) and f(S ) is bounded by: f(S) f(S ) λ X i d(ei, e i). (7) We can show that the objective functions from both examples in Section 3.4 are λ-Lipschitz for suitable kernels/distance functions: Proposition 6 Suppose that the covariance matrix of a Gaussian process is parametrized via a positive definite kernel K : V V R which is Lipschitz continuous with respect to metric d : V V R with constant L, i.e., for any triple of points x1, x2, x3 V , we have |K(x1, x3) K(x2, x3)| Ld(x1, x2). Then, the mutual information I(YS; XV ) = 1 2 log det(I+K) for the Gaussian process is λ-Lipschitz with λ = Lk3, where k is the number of elements in the selected subset S. Distributed Submodular Maximization Proposition 7 Let d : V V R be a metric on the elements of the dataset. Furthermore, let l : V V R encode the dissimilarity between elements of the underlying set V . Then for l = dα, α 1 the loss function L(S) = 1 |V | P v V mine S l(e, υ) (and hence also the corresponding submodular utility function f) is λ-Lipschitz with λ = αRα 1, where R is the diameter of the ball encompassing elements of the dataset in the metric space. In particular, for the k-medoid problem, which minimizes the loss function over all clusters with respect to l = d, we have λ = 1, and for the k-means problem, which minimizes the loss function over all clusters with respect to l = d2, we have λ = 2R. Beyond Lipschitz-continuity, many practical instances of submodular maximization can be expected to satisfy a natural density condition. Concretely, whenever we consider a representative set (i.e., optimal solution to the submodular maximization problem), we expect that any of its constituent elements has potential candidates for replacement in the ground set. For example, in our exemplar-based clustering application, we expect that cluster centers are not isolated points, but have many almost equally representative points close by. Formally, for any element v V , we define its α-neighborhood as the set of elements in V within distance α from v (i.e., α-close to v): Nα(v) = {w : d(v, w) α}. By λ-Lipschitz-continuity, it must hold that if we replace element v in set S by an αclose element v (i.e., v Nα(v)) to get a new set S of equal size, it must hold that |f(S) f(S )| αλ. As described earlier, our algorithm Gree Di partitions V into sets V1, V2, . . . Vm for parallel processing. If in addition we assume that elements are assigned uniformly at random to different machines, α-neighborhoods are sufficiently dense, and the submodular function is Lipschitz continuous, then Gree Di is guaranteed to produce a solution close to the centralized one. More formally, we have the following theorem. Theorem 8 Under the conditions that 1) elements are assigned uniformly at random to m machines, 2) for each ei Ac[k] we have |Nα(ei)| km log(k/δ1/m), and 3) f is λ-Lipschitz continuous, then with probability at least (1 δ) the following holds: f(Agd[m, κ]) (1 e κ/k)(f(Ac[k]) λαk). Note that once the above conditions are satisfied for small values of α (meaning that there is a high density of data points within a small distance from each element of the optimal solution) then the distributed solution will be close to the optimal centralized one. In particular if we let α 0, the distributed solution is guaranteed to be within a 1 eκ/k factor from the optimal centralized solution. This situation naturally corresponds to very large datasets. In the following, we discuss more thoroughly this important scenario. 4.4 Performance Guarantees for Very Large Datasets Suppose that our dataset is a finite sample V drawn i.i.d. from an underlying infinite set V, according to some (unknown) probability distribution. Let Ac[k] be an optimal solution in the infinite set, i.e., Ac[k] = arg max S V f(S), such that around each ei Ac[k], there is Mirzasoleiman, Karbasi, Sarkar and Krause a neighborhood of radius at least α where the probability density is at least β at all points (for some constants α and β). This implies that the solution consists of elements coming from reasonably dense and therefore representative regions of the dataset. Let us suppose g : R R is the growth function of the metric: g(α) is defined to be the volume of a ball of radius α centered at a point in the metric space. This means, for ei Ac[k] the probability of a random element being in Nα(ei) is at least βg(α) and the expected number of α neighbors of ei is at least E[|Nα(ei)|] = nβg(α). As a concrete example, Euclidean metrics of dimension D have g(α) = O(αD). Note that for simplicity we are assuming the metric to be homogeneous, so that the growth function is the same at every point. For heterogeneous spaces, we require g to have a uniform lower bound on the growth function at every point. In these circumstances, the following theorem guarantees that if the dataset V is sufficiently large and f is λ-Lipschitz, then Gree Di produces a solution close to the centralized one. Theorem 9 For n 8km log(k/δ1/m) λk) , where ε λk α , if the algorithm Gree Di assigns elements uniformly randomly to m processors , then with probability at least (1 δ), f(Agd[m, κ]) (1 e κ/k)(f(Ac[k]) ε). The above theorem shows that for very large datasets, Gree Di provides a solution that is within a 1 eκ/k factor of the optimal centralized solution. This result is based on the fact that for sufficiently large datasets, there is a suitably dense neighborhood around each member of the optimal solution. Thus, if the elements of the dataset are partitioned uniformly randomly to m processors, at least one partition contains a set Ac i[k] such that its elements are very close to the elements of the optimal centralized solution and provides a constant factor approximation of the optimal centralized solution. 4.5 Handling Decomposable Functions So far, we have assumed that the objective function f is given to us as a black box, which we can evaluate for any given set S independently of the dataset V . In many settings, however, the objective f depends itself on the entire dataset. In such a setting, we cannot use Gree Di as presented above, since we cannot evaluate f on the individual machines without access to the full set V . Fortunately, many such functions have a simple structure which we call decomposable. More precisely, we call a submodular function f decomposable if it can be written as a sum of submodular functions as follows (Krause and Gomes, 2010): f(S) = 1 |V | In other words, there is separate submodular function associated with every data point i V . We require that each fi can be evaluated without access to the full set V . Note that the exemplar based clustering application we discussed in Section 3.4 is an instance of this framework, among many others. Let us define the evaluation of f restricted to D V as Distributed Submodular Maximization A1gc[k] = greedy(V1, k, FV1) ABgc[k] = greedy(B, k, FU) Vm Vi V1 V2 A2gc[k] = greedy(V2, k, FV2) Aigc[k] = greedy(Vi, k, FVi) Amgc[k] = greedy(Vm, k, FVm) B = A1gc[k] U A2gc[k] U ... U Amgc[k] Amaxgc[k] = argmax{FU(A1gc[k]), ... , FU(Amgc[k])} Figure 3: Illustration of our two-round algorithm Gree Di for decomposable functions f D(S) = 1 |D| In the remaining of this section, we show that assigning each element of the dataset randomly to a machine and running Gree Di will provide a solution that is with high probability close to the optimum solution. For this, let us assume that fi s are bounded, and without loss of generality 0 fi(S) 1 for 1 i |V |, S V . Similar to Section 4.3 we assume that Gree Di performs the partition by assigning elements uniformly at random to the machines. These machines then each greedily optimize f Vi. The second stage of Gree Di optimizes f U, where U V is chosen uniformly at random with size n/m . Then, we can show the following result. First, for any fixed ϵ, m, k, let us define n0 to be the smallest integer such that for n n0 we have ln(n)/n ϵ2/(mk). Theorem 10 For n max(n0, m log(δ/4m)/ϵ2), ϵ < 1/4, and under the assumptions of Theorem 9, we have, with probability at least 1 δ, f(Agd[m, κ]) (1 e κ/k)(f(Ac[k]) 2ε). The above result demonstrates why Gree Di performs well on decomposable submodular functions with massive data even when they are evaluated locally on each machine. We will report our experimental results on exemplar-based clustering in the next section. Mirzasoleiman, Karbasi, Sarkar and Krause 4.6 Performance of Gree Di on Random Partitions Without Geometric Structure Very recently, Barbosa et al. (2015) and Mirrokni and Zadimoghaddam (2015) proved that under random partitioning of the data among m machines, the expected utility of Gree Di will be only a constant factor away from the optimum. Theorem 11 (Barbosa et al. (2015); Mirrokni and Zadimoghaddam (2015)) If elements are assigned uniformly at random to the machines, and κ = k, Gree Di gives a constant factor approximation guarantee (in the average case) to the optimum centralized solution2. E[f(Agd[m, k])] 1 1/e 2 f(Ac[k]). These results show that random partitioning of the data is sufficient to guarantee that Gree Di provides a constant factor approximation, irrespective of m and k, and without the requirement of any geometric structure. On the other hand, if geometric structure is present, the bounds from the previous sections can provide sharper approximation guarantees. 5. (Non-Monotone) Submodular Functions with General Constraints In this section we show how Gree Di can be extended to handle 1) more general constraints, and 2) non-monotone submodular functions. More precisely, we consider the following optimization setting Maximize f(S) Subject to S ζ. Here, we assume that the feasible solutions should be members of the constraint set ζ 2V . The function f( ) is submodular but may not be monotone. By overloading the notation we denote the set that achieves the above constrained optimization problem by Ac[ζ]. Throughout this section we assume that the constraint set ζ is hereditary, meaning that if A ζ then for any B A we also require that B ζ. Cardinality constraints are obviously hereditary, so are all the examples we mention below. 5.1 Matroid Constraints A matroid M is a pair (V, I) where V is a finite set (called the ground set) and I 2V is a family of subsets of V (called the independent sets) satisfying the following two properties: Heredity property: A B V and B I implies that A I, i.e. every subset of an independent set is independent. Augmentation property: If A, B I and |B| > |A|, there is an element e B \A such that A {e} I. 2. In fact, Mirrokni and Zadimoghaddam (2015) proved a 0.27-approximation guarantee which is slightly worse than (1 1/e)/2. Distributed Submodular Maximization Maximizing a submodular function subject to matroid constraints has found several applications in machine learning and data mining, ranging from content aggregation on the web (Abbassi et al., 2013) to viral marketing (Narayanam and Nanavati, 2012) and online advertising (Streeter et al., 2009). One way to approximately maximize a monotone submodular function f(S) subject to the constraint that each S is independent, i.e., S I, is to use a generalization of the greedy algorithm. This algorithm, which starts with an empty set and in each iteration picks the feasible element with maximum benefit until there is no more element e such that S {e} I, is guaranteed to provide a 1 2-approximation of the optimal solution (Fisher et al., 1978). Recently, this bound has been improved to (1 1/e) using the continuous greedy algorithm (Calinescu et al., 2011). For non-negative and non-monotone submodular functions with matroid constraints, the best known result is a 0.325-approximation based on simulated annealing (Gharan and Vondr ak, 2011). Curvature: For a submodular function f, the total curvature of f with respect to a set S is defined as: c = 1 min j V f(j|S \ j) Intuitively, the notion of curvature determines how far away f is from being modular. In other words, it measures how much the marginal gain of an element w.r.t. set S can decrease as a function of S. In general, c [0, 1], and for additive (modular) functions, c = 0, i.e., the marginal values are independent of S. In this case, the greedy algorithm returns the optimal solution to max{f(S) : S I}. In general, the greedy algorithm gives a 1 1+c-approximation to maximizing a non-decreasing submodular function with curvature c subject to a matroid constraint (Conforti and Cornu ejols, 1984). In case of the uniform matroid I = {S : |S| k}, the approximation factor is (1 e c)/c. Intersection of Matroids: A more general case is when we have p matroids M1 = (V, I1), M2 = (V, I2), ..., Mp = (V, Ip) on the same ground set V , and we want to maximize the submodular function f on the intersection of p matroids. That is, I = T i Ii consists of all subsets of V that are independent in all p matroids. This constraint arises, e.g., when optimizing over rankings (which can be modeled as intersections of two partition matroids). Another recent application considered is finding the influential set of users in viral marketing when multiple products need to be advertised and each user can tolerate only a small number of recommendations (Du et al., 2013). For p matroid constraints, the 1 p+1-approximation provided by the greedy algorithm (Fisher et al., 1978) has been improved to a (1 p ε)-approximation for p 2 by Lee et al. (2009b). For the non-monotone case, a 1/(p + 2 + 1/p + ε)-approximation based on local search is also given by Lee et al. (2009b). p-systems: p-independence systems generalize constraints given by the intersection of p matroids. Given an independence family I and a set V V , let S(V ) denote the set of maximal independent sets of I included in V , i.e., S(V ) = {A I | e V \A : A {e} / I}. Then we call (V, I) a p-system if for all nonempty V V we have max A S(V ) |A| p min A S(V ) |A|. Mirzasoleiman, Karbasi, Sarkar and Krause Similar to p matroid constraints, the greedy algorithm provides a 1 p+1-approximation guarantee for maximizing a monotone submodular function subject to a p-systems constraint (Fisher et al., 1978). For the non-monotone case, Gupta et al. (2010) provided a p/((p + 1)(3p + 3))-approximation can be achieved by combining an algorithm of Gupta et al. (2010) with the result for unconstrained submodular maximization of Buchbinder et al. (2012). This result has been recently tightened to p/((p + 1)(2p + 1)) by Mirzasoleiman et al. (2016). 5.2 Knapsack Constraints In many applications, including feature and variable selection in probabilistic models (Krause and Guestrin, 2005a) and document summarization (Lin and Bilmes, 2011), elements e V have non-uniform costs c(e) > 0, and we wish to find a collection of elements S that maximize f subject to the constraint that the total cost of elements in S does not exceed a given budget R, i.e. max S f(S) s.t. X v S c(v) R. Since the simple greedy algorithm ignores cost while iteratively adding elements with maximum marginal gains according (see Eq. 2) until |S| R, it can perform arbitrary poorly. However, it has been shown that taking the maximum over the solution returned by the greedy algorithm that works according to Eq. 2 and the solution returned by the modified greedy algorithm that optimizes the cost-benefit ratio v = arg max e V \S c(v) R c(S) f(S {e}) f(S) provides a (1 1/ e)-approximation of the optimal solution (Krause and Guestrin, 2005b). Furthermore, a more computationally expensive algorithm which starts with all feasible solutions of cardinality 3 and augments them using the cost-benefit greedy algorithm to find the set with maximum value of the objective function provides a (1 1/e)-approximation (Sviridenko, 2004). For maximizing non-monotone submodular functions subject to knapsack constraints, a (1/5 ε)-approximation algorithm based on local search was given by Lee et al. (2009a). Multiple Knapsack Constraints: In some applications such as procurement auctions (Garg et al., 2001), video-on-demand systems and e-commerce (Kulik et al., 2009), we have a d-dimensional budget vector R and a set of element e V where each element is associated with a d-dimensional cost vector. In this setting, we seek a subset of elements S V with a total cost of at most R that maximizes a non-decreasing submodular function f. Kulik et al. (2009) proposed a two-phase algorithm that provides a (1 1/e ε)-approximation for the problem by first guessing a constant number of elements of highest value, and then taking the value residual problem with respect to the guessed subset. For the non-monotone case, Lee et al. (2009a) provided a (1/5 ε)-approximation based on local search. p-system and d knapsack constraints: A more general type of constraint that has recently found interesting applications in viral marketing (Du et al., 2013) and personalized data summarization Mirzasoleiman et al. (2016) which can be cast by combining a Distributed Submodular Maximization Constraint Approximation (τ) monotone submodular functions non-monotone submodular functions Cardinality 1 1/e (Fisher et al., 1978) 0.325 (Gharan and Vondr ak, 2011) 1 matroid 1 1/e (Calinescu et al., 2011) 0.325 (Gharan and Vondr ak, 2011) p matroid 1/p ε (Lee et al., 2009b) 1/(p + 2 + 1/p + ε) (Lee et al., 2009b) 1 knapsack 1 1/e (Sviridenko, 2004) 1/5 - ε (Lee et al., 2009a) d knapsack 1 1/e ε Kulik et al. (2009) 1/5 - ε (Lee et al., 2009a) p-system 1/(p + 1) (Fisher et al., 1978) p/((p + 1)(2p + 1)) (Mirzasoleiman et al., 2016) p-system + d knapsack 1/(p+2d+1) (Badanidiyuru and Vondr ak, 2014) (1+ϵ)(p+1)(2p+2d+1)/p (Mirzasoleiman et al., 2016) Table 1: Approximation guarantees (τ) for monotone and non-monotone submodular maximization under different constraints. p-system with d knapsack constraints. For maximizing a monotone submodular function Badanidiyuru and Vondr ak (2014) proposed a modified version of the greedy algorithm that guarantees a 1/(p + 2d + 1)-approximation. By combining this algorithm with the one proposed in (Gupta et al., 2010), Mirzasoleiman et al. (2016) provided a fast algorithm for maximizing a non-monotone submodular function subject to a p-system and d knapsack constraints with (1 + ϵ)(p + 1)(2p + 2d + 1)/p-approximation. Table 1 summarizes the approximation guarantees for monotone and non-monotone submodular maximization under different constraints. 5.3 Gree Di Approximation Guarantee under More General Constraints Assume that we have a set of constraints ζ 2V that is hereditary. Further assume we have access to a black box algorithm X that gives us a constant factor approximation guarantee for maximizing a non-negative (but not necessarily monotone) submodular function f subject to ζ, i.e. X : (f, ζ) 7 AX ζ s.t. f(AX[ζ]) τ max A ζ f(A). (8) We can modify Gree Di to use any such approximation algorithm as a black box, and provide theoretical guarantees about the solution. In order to process a large dataset, it first distributes the ground set over m machines. Then instead of greedily selecting elements, each machine i in parallel separately runs the black box algorithm X on its local data in order to produce a feasible set AX i [ζ] meeting the constraints ζ. We denote by Agc max[ζ] the set with maximum value among AX i [ζ]. Next, the solutions are merged: B = m i=1AX i [ζ], and the black box algorithm is applied one more time to set B to produce a solution Agc B [ζ]. Then, the distributed solution for parameter m and constraints ζ, AXd[m, ζ], is the best among Agc max[ζ] and Agc B [ζ]. This procedure is given in more detail in Algorithm 3. The following result generalizes Theorem 4 for maximizing a submodular function subject to more general constraints. Mirzasoleiman, Karbasi, Sarkar and Krause Algorithm 3 Gree Di under General Constraints Input: Set V , #of partitions m, constraints ζ, submodular function f. Output: Set AXd[m, ζ]. 1: Partition V into m sets V1, V2, . . . , Vm. 2: In parallel: Run the approximation algorithm X on each set Vi to find a solution AX i [ζ]. 3: Find Agc max[ζ] = arg max A{F(A)|A {AX 1 [ζ], . . . , AX m[ζ]}}. 4: Merge the resulting sets: B = m i=1AX i [ζ]. 5: Run the approximation algorithm X on B to find a solution Agc B [ζ]. 6: Return AXd[m, ζ] = arg max{Agc max[ζ], Agc B [ζ]}. Theorem 12 Let f be a non-negative submodular function and X be a black box algorithm that provides a τ-approximation guarantee for submodular maximization subject to a set of hereditary constraints ζ. Then f(AXd[m, ζ])) τ min m, ρ([ζ]) f(Ac[ζ]), where f(Ac[ζ]) is the optimum centralized solution, and ρ([ζ]) = max A ζ |A|. Specifically, for submodular maximization subject to the matroid constraint M, we have ρ([A I]) = r M where r M is the rank of the matroid (i.e., the maximum size of any independent set in the system). For submodular maximization subject to the knapsack constraint R, we can bound ρ([c(A) R]) by R/ minv c(v) (i.e. the capacity of the knapsack divided by the smallest weight of any element). Performance on Datasets with Geometric Structure. When the submodular function f( ) and the constraint set ζ have more structure, then we can provide much better approximation guarantees. Assuming the elements of V are embedded in metric space with distance d : V V R+, we say that ζ is locally replaceable with respect to a set S V with parameter α > 0 if S V s.t. |S | = |S| and d (S, S ) α S ζ. Here, we define the distance d between two sets S and S of the same size k as follows. Let M be the set of all possible matchings between S and S , i.e., M = {((e1, e 1), . . . , (ek, e k)) s.t ei S and e i S for 1 i k}. Then d (S, S ) = min M maxi d(ei, e i). We require locality only with respect to Ac[ζ] to ensure that the optimum solution can be well approximated. What the locally replaceable property requires is that as elements of Ac[ζ] get replaced by nearby elements, the resulting set is also a feasible solution. Combining this property with λ-Lipschitzness will provide us with the following theorem. Theorem 13 Under the conditions that 1) elements are assigned uniformly at random to m machines, 2) for each ei Ac[ζ] we have |Nα(ei)| ρ([ζ])m log(ρ([ζ])/δ1/m), 3) f( ) is Distributed Submodular Maximization λ-Lipschitz, and 4) ζ is locally replaceable with respect to Ac[ζ] with parameter α, then with probability at least (1 δ), f(AXd[m, ζ])) τ(f(Ac[ζ]) λαρ([ζ])). The above result generalizes Theorem 8 for maximizing non-negative submodular functions subject to different constraints. Performance Guarantee for Very Large datasets. Similarly, we can generalize Theorem 9 for maximizing non-negative submodular functions subject to more general constraints. Suppose that our dataset is a finite sample V drawn i.i.d. from an underlying infinite set V, according to some (unknown) probability distribution. Let Ac[ζ] be an optimal solution in the infinite set, i.e., Ac[ζ] = arg max S V f(S), such that around each ei Ac[ζ], there is a neighborhood of radius at least α where the probability density is at least β at all points (for some constants α and β). Recall that g : R R is the growth function where g(α) measures the volume of a ball of radius α centered at a point in the metric space. Theorem 14 For n 8ρ([ζ])m log(ρ([ζ])/δ1/m) βg( ε λρ([ζ])) , where ε λρ([ζ]) α , if Gree Di assigns el- ements uniformly at random to m processors and under the conditions that f is λ-Lipschitz, and ζ is locally replaceable with respect to Ac[ζ] with parameter α , then with probability at least (1 δ), we have f(AXd[m, ζ])) τ(f(Ac[ζ]) ε). Performance Guarantee for Decomposable Functions. For the case of decomposable functions described in Section 4.5, the following generalization of Theorem 10 holds for maximizing a non-negative submodular function subject to more general constraints. Let us define n0 to be the smallest integer such that for n n0 we have ln(n)/n ϵ2/(m ρ([ζ])). Theorem 15 For n max(n0, m log(δ/4m)/ϵ2), ϵ < 1/4, and under the assumptions of Theorem 14, we have, with probability at least 1 δ, f(AXd[m, ζ])) τ(f(Ac[ζ]) 2ε). 6. Experiments In our experimental evaluation we wish to address the following questions: 1) how well does Gree Di perform compared to the centralized solution, 2) how good is the performance of Gree Di when using decomposable objective functions (see Section 4.5), and finally 3) how well does Gree Di scale in the context of massive datasets. To this end, we run Gree Di on three scenarios: exemplar based clustering, active set selection in GPs and finding the maximum cuts in graphs. We compare the performance of our Gree Di method to the following naive approaches: random/random: in the first round each machine simply outputs k randomly chosen elements from its local data points and in the second round k out of the merged mk elements, are again randomly chosen as the final output. Mirzasoleiman, Karbasi, Sarkar and Krause random/greedy: each machine outputs k randomly chosen elements from its local data points, then the standard greedy algorithm is run over mk elements to find a solution of size k. greedy/merge: in the first round k/m elements are chosen greedily from each machine and in the second round they are merged to output a solution of size k. greedy/max: in the first round each machine greedily finds a solution of size k and in the second round the solution with the maximum value is reported. For Gree Di, we let each of the m machines select a set of size αk, and select a final solution of size k among the union of the m solutions (i.e., among αkm elements). We present the performance of Gree Di for different parameters α > 0. For datasets where we are able to find the centralized solution, we report the ratio of f(Adist[k])/f(Agc[k]), where Adist[k] is the distributed solution (in particular Agd[m, αk, k] = Adist[k] for Gree Di). 6.1 Exemplar Based Clustering Our exemplar based clustering experiment involves Gree Di applied to the clustering utility f(S) (see Sec. 3.4) with d(x, x ) = x x 2. We performed our experiments on a set of 10,000 Tiny Images (Torralba et al., 2008). Each 32 by 32 RGB pixel image was represented by a 3,072 dimensional vector. We subtracted from each vector the mean value, normalized it to unit norm, and used the origin as the auxiliary exemplar. Fig. 4a compares the performance of our approach to the benchmarks with the number of exemplars set to k = 50, and varying number of partitions m. It can be seen that Gree Di significantly outperforms the benchmarks and provides a solution that is very close to the centralized one. Interestingly, even for very small α = κ/k < 1, Gree Di performs very well. Since the exemplar based clustering utility function is decomposable, we repeated the experiment for the more realistic case where the function evaluation in each machine was restricted to the local elements of the dataset in that particular machine (rather than the entire dataset). Fig 4b shows similar qualitative behavior for decomposable objective functions. Large scale experiments with Hadoop. As our first large scale experiment, we applied Gree Di to the whole dataset of 80,000,000 Tiny Images (Torralba et al., 2008) in order to select a set of 64 exemplars. Our experimental infrastructure was a cluster of 10 quad-core machines running Hadoop with the number of reducers set to m = 8000. Hereby, each machine carried out a set of reduce tasks in sequence. We first partitioned the images uniformly at random to reducers. Each reducer separately performed the lazy greedy algorithm on its own set of 10,000 images ( 123MB) to extract 64 images with the highest marginal gains w.r.t. the local elements of the dataset in that particular partition. We then merged the results and performed another round of lazy greedy selection on the merged results to extract the final 64 exemplars. Function evaluation in the second stage was performed w.r.t a randomly selected subset of 10,000 images from the entire dataset. The maximum running time per reduce task was 2.5 hours. As Fig. 5a shows, Gree Di highly outperforms the other distributed benchmarks and can scale well to very large datasets. Fig. 5b shows a set of cluster exemplars discovered by Gree Di where Fig. 5c and Fig. 5d show 100 nearest images to exemplars 26 and 63 (shown with red borders) in Fig. 5b. Distributed Submodular Maximization 2 4 6 8 10 0.8 Distributed/Centralized Random Random/ Gree DI (α=1) α=4/m (a) Global objective function 2 4 6 8 10 0.8 Distributed/Centralized Gree DI (α=1) α=4/m Random Random/ (b) Local objective function 20 40 60 80 100 Distributed/Centralized α=4/m α=2/m Gree DI (α=1) (c) Global objective function 20 40 60 80 100 Distributed/Centralized α=4/m α=2/m Gree DI (α=1) (d) Local objective function Figure 4: Performance of Gree Di compared to the other benchmarks. a) and b) show the mean and standard deviation of the ratio of distributed vs. centralized solution for global and local objective functions with budget k = 50 and varying the number m of partitions. c) and d) show the same ratio for global and local objective functions for m = 5 partitions and varying budget k, for a set of 10,000 Tiny Images. 6.2 Active Set Selection Our active set selection experiment involves Gree Di applied to the information gain f(S) (see Sec. 3.4) with Gaussian kernel, h = 0.75 and σ = 1. We used the Parkinsons Telemonitoring dataset (Tsanas et al., 2010) consisting of 5,875 bio-medical voice measurements with 22 attributes from people with early-stage Parkinson s disease. We normalized the vectors to zero mean and unit norm. Fig. 6b compares the performance Gree Di to the Mirzasoleiman, Karbasi, Sarkar and Krause 10 20 30 40 50 60 1.75 Distributed α=4/m α=2/m Max Greedy/ Gree DI (α=1) (a) Tiny Images 80M Figure 5: Performance of Gree Di compared to the other benchmarks. a) shows the distributed solution with m = 8000 and varying k for local objective functions on the whole dataset of 80,000,000 Tiny Images. b) shows a set of cluster exemplars discovered by Gree Di, and each column in c) shows 100 images nearest to exemplars 26 and d) shows 100 images nearest to exemplars 63 in b). benchmarks with fixed k = 50 and varying number of partitions m. Similarly, Fig 6a shows the results for fixed m = 10 and varying k. We find that Gree Di significantly outperforms the benchmarks. Large scale experiments with Hadoop. Our second large scale experiment consists of 45,811,883 user visits from the Featured Tab of the Today Module on Yahoo! Front Page (web, 2012). For each visit, both the user and each of the candidate articles are associated with a feature vector of dimension 6. Here, we used the normalized user features. Our experimental setup was a cluster of 8 quad-core machines running Spark with the number Distributed Submodular Maximization 1 2 3 4 5 6 7 8 9 10 0.65 Distributed/Centralized α=4/m α=2/m Gree DI (α=1) (a) Parkinsons Telemonitoring 20 40 60 80 100 0.65 Distributed/Centralized Gree DI (α=1) α=4/m α=2/m Greedy/ Greedy Greedy/ (b) Parkinsons Telemonitoring Figure 6: Performance of Gree Di compared to the other benchmarks. a) shows the ratio of distributed vs. centralized solution with k = 50 and varying m for Parkinsons Telemonitoring. b) shows the same ratio with m = 10 and varying k on the same dataset. of reducers set to m = 32. Each reducer performed the lazy greedy algorithm on its own set of 1,431,621 vectors ( 34MB) in order to extract 256 elements with the highest marginal gains w.r.t the local elements of the dataset in that particular partition. We then merged the results and performed another round of lazy greedy selection on the merged results to extract the final active set of size 256. The maximum running time per reduce task was 12 minutes for selecting 128 elements and 48 minutes for selecting 256 elements. Fig. 7 shows the performance of Gree Di compared to the benchmarks. We note again that Gree Di significantly outperforms the other distributed benchmarks and can scale well to very large datasets. Performance Comparison. Fig. 8 shows the speedup of Gree Di compared to the centralized greedy benchmark for different values of k and varying number of partitions m. As Fig. 8a shows, for small values of m, the speedup is almost linear in the number of machines. However, for large values of m the running time of the second stage of Gree Di increases and ultimately dominates the whole running time. Hence, we do not observe a linear speedup anymore. This effect can be observed in Fig. 8b. For larger values of k, the speedup is higher on fewer machines, but decreases more quickly by increasing m, as the second stage takes longer to complete. 6.3 Non-Monotone Submodular Function (Finding Maximum Cuts) We also applied Gree Di to the problem of finding maximum cuts in graphs. In our setting we used a Facebook-like social network (Opsahl and Panzarasa, 2009). This dataset includes the users that have sent or received at least one message in an online student community at University of California, Irvine and consists of 1,899 users and 20,296 directed ties. Fig. 9a and 9b show the performance of Gree Di applied to the cut function on graphs. We evaluated the objective function locally on each partition. Thus, the links Mirzasoleiman, Karbasi, Sarkar and Krause k 0 50 100 150 200 250 Distributed/Centralized Gree DI (,=1) Greedy/ ,=4/m ,=2/m Figure 7: Performance of Gree Di with m = 32 and varying budget k compared to the other benchmarks on Yahoo! Webscope data. m 5 10 15 20 25 30 k = 256 k = 128 k = 64 (a) Yahoo! front page m 100 200 300 400 500 k = 256 k = 128 k = 64 (b) Yahoo! front page Figure 8: Running time of Gree Di compared to the centralized greedy algorithm. a) shows the ratio of centralized vs. distributed solution with k = 64, 128, 256 and up to m = 32 machines for Yahoo Webscope data. b) shows the same ratio with k = 64, 128, 256 and up to m = 512 machines on the same dataset. Both experiments are performed on a cluster of 8 quad core machines. between the partitions are disconnected. Since the problem of finding the maximum cut in a graph is non-monotone submodular, we applied the Random Greedy algorithm proposed by Buchbinder et al. (2014) to find the near optimal solution in each partition. Distributed Submodular Maximization 5 10 15 20 0 Distributed/Centralized α=2/m α=4/m Gree DI (α=1) (a) Facebook-like social network 20 40 60 80 100 0 Distributed/Centralized α=2/m α=4/m Gree DI (α=1) (b) Facebook-like social network Figure 9: Performance of Gree Di compared to the other benchmarks. a) shows the mean and standard deviation of the ratio of distributed to centralized solution for budget k = 20 with varying number of machines m and b) shows the same ratio for varying budget k with m = 10 on Facebook-like social network. Although the cut function does not decompose additively over individual data points, perhaps surprisingly, Gree Di still performs very well, and significantly outperforms the benchmarks. This suggests that our approach is quite robust, and may be more generally applicable. 6.4 Comparision with Greedy Scaling. Kumar et al. (2013) recently proposed an alternative approach Greedy Scaling for parallel maximization of submodular functions. Greedy Scaling is a randomized algorithm that carries out a number (typically less than k) rounds of Map Reduce computations. We applied Gree Di to the submodular coverage problem in which given a collection V of sets, we would like to pick at most k sets from V in order to maximize the size of their union. We compared the performance of our Gree Di algorithm to the reported performance of Greedy Scaling on the same datasets, namely Accidents (Geurts et al., 2003) and Kosarak (Bodon, 2012). As Fig 10a and 10b shows, Gree Di outperforms Greedy Scaling on the Accidents dataset and its performance is comparable to that of Greedy Scaling in the Kosarak dataset. 7. Conclusion We have developed an efficient distributed protocol Gree Di, for constrained submodular maximization. We have theoretically analyzed the performance of our method and showed that under certain natural conditions it performs very close to the centralized (albeit impractical in massive datasets) solution. We have also demonstrated the effectiveness of our approach through extensive experiments, including active set selection in GPs on a dataset Mirzasoleiman, Karbasi, Sarkar and Krause 0 20 40 60 80 100 0.94 Distributed/Centralized Gree Di Greedy Scaling (a) Accidents 0 100 200 300 400 500 0.9997 Distributed/Centralized Gree Di Greedy Scaling (b) Kosarak Figure 10: Performance of Gree Di compared to the Greedy Scaling algorithm of Kumar et al. (2013) (as reported in their paper). a) shows the ratio of distributed to centralized solution on Accidents dataset with 340,183 elements and b) shows the same ratio for Kosarak dataset with 990,002 elements. The results are reported for varying budget k and varying number of machines m = n/µ where µ = O(knδ log n) and n is the size of the dataset. The results are reported for δ = 1/2. Note that the results presented by Kumar et al. (2013) indicate that Greedy Scaling generally requires a substantially larger number of Map Reduce rounds compared to Gree Di. of 45 million examples, and exemplar based summarization of a collection of 80 million images using Hadoop. We believe our results provide an important step towards solving submodular optimization problems in very large scale, real applications. Acknowledgments This research was supported by SNF 200021-137971, DARPA MSEE FA8650-11-1-7156, ERC St G 307036, a Microsoft Faculty Fellowship, an ETH Fellowship, Google Research Faculty Award, and a Scottish Informatics and Computer Science Alliance. Appendix A. Proofs This section presents the complete proofs of theorems presented in the article. Distributed Submodular Maximization A.1 Proof of Theorem 3 direction: The proof easily follows from the following lemmas. Lemma 16 max i f(Ac i[k]) 1 Proof Let Bi be the elements in Vi that are contained in the optimal solution, Bi = Ac[k] Vi. Then we have: f(Ac[k]) = f(B1 . . . Bm) = f(B1) + f(B2|B1) + . . . + f(Bm|Bm 1, . . . , B1). Using submodularity of f, for each i {1 . . . m}, we have f(Bi|Bi 1 . . . B1) f(Bi), and thus, f(Ac[k]) f(B1) + . . . + f(Bm). Since, f(Ac i[k]) f(Bi), we have f(Ac[k]) f(Ac 1[k]) + . . . + f(Ac m[k]). Therefore, f(Ac[k]) m max i f(Ac i[k]). Lemma 17 max i f(Ac i[k]) 1 Proof Let f(Ac[k]) = f({u1, . . . uk}). Using submodularity of f, we have Thus, f(Ac[k]) kf(u ) where u = arg maxif(ui). Suppose that the element with highest marginal gain (i.e., u ) is in Vj. Then the maximum value of f on Vj would be greater or equal to the marginal gain of u , i.e., f(Ac j[k]) f(u ) and since f(maxi f(Ac i[k])) f(Ac j[k]), we can conclude that f(max i f(Ac i[k])) f(u ) 1 Since f(Ad[m, k]) maxi f(Ac i[k]); from Lemma 16 and 17 we have f(Ad[m, k]) 1 min(m, k)f(Ac[k]). Mirzasoleiman, Karbasi, Sarkar and Krause direction: Let us consider a set of unbiased and independent Bernoulli random variables Xi,j for i {1, . . . , m} and j {1, . . . , k}, i.e., Pr(Xi,j = 1) = Pr(Xi,j = 0) = 1/2 and (Xi,j Xi ,j ) if i = i or j = j . Let us also define Yi = (Xi,1, . . . , Xi,k) for i {1, . . . , m}. Now assume that Vi = {Xi,1, . . . , Xi,k, Yi}, V = Sm i=1 Vi and f(S) = H(S), where H is the entropy of the subset S of random variables. Note that H is a monotone submodular function. It is easy to see that Ac i[k] = {Xi,1, . . . , Xi,k} or Ac i[k] = Yi as in both cases H(Ac i[k]) = k. If we assume Ac i[k] = {Xi,1, . . . , Xi,k}, then B = {Xi,j|1 i m, 1 j k}. Hence, by selecting at most k elements from B, we have H(Ad[m, k]) = k. On the other hand, the set of k elements that maximizes the entropy is {Y1, . . . , Ym}. Note that H(Yi) = k and Yi Yj for i = j. Hence, H(Ac) = k m if m k or otherwise H(Ac[k]) = k2. A.2 Proof of Theorem 4 Let us first mention a slight generalization over the performance of the standard greedy algorithm. It follows easily from the argument in (Nemhauser et al., 1978). Lemma 18 Let f be a non-negative submodular function, and let Agc[q] of cardinality q be the greedy selected set by the standard greedy algorithm. Then, f(Agc[q]) 1 e q k f(Ac[k]). By Lemma 18 we know that f(Agc i [κ]) (1 exp( κ/k))f(Ac i[k]). Now, let us define Bgc = m i=1Agc i [κ], Agc max[κ] = max i f(Agc i [κ]), A[κ] = arg max S Bgc&|S| κf(S). Then by using Lemma 18 again, we obtain f(Agd[m, κ]) max f(Agc max[κ]), (1 exp( κ/κ))f( A[κ]) (1 exp( κ/k)) min(m, k) f(Ac[k]). A.3 Proof of Proposition 6 Let K be a positive definite kernel matrix defined in section 3.4.1. If we replace a point ei S with another point e i V \ S, the corresponding row and column i in the modified kernel matrix K will be changed. W.l.o.g assume that we replace the first element e1 S with another element e 1 V \ S, i.e., K = K K has the following form with non-zero entries only on the first row and first column, a1 a2 ak a2 0 0 ... ... ... ... ak 0 0 Distributed Submodular Maximization Note that kernel is Lipschitz continuous with constant L, hence we have |ai| Ld(e1, e 1) for 1 i k. Then the absolute value of the change in the objective function would be f(S) f(S ) = 1 2 log det(I + K ) 1 2 log det(I + K) log det(I + K ) log det(I + K + K) log[det(I + K + K). det(I + K) 1] log det(I + K(I + K) 1) . (9) Note that since K is positive-definite, I+K is an invertible matrix. Furthermore, since K and K are symmetric matrices they both have k real eigenvalues. Therefore, (I + K) 1 has k eigenvalues λi = 1 1+λ i 1, for 1 i k, where λ 1 λ k are (non-negative) eigenvalues of kernel matrix K. Now, we bound the maximum eigenvalues of K and K(I + K) 1 respectively. Consider vectors x, x Rn, such that ||x||2 = ||x ||2 = 1. We have, x1 x2 ... xk a1 a2 ak a2 0 0 ... ... ... ... ak 0 0 x 1 x 2... x k x1 x2 ... xk Pk i=1 aix i a2x 1 ... akx 1 i=1 aix i + x 1 i=1 |aix i| + |x 1|. 2k Ld(e1, e 1), (10) where we used the following facts to derive the last inequality: 1) the Lipschitz continuity of the kernel gives us an upperbound on the values of |ai|, i.e., |ai| Ld(e1, e 1) for 1 i k; and 2) since ||x||2 = ||x ||2 = 1, the absolute value of the elements in vectors x and x cannot Mirzasoleiman, Karbasi, Sarkar and Krause be greater than 1, i.e., |xi| 1, |x i| 1, for 1 i k. Therefore, λmax( K) = max x: ||x||2=1 |x T Kx| 2k Ld(e1, e 1). Now, let v1, vk Rn be the k eigenvectors of matrix (I+K) 1. Note that {v1, vk} is an orthonormal system and thus for any x Rn we can write it as x = Pk i=1 civi, and we have ||x||2 2 = Pk i=1 c2 i . In order to bound the largest eigenvalue of K(I + K), we write x T K (I + K) 1x = x T K (I + K) 1 k X i,j=1 λicicjv T j Kvi (a) 2k Ld(e1, e 1) i,j=1 |ci||cj| = 2k Ld(e1, e 1) where in (a) we used Eq. 10 and the fact that λi 1 for 1 i k. Using Cauchy-Schwarz inequality k X and the assumption ||x||2 = 1, we conclude x T K (I + K) 1x 2k2Ld(e1, e 1) 2k2||x||2 2Ld(e1, e 1) 2k2Ld(e1, e 1). λmax K(I + K) 1 = max x: ||x||2=1 x T K (I + K) 1x 2k2Ld(e1, e 1). (11) Distributed Submodular Maximization Finally, we can write the determinant of a matrix as the product of its eigenvalues, i.e. det(I + K(I + K) 1) (1 + 2k2Ld(e1, e 1))k. (12) By substituting Eq. 11 and Eq. 12 into Eq. 9 we obtain f(S) f(S ) 1 log(1 + 2k2Ld(e1, e 1))k log(1 + 2k2Ld(e1, e 1)) k3Ld(e1, e 1), where in the last inequality we used log(1 + x) x, for x 0. Replacing all the k points in set S with another set S of the same size, we get f(S) f(S ) k3L i=1 d(ei, e i). Hence, the differential entropy of the Gaussian process is λ-Lipschitz with λ = Lk3. A.4 Proof of Proposition 7 Assume we have a set S of k exemplars, i.e., S0 = {e1, , ek}, and each element of the dataset v V is assigned to its closest exemplar. Now, if we replace set S with another set S of the same size, the loss associated with every element v V may be changed. W.l.o.g, assume we swap one exemplar at a time, i.e., in step i, 1 i k, we have Si = {e 1, , e i, ei+1, , ek}. Swapping the ith exemplar ei Si 1 with another element e i S , 4 cases may happen: 1) element v was not assigned to ei before and doesn t get assigned to e i, 2) element v was assigned to ei before and gets assigned to e i, 3) element v was not assigned to ei before and gets assigned to e i, 4) element v was assigned to ei before and gets assigned to another exemplar ex Si \ {e i}. For any element v V , we look into the four cases and show that in each case |l(e i, v) l(ei, v)| d(ei, e i) αRα 1. Case 1: In this case, element v was assigned to another exemplar ex Si \ ei and the assignment doesn t change. Therefore, there is no change in the value of the loss function. Case 2: In this case, element v was assigned to ei before and gets assigned to e i. let a = d(ei, v) and b = d(e i, v). Then we can write |l(e i, v) l(ei, v)| = |aα bα| = |(a b)|(aα 1 + aα 2b + + abα 2 + bα 1) d(ei, e i) αRα 1, (13) where in the last step we used triangle inequality |d(e t, v) d(et, v)| d(et, e t) and the fact that data points are in a ball of diameter R in the metric space. Mirzasoleiman, Karbasi, Sarkar and Krause Case 3: In this case, v was assigned to another exemplar ex Si 1 \ {ei} and gets assigned to e i, which implies that |l(e i, v) l(ex, v)| |l(ei, v) l(e i, v)|, since otherwise e would have been assigned to et before. Case 4: In the last case, element v was assigned to ei before and gets assigned to another exemplar ex Si \ {e i}. Thus, we have |l(ex, v) l(ei, v)| |l(e i, v) l(ei, v)| since otherwise v would have been assigned to ex before. Hence, in all four cases the following inequality holds: | min e Si 1 l(e, υ) min e Si l(e, υ)| |l(e i, v) l(ei, v)| d(ei, e i) αRα 1. By using Eq. 13 and averaging over all elements v V , we have |L(Si 1) L(Si)| = 1 |V | v V | min e Si 1 l(e, υ) min e Si l(e, υ)| αRα 1d(ei, e i). Thus, for any point e0 that satisfies max v V l(v, v ) l(v, e0), v V \ S, we have L({e0 S}) = L({S}) and thus |f(Si 1) f(Si)| = |L({e0}) L({e0 Si 1}) L({e0}) + L({e0 Si})| αRα 1d(ei, e i). Now, if we replace all the k points in set S with another set S of the same size, we get f(S) f(S ) = i=1 f(Si 1) f(Si) i=1 |f(Si 1) f(Si)| i=1 d(ei, e i). Therefore, for l = dα, the loss function is λ-Lipschitz with λ = αRα 1. A.5 Proof of Theorem 8 In the following, we say that sets S and S are γ-close if |f(S) f(S )| γ. First, we need the following lemma. Lemma 19 If for each ei Ac[k], |Nα(ei)| km log (k/δ1/m), and if V is partitioned into sets V1, V2, . . . Vm, where each element is randomly assigned to one set with equal probabilities, then there is at least one partition with a subset Ac i[k] such that |f(Ac[k]) f(Ac i[k])| λαk with probability at least (1 δ). Distributed Submodular Maximization Proof By the hypothesis, the α neighborhood of each element in Ac[k] contains at least km log (k/δ1/m) elements. For each ei Ac[k], let us take a set of m log (k/δ1/m) elements from its α-neighborhood. These sets can be constructed to be mutually disjoint, since each α-neighborhood contains m log (k/δ1/m) elements. We wish to show that at least one of the m partitions of V contains elements from α-neighborhoods of each element. Each of the m log (k/δ1/m) elements goes into a particular Vj with a probability 1/m. The probability that a particular Vj does not contain an element α-close to ei Ac[k] is δ1/m k . The probability that Vj does not contain elements α-close to one or more of the k elements is at most δ1/m (by union bound). The probability that each V1, V2, . . . Vm does not contain elements from the α-neighborhood of one or more of the k elements is at most δ. Thus, with high probability of at least (1 δ), at least one of V1, V2, . . . Vm contains an Ac i[k] that is λαk-close to Ac[k]. By lemma 19, for some Vi, |f(Ac[k]) f(Ac i[k]|) λαk with the given probability. Furthermore, f(Agc i [κ]) (1 e κ/k)f(Ac i[k]) by Lemma 18. Therefore, the result follows using arguments analogous to the proof of Theorem 4. A.6 Proof of Theorem 9 The following lemma says that in a sample drawn from distribution over an infinite dataset, a sufficiently large sample size guarantees a dense neighborhood near each element of Ac[k] when the elements are from representative regions of the data. Lemma 20 A number of elements: n 8km log (k/δ1/m) βg(α) , where α α , suffices to have at least 4km log (k/δ1/m) elements in the α-neighborhood of each ei Ac[k] with probability at least (1 δ), for small values of δ. Proof The expected number of α-neighbors of an ei Ac[k], is E[|Nα(ei)|] 8km log (k/δ1/m). We now show that in a random set of samples, at least a half of this number of neighbors is realized with high probability near each element of Ac[k]. This follows from a Chernoffbound: P[|Nα(ei)| 4km log (k/δ1/m)] e km log (k/δ1/m) (δ1/m/k)km. Therefore, the probability that some ei Ac[k] does not have a suitable sized neighborhood is at most k(δ1/m/k)km. For δ 1/k, kδkm δm. Therefore, with probability at least (1 δ), the α-neighborhood of each element ei Ac[k] contains at least 4km log (1/δ) elements. Lemma 21 For n 8km log(k/δ1/m) λk) , where ε λk α , if V is partitioned into sets V1, V2, . . . Vm, where each element is randomly assigned to one set with equal probabilities, then for sufficiently small values of δ, there is at least one partition with a subset Ac i[k] such that |f(Ac[k]) f(Ac i[k])| ε with probability at least (1 δ). Mirzasoleiman, Karbasi, Sarkar and Krause Proof Follows directly by combining Lemma 20 and Lemma 19. The probability that some element does not have a sufficiently dense ε/λk-neighborhood with km log(2k/δ1/m) elements is at most (δ/2) for sufficiently small δ, and the probability that some partition does not contain elements from the one or more of the dense neighborhoods is at most (δ/2). Therefore, the result holds with probability at least (1 δ). By Lemma 21, there is at least one Vi such that |f(Ac[k]) f(Ac i[k])| ε with the given probability. And f(Agd i [κ]) (1 e κ/k)f(Ac i[k]) using Lemma 18. The result follows using arguments analogous to the proof of Theorem 4. A.7 Proof of Theorem 10 Note that each machine has on the average n/m elements. Let us define Πi the event that n/2m < |Vi| < 2n/m. Then based on the Chernoffbound we know that Pr( Πi) 2 exp( n/8m). Let us also define ξi(S) the event that |f Vi(S) f(S)| < ϵ, for some fixed ϵ < 1 and a fixed set S with |S| k. Note that ξi(S) denotes the event that the empirical mean is close to the true mean. Based on the Hoeffding inequality (without replacement) we have Pr( = ξi S| 2 exp( 2nϵ2/m). Hence, Pr(ξi(S) Πi) 1 2 exp( 2nϵ2/m) 2 exp( n/8m). Let ξi be an event that |f Vi(S) f(S)| < ϵ, for any S such that |S| κ. Note that there are at most nκ sets of size at most κ. Hence, Pr(ξi Πi) 1 2nκ(exp( 2nϵ2/m) exp( n/8m)). (14) As a result, for ϵ < 1/4 we have Pr(ξi Πi) 1 4nκ exp( 2nϵ2/m). Since there are m machines, by the union bound we can conclude that Pr((ξi Πi) on all machines) 1 4mnκ exp( 2nϵ2/m). The above calculation implies that we need to choose δ 4mnκ exp( 2nϵ2/m). Let n0 be chosen in a way that for any n n0 we have ln(n)/n ϵ2/(mk). Then, we need to choose n as follows: n = max n0, m log(δ/4m) Hence for the above choice of n, there is at least one Vi such that |f(Ac[k]) f(Ac i[κ])| ε with probability 1 δ. Hence the solution is ϵ away from the optimum solution with probability 1 δ. Now if we confine the evaluation of f(Ac i) to data point only in machine i then under the assumption of Theorem 9 we lose another ϵ. Formally, the result at this point simply follows by combining Theorem 4 and Theorem 9. Distributed Submodular Maximization A.8 Proof of Theorem 12 The proof is similar to the proof of Theorem 3 and Theorem 4 and follows from the following lemmas. Lemma 22 maxi f(Ac i[ζ]) 1 Proof Let Bi be the elements in Vi that are contained in the optimal solution, Bi = Ac[ζ] Vi. Since Ac[ζ] ζ and ζ is a set of hereditary constraints, we must have Bi ζ as well. Using submodularity of f and by the same argument as in the proof of Lemma 16, we have f(Ac[ζ]) = f(B1 Bm) = f(B1) + f(B2|B1) + + f(Bm|Bm 1, , B1) f(B1) + + f(Bm). Since f(Ac i[ζ]) f(Bi) we get f(Ac[ζ]) f(Ac 1[ζ]) + + f(Ac m[ζ]) m max i f(Ac i[ζ]). Lemma 23 maxi f(Ac i[ζ]) 1 Proof The proof follows the outline of the proof of Lemma 17. Let f(Ac[ζ]) = f({u1, , uρ([ζ])}). Since Ac[ζ] ζ and ζ is a set of hereditary constraints, we have ui ζ. Using submodularity of f, we have i=1 f(ui) ρ([ζ])f(u ). where u = arg maxi f(ui). Suppose that u Vj, we get f(max i f(Ac i[ζ])) f(Ac j[ζ]) f(u ) 1 ρ([ζ])f(Ac[ζ]). Since f(Ad[m, ρ([ζ])]) maxi f(Ac i[ζ]); from Lemma 23 and 22 we have f(Ad[m, ρ([ζ])]) 1 min(m, ρ([ζ]))f(Ac[ζ]). (15) For the black box algorithm X with a τ-approximation guarantee, we have f(AX i [ζ]) τf(Ac i[ζ]). Now, we generalize the definitions used in the proof of Theorem 4 Bgc = m i=1Agc i [ζ], Agc max[ζ] = max i f(Agc i [ζ]), A[ζ] = arg max S Bgc&|S| ρ([ζ])f(S). Mirzasoleiman, Karbasi, Sarkar and Krause Then using Eq. 15 again, we obtain f(Agd[m, ζ]) max f(Agc max[ζ]), τf( A[ζ]) τ min(m, ρ([ζ]))f(Ac[ζ]). Note that since we do not use monotonicity of the submodular function in any of the proofs, the results hold in general for constrained maximization of any non-negative submodular function. A.9 Proof of Theorem 13 Lemma 24 If for each ei AX[ζ], |Nα(ei)| ρ([ζ])m log (ρ([ζ])/δ1/m), and if V is partitioned into sets V1, V2, . . . Vm, where each element is randomly assigned to one set with equal probabilities, then there is at least one partition with a subset AX i [ζ] ζ such that f(Ac[ζ]) f(AX i [ζ]) λαρ([ζ]) with probability at least (1 δ). The proof is similar to the proof of Lemma 19 by taking disjoint sets of size m log (ρ([ζ])/δ1/m) in an α-neighborhood of each ei Ac[ζ] and showing that with high probability, at least one of the m partitions of V contains elements from α-neighborhoods of each element in the optimal solution. Note that now the size of the optimal solution is at most ρ([ζ]). Since ζ is locally replaceable with parameter α, as elements of Ac[ζ] gets replaced by nearby elements in their α-neighborhood, the resulting set is also a feasible solution. By Lemma 24, for some Vi, f(Ac[ζ]) f(AX i [ζ]) λαρ([ζ]) with the given probability. On the other hand, for the black box algorithm X, we have f(AX i [ζ]) τf(Ac i[ζ]). Therefore, the result follows using arguments analogous to the proof of Theorem 12. A.10 Proof of Theorem 14 We use the following Lemmas to show that in a sample drawn from a ddistribution over an infinite dataset, a sufficiently large sample size guarantees a dense neighborhood near each element of the optimal solution. Lemma 25 A number of elements: n 8ρ([ζ])m log (ρ([ζ])/δ1/m) βg(α) , where α α , suffices to have at least 4ρ([ζ])m log (ρ([ζ])/δ1/m) elements in the α-neighborhood of each ei Ac[ζ] with probability at least (1 δ), for small values of δ. Lemma 26 For n 8ρ([ζ])m log(ρ([ζ])/δ1/m) βg( ε λρ([ζ])) , where ε λρ([ζ]) α , if V is partitioned into sets V1, V2, . . . Vm, where each element is randomly assigned to one set with equal probabilities, then for sufficiently small values of δ, there is at least one partition with a subset Ac i[ζ] such that |f(Ac[ζ]) f(Ac i[ζ])| ε with probability at least (1 δ). The proofs follows the same arguments as in the proof of Lemma 20 and 21. Recall that, by assumption ζ is locally replaceable with parameter α. Hence, for ε αλρ([ζ]), any set ε-close to the optimal solution is also a feasible solution. Distributed Submodular Maximization By Lemma 26, there is at least one Vi such that |f(Ac[ζ]) f(Ac i[ζ])| ε with the given probability. Furthermore, for the black box algorithm X, we have f(Agd i [ζ]) τf(Ac i[ζ]). Thus the result follows using arguments analogous to the proof of Theorem 12. A.11 Proof of Theorem 15 Again the proof follows the same line of reasoning as the proof of Theorem 10, except that for a constraint set ζ with ρ([ζ]) = max S ζ |S|, there are at most nρ([ζ]) feasible solutions. Using the same definitions for Πi and Ei as in the proof of Theorem 10, instead of Eq. 14 we get Pr(ξi Πi) 1 2nρ([ζ])(exp( 2nϵ2/m) exp( n/8m)). As a result, for ϵ < 1/4 and using union bound we conclude that Pr((ξi Πi) on all machines) 1 4mnρ([ζ]) exp( 2nϵ2/m). which implies that we need to choose δ 4mnρ([ζ]) exp( 2nϵ2/m). Now if n0 be chosen in a way that for any n n0 we have ln(n)/n ϵ2/(mk), we get n max(n0, m log(δ/4m)/ϵ2). Bearing in mind that ζ is locally replaceable, there is at least one Vi such that the solution Ac i[ζ] is feasible and ϵ away from the optimum solution with probability 1 δ. Now under the assumption of Theorem 14, if we evaluate f(Ac i) only on machine i, then we lose another ϵ. Now by combining Theorem 12 and Theorem 14 we get the desired result. Yahoo! academic relations. r6a, yahoo! front page today module user click log dataset, version 1.0, 2012. URL http://Webscope.sandbox.yahoo.com. Zeinab Abbassi, Vahab S Mirrokni, and Mayur Thakur. Diversity maximization under matroid constraints. In Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 32 40. ACM, 2013. Mahmoudreza Babaei, Baharan Mirzasoleiman, Mahdi Jalili, and Mohammad Ali Safari. Revenue maximization in social networks through discounting. Social Network Analysis and Mining, 3(4):1249 1262, 2013. Ashwinkumar Badanidiyuru and Jan Vondr ak. Fast algorithms for maximizing submodular functions. In Proceedings of the Twenty-Fifth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1497 1514. SIAM, 2014. Ashwinkumar Badanidiyuru, Baharan Mirzasoleiman, Amin Karbasi, and Andreas Krause. Streaming submodular maximization: Massive data summarization on the fly. In Proceedings of the 20th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 671 680. ACM, 2014. Rafael Barbosa, Alina Ene, Huy Nguyen, and Justin Ward. The power of randomization: Distributed submodular maximization on massive datasets. Proceedings of The 32nd International Conference on Machine Learning, pages 1236 1244, 2015. Mirzasoleiman, Karbasi, Sarkar and Krause Guy E Blelloch, Richard Peng, and Kanat Tangwongsan. Linear-work greedy parallel approximate set cover and variants. In Proceedings of the Twenty-Third Annual ACM Symposium on Parallelism in Algorithms and Architectures, pages 23 32. ACM, 2011. Ferenc Bodon. Kosarak dataset, 2012. URL http://fimi.ua.ac.be/data/. Niv Buchbinder, Michael Feldman, Joseph Naor, and Roy Schwartz. A tight linear time (1/2)-approximation for unconstrained submodular maximization. In 53rd Annual Symposium on Foundations of Computer Science (FOCS), pages 649 658. IEEE, 2012. Niv Buchbinder, Moran Feldman, Joseph SeffiNaor, and Roy Schwartz. Submodular maximization with cardinality constraints. In Proceedings of the Twenty-Fifth Annual ACMSIAM Symposium on Discrete Algorithms, pages 1433 1452. SIAM, 2014. Gruia Calinescu, Chandra Chekuri, Martin P al, and Jan Vondr ak. Maximizing a monotone submodular function subject to a matroid constraint. SIAM Journal on Computing, 40 (6):1740 1766, 2011. Flavio Chierichetti, Ravi Kumar, and Andrew Tomkins. Max-cover in map-reduce. In Proceedings of the 19th International Conference on World Wide Web, pages 231 240. ACM, 2010. Cheng Chu, Sang Kyun Kim, Yi-An Lin, Yuan Yuan Yu, Gary Bradski, Andrew Y Ng, and Kunle Olukotun. Map-reduce for machine learning on multicore. Advances in Neural Information Processing Systems, 19:281, 2007. Michele Conforti and G erard Cornu ejols. Submodular set functions, matroids and the greedy algorithm: tight worst-case bounds and some generalizations of the rado-edmonds theorem. Discrete Applied Mathematics, 7(3):251 274, 1984. Graham Cormode, Howard Karloff, and Anthony Wirth. Set cover algorithms for very large datasets. In Proceedings of the 19th ACM International Conference on Information and Knowledge Management, pages 479 488. ACM, 2010. Sven De Vries and Rakesh V Vohra. Combinatorial auctions: A survey. INFORMS Journal on Computing, 15(3):284 309, 2003. Jeffrey Dean and Sanjay Ghemawat. Mapreduce: simplified data processing on large clusters. Communications of the ACM, 51(1):107 113, 2008. Nan Du, Yingyu Liang, Maria Florina Balcan, and Le Song. Budgeted influence maximization for multiple products. ar Xiv preprint ar Xiv:1312.2164, 2013. Delbert Dueck and Brendan J Frey. Non-metric affinity propagation for unsupervised image categorization. In IEEE 11th International Conference on Computer Vision (ICCV), pages 1 8. IEEE, 2007. Jaliya Ekanayake, Shrideep Pallickara, and Geoffrey Fox. Mapreduce for data intensive scientific analyses. In IEEE Fourth International Conference on e Science, pages 277 284. IEEE, 2008. Distributed Submodular Maximization Uriel Feige. A threshold of ln n for approximating set cover. Journal of the ACM (JACM), 45(4):634 652, 1998. Marshall L. Fisher, George L. Nemhauser, and Laurence A. Wolsey. An analysis of approximations for maximizing submodular set functions - II. Mathematical Programming Study, (8):73 87, 1978. Rahul Garg, Vijay Kumar, and Vinayaka Pandit. Approximation algorithms for budgetconstrained auctions. In Approximation, Randomization, and Combinatorial Optimization: Algorithms and Techniques, pages 102 113. Springer, 2001. Karolien Geurts, Geert Wets, Tom Brijs, and Koen Vanhoof. Profiling of high-frequency accident locations by use of association rules. Transportation Research Record: Journal of the Transportation Research Board, 1840(1):123 130, 2003. Shayan Oveis Gharan and Jan Vondr ak. Submodular maximization by simulated annealing. In Proceedings of the Twenty-Second Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1098 1116. SIAM, 2011. Daniel Golovin and Andreas Krause. Adaptive submodularity: Theory and applications in active learning and stochastic optimization. Journal of Artificial Intelligence Research, pages 427 486, 2011. Daniel Golovin, Matthew Faulkner, and Andreas Krause. Online distributed sensor selection. In Proceedings of the 9th ACM/IEEE International Conference on Information Processing in Sensor Networks, pages 220 231. ACM, 2010. Manuel Gomez Rodriguez, Jure Leskovec, and Andreas Krause. Inferring networks of diffusion and influence. In Proceedings of the 16th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 1019 1028. ACM, 2010. Andrew Guillory and JeffBilmes. Active semi-supervised learning using submodular functions. In Proceedings of the Twenty-Seventh Conference Annual Conference on Uncertainty in Artificial Intelligence (UAI-11), pages 274 282. AUAI, 2011. Anupam Gupta, Aaron Roth, Grant Schoenebeck, and Kunal Talwar. Constrained nonmonotone submodular maximization: Offline and secretary algorithms. In Internet and Network Economics, pages 246 257. Springer, 2010. Jason Hartline, Vahab Mirrokni, and Mukund Sundararajan. Optimal marketing strategies over social networks. In Proceedings of the 17th International Conference on World Wide Web, pages 189 198. ACM, 2008. Howard Karloff, Siddharth Suri, and Sergei Vassilvitskii. A model of computation for mapreduce. In Proceedings of the Twenty-First Annual ACM-SIAM Symposium on Discrete Algorithms, pages 938 948. Society for Industrial and Applied Mathematics, 2010. Leonard Kaufman and Peter J Rousseeuw. Finding groups in data: an introduction to cluster analysis, volume 344. John Wiley & Sons, 2009. Mirzasoleiman, Karbasi, Sarkar and Krause David Kempe, Jon Kleinberg, and Eva Tardos. Maximizing the spread of influence through a social network. In Proceedings of the Ninth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 137 146. ACM, 2003. Chun-Wa Ko, Jon Lee, and Maurice Queyranne. An exact algorithm for maximum entropy sampling. Operations Research, 43(4):684 691, 1995. Andreas Krause and Daniel Golovin. Submodular function maximization. Tractability: Practical Approaches to Hard Problems, 3:19, 2012. Andreas Krause and Ryan G Gomes. Budgeted nonparametric learning from data streams. In Proceedings of the 27th International Conference on Machine Learning (ICML-10), pages 391 398, 2010. Andreas Krause and Carlos Guestrin. Near-optimal nonmyopic value of information in graphical models. In Proceedings of Uncertainty in Artificial Intelligence (UAI), page 5, 2005a. Andreas Krause and Carlos Guestrin. A note on the budgeted maximization on submodular functions. Technical Report CMU-CALD-05-103, Carnegie Mellon University, 2005b. Andreas Krause and Carlos Guestrin. Submodularity and its applications in optimized information gathering. ACM Transactions on Intelligent Systems and Technology (TIST), 2(4):32, 2011. Alex Kulesza. Determinantal point processes for machine learning. Machine Learning, 5 (2-3):123 286, 2012. Ariel Kulik, Hadas Shachnai, and Tami Tamir. Maximizing submodular set functions subject to multiple linear constraints. In Proceedings of the Twentieth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 545 554. Society for Industrial and Applied Mathematics, 2009. Ravi Kumar, Benjamin Moseley, Sergei Vassilvitskii, and Andrea Vattani. Fast greedy algorithms in mapreduce and streaming. In Proceedings of the 25th ACM Symposium on Parallelism in Algorithms and Architectures, pages 1 10. ACM, 2013. Silvio Lattanzi, Benjamin Moseley, Siddharth Suri, and Sergei Vassilvitskii. Filtering: a method for solving graph problems in mapreduce. In Proceedings of the Twenty-Third Annual ACM Symposium on Parallelism in Algorithms and Architectures, pages 85 94. ACM, 2011. Jon Lee, Vahab S Mirrokni, Viswanath Nagarajan, and Maxim Sviridenko. Non-monotone submodular maximization under matroid and knapsack constraints. In Proceedings of the 41st Annual ACM Symposium on Theory of Computing, pages 323 332. ACM, 2009a. Jon Lee, Maxim Sviridenko, and Jan Vondr ak. Submodular maximization over multiple matroids via generalized exchange properties. In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages 244 257. Springer, 2009b. Distributed Submodular Maximization Jure Leskovec, Andreas Krause, Carlos Guestrin, Christos Faloutsos, Jeanne Van Briesen, and Natalie Glance. Cost-effective outbreak detection in networks. In KDD 07: Proceedings of the 13th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 420 429. ACM, 2007. Hui Lin and JeffBilmes. A class of submodular functions for document summarization. In Proceedings of the 49th Annual Meeting of the Association for Computational Linguistics: Human Language Technologies-Volume 1, pages 510 520. Association for Computational Linguistics, 2011. Odile Macchi. The coincidence approach to stochastic point processes. Advances in Applied Probability, pages 83 122, 1975. Michel Minoux. Accelerated greedy algorithms for maximizing submodular set functions. In Optimization Techniques, pages 234 243. Springer, 1978. Vahab Mirrokni and Morteza Zadimoghaddam. Randomized composable core-sets for distributed submodular maximization. In Proceedings of the Forty-Seventh Annual ACM on Symposium on Theory of Computing, STOC 15, pages 153 162. ACM, 2015. Baharan Mirzasoleiman, Mahmoudreza Babaei, and Mahdi Jalili. Immunizing complex networks with limited budget. EPL (Europhysics Letters), 98(3):38004, 2012. Baharan Mirzasoleiman, Amin Karbasi, Rik Sarkar, and Andreas Krause. Distributed submodular maximization: Identifying representative elements in massive data. In Advances in Neural Information Processing Systems, pages 2049 2057, 2013. Baharan Mirzasoleiman, Ashwinkumar Badanidiyuru, Amin Karbasi, Jan Vondrak, and Andreas Krause. Lazier than lazy greedy. In Twenty-Ninth AAAI Conference on Artificial Intelligence, 2015a. Baharan Mirzasoleiman, Amin Karbasi, Ashwinkumar Badanidiyuru, and Andreas Krause. Distributed submodular cover: Succinctly summarizing massive data. In Advances in Neural Information Processing Systems, 2015b. Baharan Mirzasoleiman, Ashwin Badanidiyuru, and Amin Karbasi. Fast constrained submodular maximization: Personalized data summarization. In Proceedings of The 33nd International Conference on Machine Learning, 2016. Ramasuri Narayanam and Amit A Nanavati. Viral marketing for product cross-sell through social networks. In Machine Learning and Knowledge Discovery in Databases, pages 581 596. Springer, 2012. George L Nemhauser and Leonard A Wolsey. Best algorithms for approximating the maximum of a submodular set function. Mathematics of Operations Research, 3(3):177 188, 1978. George L Nemhauser, Laurence A Wolsey, and Marshall L Fisher. An analysis of approximations for maximizing submodular set functionsi. Mathematical Programming, 14(1): 265 294, 1978. Mirzasoleiman, Karbasi, Sarkar and Krause Tore Opsahl and Pietro Panzarasa. Clustering in weighted networks. Social Networks, 31 (2):155 163, 2009. Carl Edward Rasmussen. Gaussian processes in machine learning. In Advanced Lectures on Machine Learning, pages 63 71. Springer, 2004. Matthew Streeter, Daniel Golovin, and Andreas Krause. Online learning of assignments. In Advances in Neural Information Processing Systems, pages 1794 1802, 2009. Maxim Sviridenko. A note on maximizing a submodular set function subject to knapsack constraint. Operations Research Letters, 32, 2004. Antonio Torralba, Rob Fergus, and William T Freeman. 80 million tiny images: A large data set for nonparametric object and scene recognition. IEEE Transactions on Pattern Analysis and Machine Intelligence, 30(11):1958 1970, 2008. Athanasios Tsanas, Max Little, Patrick E Mc Sharry, Lorraine O Ramig, et al. Enhanced classical dysphonia measures and sparse regression for telemonitoring of parkinson s disease progression. In IEEE International Conference on Acoustics Speech and Signal Processing (ICASSP), pages 594 597. IEEE, 2010. Kai Wei, Yuzong Liu, Katrin Kirchhoff, and JeffBilmes. Using document summarization techniques for speech data subset selection. In Proceedings of NAACL-HLT, pages 721 726, 2013. Kai Wei, Rishabh Iyer, and JeffBilmes. Fast multi-stage submodular maximization. In Proceedings of the 31st International Conference on Machine Learning (ICML-14), pages 1494 1502, 2014.