arXiv:1612.00960v1 [] 3 Dec 2016

Non-monotone DR-Submodular Function Maximization
(Full version)
Tasuku Soma

Yuichi Yoshida

The University of Tokyo
tasuku soma@mist.i.u-tokyo.ac.jp

National Institute of Informatics, and
Preferred Infrastructure, Inc.
yyoshida@nii.ac.jp

Abstract
We consider non-monotone DR-submodular function maximization, where DR-submodularity (diminishing return submodularity) is an extension of submodularity for functions
over the integer lattice based on the concept of the diminishing return property. Maximizing non-monotone DRsubmodular functions has many applications in machine
learning that cannot be captured by submodular set functions.
1
In this paper, we present a 2+
-approximation algorithm with
a running time of roughly O( n log2 B), where n is the size
of the ground set, B is the maximum value of a coordinate,
and  > 0 is a parameter. The approximation ratio is almost
tight and the dependency of running time on B is exponentially smaller than the naive greedy algorithm. Experiments
on synthetic and real-world datasets demonstrate that our algorithm outputs almost the best solution compared to other
baseline algorithms, whereas its running time is several orders of magnitude faster.

1

Introduction

Submodular functions have played a key role in various
tasks in machine learning, statistics, social science, and economics. A set function f : 2E → R with a ground set E is
submodular if
f (X ∪ {e}) − f (X) ≥ f (Y ∪ {e}) − f (Y )
for arbitrary sets X, Y ⊆ E with X ⊆ Y , and an element e ∈ E \ Y . The importance and usefulness of submodularity in these areas are due to the fact that submodular functions naturally capture the diminishing return property. Various important functions in these areas such as the
entropy function, coverage function, and utility functions
satisfy this property. See, e.g., (Krause and Golovin 2014;
Fujishige 2005).
Recently, maximizing (non-monotone) submodular functions has attracted particular interest in the machine learning community. In contrast to minimizing submodular functions, which can be done in polynomial time, maximizing
submodular functions is NP-hard in general. However, we
can achieve a constant factor approximation for various settings. Notably, (Buchbinder et al. 2012) presented a very elegant double greedy algorithm for (unconstrained) submodular function maximization, which was the first algorithm
achieving 12 -approximation, and this approximation ratio is

tight (Feige, Mirrokni, and Vondrak 2011). Applications of
non-monotone submodular function maximization include
efficient sensor placement (Krause, Singh, and Guestrin
2008), privacy in online services (Krause and Horvitz 2008),
and maximum entropy sampling (Ko, Lee, and Queyranne
1995).
The models and applications mentioned so far are built
upon submodular set functions. Although set functions are
fairly powerful for describing problems such as variable selection, we sometimes face difficult situations that cannot be
cast with set functions. For example, in the budget allocation
problem (Alon, Gamzu, and Tennenholtz 2012), we would
like to decide how much budget should be set aside for each
ad source, rather than whether we use the ad source or not.
A similar issue arises when we consider models allowing
multiple choices of an element in the ground set.
To deal with such situations, several generalizations of
submodularity have been proposed. (Soma et al. 2014) devised a general framework for maximizing monotone submodular functions on the integer lattice, and showed that
the budget allocation problem and its variant fall into this
framework. In their framework, functions are defined over
the integer lattice ZE
+ and therefore effectively represent discrete allocations of budget. Regarding the original motivation for the diminishing return property, one can naturally
consider its generalization to the integer lattice: a function
f : ZE
+ → R satisfying
f (x + χe ) − f (x) ≥ f (y + χe ) − f (y)
for x ≤ y and e ∈ E, where χe ∈ RE is the vector with
χe (e) = 1 and χe (a) = 0 for every a 6= e. Such functions
are called diminishing return submodular (DR-submodular)
functions (Soma and Yoshida 2015) or coordinate-wise concave submodular functions (Milgrom and Strulovici 2009).
DR-submodular functions have found various applications
in generalized sensor placement (Soma and Yoshida 2015)
and (a natural special case of) the budget allocation problem (Soma et al. 2014).
As a related notion, a function is said to be lattice submodular if
f (x) + f (y) ≥ f (x ∨ y) + f (x ∧ y)
for arbitrary x and y, where ∨ and ∧ are coordinate-wise
max and min, respectively. Note that DR-submodularity is

stronger than lattice submodularity in general (see, e.g.,
(Soma et al. 2014)). Nevertheless, we consider the DRsubmodularity to be a “natural definition” of submodularity, at least for the applications mentioned so far, because
the diminishing return property is crucial in these real-world
scenarios.

Our contributions
We design a novel polynomial-time approximation algorithm for maximizing (non-monotone) DR-submodular
functions. More precisely, we consider the optimization
problem
maximize f (x)
subject to 0 ≤ x ≤ B,

(1)

where f : ZE
+ → R+ is a non-negative DR-submodular
function, 0 is the zero vector, and B ∈ ZE
+ is a vector representing the maximum values for each coordinate. When
B is the all-ones vector, this is equivalent to the original
(unconstrained) submodular function maximization. We assume that f is given as an evaluation oracle; when we specify x ∈ ZE
+ , the oracle returns the value of f (x).
1
Our algorithm achieves 2+
-approximation for any con∆
stant  > 0 in O( |E|
 ·log( δ ) log B ·(θ+log B)) time, where
δ and ∆ are the minimum positive marginal gain and maximum positive values, respectively, of f , B = kBk∞ :=
maxe∈E B(e), and θ is the running time of evaluating (the
oracle for) f . To our knowledge, this is the first polynomialtime algorithm achieving (roughly) 12 -approximation.
We also conduct numerical experiments on the revenue
maximization problem using real-world networks. The experimental results show that the solution quality of our algorithm is comparable to other algorithms. Furthermore, our
algorithm runs several orders of magnitude faster than other
algorithms when B is large.
DR-submodularity is necessary for obtaining polynomialtime algorithms with a meaningful approximation guarantee; if f is only lattice submodular, then we cannot obtain
constant-approximation in polynomial time. To see this, it
suffices to observe that an arbitrary univariate function is
lattice submodular, and therefore finding an (approximate)
maximum value must invoke O(B) queries. We note that
representing an integer B requires dlog2 Be bits. Hence, the
running time of O(B) is pseudopolynomial rather than polynomial.

Fast simulation of the double greedy algorithm
Naturally, one can reduce the problem (1) to maximization
of a submodular set function by simply duplicating each
element e in the ground set into B(e) distinct copies and
defining a set function over the set of all the copies. One
can then run the double greedy algorithm (Buchbinder et al.
2012) to obtain 12 -approximation. This reduction is simple
butP
has one large drawback; the size of the new ground set
is e∈E B(e), a pseudopolynomial in B. Therefore, this
naive double greedy algorithm does not scale to a situation
where B is large.

For scalability, we need an additional trick that reduces
the pseudo-polynomial running time to a polynomial one. In
monotone submodular function maximization on the integer
lattice, (Soma and Yoshida 2015; Soma and Yoshida 2016)
provide such a speedup trick, which effectively combines the
decreasing threshold technique (Badanidiyuru and Vondrák
2014) with binary search. However, a similar technique does
not apply to our setting, because the double greedy algorithm works differently from (single) greedy algorithms for
monotone submodular function maximization. The double
greedy algorithm examines each element in a fixed order and
marginal gains are used to decide whether to include the element or not. In contrast, the greedy algorithm chooses each
element in decreasing order of marginal gains, and this property is crucial for the decreasing threshold technique.
We resolve this issue by splitting the set of all marginal
gains into polynomially many small intervals. For each interval, we approximately execute multiple steps of the double greedy algorithm at once, as long as the marginal gains
remain in the interval. Because the marginal gains do not
change (much) within the interval, this simulation can be
done with polynomially many queries and polynomial-time
overhead. To our knowledge, this speedup technique is not
known in the literature and is therefore of more general interest.
Very recently, (Ene and Nguyen 2016) pointed out that
a DR-submodular function f : {0, 1, . . . , B}E → R+
can be expressed as a submodular set function g over a
polynomial-sized ground set, which turns out to be E ×
{0, 1, . . . , k − 1}, where k = dlog2 (B + 1)e. Their idea
is representing x(e) in binary form for each e ∈ E, and
bits in the binary representations form the new ground set.
We may want to apply the double greedy algorithm to g
in order to get a polynomial-time approximation algorithm.
However, this strategy has the following two drawbacks:
(i) The value of g(E × {0, 1, . . . , k − 1}) is defined as
f (x), where x(e) = 2k − 1 for every e ∈ E. This means
that we have to extend the domain of f . (ii) More crucially, the double greedy algorithm on g may return a large
set such as E × {0, 1, . . . , k − 1} whose corresponding
vector x ∈ ZE
+ may violate the constraint x ≤ B. Although we can resolve these issues by introducing a knapsack constraint, it is not a practical solution because existing algorithms for knapsack constraints (Lee et al. 2009;
Chekuri, Vondrák, and Zenklusen 2014) are slow and have
worse approximation ratios than 1/2.
Notations For an integer n ∈ N, [n] denotes the set
{1, . . . , n}. For vectors x, y ∈ ZE , we define f (x | y) :=
f (x + y) − f (y). The `1 -norm and `∞P
-norm of a vector x ∈ ZE are defined as kxk1 :=
e∈E |x(e)| and
kxk∞ := maxe∈E |x(e)|, respectively.

2

Related work

As mentioned above, there have been many efforts to maximize submodular functions on the integer lattice. Perhaps
the work most related to our interest is (Gottschalk and
Peis 2015), in which the authors considered maximizing

lattice submodular functions over the bounded integer lattice and designed 13 -approximation pseudopolynomial-time
algorithm. Their algorithm was also based on the double
greedy algorithm, but does not include a speeding up technique, as proposed in this paper.
In addition there are several studies on the constrained
maximization of submodular functions (Feige, Mirrokni,
and Vondrak 2011; Buchbinder et al. 2014; Buchbinder and
Feldman 2016), although we focus on the unconstrained
case. Many algorithms for maximizing submodular functions are randomized, but a very recent work (Buchbinder
and Feldman 2016) devised a derandomized version of the
double greedy algorithm. (Gotovos, Karbasi, and Krause
2015) considered maximizing non-monotone submodular
functions in the adaptive setting, a concept introduced
in (Golovin and Krause 2011).
A continuous analogue of DR-submodular functions is
considered in (Bian et al. 2016).

3

Algorithms

In this section, we present a polynomial-time approximation
algorithm for maximizing (non-monotone) DR-submodular
functions. We first explain a simple adaption of the double
greedy algorithm for (set) submodular functions to our setting, which runs in pseudopolynomial time. Then, we show
how to achieve a polynomial number of oracle calls. Finally,
we provide an algorithm with a polynomial running time
(details are placed in Appendix A).

Pseudopolynomial-time algorithm
Algorithm 1 is an immediate extension of the double
greedy algorithm for maximizing submodular (set) functions (Buchbinder et al. 2012) to our setting. We start with
x = 0 and y = B, and then for each e ∈ E, we tighten the
gap between x(e) and y(e) until they become exactly the
same. Let α = f (χe | x) and β = f (−χe | y). We note
that
α + β = f (x + χe ) − f (x) − (f (y) − f (y − χe )) ≥ 0
holds from the DR-submodularity of f . Hence, if β < 0,
then α > 0 must hold, and we increase x(e) by one. Similarly, if α < 0, then β > 0 must hold, and we decrease y(e)
by one. When both of them are non-negative, we increase
α
x(e) by one with probability α+β
, or decrease y(e) by one
with the complement probability

β
α+β .

Theorem 1. Algorithm 1 is a 12 -approximation algorithm
for (1) with time complexity O(kBk1 · θ + kBk1 ), where θ
is the running time of evaluating f .
We omit the proof as it is a simple modification of the
analysis of the original algorithm.

Algorithm with polynomially many oracle calls
In this section, we present an algorithm with a polynomial
number of oracle calls.
Our strategy is to simulate Algorithm 1 without evaluating
the input function f many times. A key observation is that,
at Line 4 of Algorithm 1, we do not need to know the exact

Algorithm 1 Pseudopolynomial-time algorithm
E
Input: f : ZE
+ → R+ , B ∈ Z+ .
Assumption: f is DR-submodular.
1: x ← 0, y ← B.
2: for e ∈ E do
3:
while x(e) < y(e) do
4:
α ← f (χe | x) and β ← f (−χe | y).
5:
if β < 0 then
6:
x(e) ← x(e) + 1.
7:
else if α < 0 then
8:
y(e) ← y(e) − 1.
9:
else
α
10:
x(e) ← x(e) + 1 with probability α+β
and
y(e) ← y(e) − 1 with the complement probaβ
α
bility α+β
. If α = β = 0, we assume α+β
= 1.
11:
end if
12:
end while
13: end for
14: return x.

value of f (χe | x) and f (−χe | y); good approximations
to them are sufficient to achieve an approximation guarantee
close to 21 . To exploit this observation, we first design an
algorithm that outputs (sketches of) approximations to the
functions g(b) := f (χe | x + bχe ) and h(b) := f (−χe |
y − bχe ). Note that g and h are non-increasing functions in
b because of the DR-submodularity of f .
To illustrate this idea, let us consider a non-increasing
function φ : {0, 1, . . . , B − 1} → R and suppose that φ
is non-negative (φ is either g or h later on). Let δ and ∆
be the minimum and the maximum positive values of φ, respectively. Then, for each δ ≤ τ ≤ ∆ of the form δ(1 + )k ,
we find the minimum bτ such that φ(bτ ) < τ (we regard
φ(B) = −∞). From the non-increasing property of φ, we
then have φ(b) ≥ τ for any b < bτ . Using the set of pairs
{(τ, bτ )}τ , we can obtain a good approximation to φ. The
details are provided in Algorithm 2.
Lemma 2. For any φ : {0, 1, . . . , B − 1} → R and  > 0,
Algorithm 2 outputs a set of pairs {(bτ , τ )}τ from which, for
any b ∈ {0, 1, . . . , B − 1}, we can reconstruct a value v in
O(log B) time such that v ≤ φ(b) < (1 + )v if φ(b) > 0
and v = 0 otherwise. The time complexity of Algorithm 2
is O( 1 log( ∆
δ ) log B · θ) if φ has a positive value, where δ
and ∆ are the minimum and maximum positive values of φ,
respectively, and θ is the running time of evaluating φ, and
is O(log B · θ) otherwise.
Proof. Let S = {(bτ , τ )}τ be the set of pairs output by Algorithm 2. Our reconstruction algorithm is as follows: Given
b ∈ {0, 1, . . . , B − 1}, let (bτ ∗ , τ ∗ ) be the pair with the minimum bτ ∗ , where b < bτ ∗ . Note that such a bτ ∗ always exists
because a pair of the form (B, ·) is always added to S. We
then output τ ∗ . The time complexity of this reconstruction
algorithm is clearly O(log B).
We now show the correctness of the reconstruction algorithm. If φ(b) > 0, then, in particular, we have φ(b) ≥ δ.
Then, τ ∗ is the maximum value of the form δ(1+)k at most

Algorithm 2 Sketching subroutine for Algorithm 3

Algorithm 3 Algorithm with polynomially many queries

Input: φ : {0, 1, . . . , B − 1} → R,  > 0.
Assumption: φ is non-increasing.
1: S ← ∅. We regard φ(B) = −∞.
2: Find the minimum b0 ∈ {0, 1, . . . , B} with φ(b0 ) ≤ 0
by binary search.
3: if b0 ≥ 1 then
# φ has a positive value.
4:
∆ ← φ(0) and δ ← φ(b0 − 1).
5:
for (τ ← δ; τ ≤ ∆; τ ← (1 + )τ ) do
6:
Find the minimum bτ ∈ {0, 1, . . . , B} with
φ(bτ ) < τ by binary search.
7:
S ← S ∪ {(bτ , τ )}
8:
end for
9:
if bδ 6= B then
10:
S ← S ∪ {(B, 0)}.
11:
end if
12: else
# φ is non-positive.
13:
S ← S ∪ {(B, 0)}.
14: end if
15: return S.

E
Input: f : ZE
+ → R+ , B ∈ Z+ ,  > 0.
Assumption: f is DR-submodular.
1: x ← 0, y ← B.
2: for e ∈ E do
3:
Define g, h : {0, 1, . . . , B−1} → R as g(b) = f (χe |
x + bχe ) and h(b) = f (−χe | y − bχe ).
4:
Let g̃ and h̃ be approximations to g and h, respectively, given by Algorithm 2.
5:
while x(e) < y(e) do
6:
α ← g̃(x(e)) and β ← h̃(B(e) − y(e)).
α
7:
x(e) ← x(e)+1 with probability α+β
and y(e) ←

y(e) − 1 with the complement probability
α
α = β = 0, we assume α+β
= 1.
8:
end while
9: end for
10: return x.

β
α+β .

If

the kth iteration. We then have
φ(b). Hence, we have τ ∗ ≤ φ(b) < (1 + )τ ∗ . If φ(b) ≤ 0,
(bτ ∗ , τ ∗ ) = (B, 0) and we output zero.
Finally, we analyze the time complexity of Algorithm 2.
Each binary search requires O(log B) time. The number of
1
∆
binary searches performed is O(log1+ ∆
δ ) = O(  log δ )
when φ has a positive value and 1 when φ is non-positive.
Hence, we have the desired time complexity.
We can regard Algorithm 2 as providing a value oracle for
a function φ̃ : {0, 1, . . . , B − 1} → R+ that is an approximation to the input function φ : {0, 1, . . . , B − 1} → R.
We now describe our algorithm for maximizing DRsubmodular functions. The basic idea is similar to Algorithm 1, but when we need f (χe | x) and f (−χe | y), we
use approximations to them instead. Let α and β be approximations to f (χe | x) and f (−χe | y), respectively, obtained by Algorithm 2. Then, we increase x(e) by one with
α
probability α+β
and decrease y(e) by one with the compleβ
ment probability α+β
. The details are given in Algorithm 3.
We now analyze Algorithm 3. An iteration refers to an
iteration in the while loop from Line 5. We have kBk1 iterations in total. For k ∈ {1, . . . , kBk1 }, let xk and yk be
x and y, respectively, right after the kth iteration. Note that
xkBk1 = ykBk1 is the output of the algorithm. We define
x0 = 0 and y0 = B for convenience.
Let o be an optimal solution. For k ∈ {0, 1, . . . , kBk1 },
we then define ok = (o ∨ xk ) ∧ yk . Note that o0 = o holds
and okBk1 equals the output of the algorithm. We have the
following key lemma.
Lemma 3. For every k ∈ [kBk1 ], we have

E[f (ok−1 ) − f (ok )]
1+
≤
E[f (xk ) − f (xk−1 ) + f (yk ) − f (yk−1 )] (2)
2
Proof. Fix k ∈ [kBk1 ] and let e be the element of interest
in the kth iteration. Let α and β be the values in Line 6 in

E[f (xk ) − f (xk−1 ) + f (yk ) − f (yk−1 )]
α
β
=
f (χe | xk−1 ) +
f (−χe | yk−1 )
α+β
α+β
α
β
α2 + β 2
≥
α+
β=
,
(3)
α+β
α+β
α+β
where we use the guarantee in Lemma 2 in the inequality.
We next establish an upper bound of E[f (ok−1 )−f (ok )].
As ok = (o ∨ xk ) ∧ yk , conditioned on a fixed ok−1 , we
obtain
E[f (ok−1 ) − f (ok )]

α 
=
f (ok−1 ) − f (ok−1 ∨ xk (e)χe )
α+β

β 
+
f (ok−1 ) − f (ok−1 ∧ yk (e)χe ) .
α+β
Claim 4. (4) ≤

(4)

(1+)αβ
α+β .

Proof. We show this claim by considering the following
three cases.
If xk (e) ≤ ok−1 (e) ≤ yk (e), then (4) is zero.
If ok−1 (e) < xk (e), then ok (e) = ok−1 (e) + 1, and the
first term of (4) is
f (ok−1 ) − f (ok−1 ∨ xk (e)χe )
= f (ok−1 ) − f (ok−1 + χe )
≤ f (yk−1 − χe ) − f (yk−1 )
= f (−χe | yk−1 )
≤ (1 + )β.
Here, the first inequality uses the DR-submodularity of f
and the fact that ok−1 ≤ yk−1 − χe , and the second inequality uses the guarantee in Lemma 2. The second term
of (4) is zero, and hence we have (4) ≤ (1+)αβ
α+β .
If yk (e) < ok−1 (e), then by a similar argument, we have
(4) ≤ (1+)αβ
α+β .

Polynomial-time algorithm

We now return to proving Lemma 3. By Claim 4,
(4) ≤

In many applications, the running time needed to evaluate
the input function is a bottleneck, and hence Algorithm 3 is
already satisfactory. However, it is theoretically interesting
to reduce the total running time to a polynomial, and we
show the following. The proof is deferred to Appendix A.

(1 + )αβ
1 +   α2 + β 2  1 + 
≤
≤
· (3),
α+β
2
α+β
2

which indicates the desired result.
1
2+ -approximation algorithm
∆
O( |E|
 · log( δ ) log kBk∞ · θ +

Theorem 5. Algorithm 3 is a

for (1) with time complexity
kBk1 log kBk∞ ), where δ and ∆ are the minimum positive
marginal gain and the maximum positive value, respectively,
of f and θ is the running time of evaluating f .

4

Proof. Summing up (2) for k ∈ [kBk1 ], we get

Experiments

In this section, we show our experimental results and the
superiority of our algorithm with respect to other baseline
algorithms.

kBk1

X

1
Theorem 7. There exists a 2+2
-approximation algo|E|
e
rithm with time complexity O(  log( ∆
δ ) log kBk∞ · (θ +
log kBk∞ )), where δ and ∆ are the minimum positive
marginal gain and the maximum positive value, respectively
e )
of f and θ is the running time of evaluating f . Here O(T
means O(T logc T ) for some c ∈ N.

E[f (ok−1 ) − f (ok )]

k=1
kBk1
1+ X
≤
E[f (xk ) − f (xk−1 ) + f (yk ) − f (yk−1 )].
2
k=1

The above sum is telescopic, and hence we obtain
E[f (o0 ) − f (okBk1 )]
1+
≤
E[f (xkBk1 ) − f (x0 ) + f (ykBk1 ) − f (y0 )]
2
1+
≤
E[f (xkBk1 ) + f (ykBk1 )]
2
= (1 + ) E[f (xkBk1 )].
The second inequality uses the fact that f is non-negative,
and the last equality uses ykBk1 = xkBk1 . Because
E[f (o0 ) − f (okBk1 )] = f (o) − E[f (xkBk1 )], we obtain
1
that E[f (xkBk1 )] ≥ 2+
f (o).
We now analyze the time complexity. We only query the
input function f inside of Algorithm 2, and the number of
∆
oracle calls is O( |E|
 log( δ ) log B) by Lemma 2. Note that
we invoke Algorithm 2 with g and h, and the minimum positive values of g and h are at least the minimum positive
marginal gain δ of f . The number of iterations is kBk1 , and
we need O(log B) time to access g̃ and h̃. Hence, the total
time complexity is as stated.
Remark 6. We note that even if f is not a non-negative function, the proof of Theorem 5 works as long as f (x0 ) ≥ 0 and
f (y0 ) ≥ 0, that is, f (0) ≥ 0 and f (B) ≥ 0. Hence, given
E
a DR-submodular function f : ZE
+ → R and B ∈ Z+ , we
1
can obtain a 2+ -approximation algorithm for the following
problem:
maximize f (x) − min{f (0), f (B)}
subject to 0 ≤ x ≤ B,

(5)

This observation is useful, as the objective function often
takes negative values in real-world applications.

Experimental setting
We conducted experiments on a Linux server with an Intel
Xeon E5-2690 (2.90 GHz) processor and 256 GB of main
memory. All the algorithms were implemented in C# and
were run using Mono 4.2.3.
We compared the following four algorithms:
• Single Greedy (SG): We start with x = 0. For each element e ∈ E, as long as the marginal gain of adding
χe to the current solution x is positive, we add it to x.
The reason that we do not choose the element with the
maximum marginal gain is to reduce the number of oracle calls, and our preliminary experiments showed that
such a tweak does not improve the solution quality.
• Double Greedy (DG, Algorithm 1).
• Lattice Double Greedy (Lattice-DG): The 1/3approximation algorithm for maximizing non-monotone
lattice submodular functions (Gottschalk and Peis 2015).
• Double Greedy with a polynomial number of oracle calls
with error parameter  > 0 (Fast-DG , Algorithm 3).
We measure the efficiency of an algorithm by the number
of oracle calls instead of the total running time. Indeed, the
running time for evaluating the input function is the dominant factor of the total running, because objective functions
in typical machine learning tasks contain sums over all data
points, which is time consuming. Therefore, we do not consider the polynomial-time algorithm (Theorem 7) here.

Revenue maximization
In this application, we consider revenue maximization on
an (undirected) social network G = (V, W ), where W =
(wij )i,j∈V represents the weights of edges. The goal is to
offer for free or advertise a product to users so that the revenue increases through their word-of-mouth effect on others.
If we invest x units of cost on a user i ∈ V , the user becomes
an advocate of the product (independently from other users)
with probability 1 − (1 − p)x , where p ∈ (0, 1) is a parameter. This means that, for investing a unit cost to i, we
have an extra chance that the user i becomes an advocate

109

1010
Number of oracle calls

Number of oracle calls

1010

1011
SG
DG
Lattice-DG
Fast-DG0.5
Fast-DG0.05
Fast-DG0.005

108
107
106

109

1011
SG
DG
Lattice-DG
Fast-DG0.5
Fast-DG0.05
Fast-DG0.005

108
107
106

105 2
10

103

104
B

105

106

SG
DG
Lattice-DG
Fast-DG0.5
Fast-DG0.05
Fast-DG0.005

1010
Number of oracle calls

1011

109
108
107
106

105 2
10

103

(a) Adolescent health

104
B

105

105 2
10

106

(b) Advogato

103

104
B

105

106

(c) Twitter lists

Figure 1: Number of oracle calls
Table 1: Objective values (Our methods are highlighted in
gray.)

SG
DG
Lattice-DG
Fast-DG0.5
Fast-DG0.05
Fast-DG0.005

B = 102
280.55
280.55
215.39
280.55
280.55
280.55

SG
DG
Lattice-DG
Fast-DG0.5
Fast-DG0.05
Fast-DG0.005

B = 102
993.15
993.15
753.93
993.15
993.15
993.15

SG
DG
Lattice-DG
Fast-DG0.5
Fast-DG0.05
Fast-DG0.005

B = 102
882.43
882.43
675.67
882.43
882.43
882.43

Adolescent health
103
104
105
2452.16
7093.73
7331.42
2452.16
7124.90
7332.96
1699.66
6808.97
6709.11
2452.16
7101.14
7331.36
2452.16
7100.86
7331.36
2452.16
7100.83
7331.36
Advogato
103
104
105
8680.87 25516.05 27325.78
8680.87 25330.91 27329.39
6123.39 24289.09 24878.94
8680.87 25520.83 27325.75
8680.87 25520.52 27325.75
8680.87 25520.47 27325.75
Twitter lists
103
104
105
7713.07 22452.61 25743.26
7713.07 22455.97 25751.42
5263.87 20918.89 20847.48
7713.07 22664.65 25743.06
7713.07 22658.58 25743.06
7713.07 22658.07 25743.06

106
7331.50
7331.50
5734.30
7331.48
7331.48
7331.48
106
27326.01
27326.01
21674.35
27325.98
27325.98
27325.98
106
25744.02
25744.02
15001.19
25743.88
25743.88
25743.88

with probability p. Let S ⊆ V be a set of users who advocate the product. Note that S is a random set. Following
a simplified version of the model introduced by (Hartline,
Mirrokni,
the revenue is defined as
P
P and Sundararajan 2008),
E
w
.
Let
f
:
Z
→
R be the expected revij
+
i∈S
j∈V \S
enue obtained in this model, that is,
hX X
i
f (x) = E
wij
S

=

i∈S j∈V \S

X X

wij (1 − (1 − p)x(i) )(1 − p)x(j) .

i∈S j∈V \S

It is not hard to show that f is non-monotone DRsubmodular function (see Appendix B for the proof).
In our experiment, we used three networks, Adolescent
health (2,539 vertices and 12,969 edges), Advogato (6,541
vertices and 61,127 edges), and Twitter lists (23,370 vertices
and 33,101 edges), all taken from (Kunegis 2013). We regard

all the networks as undirected. We set p = 0.0001, and set
wij = 1 when an edge exists between i and j and wij = 0
otherwise. We imposed the constraint 0 ≤ x(e) ≤ B for
every e ∈ E, where B is chosen from {102 , . . . , 106 }.
Table 1 shows the objective values obtained by each
method. As can be seen, except for Lattice-DG, which is
clearly the worst, the choice of a method does not much affect the obtained objective value for all the networks. Notably, even when  is as large as 0.5, the objective values
obtained by Fast-DG are almost the same as SG and DG.
Figure 1 illustrates the number of oracle calls of each
method. The number of oracle calls of DG and Lattice-DG
is linear in B, whereas that of Fast-DG slowly grows. Although the number of oracle calls of SG also slowly grows,
it is always orders of magnitude larger than that of Fast-DG
with  = 0.5 or  = 0.05.
In summary, Fast-DG0.5 achieves almost the best objective value, whereas the number of oracle calls is two or three
orders of magnitude smaller than those of the other methods
when B is large.

5

Conclusions

1
In this paper, we proposed a polynomial-time 2+
approximation algorithm for non-monotone DR-submodular
function maximization. Our experimental results on the
revenu maximization problem showed the superiority of our
method against other baseline algorithms.
Maximizing a submodular set function under constraints
is well studied (Lee et al. 2009; Gupta et al. 2010;
Chekuri, Vondrák, and Zenklusen 2014; Mirzasoleiman et
al. 2016). An intriguing open question is whether we can
obtain polynomial-time algorithms for maximizing DRsubmodular functions under constraints such as cardinality constraints, polymatroid constraints, and knapsack constraints.

Acknowledgments T. S. is supported by JSPS Grantin-Aid for Research Activity Start-up. Y. Y. is supported by JSPS Grant-in-Aid for Young Scientists (B)
(No. 26730009), MEXT Grant-in-Aid for Scientific Research on Innovative Areas (No. 24106003), and JST, ERATO, Kawarabayashi Large Graph Project.

References
[Alon, Gamzu, and Tennenholtz 2012] Alon, N.; Gamzu, I.;
and Tennenholtz, M. 2012. Optimizing budget allocation
among channels and influencers. In WWW, 381–388.
[Badanidiyuru and Vondrák 2014] Badanidiyuru, A., and
Vondrák, J.
2014.
Fast algorithms for maximizing
submodular functions. In SODA, 1497–1514.
[Bian et al. 2016] Bian, Y.; Mirzasoleiman, B.; Buhmann,
J. M.; and Krause, A. 2016. Guaranteed non-convex optimization: Submodular maximization over continuous domains. CoRR abs/1606.05615.
[Buchbinder and Feldman 2016] Buchbinder, N., and Feldman, M. 2016. Deterministic algorithms for submodular
maximization problems. In SODA, 392–403.
[Buchbinder et al. 2012] Buchbinder, N.; Feldman, M.;
Naor, J. S.; and Schwartz, R. 2012. A tight linear
time (1/2)-approximation for unconstrained submodular
maximization. In FOCS, 649–658.
[Buchbinder et al. 2014] Buchbinder, N.; Feldman, M.;
Naor, J. S.; and Schwartz, R. 2014. Submodular maximization with cardinality constraints. In SODA, 1433–1452.
[Chekuri, Vondrák, and Zenklusen 2014] Chekuri,
C.;
Vondrák, J.; and Zenklusen, R. 2014. Submodular function
maximization via the multilinear relaxation and contention resolution schemes. SIAM Journal on Computing
43(6):1831–1879.
[Devroye 1986] Devroye, L. 1986. Non-Uniform Random
Variate Generation. Springer.
[Ene and Nguyen 2016] Ene, A., and Nguyen, H. L. 2016. A
reduction for optimizing lattice submodular functions with
diminishing returns. CoRR abs/1606.08362.
[Feige, Mirrokni, and Vondrak 2011] Feige, U.; Mirrokni,
V. S.; and Vondrak, J. 2011. Maximizing non-monotone submodular functions. SIAM J. on Comput. 40(4):1133–1153.
[Fujishige 2005] Fujishige, S. 2005. Submodular Functions
and Optimization. Elsevier, 2nd edition.
[Golovin and Krause 2011] Golovin, D., and Krause, A.
2011. Adaptive submodularity: Theory and applications in
active learning and stochastic optimization. Journal of Artificial Intelligence Research 427–486.
[Gotovos, Karbasi, and Krause 2015] Gotovos, A.; Karbasi,
A.; and Krause, A. 2015. Non-monotone adaptive submodular maximization. In AAAI, 1996–2003.
[Gottschalk and Peis 2015] Gottschalk, C., and Peis, B.
2015. Submodular function maximization on the bounded
integer lattice. In WAOA, 133–144.
[Gupta et al. 2010] Gupta, A.; Roth, A.; Schoenebeck, G.;
and Talwar, K. 2010. Constrained non-monotone submodular maximization: offline and secretary algorithms. In WINE,
246–257.
[Hartline, Mirrokni, and Sundararajan 2008] Hartline,
J.;
Mirrokni, V.; and Sundararajan, M. 2008. Optimal
marketing strategies over social networks.
In WWW,
189–198.

[Ko, Lee, and Queyranne 1995] Ko, C. W.; Lee, J.; and
Queyranne, M. 1995. An exact algorithm for maximum
entropy sampling. Operations Research 684–691.
[Krause and Golovin 2014] Krause, A., and Golovin, D.
2014. Submodular function maximization. In Tractability:
Practical Approaches to Hard Problems. Cambridge University Press. 71–104.
[Krause and Horvitz 2008] Krause, A., and Horvitz, E. 2008.
A utility-theoretic approach to privacy and personalization.
In AAAI, 1181–1188.
[Krause, Singh, and Guestrin 2008] Krause, A.; Singh, A.;
and Guestrin, C. 2008. Near-optimal sensor placements
in gaussian processes: theory, efficient algorithms and empirical studies. The Journal of Machine Learning Research
9:235–284.
[Kunegis 2013] Kunegis, J. 2013. Konect - the koblenz network collection. In WWW Companion, 1343–1350.
[Lee et al. 2009] Lee, J.; Mirrokni, V. S.; Nagarajan, V.; and
Sviridenko, M. 2009. Non-monotone submodular maximization under matroid and knapsack constraints. In STOC,
323–332.
[Milgrom and Strulovici 2009] Milgrom, P., and Strulovici,
B. 2009. Substitute goods, auctions, and equilibrium. Journal of Economic Theory 144(1):212–247.
[Mirzasoleiman et al. 2016] Mirzasoleiman, B.; CH, E.; Karbasi, A.; and EDU, Y. 2016. Fast constrained submodular
maximization: Personalized data summarization. In ICLM
’16: Proceedings of the 33rd International Conference on
Machine Learning (ICML).
[Soma and Yoshida 2015] Soma, T., and Yoshida, Y. 2015.
A generalization of submodular cover via the diminishing
return property on the integer lattice. In NIPS, 847–855.
[Soma and Yoshida 2016] Soma, T., and Yoshida, Y. 2016.
Maximizing submodular functions with the diminishing return property over the integer lattice. In IPCO, 325–336.
[Soma et al. 2014] Soma, T.; Kakimura, N.; Inaba, K.; and
Kawarabayashi, K. 2014. Optimal budget allocation: Theoretical guarantee and efficient algorithm. In ICML, 351–359.

A

Proof of Theorem 7

A key observation to obtain an approximation algorithm
with a polynomial time complexity is that the approximate
functions g̃ and h̃ used in Algorithm 3 are piecewise constant
functions. Hence, while x(e) and y(e) lie on the intervals
for which g̃ and h̃, respectively, are constant, the values of
α and β do not change. This means that we repeat the same
random process in the while loop of Algorithm 3 as long as
x(e) and y(e) lie on the intervals. We will show that we can
simulate the entire random process in polynomial time. Because the number of possible values of g̃ and h̃ is bounded
by O( 1 log ∆
δ ), we obtain a polynomial-time algorithm.
As the model of computation, we assume that we can perform an elementary arithmetic operation on real numbers in
constant time, and that we can sample a uniform [0, 1] random variable.

Algorithm 4 Subroutine to simulate random processes for
Algorithm 5
Input: p ∈ [0, 1], `a , `b , `a+b ∈ Z+ , η ∈ (0, 1).
1: q ← 1 −p.

2: N ← O log η1 log(`a + `b + `a+b ) .
3: a ← 0, b ← 0.
4: while a < `a , b < `b , and a + b < `a+b do
5:
while `a − a ≤ N or `a+b − (a + b) < N do
6:
s ← a value sampled from G(q).
7:
if b + s ≤ `b and a + b + s ≤ `a+b then
8:
a ← a + 1, b ← b + s.
9:
else
10:
return (a, min{`b , `a+b − a}).
11:
end if
12:
end while
13:
while `b − b ≤ N do
14:
s ← a value sampled from G(p).
15:
if a + s ≤ `a and a + b + s ≤ `a+b then
16:
a ← a + s, b ← b + 1.
17:
else
18:
return (min{`a , `a+b − b}, b).
19:
end if
20:
end while
21:
n ← min{b `ap−a c, b `bq−b c, `a+b − (a + b)}, m ←
bn/2c.
22:
s ← a value sampled from B(m, p).
23:
a ← a + s, b ← b + m − s.
24:
if a > `a , b > `b , or a + b > `a+b then
25:
Fail.
26:
end if
27: end while
28: return (a, b).

The first two ingredients for simulating the random process are sampling procedures for a binomial distribution and
a geometric distribution. For n ∈ N and p ∈ [0, 1], let
B(n, p) be the binomial distribution with mean np and variance np(1 − p). For p ∈ [0, 1], let G(p) be the geometric
distribution with mean 1/p, that is, PrX∼G(p) [X = k] =
(1 − p)k−1 p for k ≥ 1. We then have the following:
Lemma 8 (See, e.g., (Devroye 1986)). For any n ∈ N, and
p ∈ [0, 1], we can sample a value from the binomial distribution B(n, p) in O(log n) time.
Lemma 9 (See, e.g., (Devroye 1986)). For any p ∈ [0, 1],
we can sample a value from the geometric distribution G(p)
in O(1) time.
We consider the following random process parameterized
by p ∈ [0, 1] and integers `a , `b , `a+b ∈ Z+ , which we denote by P(p, `a , `b , `a+b ): We start with a = b = 0. While
a < `a , b < `b , and a + b < `a+b , we increment a with
probability p, and b with the complement probability 1 − p.
Note that, in the end of the process, we have a = `a , b = `b ,
or a + b = `a+b . Let D(p, `a , `b , `a+b ) be the distribution of
the pair (a, b) generated by P(p, `a , `b , `a+b ).
We introduce an efficient procedure (Algorithm 4) that

succeeds in simulating the process P(p, `a , `b , `a+b ) with
high probability. To prove the correctness of Algorithm 4,
we use the following form of Chernoff’s bound.
Lemma 10 (Chernoff’s bound). Let X1 , . . . , Xn be independent
random variables taking values in {0, 1}. Let X =
Pn
X
i=1 i and µ = E[X]. Then, for any δ > 1, we have
Pr[X ≥ (1 + δ)µ] ≤ exp(−δµ/3).
Lemma 11. We have the following:
• Algorithm 4 succeeds to return a pair (a, b) with probability at least 1 − η.
• The distribution of the pair (a, b) output by Algorithm 4 is
distributed according to D(p, `a , `b , `a+b ).
• The pair (a, b) output by Algorithm 4 satisfies at least one
of the following: a = `a , b = `b , or a + b = `a+b .
• The time complexity of Algorithm 4 is O log( η1 log `) ·

log ` , where ` = `a + `b + `a+b .
Proof. We first note that once (i) `a −a ≤ N , (ii) `b −b ≤ N ,
or (iii) `a+b − (a + b) ≤ N holds, then we enter the while
loop from Line 5 or from 13 until the end of the algorithm.
We check the first claim. Suppose that none of (i), (ii),
and (iii) holds. Then, we reach Line 21. Here, we intend to
simulate the process P(p, `a , `b , `a+b ) to a point where we
will increment a and b m = bn/2c times in total. By union
bound and Chernoff’s bound, the probability that we fail can
be bounded by
h
i
Pr s ≥ `a − a ∨ s ≥ `b − b ∨ s ≥ `a+b − (a + b)
≤ Pr[s ≥ `a − a] + Pr[s ≥ `b − b] + Pr[s ≥ `a+b − (a + b)]
 ` − b − qm 
 ` − a − pm 
b
a
+ exp −
≤ exp −
3
3
+ exp(−`a+b − (a + b) − m)
 ` − a
 ` − b
a
b
≤ exp −
+ exp −
6
6

 `
a+b − (a + b)
+ exp −
2
≤3 exp(−O(N )).
When we do not fail, at least one of the following three
values shrinks by half: `a − a, `b − b, and `a+b − (a + b).
Hence, after O(log(`a + `b + `a+b )) iterations of the while
loop (from Line (4)), at least one of (i), (ii), or (iii) is satisfied. Once this happens, we do not fail and output a pair
(a, b). By union bound, the failure probability is at most
3 exp(−O(N )) · O(log(`a + `b + `a+b )) ≤ δ
By choosing the hidden constant in N large enough.
Next, we check the second claim. From the argument
above, as long as none of (i), (ii), or (iii) is satisfied, we
exactly simulate the process P(p, `a , `b , `a+b ). Hence, suppose that (i) is satisfied. Then, what we does is until a
reaches `a , we sample s from the geometric distribution
G(q) = G(1 − p), and if b + s ≤ `b and a + b + s ≤ `a+b ,
then we update a by a + 1 and b by b + s, and otherwise
we output the pair (a, b = min{`b , `a+b − a}). If a reaches

Algorithm 5 Polynomial time approximation algorithm
E
Input: f : ZE
+ → R+ , B ∈ Z ,  > 0.
Assumption: f is DR-submodular
1: x ← 0, y ← B.
2: for e ∈ E do
3:
Define g, h : {0, 1, . . . , B−1} → R as g(b) = f (χe |
x + bχe ) and h(b) = f (−χe | y − bχe ).
4:
Let g̃ and h̃ be approximations to g and h, respectively, given by Algorithm 2.
5:
while x(e) < y(e) do
6:
α ← g̃(x(e)) and β ← h̃(B(e) − y(e)).
7:
if β = 0 then
8:
x(e) ← y(e) and break.
9:
else if α < 0 then
10:
y(e) ← x(e) and break.
11:
else
12:
s ← max{b | g̃(b) = α}, t ← B(e) − max{b |
h̃(b) = β}.
α
13:
Call Algorithm 4 with p = α+β
, `a = s − x(e),
`b = y(e) − t, `a+b = y(e) − x(e), and η =
2
O( (2+)|E| log(∆/δ) ).
14:
if Algorithm 4 returned a pair (a, b) then
15:
x(e) ← x(e) + a, y(e) ← y(e) + b.
16:
else
17:
return 0.
18:
end if
19:
end if
20:
end while
21: end for
22: return x.

`a , then we output the pair (a = `a , b). This can be seen as
an efficient simulation of the process P(p, `a , `b , `a+b ). The
case (ii) or (iii) is satisfied can be analyzed similarly, and the
second claim holds.
The third and fourth claims are obvious from the definition of the algorithm.
Our idea for simulating Algorithm 3 efficiently is as follows. Suppose we have α = g̃(x(e)) and β = h̃(B(e) −
y(e)) for the current x and y. Let s = max{b | g̃(b) = α}
and t = B(e) − max{b | h̃(b) = β}. Then, g̃ and
h̃ are constant in the intervals [x(e), . . . , s] and [B(e) −
t, . . . , B(e) − y(e)], respectively. By running Algorithm 4
with p = α/(α + β), `a = s − x(e), `b = t − y(e), and
`a+b = y(e) − x(e), we can simulate Algorithm 3 to a point
where at least one of the following happens: x(e) reaches s,
y(e) reaches t, or x(e) is equal to y(e). When Algorithm 4
failed to output a pair, we output an arbitrary feasible solution, say, the zero vector 0. Algorithm 5 presents a formal
description of the algorithm.
Proof of Theorem 7. We first analyze the failure probability. Since the number of possible values of g̃ and h̃ is
bounded by O( 1 log ∆
δ ) for each e ∈ E, we call Algorithm 4
|E|
∆
O(  log δ ) times by the third claim of Lemma 11. Hence,

by the first claim of Lemma 11 and union bound, the fail
ure probability is at most 2+2
if the hidden constant in η at
Line 13 is chosen to be small enough.
Let D and D0 be the distributions of outputs from Algorithms 3 and 5, respectively. Conditioned on the event that
Algorithm 4 does not fail (and hence we output f (0)), D
exactly matches D0 by the second claim of Lemma 11.
By letting o be the optimal solution, we have, by Theorem 5, that
2+
1
· E f (x) ≥
f (o).
E f (x) ≥
2 + 2 x∼D
2 + 2
x∼D 0
1
Hence, we have the approximation factor of 2+
.
The number of oracle calls is exactly the same as Algorithm 3. The total time spent inside Algorithm 4 is
  |E|
 |E|


∆
∆
· O log
O
log
log log kBk∞ log kBk∞

δ

δ
 |E|

∆
e
=O
log log kBk∞

δ
by the fourth claim of Lemma 11. Because evaluating g̃ and
h̃ takes O(log kBk∞ ) time and computing s and t also takes
O(log kBk∞ ) time, we need
 |E|

 |E|
∆
∆
log
·O(log2 kBk∞ ) = O
log log2 kBk∞
O

δ

δ
time. Summing them up, we obtain the stated time complexity.
We again note that, even if the given DR-submodular
function f : ZE
+ → R is not non-negative, we can obtain
1
a 2+2
-approximation algorithm for (5), as stated in Remark 6.

B

DR-submodularity of functions used in
experiments

In this section, we will see that the objective function used
in Section 4 is indeed DR-submodular. Recall that our objective function of revenue maximization is as follows:
X X
f (x) =
wij (1 − (1 − p)x(i) )(1 − p)x(j) ,
i∈S j∈V \S

where wij is a nonnegative weight and p ∈ [0, 1] is a parameter. Since DR-submodular functions are closed under
nonnegative linear combination, it suffices to check that
g(x) = (1 − q x(i) )q x(j)
is DR-submodular, where q = 1 − p. To see the DRsubmodularity of g, we need to check that
g(χi | x + χj ) ≤ g(χi | x) (x ∈ RE
+ , i, j ∈ E).
Note that i and j may be identical. By direct algebra,
g(χi | x) = (1 − q x(i)+1 )q x(j) − (1 − q x(i) )q x(j)
= q x(i)+x(j) (1 − q),
g(χi | x + χj ) = q x(i)+x(j)+1 (1 − q)
= qg(χi | x).
Since q ∈ [0, 1], we obtain g(χi | x + χj ) ≤ g(χi | x).

