Near Optimal Algorithms for Hard Submodular Programs with Discounted Cooperative Costs
Rishabh Iyer Jeff Bilmes


[
Near Optimal Algorithms for Hard Submodular Programs with Discounted Cooperative Costs

 Rishabh Iyer Jeff Bilmes 

 Microsoft Corporation University of Washington, Seattle
]



  In this paper, we investigate a class of submodular problems which in general are very hard. These include minimizing a
  submodular cost function under combinatorial constraints, which
  include cuts, matchings, paths, etc., optimizing a submodular
  function under submodular cover and submodular knapsack
  constraints, and minimizing a ratio of submodular functions. All these problems appear in several real world problems but have hardness factors of
  Ω(√(n)) for general submodular cost functions. We show
  how we can achieve constant approximation factors when we restrict
  the cost functions to low rank sums of concave over modular
  functions. A wide variety of machine
  learning applications are very naturally modeled via this subclass
  of submodular functions. Our work therefore provides a tighter connection
  between theory and practice by enabling theoretically satisfying
  guarantees for a rich class of expressible, natural, and useful submodular
  cost models. We empirically demonstrate the utility of our models on
  real world problems of cooperative image matching and sensor
  placement with cooperative costs.-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. A set function f: 2^V →ℝ over
a finite set V = {1, 2, …, n} is submodular 
<cit.> 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 perhaps more intuitive
characterization of submodularity is as follows:
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.

In this paper, we address the following a family of hard submodular optimization problems. The first one is constrained submodular minimization <cit.>:

    min{f(X) | X ∈𝒞}

where the function f is monotone submodular, and 𝒞 is a combinatorial constraint, which could represent a cardinality lower bound constraint, or more complicated ones like cuts, matchings, trees, or paths in a graph. With cut constraints, this problem becomes cooperative cuts <cit.>, and with matching constraints, we call this cooperative matchings, which we introduce and utilize in this paper.-1

The second problem asks for minimizing a monotone submodular cost function f, while simultaneously maximizing a monotone submodular coverage function g. A natural way to model this bi-optimization problem is to introduce one of f and g as a constraint <cit.>. In particular, we obtain two optimization problems:

    min{f(X)   |  g(X) ≥ c}, 
       max{ g(X)  |  f(X) ≤ b},


The fourth problem considered in this paper is minimizing the ratio of submodular functions <cit.>.

    min{f(X)/g(X)  |  ∅⊂ X ⊂ V}


A key assumption in this paper is that the functions f and g in Problems 1-4, are monotone submodular – an assumption that, as we shall see, is natural in many applications. Problem 2 is a special case of Problem 1, with 𝒞 = {X: g(X) ≥ c}. Furthermore, Problem 2 and Problem 3 are closely related and, loosely speaking, duals of each other <cit.>. Similarly, Problem 4 is closely related to Problems 2 and 3, in that given an approximation algorithm for Problems 2 or 3, we can obtain an approximation with similar guarantees for Problem 4 <cit.> (also considered in <cit.> with general monotone set functions). Problem 1 is constrained submodular minimization, while Problems 2, 3 and 4 try to simultaneously minimize one submodular function while maximizing another.

Problems 1-4 appear naturally in several machine learning applications. However, in the worst case all four problems have polynomial hardness factors of Ω(√(n)) <cit.>. An important observation is that the
polynomial hardness of problems 1 - 3, comes up mainly due to the
submodular cost function f – they do not depend as much on the
constraints 𝒞 or the submodular function
g <cit.>. In the case of Problem 4, the hardness depends on both f and g <cit.>.



On the other hand, these problems come up as models in many machine
learning applications. These lower bounds are specific to rather
contrived classes of functions, whereas much better results can be
achieved for many practically relevant cases. The pessimistic worst
case results are somewhat discouraging, d begging the need to
quantify sub-classes of submodular functions that are more amenable to these optimization problems. Only limited past work has focused on investigating these problems with potentially a subclass of submodular functions. <cit.> provide bounds for Problems 1-4 based on the notion of curvature, and argue how several submodular functions (e.g. clustered concave over modular functions) have bounded curvature. Their curvature bounds depend on the choice of the submodular functions, and in certain cases yield no improvement over the worst case bounds. For classes of functions with bounded curvature, their bounds yield improved results.

In this paper, we
focus on a tractable yet expressive subclass of submodular cost
functions f, namely low rank sums of concave over modular functions. 

Low rank sums of concave over modular functions are the class of functions representable as f(X) = ∑_i = 1^k ψ_i(w_i(X)), where ψ_is are monotone concave, and k is constant or O(log n). 

Low Rank in this context means that the number of components in
the sum is small (i.e., k is small).  Our use of the terminology “low
rank” is identical to that used in <cit.>.  We argue
how this subclass naturally models many interesting applications










of problems 1 - 4 in machine learning. We do not
need to consider the entire class of submodular functions (which
includes rather contrived instances), but only this subclass. This
observation helps us in providing better connections between theory
and practice. The main specialty of this subclass is that these
functions effectively model cooperation between objects via discounts
provided by concave functions. Moreover, we show that this subclass
admits fully polynomial time approximation schemes for Problem 1, and
constant factor approximation guarantees for Problems 2 and
3.  Similarly, we achieve constant factor approximation guarantees for Problem 4, when f is a low rank sum of concave over Modular functions, and g is an arbitrary submodular function, a significant improvement over <cit.>. The bounds we obtain are significantly better than the worst case bounds, and also an improvement over the bounds achieved using the curvature <cit.>.
-1

Low rank sums of concave over modular functions in Problems 1 - 4, fit as natural models in several machine learning problems. Below, we summarize some of these.


Image segmentation (Cooperative Cuts): Markov random fields
with pairwise attractive potentials occur naturally in modeling image
segmentation and related
applications <cit.>. While models are tractably
solved using graph-cuts, they suffer from the shrinking bias problem,
and images with elongated edges are not segmented properly.
When modeled via a submodular
function, however, the cost of a cut is not just the sum of the edge weights,
but a richer function that allows cooperation between edges, and
yields superior results on many challenging
tasks
(see,
for example, the results of the image
segmentations in <cit.>).
This was
achieved in <cit.> by partitioning the
set of edges ℰ of the grid graph into groups of similar
edges (or types) ℰ_1, ⋯, ℰ_k, and defining a
function f(S) = ∑_i = 1^k ψ_i(w(S ∩ℰ_i)), S
⊆ℰ, where ψ_is are concave functions and w
encodes the edge potentials. This ensures that we offer a
discount to edges of the same type. Moreover, the number of
types of edges are typically much smaller than the number of pixels,
so this is a low-rank sum of concave
functions. 












 Image Correspondence (Cooperative Matchings):  The simplest model for matching key-points in pairs of images (which is also called the correspondence problem) can be posed as a bipartite matching. These models, however, do not capture interaction between the pixels. We illustrate the difficulty of this in Figure <ref>.  One kind of desirable interaction is that similar or neighboring pixels be matched together. We can achieve this as follows. First we cluster the key-points in the two images into k groups (this is illustrated in Figure <ref>-left via green, blue and red key-points).  This induces a clustering of edges that can be given a discount via a submodular function (details are given in Section <ref>).


















In practice, the number of groups (k) can be much smaller than n
and this is a low-rank sum of concave over modular
functions. Figure <ref>-right shows how the
submodular matchings improves over the simple bipartite matching. In
particular, the minimum matching approach produces many spurious
matches between clusters (shown in red) that are avoided via the
cooperation described above.

Sensor Placement or Feature Selection: Often, the problem of
choosing sensor locations A from a given set of possible locations
V can be modeled <cit.> by maximizing
the mutual information between the chosen variables A and the
unchosen set V \ A (i.e.,
g(A) = I(X_A; X_V \ A)). Alternatively, we may wish to
maximize the mutual information between a set of chosen sensors X_A
and a quantity of interest C (i.e., g(A) = I(X_A ; C)) assuming
that the set of features X_A are conditionally independent given
C <cit.>. Both these functions are submodular. Since
there are costs involved, we want to simultaneously minimize the cost
f(A). Often this cost is
submodular <cit.>, since there is
typically a discount when purchasing sensors in bulk (or computing
features), and we can express this via Problems 2 and 3. For example,
there may be diminished cost for placing a sensor in a particular
location given placement in certain other locations. Similarly,
certain features might be cheaper to use given that others are already
being computed (e.g., those that use an FFT). A natural cost model in
such cases is f(A) = ∑_i=1^k ψ_i( m(A ∩ S_i) ) where
ψ_i's are concave, m(j) is the cost of sensor (or feature) j
and S_1, ⋯, S_k are groups of similar sensors or features.
Typically, k is much smaller than n and this can be expressed as
low rank sum of concave over modular functions.



§ BACKGROUND & EXISTING ALGORITHMS



The basic idea for most combinatorial algorithms solving Problems 1 - 4, are based on approximating the cost function f with a tractable surrogate function f̂ <cit.>. Moreover, all four problems have similar guarantees. We characterize the quality of the solution via the notion of approximation factors. In particular, we say that an algorithm achieves an approximation factor of α≥ 1 for Problem 1, if we can obtain a set X̂ such that f(X̂) ≤α f(X^*), where X^* is the optimizer of Problem 1. For Problems 2 and 3, we use the notion of bi-criterion approximation factors. An algorithm is a [σ, ρ] bi-criterion algorithm for Problem 2 if it is guaranteed to obtain a set X̂ such that f(X̂) ≤σ f(X^*) (approximate optimality) and g(X̂) ≥ρ c (approximate feasibility), where X^* is an optimizer of Problem 2.  Typically, σ≥ 1 and ρ≤ 1.  Similarly, an algorithm is a [ρ,σ] bi-criterion algorithm for Problem 3 if it is guaranteed to obtain a set X̂ such that g(X̂) ≥ρ g(X^*) and f(X̂) ≤σ b, where X^* is the optimizer of Problem 3. 
Moreover, problems 2 and 3 are very closely related <cit.>, in that an approximation algorithm for one problem can be used to obtain guarantees for the other problem. The two problems also have matching hardness factors. For Problem 4, we study an algorithm which achieve α-approximation guarantees, in that we can achieve a set X̂ such that h(X̂) ≤α h(X^*) where h(X) = f(X)/g(X) and X^* is the optimal minimizer of h. -1

Supergradient based Algorithm (SGA): One such method uses the supergradients of a submodular function <cit.> to obtain modular upper bounds in an iterative manner. In particular, define a modular upper bound:

    m^f_X(Y) ≜ f(X) - ∑_j ∈ X \ Y  f(j| V \ j) + ∑_j ∈ Y \ X f(j| X) ≥ f(Y)
 
The algorithm starts with the X^0 = ∅ and sequentially sets X^i+1 as the solution of the corresponding problem (1, 2 or 3) with a surrogate function as f̂(X) = m^f_X^i(X) <cit.>. In each case, this subproblem is much easier. For example, in the case of Problem 1, the subproblem becomes, 

    X^i+1 = min{m^f_X^i(X) | X ∈𝒞},
 
which is a linear cost problem, poly-time solvable for many constraints, like cardinality, cuts, matchings, paths etc. 

In the case of Problems 2 and 3, these subproblems are 

    X^i+1 = min{m^f_X^i(X) | g(X) ≥ c}
    
    X^i+1 = max{g(X) | m^f_X^i(X) ≤ b},

which are the submodular set cover and the submodular knapsack problems respectively <cit.>, and are constant factor approximable to a factor of 1 - 1/e. 

With Problem 4, the subproblem becomes,

    X^i+1 = min{m^f_X^i(X)/g(X)  |, ∅⊂ X ⊂ V}

This can be approximated up to a factor of e/(e-1) via a Greedy algorithm <cit.>.


Define α_f(X^*) = |X^*|/1 + (|X^*| - 1)(1 - κ̂_̂f̂(X^*))≤min{|X^*|, 1/1 - κ̂_̂f̂(X^*)}, where κ̂_̂f̂(X) = 1 - ∑_j ∈ X f(j | X \ j)/∑_j ∈ X f(j) represents the average curvature of the function f.
The supergradient based iterative algorithm (SGA) achieves an approximation factor of α_f(X^*) for Problem 1, and bicriteria factors satisfying σ = α_f(X^*) and ρ = 1 - 1/e for Problems 2 and 3. Finally, SGA achieves an approximation factor of e/(e-1)*α_f(X^*) for Problem 4.

This Lemma follows easily from the results in <cit.>. We can also achieve a non-bicriteria approximation factor for Problem 2, which is worse than the bicriteria factor by a log factor <cit.>. A key quantity which defines the approximation factor above is the average curvature κ̂_̂f̂(X^*), which in turn depends on the concave functions. If the concave function is ψ_i(x) = x^a, a ∈ (0, 1), SGA admits approximation factors of O(|X^*|^1-a) <cit.>. On the other hand, if the concave function is ψ(x) = log(1+ x), the guarantees are O(|X^*|), which is much poorer. 

The supergradient based algorithm is easy to implement, and also works well in practice <cit.>. For the general class of submodular functions, these results are close to the optimal bounds, and are, in fact, tight for some constraints. Nevertheless, the worst case guarantees seem discouraging, particularly for the class of low rank sums of concave over modular functions that we consider here, and that as mentioned above are natural 
for many applications.

Ellipsoidal Approximation based Algorithm (EA): Another generic approximation of a submodular function, introduced by Goemans et.  al <cit.>, is based on approximating the submodular polyhedron by an ellipsoid. The main result states that any polymatroid (monotone submodular) function f,
can be approximated by a function of the form √(w^f(X)) for a
certain modular weight vector w^f ∈ℝ^V, such that √(w^f(X))≤ f(X) ≤
O(√(n)logn) √(w^f(X)), ∀ X ⊆ V.
A simple trick then provides a curvature-dependent approximation <cit.>. We have the following result borrowed from <cit.>.

Define α = O(√(n)/1 + (√(n) - 1)(1 - κ_f)), where κ_f = 1 - min_j ∈ X f(j | V \ j)/min_j ∈ X f(j) represents the worst case curvature of the function f.. The Ellipsoidal Approximation based algorithm (EA) achieves an approximation factor of α  for Problem 1, and bicriteria factors satisfying σ = α and ρ = 1 - 1/e for Problems 2 and 3. Similarly EA achieves an approximation guarantee of eα/(e-1) for Problem 4.

The Ellipsoidal Approximation obtains the tightest bounds for Problems 1-4 <cit.>. This is again for the general class of submodular functions and the worst case factor of O(√(n)) is quite discouraging. This algorithm, however is very expensive computationally, and is not practical for solving machine learning applications <cit.>.



§ IMPROVED ALGORITHMS FOR LOW-RANK SUMS OF CONCAVE-MODULAR FUNCTIONS
 




Our main new results are that we can achieve a fully polynomial
time approximation scheme for Problem 1, and constant factor
approximation guarantees for Problems 2, 3 and 4 when the cost function
f is a low rank sum of concave over modular functions
(Theorem <ref>). Our
techniques build on recent methods used for minimizing quasi-concave
functions over solvable polytopes <cit.>.

Assume the concave functions ψ_i's are monotone functions, i.e., ψ_i(y) ≤ψ_i(y^'), ∀ y ≤ y^'. We also assume that for all i, ψ_i(ky) ≤ k^c ψ_i(y) for k ≥ 1, y ≥ 0 and some constant c. The second assumption holds for a number of concave functions, including ψ_i(x) = x^a, a ∈ (0, 1), ψ_i(x) = log(1+ x) and ψ_i(x) = min(x, a). -1





The main idea of this approach is to replace the concave functions ψ_i's by piece-wise linear approximations ψ_i^PL(x). We define an approximation of f(X) as f^PL(X) defined as f^PL(X) = ∑_i = 1^k ψ_i^PL(w_i(X)). We then optimize this piece-wise linear approximation function, and the approximation factor comes based on the tightness of this piece-wise linear approximation. We call this procedure the piece-wise linear approximation based algorithm (PLA).

We compute this approximation as follows. In the case of Problem 1, compute l_i = min{w_i(X) | X ∈𝒞} and u_i = max{w_i(X) | X ∈𝒞} for each i = 1, 2, ⋯, k. Both these computations are linear cost problems and are polynomial time for most constraints. In case these are NP hard for Problem 1, or in the case of Problems 2, 3 and 4, we set l_i = min{w_i(j), j ∈ V: w_i(j) > 0} and u_i = w_i(V). Then divide the range [l_i, u_i] into pieces with breakpoints b_i^1, b_i^2, ⋯, b_i^N_i such that b_i^1 = l_i, b_i^2 = l_i (1 + ϵ), b_i^3 = l_i (1 + ϵ)^2 and so on, for any ϵ > 0. It is easy to see that N_i = log_1+ϵ u_i/l_i ≈log(u_i/l_i)/ϵ. The precision ϵ defines the fineness of the points, and the quality of the approximation.

For all i = 1, 2, ⋯, k, define the piece-wise linear function ψ^PL_i, via the breakpoints b_i^1, b_i^2, ⋯, b_i^N_i. A visualization of this is shown in Figure <ref>, where the dotted lines are the piece-wise approximation, while the solid curve is the concave function ψ_i. We first show that the function f^PL approximates the function f within a factor of 1 + ϵ.


The piece-wise linear function f^PL defined with a precision ϵ^' satisfies,

    f^PL(X) ≤ f(X) ≤ (1 + ϵ^')^c f^PL(X) = (1 + ϵ) f^PL(X)

where c is a constant such that ψ_i(ky) ≤ k^c ψ_i(y) for k ≥ 1, y ≥ 0 for all i = 1, 2, ⋯, k.


By the construction of f^PL, and the concavity of the ψ_is, it is easy to see that f^PL(X) ≤ f(X), ∀ X ⊆ V. To show the upper bound, consider a region defined by breakpoints b_i^j and b_i^j+1. Due to concavity of ψ_i, there exists a tangent at some point in [b_i^j, b_i^j+1] whose slope equals that of the line connecting (b_i^j, ψ_i(b_i^j)) and (b_i^j+1, ψ_i(b_i^j+1)). This tangent line upper bounds the concave function ψ_i, and we can denote the corresponding upper bound as ψ_i^T. It then holds that ψ_i^PL(y) ≤ψ_i(y) ≤ψ_i^T(y). We now show that ψ_i^T(y) ≤ (1 + ϵ^')^c ψ_i^PL(y).  

We now focus on the region [b_i^j, b_i^j+1]. Let β be the constant difference between the two (parallel) lines, in terms of the y value. A visualization of this is shown in Figure <ref>. We would like to give a worst case bound on ψ^T_i(y)/ψ^PL_i(y), ∀ y ∈ [b_i^j, b_i^j+1]. Notice that ψ^T_i(y)/ψ^PL_i(y) = 1 + β/ψ^PL_i(y) ≤ 1 + β/ψ^PL_i(b^j_i) = 1 + β/ψ_i(b^j_i) ≤ψ_i(b^j+1_i)/ψ_i(b^j_i). The last inequality holds since ψ_i(b^j_i) + β≤ψ_i(b^j+1_i), and the second last one holds since b^j_i is a break point.

Moreover, b^j+1_i = b^j_i (1 + ϵ^') and hence ψ_i(b^j+1_i)/ψ_i(b^j_i) ≤ψ_i((1 + ϵ^')b^j_i)/ψ_i(b^j_i) ≤ (1 + ϵ^')^c = 1 + ϵ.


We now show how we can exactly solve Problems 1, 2 and 3 using the
cost function f^PL.  Let s^j_i denote the slopes of the
piece-wise linear functions – in other words, s^j_i =
[ψ_i(b^j+1_i) - ψ_i(b^j_i)]/[b^j+1_i - b^j_i]. Also, we
denote c^j_i as the corresponding intercepts. The functions
ψ_i^PL are characterized by the pairs {(s^1_i, c^1_i),
(s^2_i, c^2_i), ⋯, s^N_i_i, c^N_i_i)}, and ψ_i^PL(y)
= s^j_i.y + c^j_i, ∀ y ∈ [b^j_i, b^j+1_i]. We then 
consider the ∏_i = 1^k N_i different possibilities of the
cross-terms. Define J = [j_1, j_2, ⋯, j_k] as a vector such that
J ∈ [1, N_1] × [1, N_2] ×⋯× [1, N_k].

In the case of Problem 1, PLA solves a set of optimization problems,

    X̂_J =    {∑_i = 1^k s^j_i_i w_i(X) + c^j_i_i  |  X ∈𝒞}, 
       ∀ J ∈ [1, N_1] × [1, N_2] ×⋯× [1, N_k].

The final solution X̂ is the minimum among the ones above. For problem 2, we consider the set of problems,

    X̂_J =    {∑_i = 1^k s^j_i_i w_i(X) + c^j_i_i  |  g(X) ≥ c}, 
       ∀ J ∈ [1, N_1] × [1, N_2] ×⋯× [1, N_k],

and again set X̂ is the minimum among the X̂_J's above. Similarly, for Problem 3, we solve,

    X̂_J =    {g(X)  | ∑_i = 1^k s^j_i_i w_i(X) + c^j_i_i ≤ b}, 
       ∀ J ∈ [1, N_1] × [1, N_2] ×⋯× [1, N_k]

We set X̂ corresponding to the set with the largest value of g(X̂_J).
Finally, for Problem 4, we have:

    X̂_J =    {∑_i = 1^k s^j_i_i w_i(X) + c^j_i_i/g(X)}, 
       ∀ J ∈ [1, N_1] × [1, N_2] ×⋯× [1, N_k].


Our main result is that these simple procedures provide improved guarantees for all three problems.

PLA achieves an approximation factor of 1+ϵ for Problem 1 as long as a linear function can be exactly minimized under 𝒞. PLA also achieves a bi-criterion approximation factor satisfying σ = 1 + ϵ and ρ = 1 - 1/e for Problems 2 and 3. PLA also achieves a non bicriterion approximation factor of (1 + ϵ) log g(V) for Problem 2. PLA also achieves an approximation factor of e(1 + ϵ)/(e-1) for Problem 4. The worst case complexity of PLA is ∏_i = 1^k log(u_i/l_i) (1/ϵ)^k T = O((1/ϵ)^k T), where T is the complexity of Problems 1-4, with a linear cost function f.


We first show that PLA solves Problems 1-4 with the surrogate function f^PL. Note that with the piece-wise linear approximation, Problem 1 becomes min_X ∈𝒞∑_i = 1^k ψ^PL_i(w_i(X)) = min_X ∈𝒞∑_i = 1^k min{s^1_i w_i(X) + c^1_i, s^2_i w_i(X) + c^2_i, ⋯, s^N_i_i w_i(X) + c^N_i_i}. This holds since ψ^PL_i(y) = min{s^1_i y + c^1_i, s^2_i y + c^2_i, ⋯, s^N_i_i y + c^N_i_i}, due to the concavity of ψ_i's. We can then rewrite this as ∑_i = 1^k ψ^PL_i(w_i(X)) = min_J ∈𝒥∑_i = 1^k s^j_i_i w_i(X) + c^j_i_i, where 𝒥 = [1, N_1] × [1, N_2] ×⋯× [1, N_k]. Combining these facts, we can rewrite the problem as min_X ∈𝒞min_J ∈𝒥∑_i = 1^k s^j_i_i w_i(X) + c^j_i_i, which after interchanging the min's becomes min_J ∈𝒥min_X ∈𝒞∑_i = 1^k s^j_i_i w_i(X) + c^j_i_i. This is exactly Eq. (<ref>).

The algorithm for Problem 2 (Eq. (<ref>)) is basically the same as that of Problem 1, since it is a special case. Similarly we can write Problem 4 as min_X ∈𝒞min_J ∈𝒥∑_i = 1^k (s^j_i_i w_i(X) + c^j_i_i)/g(X) which is equivalent to min_J ∈𝒥min_X ∈𝒞∑_i = 1^k (s^j_i_i w_i(X) + c^j_i_i)/g(X), which becomes equation (<ref>). Equations (<ref>), (<ref>) and (<ref>) each become instances of Problems 1, 2 and 4 with f being modular and the approximation guarantees follow directly from <cit.>.

To deal with Problem 3, we use the fact that f^PL(X) = min_J ∈𝒥∑_i = 1^k (s^j_i_i w_i(X) + c^j_i_i), and hence we have the constraint, {min_J ∈𝒥∑_i = 1^k s^j_i_i w_i(X) + c^j_i_i ≤ b}. First we show that X_J is feasible for all ∈𝒥. This follows easily from the fact that if for any X, w_J(X) ≤ b, it holds that min_J ∈𝒥 w_J(X) ≤ b}. Next, let X^* be the optimal solution of Problem 3, and let J^* be such that ∑_i = 1^k s^j^*_i_i w_i(X^*) + c^j^*_i_i = f^PL(X^*). Note that our algorithm covers J^* and hence g(X̂) ≥ g(X̂_J^*) ≥ (1 - 1/e)g(X^*), where 1 - 1/e is the approximation factor of the submodular knapsack problem <cit.>. Note that the approximation factor of Problem 1 with f^PL is 1 assuming 𝒞 admits an exact solution with linear cost functions, while the factor for problem 2 is log g(V) for non-bicriterion algorithms, and a bicriterion factor of [1, 1 - 1/e] with a bi-criterion algorithm <cit.>.


Results similar to Theorem <ref> have been shown for a generalization of Problem 1, which asks for constrained optimization of low rank functions <cit.>. This problem in general is not a combinatorial optimization problem. However, when the functions are quasi-concave, the optimum lies on an extreme point, and hence, can be posed as a combinatorial optimization problem. Problem 1 asks for optimizing a specific subclass of concave (and hence quasi-concave) functions.  <cit.> focus on the class of low rank quasi-concave functions, while <cit.> consider the general class of low rank functions. While their algorithms apply to our class of functions as well, their approach while being more general, is also more complicated and involved.  <cit.> also consider a special case of Problem 1, with 𝒞 being the family of cuts (i.e., the cooperative cut problem). Interestingly, they suggest an algorithm that is identical to PLA when f is a (low rank) sum of truncations 
(i.e., ψ_i(x) = min(x, a)).
For general sums of low-rank concave functions, they resort to the algorithms of <cit.>. We provide a generic algorithm, which not only works for a much large class of constraints and functions, but also extends to the Problems 2, 3 and 4. Moreover, it is easy to see that our algorithms would also work for the more general problem of minimizing low rank sums of concave functions, over a solvable polytope.-1

Note that the complexity of PLA is polynomial in 1/ϵ, but exponential in k. Hence this makes sense only if k is a constant or is O(log n). If k is a constant (with respect to n), PLA is a fully polynomial time approximation scheme (FPTAS) <cit.>. If k = O(log n), then PLA is a polynomial-time approximation scheme (PTAS). This assumption is reasonable for many of the applications of Problems 1-4 (see details of this in the experiments section). Moreover, there are a number of ways one can speed up PLA. A very simple observation is that PLA is amenable to a distributive implementation via Map-Reduce. In particular, let N = O(1/ϵ^k) denote the total number of computations of PLA (i.e., this is the number of times one performs an instance of Problems 1-4 with a modular function). All these can be performed in parallel on m processing
systems. We output the best from each system to a central processor, which finds the optimal amongst these. The complexity of this distributive procedure is O(NT/m + m), (where T is the complexity of using a modular function in the place of f in Problems 1-4), which improves the overall complexity by a factor of m.check this, as originally
m was not defined. Also, T is currently not defined but should be.

In addition, we can also provide early stopping criterion and heuristics for speeding up PLA. One strategy of implementing PLA, is to start with j_i = 1, ∀ i = 1, 2, ⋯, k, and incrementally increase j_i in a coordinate ascent fashion. The following lemma gives a sufficient condition for stopping PLA.

Let J = [j_1, j_2, ⋯, j_k] be such that the corresponding solution X̂_J satisfies w_i(X_J) ∈ [b_i^j_i, b_i^j_i + 1], ∀ i. Then X̂_J is the (near) optimal solution for Problems 1, 2 and 3.

The values of w_i(X_J) also suggest the direction of the co-ordinate wise algorithm. For example, if w_i(X_J) < b_i^j_i, it suggests that the value of j_i be decreased. Similarly, if w_i(X_J) > b_i^j_i + 1, its a sign that j_i be decreased. In this manner, one can define a greedy like heuristic to implement PLA <cit.>, which picks for every coordinate, the slope which increases the objective value the most. Many of these heuristics have been considered in <cit.> in the case of cuts, and when the function class is low rank sums of truncations. These heuristics are all polynomial in k, but are not guaranteed to obtain the optimal solutions. Moreover, in certain cases (for example, the case of cuts), one can do parametric versions, thereby solving a set of related problems simultaneously <cit.>.



§ EXPERIMENTS



We next experimentally evaluate the performance of our methods. The utility of the constrained minimization algorithms for cooperative cuts have been investigated in <cit.>. In this paper, we consider the applications of cooperative image matching and sensor placement.



 §.§ Cooperative Image matching





The problem of matching key-points in images, also called image
correspondence, is an important problem in computer
vision <cit.>. The simplest model for this problem
constructs a matching with linear scores, i.e., a max bipartite
matching <cit.>, called a linear
  assignment. This model does not allow a representation of
interaction between the pixels. For example, we see many obviously
spurious matches in figure <ref>b. Many models try to
capture this, via, for example via quadratic
assignments <cit.>. Instead of just looking at the
best linear assignment, the quadratic models try to incorporate
pairwise constraints. This is also called graph matching.



We describe a new and different model here.  First, we cluster
key-points, separately in each of the two images, into k
clusters. Figure <ref>a shows a particular clustering of an
image into k = 3 groups. The clustering can be performed based on
the pixel color map, or simply the distance of the key-points. That is,
each image has k clusters. Let { V_i^(1)}_i=1^k and {
V_i^(2)}_i=1^k be the two sets of clusters.  We then compute
the linear assignment problem, letting ℳ⊆ℰ be the resulting maximum matching. We then partition the edge set ℰ = ℰ_1 ∪ℰ_2 ∪…ℰ_k ∪ℰ' where ℰ_i = ℳ∩ (V_ℓ^(1)× V_s^(2) ) for ℓ,s ∈{ 1, 2, …, k} corresponding
to the i'th largest intersection, and ℰ' = {ℰ\∪_i = 1^k ℰ_i} are the remaining edges
either that were not matched or that did not lie within a frequently
associated pair of image key-point clusters. We then
define a submodular function as follows:

    f(S) = ∑_i = 1^k ψ_i(w(S ∩ℰ_i)) + 
    w(S ∩ℰ'),

which provides an additional discount to the edges {ℰ_i}_i = 1^k corresponding to key-points that were frequently
associated in the initial pass. The problem of co-operative matching then becomes an instance of Problem 1 with the submodular function (over the edges) defined above, and a constraint that the edges form a matching. Figures <ref>b and
<ref>c shows how the submodular matchings improve over 
the simple bipartite matching, with k = 3. The minimum matching
approach obtains many spurious matches between clusters (shown in
red), while the cooperation described above reduces these spurious
matches. The cooperative matching improves the performance over the
modular method on these images by about 20%.

We also test the performance of our algorithms on the CMU House and Hotel dataset <cit.>. The house dataset has 111 images, while the hotel dataset has 101 images. We consider all possible pairs of images, with differences between the two images ranging from 0:10:90 in both cases. We consider three algorithms: PLA, SGA (both using Equation (<ref>)) and the simple modular bipartite matching as a baseline (Mod). Again, we set k = 3. 
The results are shown in Figure <ref>(d-e) where 
we observe that PLA and SGA beat Mod by about 3-5% on average. 
Moreover, we also see that PLA, in general, outperforms SGA, thus showing how superior theoretical guarantees translate into better empirical performance. In PLA, we chose ϵ such that each concave function ψ_i has four break points. We observed, moreover, that setting lower values of ϵ does not improve the objective value in this application. 
Should definitely give the ϵ used, and also report running times. Ideally, we would have multiple plots for various different ϵs and
that also show running time as well as accuracy.
We observe, moreover, that PLA also beats SGA in terms of objective value. We do not compare the ellipsoidal approximation algorithm (EA) <cit.>, mainly because it is too slow to run on real world problems. Moreover, this algorithm has been observed to perform comparably to the much simpler SGA <cit.>. While we considered the simple linear assignment as a baseline for the cooperative matching, it seems possible to embed this cooperation on more involved graph matching models as well. 




 §.§ Sensor Placement





We next consider an application of sensor placement. A number of natural models for this problem are forms of submodular maximization <cit.>. A natural model, that performs very well in practice, is to maximize the mutual information I(X_A; X_V \ A), where A refers to the set of sensors chosen. <cit.> investigate this in the setting of additive costs on the sensors. Often however, the costs are not additive in practice. In fact, very often, they are also submodular <cit.>, and a natural model is,

    f(X) = ∑_i=1^k ψ_i( c(X ∩ S_i) )

 where ψ_is are concave, c(j) is the cost of sensor j and S_1, ⋯, S_k are groups of similar sensors. This was posed as an open problem in <cit.>. We can naturally pose this as instances of Problem 3, where g(X) = I(X_A; X_V \ A) and f(X) is the cost function above. Note that we could equivalently also express this as an instance of Problem 2 with a constraint on g while minimizing f.
 
We consider real world data of placing sensors to predict the pH
values from the lake of Merced <cit.>. We also
assume that the function f is piece-wise linear, shown in
Figure <ref>a (far left).  



Figure <ref>b shows the locations (horizontal and vertical). We assume that there are three kinds of locations, shown in blue, green and red colors respectively, and the costs of placing sensors in the same kind of location is discounted. Correspondingly, we assume the cost function is an instance of function Equation (<ref>) with k = 3. For simplicity, we assume also that all three types of sensor locations have the same coverage model (though, in general, it would make sense for them to have different models for coverage, based on their type). Under this assumption, the optimal configuration would tend to be spatially diverse, yet cooperative (in the sense, that the same type of sensors would be chosen). 
TODO: This section needs to be fixed. We only have one type of sensor,
but multiple types of locations. The types of locations interact (i.e., underwater,
mountain, on a tower) in that it is cheaper to place sensors
in similar types of locations (due to shared fixed costs of purchasing
equipment, human capital costs and shared fixed finding and hiring costs of personnel, economies of scale of equipment, etc.). Also, the above right plot should
be changed to “cost” rather than budget.

We compare three algorithms: PLA, and SGA (both on Problem 3), and a simple cost agnostic greedy algorithm (AG), which ignores the cost function f, and greedily adds sensors. Figure <ref>c shows the sensors chosen by PLA (the cost sensitive one), and Figure <ref>d shows the choices of AG (the cost agnostic one). While both have the same cost budget, the cost agnostic one does not utilize the discounts of placing sensors in similar locations, and correspondingly, places fewer sensors. The cost sensitive algorithms (PLA and SGA) on the other hand, simultaneously achieve coverage, while making use of the discounts. Figure <ref>e plots the objective functions attained by the three algorithms. We see that both PLA and SGA, outperform the agnostic greedy algorithm. Moreover, PLA also performs better than SGA. Note that the function f, used in this case is piece-wise linear, and correspondingly PLA is exact in this case.



§ CONCLUSIONS



In this paper, we investigated a new class of algorithms for various
forms of constrained submodular programs, with a special subclass of
submodular cost functions. We focus on problems that for the general
class of submodular functions are hard, and yet occur naturally in
many applications. We showed that when we restrict the class of
functions to low rank sums of concave over modular functions, we can
obtain significantly improved worst case theoretical results. We also
complemented our results with experimental results in sensor placement
and image correspondence. An immediate open question is whether there
are similar algorithms for other rich and useful subclasses of
submodular functions. In particular, it would be interesting if one
can remove the low rank assumption, and provide tighter approximation
algorithms for general sums of concave over modular functions, which
would be very powerful.

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

abbrv

