A Memoization Framework for Scaling Submodular Optimization to Large Scale Problems
Rishabh Iyer Jeff Bilmes



[
A Memoization Framework for Scaling Submodular Optimization to Large Scale Problems

 Rishabh Iyer Jeff Bilmes 

 Microsoft Corporation University of Washington, Seattle
]


We are motivated by large scale submodular optimization problems, where standard algorithms that treat the submodular functions in the value oracle model do not scale. In this paper, we present a model called the precomputational complexity model, along with a unifying memoization based framework, which looks at the specific form of the given submodular function. A key ingredient in this framework is the notion of a precomputed statistic, which is maintained in the course of the algorithms. We show that we can easily integrate this idea into a large class of submodular optimization problems including constrained and unconstrained submodular maximization, minimization, difference of submodular optimization, optimization with submodular constraints and several other related optimization problems. Moreover, memoization can be integrated in both discrete and continuous relaxation flavors of algorithms for these problems. We demonstrate this idea for several commonly occurring submodular functions, and show how the precomputational model provides significant speedups compared to the value oracle model. Finally, we empirically demonstrate this for large scale machine learning problems of data subset selection and summarization.-1





§ INTRODUCTION
 
Submodular functions provide a rich class of expressible models for a variety of machine learning problems. Submodular functions occur naturally in two flavors. In minimization problems, they model notions of cooperation, attractive potentials and economies of scale, while in maximization problems, they model aspects of coverage, diversity and information. As a result, they have repeatedly appeared in several real world problems including document summarization <cit.>, image summarization <cit.>, data subset
selection and active learning <cit.>, image segmentation and denoising <cit.> and many others. A set function
f: 2^V →ℝ over a finite set V = {1, 2, …, n} is
submodular if for all subsets S, T ⊆ V, it holds that
f(S) + f(T) ≥ f(S ∪ T) + f(S ∩ T). Given a set S ⊆
V, we define the gain of an element j ∉ S in the
context S as f(j | S) ≜ f(S ∪ j) - f(S). A more intuitive characterization is the diminishing returns characterization: A function f is
submodular if it satisfies diminishing marginal returns,
namely f(j | S) ≥ f(j | T) for all S ⊆ T, j ∉ T,
and is monotone if f(j | S) ≥ 0 for all j ∉ S, S
⊆ V.

While submodular functions naturally occur in a number of real world applications, they also admit nice theoretical characterizations and algorithms. In particular, many simple iterative procedures like greedy <cit.>, local search <cit.> and majorization-minimization <cit.> yield theoretical guarantees for these problems. These algorithms are very efficient, scalable and easy to implement, and hence are being used more often in several large scale machine learning problems. The scale of machine learning problems are often massive, with dataset sizes of several hundreds of millions of examples. This has led to significant research in providing distributive, streaming, and multi-stage procedures for scaling these problems <cit.>. 

Many existing submodular optimization algorithms, treat the submodular function as a black box – a model, called the value oracle model. While this yields a useful way of quantifying the complexity of these algorithms, it does not provide an efficient way of implementing real world instances of submodular functions, which all appear in succinct representation, and can be can be stored in time and space, polynomial in the size of the ground set. A few recent papers have attempted to go beyond the value oracle model. For example, <cit.> study submodular optimization with noisy oracles. Similarly <cit.> study scenarios where we don't have the closed form expression of the submodular functions, but have a stochastic approximation (available for example through simulations). Similarly <cit.> study the problem of maximizing submodular functions drawn from specific probability distributions. 

In this paper, we take an orthogonal direction by introducing a new complexity model called the precomputational complexity model, along with a unifying memoization framework. We define the notion of a precomputed statistic, which is specific to a submodular function, and can be integrated easily into a large class of existing submodular optimization problems and algorithms.
We then show how to compute these statistics for several real world submodular functions occurring in applications, and theoretically demonstrate how the precomputational model reveals improved complexity results in comparison to the value oracle model. Finally we consider a real world application of speech corpus summarization and data subset selection, and show the scalability of our framework.-1

The idea of using precomputational statistics for speeding up submodular optimization in itself is not new. For example, <cit.> introduce the idea of incremental evaluation for the greedy algorithms for specific submodular functions such as entropy and mutual information. Similar ideas have also been discussed for Determinantal Point Processes <cit.>. In this work, we show that we can provide a unified framework via the construct of a precomputed statistics and show how this extends to both, a large class of submodular functions as well as a large class of optimization algorithms.


§ BASIC IDEAS AND BACKGROUND


We first introduce several key concepts we shall use in this paper. We start out by defining some fundamental characteristics of submodular functions.

The Submodular Polyhedron and  extension: 
For a submodular function f, the submodular polyhedron 𝒫_f and the corresponding base polytope ℬ_f are respectively defined as:  -1

    𝒫_f = { x : x(S) ≤ f(S), ∀ S ⊆ V }   
    ℬ_f = 𝒫_f ∩{ x : x(V) = f(V) }.

For a vector x ∈ℝ^V and a set X ⊆ V, we write x(X)
= ∑_j ∈ X x(j). Though 𝒫_f is defined via 2^n inequalities, its extreme point can be easily characterized <cit.>. Given any permutation σ of the ground set {1, 2, ⋯, n}, and an associated chain ∅ = S^σ_0 ⊆ S^σ_1 ⊆⋯⊆ S^σ_n = V with S^σ_i =
{σ(1), σ(2), …, σ(i) }, a vector h^f_σ satisfying,
 
    h^f_σ(σ(i) = f(S^σ_i) - f(S^σ_i-1) = f(σ(i) | S^σ_i-1), 
    ∀ i = 1, ⋯, n

forms an extreme point of 𝒫_f. Moreover, a natural convex extension of a submodular function, called the  extension <cit.> is closely related to the submodular polyhedron, and is defined as f̂(x) = max_h ∈𝒫_f⟨ h, x ⟩. Thanks to the properties of the polyhedron, f̂(x) can be efficiently computed: Denote σ_x as an ordering induced by x, such that x(σ_x(1)) ≥ x(σ_x(2)) ≥⋯ x(σ_x(n)). Then the  extension is f̂(x) = ⟨ h^f_σ, x ⟩ <cit.>.


Modular lower bounds (Subgradients):  Akin to convex functions, submodular functions have tight modular lower bounds. These bounds are related to the subdifferential ∂_f(Y) of the submodular set function f at a set Y ⊆ V, which is defined 
<cit.>
as:

    ∂_f(Y) = {y ∈ℝ^n: f(X) - y(X) ≥ f(Y) - y(Y), 
    for all  X ⊆ V}

Denote a subgradient at Y by h_Y ∈∂_f(Y). The extreme points of
∂_f(Y) may be computed in a manner similar to those of the submodular polyhedron. Let σ_Y
be a permutation of V that assigns the elements in Y to the first
|Y| positions (σ_Y(i) ∈ Y if and only if i ≤  |Y|). Then, h_Y = h^f_σ_Y (where h^f_σ_Y is as defined in eqn. (<ref>)) forms a lower bound of f, tight at Y — i.e.,
h_Y(X) = ∑_j ∈ X h_Y(j) ≤ f(X), ∀ X
⊆ V and h_Y(Y) = f(Y). Notice that the extreme points of a subdifferential are a subset of the extreme points of the submodular polyhedron.-1

Modular upper bounds (Supergradients): We can also define superdifferentials ∂^f(Y) of a submodular
function 
<cit.>
at
Y: 

    ∂^f(Y) = {y ∈ℝ^n: f(X) - y(X) ≤ f(Y) - y(Y); 
    for all  X ⊆ V}

It is possible, moreover, to provide specific supergradients <cit.> that define the following two modular upper bounds (when referring either one, we use m^f_X):


    m^f_X, 1(Y) ≜ f(X) - ∑_j ∈ X \ Y  f(j| X \ j) + ∑_j ∈ Y \ X f(j| ∅)1.3,   
    
    m^f_X, 2(Y) ≜ f(X) - ∑_j ∈ X \ Y  f(j| V \ j) + ∑_j ∈ Y \ X f(j| X).

Then m^f_X, 1(Y) ≥ f(Y) and m^f_X, 2(Y) ≥ f(Y), ∀ Y ⊆ V and m^f_X, 1(X) = m^f_X, 2(X) = f(X).





§ PRECOMPUTATIONAL COMPLEXITY MODEL

Given a specific submodular function f, we denote the precomputed statistics as p_X. This data structure depends on the specific submodular function, and stores information about a given set X. The idea here is that given p_X, computing the gains f(j | X, p_X) is often less costly compared to computing f(j | X) by evaluating f(X ∪ j) and f(X). This same idea holds for computing f(j | X \ j, p_X). 
The simplest precomputed statistic, which applies to any submodular function in the oracle model, is storing p_X = f(X). Evaluating the gains f(j | X, p_X) just requires one oracle call as opposed to two, when computing f(j | X). 
This simple trick can matter a lot since many of the algorithms require computing many gains (i.e., f(j | X) for many j ∈ V). For example, given the statistic p_X, the modular upper bounds can be easily computed since they depend on f(j | X), ∀ j ∉ X and f(j | X \ j), ∀  j ∈ X. Given a general submodular function in the value oracle model, this gives a factor 2 speedup. Many real world submodular functions, however, have richer statistics, often enabling speedups up to a factor of O(n) or even O(n^2). The same idea holds in computing modular lower bounds and extreme points of the submodular polyhedron. In this case, we are given a permutation σ, and we compute a chain of gains f(σ(i) | S^σ_i-1), ∀ i ∈ V. Often it is also much easier to update the statistic when adding an element to a set – i.e., given a set X and its statistic p_X, one can easily compute p_X ∪ j using p_X (i.e., without needing to compute it from scratch). Similarly, we can “downdate” the statistics, i.e., given the statistics p_X, we would like to compute the statistics for p_X \ j.

Most submodular function optimization algorithms either use the modular upper bounds (in which case, they compute f(j | X) repeatedly), or they compute extreme points of the submodular polyhedron (in the form of subgradients, or via the  extension), or they greedily add or remove elements. These algorithms easily admit our memoization framework.

Denote 𝐓^𝐨_𝐟 as the oracle complexity of the submodular function, a quantity which most existing algorithms use for their analysis. Define 𝐓^𝐠_𝐟 as the complexity of evaluating the gains, given the precomputed statistics. This depends in general on whether we are adding or removing the item (i.e., for f(j | X, p_X) or f(j
| X \ j, p_X)). For submodular functions we consider here, however, both complexities are the same. Also denote 𝐓^𝐮_𝐟 as the complexity of updating the precomputed statistics. Again, the complexity of updating or downdating could be different in general, but they are the same for the submodular functions we consider here. Finally, denote 𝐓^𝐩_𝐟 as the complexity of computing the precomputed statistic for a set X from scratch.

Under the value oracle model, the complexity of evaluating the modular upper bounds is O(nT^o_f), since the modular upper bound requires n oracle queries. In the precomputational model, the complexity is O(nT^g_f + T^p_f), since we would need to compute the precomputed statistic p_X (in the worst case) from scratch, and then compute f(j | X, p_X). For most submodular functions, T^p_f is roughly the same order as T^o_f, but T^g_f is often at least a factor n cheaper. Similarly, the complexity of computing the modular lower bound is O(nT^o_f) in the value oracle model, but is O(n[T^g_f + T^u_f]) in the precomputational model. Again for almost all submodular functions, T^g_f and T^u_f are at least a factor n cheaper than T^o_f. Hence the precomputational model provides significant speedups to the algorithms, which can be very important in large scale machine learning problems.

We make the ideas above more concrete in the following sections, by first explicitly defining p_X and the procedure for updating p_X, for several classes of submodular functions which occur in applications. Table 1 summarizes the precompute statistics p_f(X), the complexity T_f^o/T_f^p and T_f^u/T_f^g for different functions f.



 §.§ Graph Based Submodular Functions

Facility Location Functions: Given a similarity matrix {s_ij}_i, j ∈ V the facility location function is f(X) = ∑_i ∈ Vmax_j ∈ X s_ij. This function has successfully been used in summarization and data subset selection <cit.>. It is easy to check that the oracle complexity T_f^o = O(n^2). The precomputed statistics in this case is p_X[i] is the pair of the largest and second largest value of s_ij: i ∈ V, j ∈ X. Given these statistics however computing the gains are much easier, since f(k | X, p_X) = ∑_i ∈ Vmax{0, s_ik - max_j ∈ X s_ij}, and hence T^g_f = O(n). We can similarly efficiently compute the gains f(k | X \ k, p_X) in O(n). Moreover, updating the statistics is also easy, since max_j ∈ X ∪ k s_ij = max{s_ik, max_j ∈ X s_ij} – we can similarly update the second largest. Hence T^u_f = O(n). Moreover, computing p_X is the same complexity of computing f, and T_f^p = T_f^o. We can easily extend this idea to the top-k facility location function, where instead of taking the max, we can take the top k similarities. In that case, the precomputed statistics would be a matrix p_X[i, l], i ∈ V, l ∈{1, ⋯, k} as the lth largest value of s_ij for a given i. It is easy to see that we obtain an O(n) speedup in this case as well.

Saturated Coverage functions:  The saturated coverage function, f(X) = ∑_i ∈ Vmin{∑_j ∈ X s_ij, α_i}, has successfully been used in document summarization <cit.>. The oracle complexity of this function is also T_f^o = O(n^2). A natural choice of the precomputed statistics in this case is p_X[i] = ∑_j ∈ X s_ij, i ∈ V. Given this, it is easy to compute f(j | X, p_X) in O(n) time. Moreover, updating p_X can also be done directly in O(n), since p_X ∪ k[i] = p_X[i] + s_ik. Hence in this case also, T^g_f = T^u_f = O(n) and T^p_f = O(n^2). 

Graph Cut like functions: This class of functions have been used extensively in both summarization problems, while modeling coverage and diversity <cit.>, as well in image segmentation and denoising, by capturing cooperation <cit.>. We can denote the general class as f(X) = λ∑_i ∈ V∑_j ∈ X s_ij - ∑_i, j ∈ X s_ij. Setting λ = 1 is the standard graph cut, while λ = 0 gives the redundancy penalty <cit.>. Again, the natural choice of p_X[i] = ∑_j ∈ X s_ij. Though T_f^o = O(n^2), T_f^g = T^u_f are both O(n). Similarly, T_f^p = O(n^2).-1



 §.§ Coverage Functions

Set Covers and Neighborhood Functions: This is another important function,
capturing the notion of coverage <cit.> in maximization problems. This function also captures notions of complexity (like the size of the vocabulary in a speech corpus) in minimization problems <cit.>. Given a
set of sets {𝒮_1, ⋯, S_n} and the universe 𝒰 = ∪_i 𝒮_i, define f(X) = w(∪_i ∈ X𝒮_i), where w_j denotes the weight of item j ∈𝒰. This
setup can alternatively be expressed via a neighborhood function
Γ: 2^V → 2^𝒰 such that Γ(X) =
∪_i ∈ X𝒮_i. The oracle complexity of the function in this case is T^o_f = O(|𝒰|). The precomputed statistics here is a vector p_X ∈ℝ^𝒰 which stores the number of times u ∈𝒰 appears in the sets {𝒮_i}_i ∈ X. Then, f(j | X, p_X)  = |{u ∈𝒮_j: p_X[u] = 0}|, and f(j | X \ j, p_X) = |{u ∈𝒮_j: p_X[u] = 1}|. Similarly, p_X ∪ j can be updated as p_X ∪ j[u] = p_X[u] + 1, for u ∈𝒮_j. It is easy to see that both T_f^u = T_f^g = O(1), and T^p_f = T^u_f (assuming |S_j| is a constant, which is often the case).

Clustered Set Cover:  We can generalize this idea to clustered set cover functions, often used as confusability functions in corpus selection <cit.>. This is defined as f(X) = ∑_i = 1^k w(Γ(X) ∩ C_i), where C_1, C_2, ⋯ C_k ⊆𝒰 are clusters. The precomputed statistics here is a vector of sets. For each j = 1, ⋯, k, p_X[j] = Γ(X) ∩ C_j.

Probabilistic Coverage Functions: Another generalization of the set cover function, which has been used in a
number of models for summarization problems <cit.>. This
provides a probabilistic notion to the set cover function, and is defined as f(X) = ∑_i ∈𝒰 w_i [1 - ∏_j ∈ X (1
- p_ij)]. The complexity of evaluating this function T^o_f = O(n|𝒰|). The precomputed statistics in this case is p_X[i] = ∏_j ∈ X (1
- p_ij). Note that here, both T^u_f = T^g_f = O(|𝒰|), thereby providing a factor n speedup. Similarly, T^p_f = T^u_f.



 §.§ Feature Based Functions and Clustered Concave over Modular Functions

Another class of submodular functions are sums of concave over modular functions. They appear in maximization problems as feature based functions, defined as f(X) = ∑_e ∈ℱψ(m_e(X)), and have been used in data subset selection applications <cit.>.  m_e(j) captures how much item j covers feature e ∈ℱ. Another related function is f(X) = ∑_j = 1^k ψ(m_j(X ∩ C_j)), where C_1, C_2, ⋯, C_k are clusters of similar items in the ground set V. This function simultaneously captures diversity in maximization problems <cit.>, and notions of cooperation in minimization problems <cit.>. The complexity of evaluating these functions is O(n |ℱ|) and O(nk) respectively. A natural choice of the precomputed statistics is p_X[e] = m_e(X), in the case of feature based functions, and p_X[j] = m_j(X ∩ C_j) for clustered concave over modular functions. Again, it is easy to see that T^u_f = T^g_f are O(|ℱ|) and O(k) respectively, thus saving a factor of n. Moreover, in this case we also have that T^p_f = T^o_f.-1



 §.§ Spectral Submodular Functions

Another rich class of submodular functions, defined as f(X) = ∑_i = 1^|X|ψ(λ_i(S_X)), where S is a PSD matrix, S_X represents the principal submatrix formed by the rows and columns corresponding to X, and ψ is a concave function. This function is submodular for a large class of concave functions <cit.>. This class of spectral regularizers has been shown to promote diversity <cit.>, and includes as special cases the log det function, which occurs in the context of determinantal point processes <cit.>, since f(X) = log(S_X) = ∑_i = 1^|X|logλ_i(S_X). Another example of this function is f(X) = ∑_i = 1^|X|√(λ_i(S_X)). Evaluating this class of functions is T^o_f = O(n^3), since we need to perform the eigenvalue decomposition. A natural choice of the precomputed p_X statistic here is the eigenvalue decomposition of S_X. Using the result from <cit.>, it is possible to update (or downdate) the eigenvalue decomposition of S_X to S_X ∪ j (or S_X \ J) in O(n^2) computations given the eigenvalue decomposition of S_X (note that the result of <cit.> is in terms adding a single row and column to S_X. However, converting S_X to S_X ∪ j is the same as adding a row, followed by a column, and hence two updates.)-1



 §.§ Dispersion Functions

Denote d_ij as a distance measure between objects i and j. Define
the `Dispersion Min function as f(X) = min_i, j ∈ X d_ij. This function is not submodular, but can be efficiently optimized via a greedy algorithm <cit.>. It is easy to see that maximizing this function involves obtaining a subset with maximal minimum pairwise distance, thereby ensuring a diverse subset. Similarly, we can define two more variants. One is
the supermodular Dispersion Sum, defined as f(X) = ∑_i, j ∈ X d_ij. Another
is Dispersion Min-Sum, a combination of two forms, defined as f(X) = ∑_i ∈ Xmin_j ∈ X d_ij. This function is submodular <cit.>. 



 §.§ Mutual Information and Entropy

This is another class of functions, used often for feature subset
selection <cit.>.  The entropy function f(A) = H(X_A)
is submodular, and while mutual information f(A) = I(X_A; C) is
always a difference of submodular functions
<cit.>, it is also sometimes submodular
under some assumptions.  While both these functions require
exponential complexity to evaluate, they can be estimated easily from
data via a single sweep, and using techniques like Laplacian
smoothing <cit.> – the computational complexity being
O(n|𝒟|), where |𝒟| is the size of the training
data. The mutual information and entropy estimates also amenable to
precomputation, since we precompute the data tables for the given set
of feature X_A. Adding and removing features to this corresponds to
further dividing tables, which can be done in O(|𝒟|)
complexity.



 §.§ Influence Maximization

A number of models for influence maximization have been shown to be related to submodular maximization <cit.>. Evaluating the objective function, however, requires MCMC simulations, which is quite expensive. We can define precomputed statistics however, which can significantly speedup the greedy algorithms <cit.>. In particular, they avoid the expensive MCMC by enumerating simple paths, and relying on memoization and look ahead optimization <cit.>.



 §.§ Mixtures of Submodular Functions

Often it is desirable to consider not just one submodular function, but a mixture of many submodular bases functions <cit.>. In particular, we often express f(X) = ∑_i = 1^m w_i f_i(X), where f_i's are bases submodular functions (like one of the submodular functions above), and w_i's are weights. Assuming the submodular functions f_i's have precomputed statistics, p^f_i_X, the gain f(j | X, p^f_X) is exactly ∑_i = 1^m w_i f_i(j | X, p^f_i_X). Moreover, to update the statistics p^f_X, we simply update the individual p^f_i_X.-1




 §.§ Deep Submodular Functions

Almost all of the above classes of submodular functions are subsumed
by Deep Submodular Functions <cit.>, defined as nested
sums of concave over fewer-layer deep submodular functions. The
general form of Deep SFs are:

    f(X) = ∑_i_1 ∈ℱ_1 w^1_i_1ψ_1(⋯ψ_k-1(∑_i_k-1∈ℱ_k w^k_i^kψ_k (m^k_i_k(X))))

For simplicity we consider the case with k = 2, i.e., a two layer function. The complexity of evaluating f is |ℱ_1||ℱ_2|n. Similar to a feature based function, the precomputed statistics here is p_X = m^2_i_2(X), ∀ i_2 ∈ℱ_2. The complexity of updating the precomputed statistics is T_f^u = |ℱ_2| while the complexity of computing the gain T_f^g = |ℱ_1||ℱ_2|. Both these quantities are a factor n less expensive compared to T_f^o.






§ ALGORITHMS FOR SUBMODULAR OPTIMIZATION

We now investigate several known algorithms for submodular optimization problems, and show how we can easily integrate the precomputation idea into them. In almost all cases, we shall see that this entails only a few additional lines of code, while providing significant speedups in applications. Table 2 summarizes the complexity of various submodular optimization algorithms with the precomputational model and value oracle model.-1



 §.§ Computing Subgradients, Supergradients, and extreme points of Submodular Polyhedron

Most submodular optimization algorithms either rely on computing subgradients, supergradients or some extreme points of the submodular polyhedron. So we first compare the complexity of computing these quantities under the Precomputational Model and the Value Oracle Model. Recall that computing a supergradient of X requires computing f(j | X) for every j ∉ X (or equivalently computing f(j | X \ j) for j ∈ X). The complexity of doing this in the Value Oracle Model is O(nT_f^o). Under the Precomputational Model, the complexity if O(T_f^p + nT_f^g) since we can compute the precompute statistics p_X and using that, evaluate the gains f(j | X, p_X). This is a factor n speedup since in all cases, T_f^g is a factor n cheaper compared to T_f^0 (see Table 1). Similarly, computing a subgradient (or equivalently computing an extreme point of the submodular polyhedron) requires forming a chain of sets ∅ = X_0 ⊆ X_1 ⊆⋯⊆ X_n = V, and computing f(x_i | X_i) where x_i = X_i+1\ X_i. Computing the subgradient (or extreme point) using the Value Oracle Model is O(nT_f^o). Using the precompute statistics, we can at every step use the precompute statistics from X_i to compute f(x_i | X_i) and then update the precompute statistics. Since this is a greedy algorithm, the complexity of this is O(n[T_f^g + T_f^u]. Again, from Table 1, it is evident that we can achieve a speedup at least of a factor of n.
 


 §.§ Submodular Maximization

Submodular maximization is particularly important in applications like summarization <cit.>, data subset selection <cit.> etc. where we want to find diverse and relevant subsets.
A large class of existing submodular maximization algorithms can be expressed via a common minorization-maximization framework <cit.> –  an iterative procedure which optimizes the modular lower bound h^f_σ^t(X) (which is tight at X^t). This algorithm essentially chooses a sequence of orderings σ^t, each of which is tight with respect to the set X^t, and different known algorithms use different orderings. The simplest algorithm is to just choose a random subgradient (or ordering) σ^t at every iteration. We can compute the subgradients h^f_σ^t using memoization – start with the empty set, and compute f(σ^t(i) | S^σ^t_i-1, p_S^σ^t_i-1), and update p_S^σ^t_i, for i = 1, 2, ⋯, n. The complexity of this algorithm using the precompute statistics is Õ(n[T_f^g + T_f^u]) – where Õ hides the complexity of the outer loop which is weakly polynomial <cit.> (in practice, it is a constant). With the value oracle model, the complexity of minorize-maximize is Õ(n[T_f^o]).

While this simple algorithm provides guarantees, one gets much tighter bounds with more intelligent choices of subgradients. Below we look at a few such algorithms for various variants of submodular maximization.-1

Greedy and Lazy Greedy Algorithm:  This is a common algorithmic idea, which provides constant factor guarantees for a large class of monotone submodular maximization problems, under cardinality, knapsack and matroid constraints <cit.>. Starting with X^0 = ∅, we sequentially update X^t+1 = _j ∈ V \ X^t f(j | X^t). The complexity of this algorithm is O(nk T_f^o). We can easily integrate precomputation into this, by setting the update rule to X^t+1 = _j ∈ V \ X^t f(j | X^t, p_X^t), and updating p_X^t+1. The complexity then is essentially O(nk [T_f^u + T_f^g]). Thanks to submodularity, however, we can significantly accelerate this algorithm, to what is known as the lazy greedy algorithm <cit.>. The idea is that instead of recomputing f(j | X^t), ∀ j ∉ ^t, we maintain a priority queue of sorted gains ρ(j), ∀ j ∈ V. Initially ρ(j) is set to f(j), ∀ j ∈ V. The algorithm selects an element j ∉ X^t, if ρ(j) ≥ f(j | X^t), we add j to X^t (thank to submodularity). If ρ(j) ≤ f(j | X^t), we update ρ(j) to f(j | X^t) and resort the priority queue. The complexity of this algorithm is roughly  O(k n_R T_f^o), where n_R is the average number of resorts in each iteration. Note that n_R ≤ n, while in practice, it is a constant thus offering almost a factor n speedup compared to the simple greedy algorithm. We can use the notion of precompute here too, and use f(j | X^t, p_X^t) in place of the gain f(j | X^t). Note that additionally, whenever we add an element j to X^t, we also need to update the precomputed statistic. The resulting complexity in precomputational model is O(k n_R [T^g_f + T^u_f]), again providing a factor n speedup.-1

Lazier than Lazy greedy Algorithm:  While the lazy greedy algorithm above runs much faster in practice, in the worst case, the complexity is the same as the naïve greedy algorithm. The Lazier than Lazy greedy algorithm <cit.> attempts to obtain an approximation guarantee of 1 - 1/e - ϵ in O(n log(1/ϵ) function evaluations rather than O(nk) from the greedy (or lazy greedy) algorithm. The idea is to select a random set R of size n/k log(1/ϵ) at select the element with the largest gain of adding the element to the current set. We run this until our chosen subset has k elements (in the cardinality constrained case). Since we evaluate the gains f(j | X), j ∈ R ⊆ V \ X. We start with the X = ∅, and at every step we can compute f(j | X, p_X) using the precompute statistics, and after choosing the best element, add that element to X and update p_X. The complexity of this algorithm under the precompute model is O(n log(1/ϵ) [T_f^g + T_f^u]), while using the value oracle model, the complexity is O[n log(1/ϵ) T_f^o].

Distributed Greedy Algorithm:  The Distributed Greedy algorithm <cit.> attempts to extend the Greedy algorithm to a setting where all the data cannot fit into memory. The basic idea is we have m machines, where we equally partition the data into, i.e., V_1, ⋯, V_m. The distributed greedy algorithm then runs a greedy algorithm on each of the partitions to obtain k elements, followed by a second round of greedy on the mk elements to obtain k elements. It is easy to see that the complexity of this in the value oracle model is O([nk/m + mk^2]T_f^o). Moreover, since we run a two-stage greedy algorithm, the memoization can be used exactly like the memoization discussed above.

Sieve Streaming Algorithm:  The Sieve streaming algorithm basically performs submodular maximization in the streaming setting <cit.>. The idea of the algorithm is to maintain a set of thresholds in a set O, and corresponding to each threshold, maintain a set S. It simultaneously grows these sets depending on constraints (see <cit.> for details on the algorithm). The way we can incorporate memoization is by storing |O| copies of submodular functions, with their own precompute statistics, which are updated through the process of the algorithm. Given these precompute statistics, we can easily compute the gains f(e_i | S_i) for each of the sets maintained by the algorithms, and once the elements are added to the sets, we can update the precompute statistics. The complexity analysis (shown in Table 1) follows directly from the results in <cit.> and the precompute statistics (note that each subset is simultaneously updated through the course of the algorithm in a greedy manner).  

Local Search Algorithm:  This algorithm for unconstrained submodular maximization (USM), essentially runs multiple rounds of the greedy algorithm <cit.>, and provides a 1/3 approximation for USM. In particular, we start with X^0 = ∅, and run the forward greedy algorithm until we can no longer add elements, followed by the reverse greedy of removing elements. We continue this procedure until we converge to a local optimum. The forward greedy algorithm is essentially the same algorithm as above, while in the reverse greedy case, we remove elements with the smallest value of f(j | X^t \ j). We can use the precompute ideas and use f(j | X^t \ j, p_X^t) (in this case, we downdate the statistics after removing the elements).  The complexity of the algorithm follows from the complexity of the local search from <cit.> – see Table 2 for the complete expression. Again, we simply need to update and downdate the precompute statistics, and given these statistics, compute the gains of adding and removing the elements. For all the functions we consider, the complexity of adding/removing elements (correspondingly updating the downdating the precompute statistics) is the same, the complexity analysis follows.

Bidirectional Greedy Algorithm:  The bidirectional greedy algorithm <cit.> provides the tight 1/2 approximation for USM. Surprisingly, this is a very simple linear time algorithm. This algorithm maintains two sets A and B (initially set to A = ∅ and B = V) which increases and decreases respectively in the course of the algorithm, and depends on an initial ordering π. Then from i = 1, 2, ⋯, n, we either add π(i) to A or remove π(i) from B, depending on which of the gains is larger. The complexity of this algorithm is O(n T_f^o). We can use precomputation here, by storing two sets of statistics (one for the set A, and another for B) – in practice, we can achieve this by maintaining two submodular functions which each store their statistics. As A grows, we update its statistics and similarly downdate B's statistics as it shrinks. The complexity with precomputation is O(n [T^g_f + T^u_f]), which is in practice a factor n faster.-1


Randomized Greedy Algorithm:  
The randomized greedy algorithm provides an efficient algorithmic framework for cardinality constrained non-monotone submodular functions. The idea of this algorithm is very similar to the bidirectional greedy above, except that instead of choosing the best gain f(j | X^t), we choose at random, one of the top k gains (where k is the given cardinality constraint). The complexity of this algorithm is O(nk T^o_f). Similar to the greedy algorithm, we can incorporate the precomputation by using f(j | X^t, p_X^t), and updating p_X^t when we add the new element. The complexity of this is O(nk [T^u_f + T^g_f]).-1



 §.§ Submodular Minimization

Submodular minimization comes up in applications where we want to minimize cooperative costs and complexity measures, like image segmentation <cit.>, and limited vocabulary corpus selection <cit.>.
We investigate three important submodular minimization algorithms, and show how the notions of precomputation yields significant computational gains.

The minimum norm point algorithm:  This is one of the most practical algorithm available for general purpose submodular minimization <cit.>. This algorithm solves quadratic programming problem, which is equivalent to of the discrete minimization problem <cit.>, and uses the Wolfe algorithm <cit.>. One of the most important step in this algorithm (and the only step which requires oracle access to the submodular function) is the greedy algorithm for solving a linear programming problem max_x ∈𝒫_f⟨ x, x̂⟩.  As discussed in section <ref>, solving this problem requires computing the subgradient according to the ordering σ_x̂, and we can use precomputations to efficiently find these subgradients, thereby providing significant speedups in practice. While in practice the Minimum Norm Point algorithm is the fastest (compared to the other combinatorial algorithms), the worst case complexity is still high order polynomial <cit.>. The worst case complexity with and without memoization is in Table 2. Having said that, we demonstrate in our experimental section, that memoization can provide significant compute gains.

 extension based:  Another class of algorithms <cit.>, is based on relaxing the discrete optimization problem to a continuous one, via the  extension. This procedure, moreover, works for a large class of constrained problems, and uses convex optimization techniques <cit.>. The precomputational ideas apply in these cases too, since we can both compute the  extension (which requires solving the linear program over the submodular polyhedron) and its subgradient (which is the same as the subgradient of the submodular function) very efficiently.

Majorization-Minimization:  This is a discrete gradient based framework which applies to a large class of constrained submodular minimization problems <cit.>. This is an iterative procedure, which starts with X^0 = ∅, and minimizes the modular upper bounds m^f_X^t as a proxy to f. Each step of this algorithm is a linear cost problem, which is easy for many combinatorial constraints.  Thanks to the nature of this algorithm, we are guaranteed improvement at every iteration. Moreover, this algorithm also admits guarantees, and works very well in practice <cit.>. Moreover, as we saw earlier, we can efficiently compute the modular upper bounds using precomputations. The complexity under this model is O(nT^g_f + T^p_f), which is in general a factor n faster than the oracle model O(n T^o_f).



 §.§ Optimization with Submodular Constraints

Another class of optimization problems related to submodular functions, are ones where submodular functions appear as upper or lower bound constraints. Two general problem classes here are <cit.> (a) SCSC: min{f(X) | g(X) ≥ c}, and (b) SCSK: max{g(X) | f(X) ≤ b}. This class of problems comes in applications where we simultaneously want to maximize one submodular function, while minimizing another. 
This generalizes a number of useful problems, including for example, the submodular set cover <cit.> and the submodular knapsack <cit.>, which are instances of SCSC and SCSK respectively, when f(X) = w(X) is a modular function and g(X) is submodular function. The same lazy greedy algorithm actually works for the submodular set cover and the submodular knapsack problem, and hence precomputation directly applies in this case. The general case of SCSC and SCSK, can be handled by replacing f by its modular upper bound and iteratively solving submodular set cover and submodular knapsack respectively <cit.>. Moreover, since the modular upper bound computation and the greedy algorithm can be both done efficiently via precomputations, we can obtain significant speedups for these problems in large scale problems.



 §.§ Difference of Submodular Functions

A very common and general optimization problem involves minimizing the difference between submodular functions min_X ⊆ V f(X) - g(X), and comes up in several machine learning applications including feature subset selection, and graphical model inference <cit.>. A common class of heuristics for this problem is the submodular-supermodular procedure, and its variants <cit.>. They are essentially majorization-minimization based iterative procedures starting with X^0 = ∅, and iteratively to replace either f by its modular upper bound m^f_X^t, and g by its modular lower bound h^g_X^t, or both. At every iteration, the resulting problem is either submodular minimization, maximization or modular minimization <cit.>. Since both the upper and lower bounds can be efficiently computed, and the algorithms for submodular minimization and maximization are efficient, thanks to precomputations, we can achieve substantial speedups in applications.



§ EXPERIMENTS
 
In this section, we compare the performance of different submodular optimization problems (functions and algorithms). We study two specific applications: 1) speech data subset selection <cit.>, and 2) low complexity speech corpus creation <cit.>.  We divide the sections below by the different submodular optimization problems. We denote the precompute model as PM and the value oracle model as VO. 

Computing Sub- and Super-gradients:
We first compare the running time of computing the sub and super gradients for the speech data subset selection problem <cit.> on TIMIT. In this case, n = 4620. We compare the running time of three classes of submodular functions: Facility Location, Feature Based, and the Complexity function from <cit.> (which is a set cover function). The Results (Table 3) shows substantial gains using memoization compared to the VO model for all three classes of functions. We also compute the sub and supergradients for large scale problems (|V| = 200000 from Switchboard <cit.>. In this case, the subgradient computation takes 2.8 seconds for the Feature based function and 0.7 seconds for the Complexity function, while the supergradient compute takes 3.09 seconds for the Feat based and 1.14 sec for the Complexity Function (under the PM model). With the VO models, it would have taken around 7 days to compute these. 




Submodular Maximization
Next we compare the different functions on submodular maximization with the lazy greedy algorithm. Again, we see substantial speedups using the precompute model across the board for different submodular functions (Table 4). We also compare the Feature Based and Set Cover function for large scale submodular optimization |V| = 200000. The FB function takes 16 seconds, while the Set Cover function takes about 27 seconds for running the Lazy greedy algorithm. 

Submodular Minimization
We next compare the complexity of unconstrained and constrained submodular minimization. We define f(X) = c(X) - λ |X|, where c(X) is the complexity function from <cit.>. We use the TIMIT dataset (|V| = 4620. For unconstrained minimization, the minimum norm point algorithm takes about 7.2 seconds with PM, while under the VO Model, it takes around 2000 seconds. For larger scale problem (n = 200000), the MN algorithm takes around 88 seconds under the PM Model (the VO model will take several days to complete). Consider next the problem of constrained submodular minimization under cardinality constraints (minimizing c(X) subject to a cardinality constraint). We use the Majorization-Minimization algorithm here <cit.>. Under PM, majorization-minimization takes around 0.13 seconds for a budget of 10%, while the value oracle model, the time is 163 seconds. For the large scale version of this problem with Switchboard (n = 200000), MMin takes around 23 seconds with the PM Model.



§ CONCLUSIONS

This paper introduces the idea of precompute Statistics and Memoization for Submodular Optimization. We show how several real world submodular functions admit natural precompute statistics, and how we can integrate this idea into a large family of algorithms for submodular maximization, minimization and other forms of constrained submodular programs. We empirically demonstrate the utility of our Memoization framework on several large scale problems.

This material is based upon work supported by the National Science
Foundation under Grant No. (IIS-1162606), and a Google and a Microsoft
award. This work was supported in part by the CONIX Research Center,
one of six centers in JUMP, a Semiconductor Research Corporation (SRC)
program sponsored by DARPA.

abbrv

