Maximizing Expected Utility for Stochastic Combinatorial
Optimization Problems ∗
Jian Li†1 and Amol Deshpande‡2

arXiv:1012.3189v7 [] 10 Aug 2016

1

Institute for Interdisciplinary Information Sciences, Tsinghua University, Beijing, P.R.China
2
Department of Computer Science, University of Maryland, College Park, USA

Abstract
We study the stochastic versions of a broad class of combinatorial problems where the weights of the
elements in the input dataset are uncertain. The class of problems that we study includes shortest paths,
minimum weight spanning trees, and minimum weight matchings, and other combinatorial problems like
knapsack. We observe that the expected value is inadequate in capturing different types of risk-averse
or risk-prone behaviors, and instead we consider a more general objective which is to maximize the
expected utility of the solution for some given utility function, rather than the expected weight (expected
weight becomes a special case). Under the assumption that there is a pseudopolynomial time algorithm
for the exact version of the problem (This is true for the problems mentioned above), 1 we can obtain the
following approximation results for several important classes of utility functions:
1. If the utility function µ is continuous, upper-bounded by a constant and limx→+∞ µ(x) = 0, we
show that we can obtain a polynomial time approximation algorithm with an additive error  for
any constant  > 0.
2. If the utility function µ is a concave increasing function, we can obtain a polynomial time approximation scheme (PTAS).
3. If the utility function µ is increasing and has a bounded derivative, we can obtain a polynomial
time approximation scheme.
Our results recover or generalize several prior results on stochastic shortest path, stochastic spanning
tree, and stochastic knapsack. Our algorithm for utility maximization makes use of the separability
of exponential utility and a technique to decompose a general utility function into exponential utility
functions, which may be useful in other stochastic optimization problems.

1

Introduction

The most common approach to deal with optimization problems in presence of uncertainty is to optimize
the expected value of the solution. However, expected value is inadequate in expressing diverse people’s
∗

A preliminary version of the paper appeared in the Proceedings of the 52nd Annual IEEE Symposium on Foundations of
Computer Science (FOCS), 2011.
†
lijian83@mail.tsinghua.edu.cn
‡
amol@cs.umd.edu
1
Following the literature [55], we differentiate between exact version and deterministic version of a problem; in the exact version
of the problem, we are given a target value and asked to find a solution (e.g., a path) with exactly that value (i.e., path length).

1

preferences towards decision-making under uncertain scenarios. In particular, it fails at capturing different
risk-averse or risk-prone behaviors that are commonly observed. Consider the following simple example
where we have two lotteries L1 and L2 . In L1 , the player could win 1000 dollars with probability 1.0, while
in L2 the player could win 2000 dollars with probability 0.5 and 0 dollars otherwise. It is easy to see that
both have the same expected payoff of 1000 dollars. However, many, if not most, people would treat L1
and L2 as two completely different choices. Specifically, a risk-averse player is likely to choose L1 and a
risk-prone player may prefer L2 (Consider a gambler who would like to spend 1000 dollars to play doubleor-nothing). A more involved but also more surprising example is the St. Petersburg paradox (see e.g., [45])
which has been widely used in the economics literature as a criticism of expected value. The paradox is
named from Daniel Bernoulli’s presentation of the problem, published in 1738 in the Commentaries of the
Imperial Academy of Science of Saint Petersburg. Consider the following game: you pay a fixed fee X to
enter the game. In the game, a fair coin is tossed repeatedly until a tail appears ending the game. The payoff
of the game is 2k where k is the number of heads that appear, i.e., you win 1 dollar if a tail appears on the
first toss, 2 dollars if a head appears on the first toss and a tail on the second, 4 dollars if a head appears on
the first two tosses and a tail on the third and so on. The question is what would be a fair fee X to enter the
game? First, it is easy to see that the expected payoff is
∞

E[payoff] =

X1
1
1
1
1
1 1 1 1
·1+ ·2+ ·4+
· 8 + ··· = + + + + ··· =
=∞
2
4
8
16
2 2 2 2
2
k=1

If we use the expected payoff as a criterion for decision making, we should therefore play the game at any
finite price X (no matter how large X is) since the expected payoff is always larger. However, researchers
have done extensive survey and found that not many people would pay even 25 dollars to play the game
[45], which significantly deviates from what the expected value criterion predicts. In fact, the paradox can
be resolved by expected utility theory with a logarithmic utility function, suggested by Bernoulli himself [7].
We refer interested reader to [59, 45] for more information. These observations and criticisms have led
researchers, especially in Economics, to study the problem from a more fundamental perspective and to
directly maximize user satisfaction, often called utility. The uncertainty present in the problem instance
naturally leads us to optimize the expected utility.
Let F be the set of feasible solutions to an optimization problem. Each solution S ∈ F is associated
with a random weight w(S). For instance, F could be a set of lotteries and w(S) is the (random) payoff of
lottery S. We model the risk awareness of a user by a utility function µ : R → R: the user obtains µ(x)
units of utility if the outcome is x, i.e., w(S) = x. Formally, the expected utility maximization principle is
simply stated as follows: the most desirable solution S is the one that maximizes the expected utility, i.e.,
S = arg max
E[µ(w(S 0 ))]
0
S ∈F

Indeed, expected utility theory is a branch of utility theory that studies “betting preferences” of people with
regard to uncertain outcomes (gambles). The theory was formally initiated by von Neumann and Morgenstern in 1940s [65, 24] 2 who gave an axiomatization of the theory (known as von Neumann-Morgenstern
expected utility theorem). The theory is well known to be versatile in expressing diverse risk-averse or
risk-prone behaviors.
In this paper, we focus on the following broad class of combinatorial optimization problems. The deterministic version of the problem has the following form: we are given a ground set of elements U =
2

Daniel Bernoulli also developed many ideas, such as risk aversion and utility, in his work Specimen theoriae novae de mensura
sortis (Exposition of a New Theory on the Measurement of Risk) in 1738 [8].

2

{ei }i=1...n ; each element e is associated with a weight we ; each feasible solution is a subset of the elements
satisfying some property. Let F denote the set of feasible solutions. The objective for the deterministic
P
problem is to find a feasible solution S with the minimum (or maximum) total weight w(S) = e∈S we .
We can see that many combinatorial problems such as shortest path, minimum spanning tree, and minimum
weight matching belong to this class. In the stochastic version of the problem, the weight we of each element e is a nonnegative random variable. We assume all we s are independent of each other. We use pe (.)
to denote the probability density function for we (or probability mass function in the discrete case). We are
also given a utility function µ : R+ → R+ which maps a weight value to a utility value. By the expected
utility maximization principle, our goal here is to find a feasible solution S ∈ F that maximizes the expected
utility, i.e., E[µ(w(S))]. We call this problem the expected utility maximization (EUM) problem.
Let us use the following toy example to illustrate the rationale behind EUM. There is a graph with two
nodes s and t and two parallel links e1 and e2 . Edge e1 has a fixed length 1 while the length of e2 is 0.9 with
probability 0.9 and 1.9 with probability 0.1 (the expected value is also 1). We want to choose one edge to
connect s and t. It is not hard to imagine that a risk-averse user would choose e1 since e2 may turn out to
be a much larger value with a nontrivial probability. We can capture such behavior using the utility function
(1) (defined in Section 1.1). Similarly, we can capture the risk-prone behavior by using, for example, the
1
utility function µ(x) = x+1
. It is easy to see that e1 maximizes the expected utility in the former case, and
e2 in the latter.

1.1

Our Contributions

In order to state our contribution, we first recall some standard terminologies. A polynomial time approximation scheme (PTAS) is an algorithm which takes an instance of a minimization problem (a maximization
problem resp.) and a parameter  > 0 and produces a solution whose cost is at most (1 + )OPT (at least
(1 − )OPT resp.), and the running time, for any fixed constant  > 0, is polynomial in the size of the
input, where OPT is the optimal solution. We use A to denote the deterministic combinatorial optimization
problem under consideration, and EUM(A) the corresponding expected utility maximization problem. The
exact version of A asks the question whether there is a feasible solution of A with weight exactly equal to
a given integer K. We say an algorithm runs in pseudopolynomial time for the exact version of A if the
running time is polynomial in n and K. For many combinatorial problems, a pseudopolynomial algorithm
for the exact version is known. Examples include shortest path, spanning tree, matching and knapsack.
We discuss in detail our results for EUM. We start with a theorem which underpins our other results. We
denote kµk∞ = supx≥0 |µ(x)|. We say a function µ
e(x) is an -approximation of µ(x) if |e
µ(x) − µ(x)| ≤
kµk∞ for all x ≥ 0. We allow µ
e(x) to be a complex function and |e
µ(x)| denote its absolute value (as we
will see shortly, µ
e(x) takes the form of a finite sum of complex exponentials). 3
Theorem 1 Assume that there is a pseudopolynomial algorithm for the exact version of A. Further assume
e(x) =
PL that xgiven any constant  > 0, we can find an -approximation of the utility function µ as µ
k=1 ck φk , where |φk | ≤ 1 for all 1 ≤ k ≤ L (φk may be complex numbers). Let τ = maxk |ck |/kµk∞ .
Then, there is an algorithm that runs in time (nτ /)O(L) and finds a feasible solution S ∈ F such that
E[µ(w(S))] ≥ OPT − kµk∞ .
3

In practice, the user only needs to specify a real utility function µ(x). The complex function µ
e(x) is used to approximate the
real utility function µ(x).

3

From the above theorem, we can see that if we can -approximate the utility function µ by a short sum
of exponentials, we can obtain good approximation algorithms for EUM. In this paper, we consider three
important classes of utility functions.
1. (Class Cbounded ) Consider the deterministic problem which A is a minimization problem, i.e., we
would like the cost of our solution to be as small as possible. In the corresponding stochastic version
of A, we assume that any utility function µ(x) ∈ Cbounded is nonnegative, bounded, continuous and
limx→∞ µ(x) = 0 (please see below for the detailed technical assumptions). The last condition
captures the fact that if the cost of solution is too large, it becomes almost useless for us. We denote
the class of such utility functions by Cbounded .
2. (Class Cconcave ) Consider the deterministic problem A which is a maximization problem. In other
words, we want the value of our solution to be as large as possible. In the corresponding stochastic
version of A, we assume that µ(x) is a nonnegative, monotone nondecreasing and concave function.
Note that concave functions are extensively used to model risk-averse behaviors in the economics
literature. We denote the class of such utility functions by Cconcave .
3. (Class Cincreasing ) Consider a deterministic maximization problem A. In the corresponding stochastic
version of A, we assume that µ(x) is a nonnegative, differentiable and increasing function. We assume
d
dx µ(x) ∈ [L, U] for x ≥ 0, where L, U > 0 are constants. We denote the class of such utility
functions by Cincreasing . We can see that functions in Cincreasing can be concave, nonconcave, convex
or nonconvex. Convex functions are often associated with risk-prone behaviors, while nonconvexnonconcave utility functions have been also observed in various settings [36, 23].
Now, we state in details our assumptions and results for the above classes of utility functions.
Class Cbounded : Since µ is bounded, by scaling, without loss of generality, we can assume kµk∞ = 1. Since
limx→∞ µ(x) = 0, for any  > 0, there exist a point T such that µ(x) ≤  for x > T . We assume that
T is a constant only depending on . We further assume that the continuous utility function µ satisfies the
α-Hölder condition, i.e., |µ(x) − µ(y)| ≤ C |x − y|α , for some constant C and some constant α > 1/2. We
say f is C-Lipschitz if f satisfies 1-Hölder condition with coefficient C. Under the above conditions, we
can prove Theorem 2.
Theorem
P2 If the xutility function µ belongs to Cbounded , then, for any  > 0, we can obtain a function
µ(x) − µ(x)| ≤ , for x ≥ 0,where
µ
e(x) = L
k=1 ck φk , such that |e
L = 2O(T ) poly(1/),

|ck | ≤ 2O(T ) poly(1/),

|φk | ≤ 1 for all k = 1, . . . , L,

To show the above theorem, we use the Fourier series technique. However, the technique cannot be
used directly since it works only for periodic functions with bounded periodicities. In order to get a good
approximation for x ∈ [0, ∞), we leverage the fact that limx→∞ µ(x) = 0 and develop a general framework
that uses the Fourier series decomposition as a subroutine.
Now, we state some implications of the above results. Consider the utility function

x ∈ [0, 1]
 1
x
1
− + δ + 1 x ∈ [1, 1 + δ]
χ
e(x) =
(1)
 δ
0
x>1+δ

4

Figure 1: (1) The utility function χ
e(x), a continuous variant of the threshold function χ(x); (2) A smoother
variant of χ(x); (3) The utility function χ
e2 (x), a continuous variant of the 2-d threshold function χ2 (x).
where δ > 0 is a small constant (See Figure 1(1)). It is easy to verify that χ
e is 1/δ-Lipschitz and T = 2
for any δ < 1. Therefore, Theorem 2 is applicable. This example is interesting since χ
e can be viewed as a
continuous variant of the threshold function

1 x ∈ [0, 1]
χ(x) =
,
(2)
0 x>1
for which maximizing the expected utility is equivalent to maximizing Pr(w(S) ≤ 1). We first note that
even the problem of computing the probability Pr(w(S) ≤ 1) exactly for a fixed set S is #P-hard [38] and
there is an FPTAS [42]. Designing approximation algorithms for such special case has been considered
several times in the literature for various combinatorial problems including stochastic shortest path [52],
stochastic spanning tree [35, 26], stochastic knapsack [27] and some other stochastic problems [2, 50].
It is interesting to compare our result with the result for the stochastic shortest path problem considered
by Nikolova et al. [52, 50]. In [52], they show that there is an exact O(nlog n ) time algorithm for maximizing
the probability that the length of the path is at most 1, i.e., Pr(w(S) ≤ 1), assuming all edges are normally
distributed and there is a path with its mean at most 1. Later, Nikolova [50] extends the result to an FPTAS for
any problem under the same assumptions, if the deterministic version of the problem has a polynomial time
exact algorithm. We can see that under such assumptions, the optimal probability is at least 1/2. 4 Therefore,
provided the same assumption and further assuming that Pr(we < 0) is miniscule, 5 our algorithm is a PTAS
for maximizing E[e
χ(w(S))], which can be thought as a variant of the problem of maximizing E[χ(w(S))].
Indeed, we can translate this result to a bi-criterion approximation result of the following form: for any fixed
constants δ,  > 0, we can find in polynomial time a solution S such that
Pr(w(S) ≤ 1 + δ) ≥ (1 − ) Pr(w(S ∗ ) ≤ 1).
where S ∗ is the optimal solution (Corollary 2). We note that such a bi-criterion approximation was only
known for exponentially distributed edges before [52].
Let us consider another application of our results to the stochastic knapsack problem defined in [27].
Given a set U of independent random variables {x1 , . . . , xn }, with associated profits {v1 , . . . , vn } and an
overflow probability γ, we are asked to pick a subset S of U such that
!
X
Pr
xi ≥ 1 ≤ γ
i∈S
4

The sum of multiple Gaussians is also a Gaussian. Hence, if we assume the mean of the length of a path (which is a Gaussian)
is at most 1, the probability that the length of the path is at most 1 is at least 1/2.
5
Our technique can only handle distributions with positive supports. Thus, we have to assume that the probability that a negative
value appears is miniscule (e.g., less than 1/n2 ) and can be safely ignored (because the probability that there is any realized negative
value is at most 1/n).

5

P
and the total profit i∈S vi is maximized. Goel and Indyk [27] showed that, for any constant  > 0,
there
Pis a polynomial time algorithm that can find a solution S with the profit as least the optimum and
Pr( i∈S xi ≥ 1 + ) ≤ γ(1 + ) for exponentially distributed variables. They also gave a quasi-polynomial
time approximation scheme for Bernoulli distributed random variables. Quite recently, in parallel with
our work, Bhalgat et al. [13] obtained the same result for arbitrary distributions under the assumption that
γ = Θ(1). Their technique is based on discretizing the distributions and is quite involved. 6 Our result,
applied to stochastic knapsack, matches that of Bhalgat et al. under the same assumption. Our algorithm is
arguably simpler and has a much better running time (Theorem 7).
Equally importantly, we can extend our basic approximation scheme to handle generalizations such as
multiple utility functions and multidimensional weights. Interesting applications of these extensions include
various generalizations of stochastic knapsack, such as stochastic multiple knapsack (Theorem 10) and
stochastic multidimensional knapsack (stochastic packing) (Theorem 11).
Class Cconcave : We assume the utility function µ : [0, ∞) → [0, ∞) is a concave, monotone nondecreasing
function. This is a popular class of utility functions used to model risk-averse behaviors. For this class of
utility functions, we can obtain the following theorem in Section 5.
Theorem 3 Assume the utility function µ belongs to Cconcave , and there is a pseudopolynomial algorithm
for the exact version of A. Then, there is a PTAS for EUM(A).
Theorem 3 is also obtained by an application of Theorem 1. However, instead of approximating the
original utility function µ using a short sum of exponentials, which may not be possible in general, 7 we try
approximate a truncated version of µ. Theorem 3 recovers the recent result of [14]. Finally, we remark the
technique of [14] strongly relies on the concavity of µ, and seems difficult to extend to handle non-concave
utility functions.
Class Cincreasing : We assume the utility function µ : [0, ∞) → [0, ∞) is a positive, differentiable, and
d
increasing function. For technical reasons, we assume dx
µ(x) ∈ [L, U] for some constants L, U > 0 and all
x ≥ 0. For this class of utility functions, we can obtain the following theorem in Section 6.
Theorem 4 Assume the utility function µ belongs to Cincreasing , and there is a pseudopolynomial algorithm
for the exact version of A. Then, there is a PTAS for EUM(A).
Again, it may not be possible in general to approximate such an increasing function using a finite sum
of exponentials. Instead, we approximate a truncated version of µ, similar to the concave case. We note this
is the first such result for general increasing utility functions. Removing the bounded derivative assumption
remains an interesting open problem.
We believe our technique can be used to handle other classes of utility functions or other stochastic
optimization problems.

1.2

Related Work

In recent years stochastic optimization problems have drawn much attention from the computer science community and stochastic versions of many classical combinatorial optimization problems have been studied.
6
They also obtain several results related to stochastic knapsack, using the their discretization technique, together with other
ideas. Notably, they obtained a bi-criteria PTAS for the adaptive stochastic knapsack problem [13].
7
Suppose µ is a finite sum of exponentials. When x approaches to infinity, either |µ(x)| is periodic, or approaches to infinity, or
approaches to 0.

6

In particular, a significant portion of the efforts has been devoted to the two-stage stochastic optimization
problem. In such a problem, in a first stage, we are given probabilistic information about the input but the
cost of selecting an item is low; in a second stage, the actual input is revealed but the costs for the elements
are higher. We are asked to make decision after each stage and minimize the expected cost. Some general
techniques have been developed [31, 60]. We refer interested reader to [64] for a comprehensive survey.
Another widely studied type of problems considers designing adaptive probing policies for stochastic optimization problems where the existence or the exact weight of an element can be only known upon a probe.
There is typically a budget for the number of probes (see e.g., [30, 19]), or we require an irrevocable decision whether to include the probed element in the solution right after the probe (see e.g., [22, 17, 4, 21, 13]).
However, most of those works focus on optimizing the expected value of the solution. There is also sporadic
work on optimizing the overflow probability or some other objectives subject to the overflow probability
constraints. In particular, a few recent works have explicitly motivated such objectives as a way to capture
the risk-averse type of behaviors [2, 50, 63]. Besides those works, there has been little work on optimizing
more general utility functions for combinatorial stochastic optimization problems from an approximation
algorithms perspective.
The most related work to ours is the stochastic shortest path problem (Stoch-SP), which was also the
initial motivation for this work. The problem has been studied extensively for several special utility functions
in operation research community. Sigal et al. [61] studied the problem of finding the path with greatest
probability of being the shortest path. Loui [44] showed that Stoch-SP reduces to the shortest path (and
sometimes longest path) problem if the utility function is linear or exponential. Nikolova et al. [51] identified
more specific utility and distribution combinations that can be solved optimally in polynomial time. Much
work considered dealing with more general utility functions, such as piecewise linear or concave functions,
e.g., [48, 49, 6]. However, these algorithms are essentially heuristics and the worst case running times are
still exponential. Nikolova et al. [52] studied the problem of maximizing the probability that the length of the
chosen path is less than some given parameter. Besides the result we mentioned before, they also considered
Poisson and exponential distributions. Despite much effort on this problem, no algorithm is known to run
in polynomial time and have provable performance guarantees, especially for more general utility functions
or more general distributions. This is perhaps because the hardness comes from different sources, as also
noted in [52]: the shortest path selection per se is combinatorial; the distribution of the length of a path is
the convolution of the distributions of its edges; the objective is nonlinear; to list a few.
Kleinberg et al. [38] first considered the stochastic knapsack problem with Bernoulli-type distributions
and provided a polynomial-time O(log 1/γ) approximation where γ is the given overflow probability. In
the same paper, they noticed that even computing the overflow probability for a fixed set of items is #P-hard.
Li and Shi [42] provided an FPTAS for computing the overflow probability (or the threshold probability
for a sum of random variables). For item sizes with exponential distributions, Goel and Indyk [27] provided a bi-criterion PTAS, and for Bernoulli-distributed items they gave a quasi-polynomial approximation
scheme. Chekuri and Khanna [16] pointed out that a PTAS can be obtained for the Bernoulli case using
their techniques for the multiple knapsack problem. Goyal and Ravi [29] showed a PTAS for Gaussian distributed sizes. Bhalgat, Goel and Khanna [13] developed a general discretizaton technique that reduces the
distributions to a small number of equivalent classes which we can efficiently enumerate for both adaptive
and nonadaptive versions of stochastic knapsack. They used this technique to obtain improved results for
several variants of stochastic knapsack, notably a bi-criterion PTAS for the adaptive version of the problem.
In a recent work [43], the bi-criterion PTAS was further simplified and extended to the more general case
where the profit and size of an item can be correlated and an item can be cancelled in the middle. Dean at
al. [22] gave the first constant approximation for the adaptive version of stochastic knapsack. The adaptive

7

version of stochastic multidimensional knapsack (or equivalently stochastic packing) has been considered
in [21, 13] where constant approximations and a bi-criterion PTAS were developed.
This work is partially inspired by our prior work on top-k and other queries over probabilistic datasets [39,
41]. In fact, we can show that both the consensus answers proposed in [39] and the parameterized ranking
functions proposed in [41] follow the expected utility maximization principle where the utility functions
are materialized as distance metrics for the former and the weight functions for the latter. Our technique
for approximating the utility functions is also similar to the approximation scheme used in [41] in spirit.
However, no performance guarantees are provided in that work.
Recently, Li and Yuan [43] showed that an additive PTAS for µ ∈ Cbounded can be obtained using a
completely different approach, called the Poisson approximation technique. Roughly speaking, the Poisson
approximation technique allows us to extract a constant (depending on ) number of features from each distribution (called signature in [43]) and reduce the stochastic problem to a constant dimensional deterministic
optimization problem, which is similar to the algorithm presented in this paper. We suspect that besides this
superficial similarity, there may be deeper connections between two different techniques.
There is a large volume of work on approximating functions using short exponential sums over a
bounded domain, e.g., [54, 9, 10, 11]. Some works also consider using linear combinations of Gaussians or
other kernels to approximate functions with finite support over the entire real axis (−∞, +∞) [18]. This
is however impossible using exponentials since αx is either periodic (if |α| = 1) or approaches to infinity
when x → +∞ or x → −∞ (if |α| =
6 1).

2

An Overview of Our Approach

The high level idea of our approach is very simple and consists of the following steps:
1. We first observe that the problem is easy if the utility function is an exponential function. Specifically,
consider the exponential utility function µ(x) = φx for some complex number φ ∈ C. Fix an arbitrary
solution S. Due to independence of the elements, we can see that
hY
i Y
 P

E[φw(S) ] = E φ e∈S we = E
φwe =
E[φwe ]
e∈S

e∈S

we
Taking log on both sides, we get log E[φw(S) ] =
e∈S log E[φ ]. If φ is a positive real number
and E[φwe ] ≤ 1 (or equivalently, − log E[φwe ] ≥ 0), this reduces to the deterministic optimization
problem.

P

2. In light of P
the above observation, we -approximate the utility function µ(x) by a short exponential
x
sum, i.e., L
i=1 ci φi with L being a small value (only depending on ), where (ci and φi may be
P
w(S)
complex numbers. Hence, E[µ(w(S))] can be approximated by L
].
i=1 ci E[φi
w(S)

3. Consider the following multi-criterion version of the problem with L objectives {E[φi ]}i=1,...,L :
w(S)
given L complex numbers v1 , . . . , vL , we want to find a solution S such that E[φi ] ≈ vi for i =
1, . . . , L. We achieve this by utilizing the pseudopolynomial time algorithm for the exact version of
the problem. We argue that we only need to consider a polynomial number of v1 , . . . , vL combinations
(which we call configurations) to find out the approximate optimum.
In Section 3, we show how to solve the multi-criterion problem provided that a short exponential sum
approximation of µ is given. In particular, we prove Theorem 1. Then, we show how to approximate
8

µ ∈ Cbounded by a short exponential sum by proving Theorem 2 in Section 4.1 and Section 4.2. For µ ∈
Cconcave or µ ∈ Cincreasing , it may not be possible to approximate µ directly by an exponential sum, and
some additional ideas are required. The details are provided in Section 5 and Section 6.
We still need to show how to compute E[φwe ]. If we is a discrete random variable with a polynomial
size support, we can easily compute E[φwe ] in polynomial time. If we has an infinite discrete or continuous
support, we can not compute E[φwe ] directly and need to approximate it. We briefly discuss this issue and
its implications in Appendix A.

3

Proof of Theorem 1

Now, we prove Theorem 1. We start with some notations. We use |c| and arg(c) to denote the absolute
value and the argument of the complex number c ∈ C, respectively. In other words, c = |c| · (cos(arg(c)) +
i sin(arg(c))) = P
|c|ei arg(c) . We always require arg(c) ∈ [0, 2π) for any c ∈ C. Recall that we say the
x
exponential sum L
i=1 ci φi is an -approximation for µ(x) if the following holds:
|µ(x) −

L
X

ci φxi | ≤ kµk∞

for x ≥ 0.

i=1

We first show that if the utility function can be decomposed exactly into a short exponential sum, we can
approximate the optimal expected utility well.
P
x
Theorem 5 Assume that µ
e(x) = L
k=1 ck φk is the utility function where |φk | ≤ 1 for 1 ≤ k ≤ L. Let
τ = maxk |ck |/kµk∞ . We also assume that there is a pseudopolynomial algorithm for the exact version of
A. Then, for any  > 0, there is an algorithm that runs in time (n/)O(L) and finds a solution S such that
e
|E[e
µ(w(S))] − E[e
µ(w(S))]|
< kµk∞ ,
where Se = arg maxS 0 |E[e
µ(w(S 0 ))|.
We use the scaling and rounding technique that has been used often in multi-criterion optimization
problems (e.g., [58, 55]). Since our objective function is not additive and not monotone, the general results
for multi-criterion optimization [55, 46, 58, 1] do not directly apply here. We provide the details of the
algorithm here. We use the following parameters:

 


− ln(/Lτ )n
2πn
γ=
, J = max
,
.
Lnτ
γ
γ
Let V be the set of all 2L-dimensional integer vectors of the form v = hx1 , y1 , . . . , xL , yL i where 1 ≤ xi ≤
J and 1 ≤ yi ≤ J for i = 1, . . . , L.
For each element e ∈ U , we associate it with a 2L-dimensional integer vector
Ft(e) = hα1 (e), β1 (e), . . . , αL (e), βL (e)i,

where




e
− ln |E[φw
i ]| J
αi (e) = min
,
γ
n

9


e
arg(E[φw
i ])
βi (e) =
.
γ


and

(3)

We call Ft(e) the feature vector of e. Since
P |φi | ≤ 1, we can see that αi (e) ≥ 0 for any e ∈ U . It is easy
to see that Ft(e) ∈ V for all e ∈ U and e∈S Ft(e) ∈ V for all S ⊆ U . Intuitively, αi (e) and βi (e) can be
we
e
thought as the scaled and rounded versions of − ln |E[φw
i ]| and arg(E[φi ]), respectively.
2L
O(L)
We maintain J = (n/)
configurations (a configuration is just like a state in a dynamic program).
Each configuration Cf(v) is indexed by a 2L-dimensional vector v ∈ V and takes 0/1 value. In particular,
the value of Cf(v) for each v ∈ V is defined as follows: For each vector v ∈ V,
P
1. Cf(v) = 1 if and only if there is a feasible solution S ∈ F such that e∈S Ft(e) = v.
2. Cf(v) = 0 otherwise.
For any v = hx1 , y1 , . . . , xL , yL i, define the value of v to be
Val(v) =

L
X

ck e−xk γ+iyk γ .

k=1

Lemma 1 tells us the value of a configuration is close to the expected utility of the corresponding solution. Lemma 2 shows we can compute those configurations in polynomial time.
P
x
Lemma 1 Suppose µ
e(x) = L
k=1 ck φk , where |φk | ≤ 1 for all k = 1, . . . , L. Let τ = maxk |ck |/kµk∞ .
For any vector v = hx1 , y1 , . . . , xL , yL i ∈ V, Cf v (v) = 1 if and only if there is a feasible solution S ∈ F
such that
L
X
E[e
µ(w(S))] − Val(v) = E[e
µ(w(S))] −
ck e−xk γ+iyk γ ≤ O(kµk∞ ).
k=1

PL
P
w(S)
w(S)
]. Therefore, it suffices
] =
Proof: We first notice that E[e
µ(w(S))] = E[ L
k=1 ck E[φk
k=1 ck φk
w(S)

−x
γ+iy
γ
k
k
to
| ≤ O( Lτ ). Since Cf(v)
P show that for all k = 1, . . . , L, |E[φk ] − e
P = 1, we know that
Ft(e)
=
v
for
some
feasible
solution
S
∈
F.
In
other
words,
we
have
e∈S αk (e) = xk and
Pe∈S
β
(e)
=
y
for
all
1
≤
k
≤
L.
k
e∈S k
w(S)

Fix an arbitrary 1 ≤ k ≤ L. First, we can see that the arguments of E[φk
w(S)

arg(E[φk

]) − yk γ ≤

X

e
arg(E[φw
k ]) − βk (e)γ ≤

X

] and e−xk γ+iyk γ are close:

γ ≤ nγ =

e∈S

e∈S


,
Lτ

where we use arg(c) to denote the argument of the complex number c. Now, we show the magnitude of
w(S)
E[φk ] and e−xk γ+iyk γ are also close. We distinguish two cases:
j

k
e
e
− ln |E[φw
− ln |E[φw
i ]| J
i ]|
1. Recall that αi (e) = min
,
. If there is some e ∈ S such that
> Jn (which
γ
n
γ
implies that αk (e) = b Jn c), we know that
w(S)

− ln(|E[φk

]|) =

X

e
(− ln(|E[φw
k |)) >

e∈S

In this case, we have xk =

P

e∈S

w(S)

|E[φk

Jγ
.
n

αk (e) ≥ Jn . Thus, we have that
γd

]| − |e−xk γ | < e−Jγ/n ≤ e
10

n ln(/Lτ )
e/n
γ

<


.
Lτ

2. On the other hand, if αk (e) =
w(S)

− ln(|E[φk

j

e
− ln |E[φw
k ]|
γ

)|) − xk γ =

X

k

for all e ∈ S, we can see that

(− ln(|E[φwe |) − αk (e)γ) ≤

e∈S

X

γ ≤ nγ ≤

e∈S


.
Lτ

Since the derivative of ex is less than 1 for x < 0, we can get that
w(S)

|E[φk

]| − |e−xk γ | ≤ |e−xk γ−/Lτ − e−xk γ | ≤


.
Lτ

For any two complex numbers a, b with |a| ≤ 1 and |b| ≤ 1, if |a| − |b| <  and |∠ab| = | arg(a) −
arg(b)| < , we can see that
|a − b|2 = |a|2 + |b|2 − 2|a||b| cos(∠ab)
= (|a| − |b|)2 + 2|a||b|(1 − cos(∠ab))
≤ 2 + 2(1 − cos(∠ab)2 )
≤ 2 + 2 sin(∠ab)2
≤ 2 + 2| arg(a) − arg(b)|2 ≤ 32 .
In the third inequality, we use the fact that sin x < x for all x > 0. The proof is completed.

2

Lemma 2 Suppose there is a pseudopolynomial time algorithm for the exact version of A, which runs in
time polynomial in n and t (t is the maximum integer in the instance of A). Then, we can compute the values
O(L) .
for all configurations {Cf(v)}v∈V in time ( nτ
 )
Proof: For each vector v ∈ V, we can encode it as a nonnegative integer I(v) upper bounded by J 2L =
( n )O(L) . In particular, each coordinate of v takes the position of a specific digit in the integral representation, and the base is chosen to be J no carry can occur when we add at most n feature vectors. Then,
determining the value of a configuration Cf(v) P
is equivalent to determining whether there is a feasible solution S ∈ F such that the total weight of S (i.e., e∈S I(Ft(e))) is exactly the given value I(v). Suppose the
pseudopolynomial time algorithm for the exact version of A runs in time PA (n, t) for some polynomial PA .
Therefore, the value of each such Cf(v) can be also computed in time PA (n, I(v)) = PA (n, ( n )O(L) ) =
O(1) , the number of configuration is ( nτ )O(L) . The total running time
( n )O(L) . Since J are bounded by ( nτ
 )

O(L) × ( nτ )O(L) = ( nτ )O(L) .
is ( nτ
)
2



Now, everything is ready to prove Theorem 5.
Proof of Theorem 5: We first use the algorithm in Lemma 2 to compute the values for all configurations.
Then, we findPthe configuration Cf(hx1 , y1 , . . . , xL , yL i) that has value 1 and that maximizes the quantity
−xk γ+iyk γ |. The feasible solution S corresponding to this configuration is our final
|Val(v)| = | L
k=1 ck e
solution. It is easy to see that the theorem follows from Lemma 1.
2
Theorem 1 can be readily obtained from Theorem 5 and the fact µ
e is an -approximation of µ.
Proof ofP
Theorem 1: Suppose S is our solution and S ∗ is the optimal solution for utility function µ. Recall
x
µ
e(x) = L
k=1 ck φk . From Theorem 5, we know that
|E[e
µ(w(S))]| ≥ |E[e
µ(w(S ∗ ))]| − O(kµk∞ ).
11

Since µ
e is an -approximation of µ, we can see that
Z
Z
(µ(x) − µ
e(x))dPS (x) ≤
kµk∞ dPS (x) ≤ kµk∞ .
E[µ(w(S))] − E[e
µ(w(S))] =
for any solution S, where PS is the probability measure of w(S). Therefore, we have
|E[µ(w(S))]| ≥ |E[e
µ(w(S))]| − kµk∞ ≥ |E[e
µ(w(S ∗ ))]| − O(kµk∞ )
≥ |E[µ(w(S ∗ ))]| − O(kµk∞ ).
2

This completes the proof of Theorem 1.

4

Class Cbounded

The main goal of this section is to prove Theorem 2. In Section 4.1, we develop a generic algorithm that
takes as a subroutine an algorithm F OURIER for approximating functions in a bounded interval domain,
and approximates µ(x) ∈ Cbounded in the infinite domain [0, +∞). In the Section 4.2, we use the Fourier
series expansion as the choice of F OURIER and show that important classes of utility functions can be
approximated well.

4.1

Approximating the Utility Function

There are many works on approximating functions using short exponential sums, e.g., the Fourier decomposition approach [62], Prony’s method [54], and many others [9, 10]. However, their approximations are
done over a finite interval domain, say [−π, π] or over a finite number of discrete points. No error bound
can be guaranteed outside the domain. Our algorithm is a generic procedure that turns an algorithm that can
approximate functions over [−π, π] into one that can approximate our utility function µ over [0, +∞), by
utilizing the fact that limx→∞ µ(x) = 0.
Recall for µ ∈ Cbounded , we assume that for any constant  > 0, there exist a constant T such that
µ(x) ≤  for x > T . We also assume there is an algorithm F OURIER
any function f (under some
P that, for
x which is an -approximation
c
φ
conditions specified later), can produce an exponential sum fb(x) = L
i=1 i i
of f (x) in [−π, π] such that |φi | ≤ 1 and L depends only on  and f . In fact, we can assume w.l.o.g. that
F OURIER can approximate f (x) over [−B, B] for any B = O(1). This is because we can apply F OURIER to
the scaled version g(x) = f (x · B
π ) (which is defined on [−π, π]) and then scale the obtained approximation
π
gb(x) back to [−B, B], i.e., the final approximation is fb(x) = gb( B
· x). Scaling a function by a constant
B
factor π typically does not affect the smoothness of f in any essential way and we can still apply F OURIER.
Recall that our goal is to produce an exponential sum that is an -approximation for µ(x) in [0, +∞). We
denote this procedure by E XP S UM -A PPROX.

12

Algorithm: E XP S UM -A PPROX(µ)
1. Initially, we slightly change function µ(x) to a new function µ
b(x) as follows: We require µ
b(x) is a
“smooth ” function in [−2T , 2T ] such that µ
b(x) = µ(x) for all x ∈ [0, T ]; µ
b(x) = 0 for |x| ≥ 2T .
We choose µ
b(x) in [−2T , 0] and [T , 2T ] such that µ
b(x) is smooth. We do not specify the exact
smoothness requirements now since they may depend on the choice of F OURIER. Note that there may
be many ways to interpolate µ such that the above conditions are satisfied (see Example 1 below).
The only properties we need are: (1) µ
b is amenable to algorithm F OURIER; (2) |b
µ(x) − µ(x)| ≤  for
x ≥ 0.
2. We apply F OURIER to g(x) = η x µ
b(x) over domain [−hT , hT ] (η ≥P
1 and h ≥ 2 are constants to
x
be determined later). Suppose the resulting exponential sum gb(x) = L
g (x) −
i=1 ci φi , such that |b
g(x)| ≤  for all x ∈ [−hT , hT ].
P
φi x
3. Let µ
e(x) = L
i=1 ci ( η ) , which is our final approximation of µ(x) on [0, ∞).
Example 1 Consider the utility function µ(x) = 1/(x + 1). Let T = 1 − 1. So µ(x) <  for all x > T .
Now we create function µ
b(x) according to the first step of E XP S UM -A PPROX. If we only require µ
b(x) to be
1
continuous, then we can use, for instance, the following piecewise function: µ
b(x) = x+1 for x ∈ [0, T ];
x
2
µ
b(x) = − T +  for x ∈ [T , 2T ]; µ
b(x) = 0 for x > 2T ; µ
b(x) = µ
b(−x) for x < 0. It is easy to see that
µ
b is continuous and -approximates µ.
2
By setting η = 2 and
!
P
|c
|/)
log( L
i
i=1
h ≥ max 2,
,
T

(4)

we can show the following theorem.
Lemma 3 µ
e(x) is a 2-approximation of µ(x).
Proof: We know that |b
g (x) − g(x)| ≤  for x ∈ [0, hT ]. Therefore, we have that
|e
µ(x) − µ
b(x)| =

gb(x) g(x)

− x ≤ x ≤ .
x
η
η
η

Combining with |b
µ(x) − µ(x)| ≤ , we obtain |e
µ(x) − µ(x)| ≤ 2 for x ∈ [0, hT ]. For x > hT , we can
see that
 x
 x
L
L
L
L
X
X
φi
φi
1 X
1 X
ci
≤
ci
≤ x
|ci | ≤ hT
|e
µ(x)| =
|ci | ≤ 
η
η
2
2
i=1

i=1

i=1

i=1

2

Since µ(x) <  for x > hT , the proof is complete.

Remark: Since we do not know ci before applying F OURIER, we need to set h to be a quantity (only
depending
PL on  and T ) such that (4) is always satisfied. In particular, we need to provide an upper bound
for i=1 |ci |. In the next subsection, we use the Fourier series decomposition as the choice for F OURIER,
which allows us to provide such a bound for a large class of functions.
13

4.2

Implementing F OURIER

Now, we discuss the choice of algorithm F OURIER and the conditions that f (x) needs to satisfy so that it
is possible to approximate f (x) by a short exponential sum in a bounded interval. In fact, if we know in
advance that there is a short exponential sum that can approximate f , we can use the algorithms developed
in [10, 11] (for continuous case) and [9] (for the discrete case). However, those works do not provide an
easy characterization of the class of functions. From now on, we restrict ourselves to the classic Fourier
series technique, which has been studied extensively and allows such characterizations.
Suppose from now on that f (x) is a real periodic function defined on [−π, π]. Consider the partial sum
of the Fourier series of the function f (x):
N
X

(SN f )(x) =

ck eikx

k=−N

Rπ
1
−ikx dx. It has L = 2N + 1 terms. Since f (x) is a real
where the Fourier coefficient ck = 2π
−π f (x)e
function, we have ck = c−k and the partial sum is also real. We are interested in the question under which
conditions does the function SN f converge to f (as N increases) and what is convergence rate? Roughly
speaking, the “smoother” f is, the faster SN f converges to f . In the following, we need one classic result
about the convergence of Fourier series and show how to use it in our problem.
We need a few more definitions. We say f satisfies the α-Hölder condition if |f (x)−f (y)| ≤ C |x−y|α ,
for some constant C and α > 0 and any x and y. The constant C is called the Hölder coefficient of f , also
denoted as |f |C 0,α . We say f is C-Lipschitz if f satisfies 1-Hölder condition with coefficient C.
Example 2 It is easy to check that the utility function µ in Example 1 is 1-Lipschitz since | dµ(x)
dx | ≤ 1 for
x ≥ 0. We can also see that χ
e(x) (defined in (1)) is 1δ -Lipschitz.
We need the following classic result of Jackson.
Theorem 6 (See e.g., [56]) Suppose that f (x) is a real periodic function defined on [−π, π]. If f satisfies
the α-Hölder condition, it holds that
|f (x) − (SN f )(x)| ≤ O

 |f |

C 0,α ln N
Nα



.

We are ready to spell the details of F OURIER. Recall g(x) is obtained in step 2 in Algorithm E XP S UM A PPROX. By construction, g(−hT ) = g(hT ) = 0 for h ≥ 2. Hence, it can be considered as a periodic
function with period 2hT . Note that in Jackson’s theorem, the periodic function f is defined on [−π, π]. In
order to apply Jackson’s theorem to g(x) over [−hT , hT ], we consider the following function f , which is
the scaled version of g:
f (x) = g(xhT /π).
Then, F OURIER returns the following function gb, which is a sum of exponential functions:


xπ
gb(x) = SN f
.
hT
Now, we show that |b
g (x) − g(x)| ≤  for all x ∈ [−hT , hT ]. For the later parts of the analysis, we
need a few simple lemmas. The proofs of these lemmas are straightforward and thus omitted here.

14

Lemma 4 Suppose f : [a, c] → R is a continuous function which consists of two pieces f1 : [a, b] → R
and f2 : [b, c] → R. If both f1 and f2 satisfy the α-Hölder condition with Hölder coefficient C, then
|f |C 0,α ≤ 2C.
Lemma 5 Suppose g : [a, c] → R is a continuous function satisfying the α-Hölder condition with Hölder
coefficient C. Then, for f (x) = g(tx) for some t > 0, we have |f |C 0,α ≤ Ctα .
By Lemma 4, we know that the piecewise function µ
b (defined in step 1 in E XP S UM -A PPROX) satisfies
α-Hölder condition with coefficient 2C. Therefore, we can easily see that g(x) = µ
b(x)η x satisfies α1+2T

Hölder condition with coefficient at most 2
C on [−hT , hT ] (This is because µ
b is non-zero only in
[−2T , 2T ]). According to Lemma 5, we have |f (x)|C 0,α = |g(xhT /π)|C 0,α ≤ 21+2T (hT /π)α C. Using
Theorem 6, we obtain the following corollary.
Corollary 1 Suppose µ ∈ Cbounded satisfies the α-Hölder condition with |µ|C 0,α = O(1). For
N = 2O(T ) (h/)1+1/α ,
it holds that |g(x) − gb(x)| ≤  for x ∈ [−hT , hT ].
Proof: Applying Theorem 6 to f and plugging in the given value of N , we can see that |f (x)−(SN f )(x)| ≤
xπ
xπ
 for x ∈ [−π, π]. Hence, we have that |g(x) − gb(x)| = |f ( hT
) − SN f ( hT
)| ≤  for x ∈ [−hT , hT ]. 2


How to Choose h: Now, we discuss the issue left in Section 4.1, that is how to choose h (the value should
be independent of ci s and L) to satisfy (4), when µ satisfies the α-Hölder condition for some α > 1/2. We
need the following results about thePabsolute convergence of Fourier coefficients. If f satisfies the α-Hölder
condition for some α > 1/2, then +∞
i=−∞ |ci | ≤ |f |C 0,α · cα where cα only depends on α [62]. We can see
that in order to ensure (4), it suffices to to set value h such that
hT ≥ log


21+2T (hT /π)α Ccα
= 2T + O log hT / .


We can easily verify that the above condition can be satisfied by letting h = max(O( T1 log 1 ), 2).
Proof of Theorem 2: Everything is in place to prove Theorem 2. First, we bound L by Corollary 1:
L = 2N + 1 = 2O(T ) poly(1/).
Rπ
1
−ikx dx.
Next, we bound the magnitude of each ck . Recall ck is the Fourier coefficient: ck = 2π
−π f (x)e
where f (x) = g(xhT /π) = η xhT /π µ
b(xhT /π) for x ∈ [−π, π]. Since hT = max(O(T ), O(log 1/)),
O(T
)

we can see |f (x)| ≤ 2
poly(1/) for x ∈ [−π, π]. Therefore,
Z π
1
|ck | ≤
|f (x)|dx ≤ 2O(T ) poly(1/).
2π −π
Finally, combining Corollary 1 and Lemma 3, we complete the proof of Theorem 2.

15

2

H

H

H
µ(x)

H

=

ν(x)

H

+

+

ν(x)
x0 x1

T2

x2

x3

x0 x1

x3 x0 x1

x2

x2

x3

x0 x1

x2

x3

Figure 2: (1) The concave utility function µ(x) and ν(x) = H − µ(x) (2) The piecewise linear function
ν(x). (3)-(5) Decomposing ν(x) into three scaled copies of τ (x).

Class Cconcave

5

In this section, we handle the case where the utility function µ : [0, ∞) → [0, ∞) is a concave nondecreasing
function and our goal is to prove Theorem 3.
We use OPT to denote the optimal value of our problem EUM(A). We can assume without loss of
generality that we know OPT, modulo a multiplicative factor of (1 ± ). This can be done by guessing all
powers of (1 + ) between maxe∈U E[µ(we )] and E[µ(w(U ))], 8 and run our algorithm for each guess. For
ease of notation, we assume that our current guess is exact OPT. Let
H = OPT/2 ,

T1 = µ−1 (OPT/)

and

T2 = µ−1 (OPT/2 ).

(5)

We first make the following simplifying assumption and show how to remove it later:
S1. We assume µ(0) = 0 and µ(x) = µ(T2 ) = H for all x > T2 .
Lemma 6 If the utility function µ ∈ Cconcave
the additional assumption S1, then, for any  > 0, we
P satisfies
x , such that |e
c
φ
µ(x) − µ(x)| ≤ O(OPT) for all x > 0,
can obtain an exponential sum µ
e(x) = L
k
k=1
k
where L = poly(1/), |ck | ≤ poly(1/)H and |φk | ≤ 1 for all k = 1, . . . , L.
Proof: Consider the function ν(x) = H − µ(x). We can see ν is a nonincreasing convex function and
ν(x) = 0 for all x > T2 . We first approximate ν by a piecewise linear function ν as follows. Let N = 1/3 .
For all 0 ≤ i ≤ N , let
 


−1 iH
−1 (N − i)H
xi = µ
=ν
and xN +1 = ∞.
N
N
)−ν(xi )
for 0 ≤ i ≤ N . The piecewise linear function ν is defined by ν(xi ) = ν(xi ) for all
Let hi = ν(xxi+1
i+1 −xi
0 ≤ i ≤ N and
ν(x) = ν(xi ) + (x − xi )hi , for x ∈ [xi , xi+1 ].

It is easy to see ν is also a convex function (see Figure 2) and |ν(x) − ν(x)| ≤ H/N ≤ OPT.
Now we show ν can be written as a linear sum of N scaled copies of the following function ρ:
ρ(x) = 1 − x

for 0 ≤ x ≤ 1, and ρ(x) = 0

8

for x > 1.

We can assume every e ∈ U is in at least one feasible solution S ∈ F. Otherwise, we can simply remove those irrelevant
elements. Then, OPT is at least maxe∈U E[µ(we )]. We can test whether an item is an irrelevant element by using the pseudopolynomial time algorithm as follows: we assign the item with weight 1 and other items weight 0. We ask whether there is a feasible
solution with weight exactly 1.

16

We let ρh,a (x) = (−ha)ρ(x/a). It is easy to see that the first piece of ρh,a has slope h and ends at x = a.
Define


x
for 0 ≤ i ≤ N.
ρi (x) = ρhi −hi+1 ,xi+1 (x) = xi+1 (hi+1 − hi )ρ
xi+1
P −1
It is not hard to verify that ν(x) = N
2).
i=0 ρi (x) (see Figure
PD
x
2
By Theorem 2, we can find a function τe(x) =
k=1 dk ψk with D = poly(N/ ) = poly(1/),
9
2
|dk | = poly(1/) and |ψk | ≤ 1 for k = 1, . . . , D, such that |e
τ (x) − ρ(x)| ≤  /N for x ≥ 0. Consider
the function
νe(x) =

N
−1
X


xi+1 (hi+1 − hi )e
τ

i=0

x



xi+1

=

N
−1 X
D
X

1/xi+1 x

xi+1 (hi+1 − hi )dk (ψk

) .

i=0 k=1

Clearly, νe is the summation of N D exponentials. It is not difficult to see the magnitude of each coefficient,
|xi+1 (hi+1 − hi )dk |, is at most −xi+1 hi |dk | ≤ poly(1/)H. We can also see that
|e
ν (x) − ν(x)| ≤ |e
ν (x) − ν(x)| + |e
ν (x) − ν(x)| ≤ O(H2 ) ≤ O(OPT) for x ≥ 0.
2

Finally, letting µ
e(x) = H − νe(x) finishes the proof.

Since kµk∞ = OPT/2 , Lemma 6 implies that µ
e is an 3 -approximation of µ. Then, applying Theorem 1, we can immediately obtain a polynomial time algorithm that runs in time (n/)poly(1/) and finds a
solution S ∈ F such that OPT − E[µ(S)] ≤ 2 kµk∞ ≤ OPT, i.e., a PTAS.
Now, we show how to get rid of the assumption S1. From now on, the utility function µ is a general
increasing concave utility function with µ(0) = 0. 10 Let µH (x) = min(µ(x), H). We can see that
µH (x) satisfies S1. We say a value p ∈ R+ is huge is if p > T2 . Otherwise, we call it normal. We use
Huge to denote the set of huge values. For each element e, let wenm be the random variable which has the
same distribution
we in the normal value region, and zero probability elsewhere. For any S ⊆ U , let
P as nm
nm
w (S) = e∈S we . In the following lemma, we show µH is a good approximation for µ for normal
values.
Lemma 7 For any S ∈ F, we have that
E[µ(wnm (S))] − O()OPT ≤ E[µH (wnm (S))] ≤ E[µ(wnm (S))].
Proof: It is obvious that E[µH (wnm (S))] ≤ E[µ(wnm (S))]. So, we only need to prove the first inequality.
For any S ⊆ F, we have E[µ(w(S))] ≤ OPT. By Markov inequality, Pr[w(S) ≥ T1 ] ≤ , which implies
Pr[wnm (S) ≥ T1 ] ≤ . Now, we claim that for any integer k ≥ 1,
Pr [wnm (S) ≥ (k + 2)T1 ] ≤  Pr [wnm (S) ≥ kT1 ] .

(6)

Consider the following stochastic process. Suppose the weights of the elements in S are realized one by
one (say wenm
, . . . , wenm
). Let Zt be the sum of the first t realized values. Let t1 be the first time such that
n
1
Zt1 +1 ≥ kT1 . If this never happens, let t1 = ∞ and Zt1 = Zn . Let event E1 be t1 ≤ n and E2 be
9

It suffices to let T = 1 (i.e., ρ(x) = 0 for x ≥ 1).
The assumption that µ(0) = 0 is without loss of generality. If µ(0) > 0, we can solve the problem with the new utility function
µ(x) − µ(0). It is easy to verify a PTAS for the new problem is a PTAS for the original problem.
10

17

Zn ≥ (t + 2)T1 . Consider the random value Zn − Zt1 =
and all wenm
are nonnegative, we can see that
t
Pr[Zn − Zt1 > T1 | E1 ] = Pr

Pn

nm
t=t1 +1 wet .

As wnm (S) = Zn =

Pn

nm
t=1 wet

n
i
i
hX
nm
w
>
T
≤ .
wenm
>
T
|
t
≤
n
≤
Pr
1
1
1
et
t

n
h X

t=1

t=t1 +1

Moreover, we can see that event E1 ∧ (Zn − Zt1 > T1 ) is a necessary condition for event E2 . Hence, the
claim holds because
Pr[E2 ] ≤ Pr[E1 ] Pr[Zn − Zt1 > T1 | E1 ] ≤  Pr[E1 ].
From (6), we can see that Pr [wnm (S) ≥ 3T1 ] ≤ 2 , Pr [wnm (S) ≥ 5T1 ] ≤ 3 , . . . so on and so forth.
Furthermore, we can see that
Z ∞
Z ∞
nm
nm
nm
Pr[wnm (S) ≥ µ−1 (x)]dx
Pr[µ(w (S)) ≥ x]dx =
E[µ(w (S))] − E[µH (w (S))] =
H
H
Z
OPT ∞
nm
−1
=
Pr[w (S) ≥ µ (H + kOPT/)]dk

0
Z
OPT ∞
Pr[wnm (S) ≥ T2 + kT1 ]dk
≤

0
∞
2OPT X k
≤
 ≤ O(OPT).

k=2/

The first inequality holds due to the concavity of µ (or equivalently, the convexity of µ−1 ):
µ−1 (H + kOPT/) ≥ T2 + µ−1 (kOPT/) = T2 + kµ−1 (OPT/) = T2 + kT1 for k ≥ 0.

2

Now, we handle the contribution from huge values. Let Huge = {p | p ≥ T2 } and
X
Hg(e) =
Pr[we = p]µ(p).11
p∈Huge

Hg(e) can be thought as the expected contribution of huge values of e. We need the following observation
in [14]: the contribution of the huge values can be essentially linearized and separated from the contribution
of normal values, in the sense of the following lemma. We note that the simple insight has been used in a
variety of contexts in stochastic optimization problems (e.g., [47, 33, 34]).
Lemma 8 (The first half of Theorem 2 in [14]) For any S ∈ F, we have that


X
nm
E[µ(w(S))] ∈ (1 ± O()) E[µ(w (S))] +
Hg(e) .
e∈S
11
If we is continuously distributed, we let Hg(e) =
exactly the same as the discrete case.

R∞
T2

µ(x)pe (x)dx. The algorithm and analysis for the continuous case are

18

Now, we are ready to state our algorithm,P
which is an extension of the algorithm in Section 3. Using
x
Lemma 6, we first obtain a function µ
eH (x) = L
µH (x) − µH (x)| ≤ OPT. The feature
k=1 ck φk such that |e
vector Ft(e) is a 2L + 1-dimensional integer vector
Ft(e) = hα1 (e), β1 (e), . . . , αL (e), βL (e), bnHg(e)/OPTci,
where αi (e), βi (e) are defined as in (3) with respect to wenm . In other words, we extend the original feature vector by one more coordinate which represents the (scaled and rouned) contribution of huge values.
Similarly, each configuration Cf(v) is indexed by such a 2L + 1-dimensional vector v. The last coordinate
of v P
is at most n2 /. As before, we let Cf(v) = 1 if and only if there is a feasible solution S ∈ F such
that e∈S Ft(e) = v. We slightly modify the definition of Val(v) to incorporate the contribution of huge
values, as the following:
Val(hx1 , y1 , . . . , xL , yL , zi)| =

L
X

ck e−xk γ+iyk δ + z ·

k=1

OPT
.
n

Using the same technique as in Lemma 2 and the pseudopolynomail time algorithm for A, we can compute
the values of all configurations in time (n/)poly(1/) . Then, we return the solution for which the corresponding configuration Cf(hx1 , y1 , . . . , xL , yL , zi) that takes value 1 and maximizes |Val(v)|.
P
Proof of Theorem 3: The proof is similar to that of Theorem 1. Let any S ⊆ U , let vS = e∈S Ft(e)
Using the same proof of Lemma 1 and the fact that |e
µH (x) − µH (x)| ≤ OPT, we can see that for any
S ∈ F,
X
|Val(vS )| = E[µ(wnm (S))] +
Hg(e) ± O(OPT).
e∈S

Combining with Lemma 7 and Lemma 8, we can further see that for any S ∈ F,
|Val(vS )| = (1 ± O())E[µ(w(S))] ± O(OPT).
Suppose S is our solution and S ∗ is the optimal solution for utility function µ. From our algorithm, we
know that |Val(vS )| ≥ |Val(vS ∗ )|, which implies E[µ(w(S))] ≥ (1 − O())OPT and completes the proof.
2

6

Class Cincreasing

d
Recall that µ(x) ∈ Cincreasing is a positive, differentiable and increasing function and dx
µ(x) ∈ [L, U]
for some constants L, U > 0 and all x ≥ 0. By scaling, we can assume without loss of generality that
L ≤ 1 ≤ U. Our algorithm is almost the same as the one in Section 5 except that we use a slightly different
set of parameters:


OPT L
OPT L
−1 OPT L
H = 2 · , T1 =
· , T2 = µ
·
, and Huge = {p | p ≥ T2 }

U

U
2
U

Let µH (x) = min(µ(x), H). So, µH satisfies assumption S1. However, we can not use Lemma 6 since
it requires concavity. Nevertheless, we can still approximate µH by a short exponential sum, as in Lemma 9.
The remaining algorithm is exactly the same as the one in Section 5. To prove the performance guarantee,
we only need to prove analogues of Lemma 7 and Lemma 8. Now, we prove the aforementioned lemmas.
19

P
x
µH (x) −
Lemma 9 For any  > 0, we can obtain an exponential sum µ
eH (x) = L
k=1 ck φk , such that |e
µH (x)| ≤ O(OPT) for all x > 0, where L = poly(1/), |ck | ≤ poly(1/)H and |φk | ≤ 1 for all k.
d
Proof: Since dx
µ(x) ∈ [L, U], we can see that H/U ≤ T2 ≤ H/L. Consider ν(x) = H − µH (x). We
can see ν is a decreasing, differentiable function and ν(x) = 0 for all x > T2 . Consider the function
ν(x) = H1 ν(xH). First, let T = 1/L ≥ T2 /H and we can see ν(x) = 0 for x > T . Hence, ν ∈ Cbounded
P
x
and satisfies U-Lipschitz condition. By Theorem 2, we can compute a function νe(x) = L
k=1 dk ψk , which
2
is an  -approximation of ν, with L = poly(1/), |ψk | ≤ 1 and dk = poly(1/) for all k. Therefore,
P
1/H x
) is the desired approximation.
2
µ
eH (x) = H − H νe(x/H) = H − L
k=1 (Hdk )(ψk

The following lemma is an analogue of Lemma 7.
Lemma 10 For any S ∈ F, we have that
E[µH (wnm (S))] ∈ (1 ± O())E[µ(wnm (S))].
Proof: The proof is almost the same as that of Lemma 7, except that the last line makes use of the bounded
derivative assumption (instead of the concavity):

kOPT 
kOPT
µ−1 H +
≥ T2 +
≥ T2 + kT1 for k > 0.

U
We handle the contribution from huge values in the same way. Recall Hg(e) =
p]µ(p). The following lemma is an analogue of Lemma 8.

2
P

p∈Huge Pr[we

=

Lemma 11 For any S ∈ F, we have that


X
E[µ(w(S))] ∈ (1 ± O()) E[µ(wnm (S))] +
Hg(e) .
e∈S

Proof: We can use exactly
P the same proof of Theorem 2 in [14] to show that E[µ(w(S))] ≥ (1 −
nm
O())(E[µ(w (S))] + e∈S Hg(e)), as the proof holds even without the concavity assumption. The other
direction requires a different argument, which goes as follows. Let E0 be the event that no we is realized to
a huge value and Ee,p be the event that we is realized to value p ∈ Huge. By Markov inequality, we have
Pr[E0 ] ≥ 1 − L/U. Moreover, using the fact that e−x ≥ 1 − x, we have that

 Y


X
X
exp −
Pr[Ee,p ] =
exp −
Pr[Ee,p ]
e∈S

e∈S,p∈Huge

Y

≥

e∈S

p∈Huge

1−

X


Pr[Ee,p ]

p∈Huge

= Pr[E0 ] ≥ 1 − L/U ≥ 1 − .
P
Hence, e∈S,p∈Huge Pr[Ee,p ] ≤ −(ln(1 − )) ≤ 2 for  < 1/2.
Next, we can see that E[µ(w(S)) | E0 ] Pr[E0 ] ≤ E[µ(wnm (S))] (for each realization of {we }e∈S
satisfying E0 , there is a corresponding realization of {wenm }e∈S ). From the bounded derivative assumption,

20

we can also see that E[µ(w(S)) | Ee,p ] ≤ µ(p) + U · E[w(S)] By inclusion-exclusion, we have that
X
E[µ(w(S))] ≤ E[µ(w(S)) | E0 ] Pr[E0 ] +
E[µ(w(S)) | Ee,p ] Pr[Ee,p ]
e∈S,p∈Huge

≤ E[µ(wnm (S))] +

X

Pr[Ee,p ]µ(p) +

e∈S,p∈Huge

≤ E[µ(wnm (S))] +

X

X

Pr[Ee,p ] · U · E[w(S)]

e∈S,p∈Huge

Hg(e) + O()E[µ(w(S))].

e∈S

The last inequality holds since E[w(S)] ≤ E[µ(w(S))]/L = O(E[µ(w(S))]).

7

2

Applications

We first consider two utility functions χ(x) and χ
e(x) presented in the introduction. Note that maximizing
E[χ(w(S))] is equivalent to maximizing Pr(w(S) ≤ 1). The following lemma is straightforward.
Lemma 12 For any solution S,
Pr(w(S) ≤ 1) ≤ E[e
χ(w(S))] ≤ Pr(w(S) ≤ 1 + δ).
Corollary 2 Suppose there is a pseudopolynomial time algorithm for the exact version of A. Then, for any
fixed constants  > 0 and δ > 0, there is an algorithm that runs in time ( n )poly(1/) , and produces a solution
S ∈ F such that
Pr(w(S) ≤ 1 + δ) +  ≥ max
Pr(w(S 0 ) ≤ 1)
0
S ∈F

Proof: By Theorem 1, Theorem 2 and Lemma 12, we can easily obtain the corollary. Note that we can
choose T = 2 for any δ ∈ (0, 1) and  > 0. Thus L = poly(1/).
2
Now, let us see some applications of our general results to specific problems.
Stochastic Shortest Path: Finding a path with the exact target length (we allow non-simple paths)12 can be
easily done in pseudopolynomial time by dynamic programming.
Stochastic Spanning Tree: We are given a graph G, where the weight of each edge e is an independent,
nonnegative random variable. Our objective is to find a spanning tree T in G, such that Pr(w(T ) ≤ 1)
is maximized. Polynomial time algorithms have been developed for Gaussian distributed edges [35, 26].
To the best of our knowledge, no approximation algorithm with provable guarantee is known for other
distributions. Noticing there exists a pseudopolynomial time algorithm for the exact spanning tree problem
[5], we can directly apply Corollary 2.
Stochastic k-Median
on Trees: The problem asks for a set S of k nodes in the given probabilistic tree G
P
such that Pr( v∈V (G) dis(v, S) ≤ 1) is maximized, where dis(v, S) is the minimum distance from v to any
node in S in the tree metric. The k-median problem can be solved optimally in polynomial time on trees by
dynamic programming [37]. It is straightforward to modify the dynamic program to get a pseudopolynomial
time algorithm for the exact version.
12

The exact version of simple path is NP-hard, since it includes the Hamiltonian path problem as a special case.

21

Stochastic Knapsack with Random Sizes: We are given a set U of n items. Each item i has a random size
wi and a deterministic profit vi . We are also given a positive constant
P 0 ≤ γ ≤ 1. The goal is to find a subset
S ⊆ U such that Pr(w(S) ≤ 1) ≥ γ and the total profit v(S) = i∈S vi is maximized.
If the profits of the items are polynomially bounded integers, we can see the optimal profit is also
a polynomially bounded integer. We can first guess the optimal profit. For each guess g, we solve the
following problem: find a subset S of items such that the total profit of S is exactly g and E[e
χ(w(S))] is
maximized. The exact version of the deterministic problem is to find a solution S with a given total size
and a given total profit, which can be easily solved in pseudopolynomial time by dynamic programming.
Therefore, by Corollary 2, we can easily show that we can find in polynomial time a set S of items such that
the total profit v(S) is at least the optimum and Pr(w(S) ≤ 1 + ) ≥ (1 − )γ for any constant  and γ.
If the profits are general integers, we can use the standard scaling technique to get a (1−)-approximation
for the total profit. We first make a guess of the optimal profit, rounded down to the nearest power of (1 + ).
maxi vi
There are at most log1+ nmin
guesses. For each guess g, we solve the following problem. We discard
i vi
all items with a profit larger than g. Let ∆ = ng2 . For each item with a profit smaller than g
n , we set its new
vi
profit to be v̄i = 0. Then, we scale each of the rest profits vi to v̄i = ∆b ∆ c. Now, we define the feasible set
F(g) = {S |

X

(1 − 2)g ≤

i∈S

X

v̄i ≤ (1 + 2)g}.

i∈S

2

Since there are at most n distinct v̄ values, we can easily show that finding a solution S in F(g) with a
given total size can be solved in pseudopolynomial time by dynamic programming.
Denote the optimal solution by S ∗ and the optimal profit by OP T . Suppose g is the right guess, i.e.,
1
( 1+ )OP T ≤ g ≤ OP T . We can easily see that for any solution S, we have that
(1 −

X
X
1 X
)
vi − g ≤
v̄i ≤
vi
n
i∈S

i∈S

i∈S

where the first inequalities are due to vi ≥ g
n and we set at most g profit to zero. Therefore, we can see
∗
S ∗ ∈ F(g). Applying Corollary 2, we obtain a solution
P Pr(w(S) ≤ 1 + δ) +  ≥ Pr(w(S ) ≤
P S such that
1 + δ). Moreover, the profit of this solution v(S) = i∈S vi ≥ i∈S v̄i ≥ (1 − 2)g ≥ (1 − O())OP T.
In sum, we have obtained the following result.
Theorem 7 For any constants  > 0 and γ > 0, there is a polynomial time algorithm to compute a set S of
items such that the total profit v(S) is within a 1− factor of the optimum and Pr(w(S) ≤ 1+) ≥ (1−)γ.
poly(1/)

Bhalgat et al. [13, Theorem 8.1] obtained the same result, with a running time n2
, while our running
time is npoly(1/) .
Moreover, we can easily extend our algorithm to generalizations of the knapsack problem if the corresponding exact version has a pseudopolynomial time algorithm. For example, we can get the same result for
the partial-ordered knapsack problem with tree constraints [25, 58]. In this problem, items must be chosen
in accordance with specified precedence constraints and these precedence constraints form a partial order
and the underlining undirected graph is a tree (or forest). A pseudopolynomial algorithm for this problem is
presented in [58].
Stochastic Knapsack with Random Profits: We are given a set U of n items. Each item i has a deterministic size wi and a random profit vi . The goal is to find a subset of items that can be packed into a
knapsack with capacity 1 and the probability that the profit is at least a given threshold T is maximized.
22

Henig [32] and Carraway et al. [15] studied this problem for normally distributed profits and presented
dynamic programming and branch and bound heuristics to solve this problem optimally.
We can solve the equivalent problem of minimizing the probability that the profit is at most the given
threshold, subject to the capacity constraint. We first show that relaxing the capacity constraint is necessary.
Consider the following deterministic knapsack instance. The profit of each item is the same as its size. The
given threshold is 1. We can see that the optimal probability is 1 if and only if there is a subset of items of
total size exactly 1. Otherwise, the optimal probability is 0. However, determining whether these is a subset
of items with total size exactly 1 is NP-Complete. Therefore, it is NP-hard to approximate the original
problem within any additive error less than 1 without violating the capacity constraint.
The corresponding exact version of the deterministic problem is to find a set of items S such that w(S) ≤
1 and v(S) is equal to a given target value. In fact, there is no pseudopolynomial time algorithm for this
problem. Since otherwise we can get an  additive approximation without violating the capacity constraint,
contradicting the lower bound argument. Note that a pseudopolynomial time algorithm here should run in
time polynomial in the profit value (not the size). However, if the sizes can be encoded in O(log n) bits (we
only have a polynomial number of different sizes), we can solve the problem in time polynomial in n and
the largest profit value by standard dynamic programming.
For general sizes, we can round the size of each item down to the nearest multiple of nδ . Then, we can
solve the exact version in pseudopolynomial time poly(maxi vi , n, 1/δ) by dynamic programming. It is
easy to show that for any subset of items, its total size is at most the total rounded size plus δ. Therefore, the
total size of our solution is at most 1 + δ. We summarize the above discussion in the following theorem.
Theorem 8 If the optimal probability is Ω(1), we can find in time (n/δ)poly(1/) a subset S of items such
that Pr(v(S) > (1 − )T ) ≥ (1 − )OPT and w(S) ≤ 1 + δ, for any constant  > 0.

8

Extensions

In this section, we discuss some extensions to our basic approximation scheme. We first consider optimizing
a constant number of utility functions in Section 8.1. Then, we study the problem where the weight of each
element is a random vector in Section 8.2.

8.1

Multiple Utility Functions

The problem we study in this section contains a set U of n elements. Each element e has a random weight
we . We are also given d utility functions µ1 , . . . , µd and d positive numbers λ1 , . . . , λd . We assume d is a
constant. A feasible solution consists of d subsets of elements that satisfy some property. Our objective is
to find a feasible solution S1 , . . . , Sd such that E[µi (w(Si ))] ≥ λi for all 1 ≤ i ≤ d.
We can easily extend our basic approximation scheme to the multiple utility functions case as follows.
We decompose these utility functions into short exponential sums using E XP S UM -A PPROX as before. Then,
for each utility function, we maintain (n/)O(L) configurations. Therefore, we have (n/)O(dL) configurations in total and we would like to compute the values for these configurations. We denote the deterministic
version of the problem under consideration by A. The exact version of A asks for a feasible solution
S1 , . . . , Sd such that the total weight of Si is exactly the given number ti for all i. Following an argument
similar to Lemma 2, we can easily get the following generalization of Theorem 1.
Theorem 9 Assume that there is a pseudopolynomial algorithm for the exact version of A. Further assume
that given any  > 0, we can -approximate each utility function by an exponential sum with at most L
23

terms. Then, there is an algorithm that runs in time (n/)O(dL) and finds a feasible solution S1 , . . . , Sd such
that E[µi (w(Si )] ≥ λi −  for 1 ≤ i ≤ d, if there is a feasible solution for the original problem.
Now let us consider two simple applications of the above theorem.
Stochastic Multiple Knapsack: In this problem we are given a set U of n items, d knapsacks with capacity
1, and d constants 0 ≤ γi ≤ 1. We assume d is a constant. Each item i has a random size wi and a
deterministic profit vi . Our
Pd objective is to find d disjoint subsets S1 , . . . , Sd such that Pr(w(Si ) ≤ 1) ≥ γi
for all 1 ≤ i ≤ d and i=1 v(Si ) is maximized. The exact version of the problem is to find a packing
such that the load of each knapsack i is exactly the given value ti . It is not hard to show this problem can
be solved in pseudopolynomial time by standard dynamic programming. If the profits are general integers,
we also need the scaling technique as in stochastic knapsack with random sizes. In sum, we can get the
following generalization of Theorem 7.
Theorem 10 For any constants d ∈ N,  > 0 and 0 ≤ γi ≤ 1 for 1 ≤ i ≤ P
d, there is a polynomial time
algorithm to compute d disjoint subsets S1 , . . . , Sd such that the total profit di=1 v(Si ) is within a 1 − 
factor of the optimum and Pr(w(Si ) ≤ 1 + ) ≥ (1 − )γi for 1 ≤ i ≤ d.
Stochastic Multidimensional Knapsack: In this problem we are given a set U of n items and a constant
0 ≤ γ ≤ 1. Each item i has a deterministic profit vi and a random size which is a random d-dimensional
vector wVi = {w
Pi1 , . . . , wid }. We assume d is a constant. Our objective is to find a subset S of items such
d
that Pr( j=1 ( i∈S wij ≤ 1)) ≥ γ and total profit is maximized. This problem can be also thought as the
fixed set version of the stochastic packing problem considered in [21, 13]. We first assume the components
of each size vector are independent. The correlated case will be addressed in the next subsection.
For ease of presentation, we assume d = 2 from now on. Extension to general constant d is straightforward. We can solve the problem by casting it into a multiple utility problem as follows. For each item
i, we create two copies i1 and i2 . The copy ij has a random weight wij . A feasible solution consists of
two sets S1 and S2 such that S1 (S2 ) only contains the first (second) copies of the elements and S1 and
S2 correspond to exactly the same subset of original elements. We enumerate all such pairs (γ1 , γ2 ) such
that γ1 γ2 ≥ γ and γi ∈ [γ, 1] is a power of 1 −  for i = 1, 2. Clearly, there are a polynomial number
of such pairs.
For each pair (γ1 , γ2 ), we solve the following problem: find a feasible solution S1 , S2 such
P
that Pr( i∈Sj wij ≤ 1) ≥ γj for all j = 1, 2 and total profit is maximized. Using the scaling technique
and Theorem 9Vfor optimizing
multiple utility
P we can get a (1 − )-approximation for the optimal
P
Q2 functions,
2
profit and Pr( j=1 ( i∈Sj wij ≤ 1)) = j=1 Pr( i∈Sj wij ≤ 1) ≥ (1 − O())γ1 γ2 ≥ (1 − O())γ.
We note that the same result for independent components can be also obtained by using the discretization
technique developed for the adaptive version of the problem in [13] 13 . If the components of each size vector
are correlated, we can not decompose the problem into two 1-dimensional utilities as in the independent case.
Now, we introduce a new technique to handle the correlated case.

8.2

Multidimensional Weight

The general problem we study contains a set U of n elements. Each element e has a random weight vector
wi = (wi1 , . . . , wid ). We assume d is a constant. We are also given a utility function µ : Rd → R+ . A
feasible solution is a subset of elements satisfying some property. We use w(S) as a shorthand notation
13

With some changes of the discretization technique, the correlated case can be also handled [12].

24

P
P
for vector ( i∈S wi1 , . . . , i∈S wid ). Our objective is to find a feasible solution S such that E[µ(w(S)] is
maximized.
From now on, x and k denote d-dimensional vectors and kx (or k·x) denotes the inner product of k and x.
As before, we assume µ(x) ∈ [0, 1] for all x ≥ 0 and lim|x|→+∞ µ(x) = 0, where |x| = max(x1 , . . . , xd ),
Our algorithm is almost the same as in the one dimension case and we briefly sketch itQhere. We first notice
that expected utilities decompose for exponential utility functions, i.e., E[φk·w(S) ] = Pi∈S E[φk·wi ]. Then,
we attempt to -approximate the utility function µ(x) by a short exponential sum |k|≤N ck φkx
k (there
P
d
k·w(S)
k·w(S)
are O(N ) terms). If this can be done, E[φ
] can be approximated by |k|≤N ck E[φ
]. Using
the same argument as in Theorem 1, we can show that there is a polynomial time algorithm that can find a
feasible solution S with E[µ(w(S))] ≥ OP T − for any  > 0, provided that a pseudopolynomial algorithm
exists for the exact version of the deterministic problem.
To approximate the utility function µ(x), we need the multidimensional Fourier
series expansion of
P
d
ikx where c =
a function
f : C → C (assuming f is 2π-periodic in each axis): f (x) ∼
k
k∈Zd ck e
R
1
−ikx
f
(x)e
dx.
The
rectangular
partial
sum
is
defined
to
be
d
d
x∈[−π,π]
(2π)
X

SN f (x) =

X

...

|k1 |≤N

ck eikx .

|kd |≤N

It is known that the rectangular partial sum SN f (x) converges uniformly to f (x) in [−π, π]d for many
function classes as n tends to infinity. In fact, a generalization of Theorem 6 to [−π, π]d also holds [3]: If f
satisfies the α-Hölder condition, then
|f (x) − (SN f )(x)| ≤ O

 |f |

C 0,α

lnd N 

Nα

for x ∈ [−π, π]d .

Now, we have an algorithm F OURIER that can approximate a function in a bounded domain. It is also
straightforward to extend E XP S UM -A PPROX to the multidimensional case. Hence, we can -approximate µ
by a short exponential sum in [0, +∞)d , thereby proving the multidimensional generalization of Theorem 2.
Let us consider an application of our result.
Stochastic Multidimensional Knapsack (Revisited): We consider the case where the components of each
weight vector can be correlated. Note that the utility function χ2 corresponding to this problem is the two
dimensional threshold function: χ2 (x, y) = 1 if x ≤ 1 and y ≤ 1; χ2 (x, y) = 0 otherwise. As in the one
dimensional case, we need to consider a continuous version χ
e2 of χ2 (see Figure 1(3)). By the result in this
section and a generalization of Lemma 12 to higher dimension, we can get the following.
Theorem 11 For any constants d ∈ N,  > 0 and 0 ≤ γ ≤ 1, there is a polynomial V
time algorithm
for
P
finding a set S of items such that the total profit v(S) is 1− factor of the optimum and Pr( dj=1 ( i∈S wij ≤
1 + )) ≥ (1 − )γ.

9

A Few Remarks

Convergence of Fourier series: The convergence of the Fourier series of a function is a classic topic in harmonic analysis. Whether the Fourier series converges to the given function and the rate of the convergence
typically depends on a variety of smoothness condition of the function. We refer the readers to [62] for a
more comprehensive treatment of this topic. We note that we could obtain a smoother version of χ (e.g., see

25

Figure 1(2)), instead of the piecewise linear χ
e, and then use Theorem 6 to obtain a better bound for L. This
would result in an even better running time. Our choice is simply for the ease of presentation.
Discontinuous utility functions: If the utility function µ is discontinuous, e.g., the threshold function, then
the partial Fourier series behaves poorly around the discontinuity (this is known as the Gibbs phenomenon).
However, informally speaking, as the number of Fourier terms increases, the poorly-behaved strip around
the edge becomes narrower. Therefore, if the majority of the probability mass of our solution lies outside
the strip, we can still guarantee a good approximation of the expected utility. There are also techniques to
reduce the effects of the Gibbs phenomenon (See e.g., [28]). However, the techniques are not sufficient to
handle discontinuous functions. We note that very recently, Daskalakis et al. [20] obtained a true additive
PTAS (instead of a bi-criterion additive PTAS) for a closely related problem, called the fault tolerant storage
problem, under certain technical assumptions. 14 It is not clear how to use their technique to obtain a true
additive PTAS for our expected utility maximization problem. We leave this problem as an interesting open
problem.

10

Conclusion

We study the problem of maximizing expected utility for several stochastic combinatorial problems, such as
shortest path, spanning tree and knapsack, and several classes of utility functions. A key ingredient in our
algorithm is to decompose the utility function into a short exponential sum, using the Fourier series decomposition. Our general approximation framework may be useful for other stochastic optimization problems.
We leave the problems of obtaining a true additive PTAS, or nontrivial multiplicative approximation factors
for Cbounded as interesting open problems.

11

Acknowledgments

We would like to thank Evdokia Nikolova for providing an extended version of [52] and many helpful
discussions. We also would like to thank Chandra Chekuri for pointing to us the work [13] and Anand
Bhalgat for some clarifications of the same work.

References
[1] H. Ackermann, A. Newman, H. Röglin, and B. Vöcking. Decision making based on approximate and
smoothed pareto curves. Algorithms and Computation, pages 675–684, 2005.
[2] S. Agrawal, A. Saberi, and Y. Ye. Stochastic Combinatorial Optimization under Probabilistic Constraints. Arxiv preprint arXiv:0809.0460, 2008.
[3] S. Alimov, R. Ashurov, and A. Pulatov. Multiple fourier series and fourier integrals, in commutative
harmonic analysis. IV: Harmonic analysis in Rn . Encyclopedia of Mathematical Science, 42, 1992.
14

In the fault tolerant storage problem, we are given n real numbers
and an addition number 0 < θ < 1.
P 0 < p1 ≤ . . . ≤ pn < 1,P
n
Our goal is to partition 1 into n positive values x1 , . . . , xn (i.e., n
i=1 xi = 1), such that Pr[
i=1 Xi ≥ θ] is maximized, where
Xi is the Bernoulli random variable which takes value xi with probability pi . In order to obtained an additive PTAS, Daskalakis et
al. [20] assumed that all pi s are bounded below by a constant.

26

[4] N. Bansal, A. Gupta, J. Li, J. Mestre, V. Nagarajan, and A. Rudra. When LP is the Cure for Your
Matching Woes: Improved Bounds for Stochastic Matchings. European Symposium on Algorithms,
pages 218–229, 2010.
[5] F. Barahona and W. Pulleyblank. Exact arborescences, matchings and cycles. Discrete Applied Mathematics, 16(2):91–99, 1987.
[6] J. Bard and J. Bennett. Arc reduction and path preference in stochastic acyclic networks. Management
Science, 37(2):198–215, 1991.
[7] D. Bernoulli. Exposition of a new theory on the measurement of risk. Econometrica: Journal of
the Econometric Society, pages 23–36, 1954. Originally published in 1738; translated by Dr. Louise
Sommer.
[8] D. Bernoulli. Exposition of a new theory on the measurement of risk. Econometrica, 22(1):22–36,
1954. Originally published in 1738; translated by Dr. Lousie Sommer.
[9] G. Beylkin and L. Monzón. On Generalized Gaussian Quadratures for Exponentials and Their Applications* 1. Applied and Computational Harmonic Analysis, 12(3):332–373, 2002.
[10] G. Beylkin and L. Monzón. On approximation of functions by exponential sums. Applied and Computational Harmonic Analysis, 19(1):17–48, 2005.
[11] G. Beylkin and L. Monzón. Approximation by exponential sums revisited. Applied and Computational
Harmonic Analysis, 28(2):131–149, 2010.
[12] A. Bhalgat, 2011. Personal Communication.
[13] A. Bhalgat, A. Goel, and S. Khanna. Improved approximation results for stochastic knapsack problems. In ACM-SIAM symposium on Discrete algorithms, 2011.
[14] A. Bhalgat and S. Khanna. A utility equivalence theorem for concave functions. In Integer Programming and Combinatorial Optimization, pages 126–137. Springer, 2014.
[15] R. Carraway, R. Schmidt, and L. Weatherford. An algorithm for maximizing target achievement in the
stochastic knapsack problem with normal returns. Naval research logistics, 40(2):161–173, 1993.
[16] C. Chekuri and S. Khanna. A PTAS for the multiple knapsack problem. In ACM-SIAM symposium on
Discrete algorithms, pages 213–222, 2000.
[17] N. Chen, N. Immorlica, A. Karlin, M. Mahdian, and A. Rudra. Approximating matches made in
heaven. International Colloquium on Automata, Languages and Programming, pages 266–278, 2009.
[18] W. Cheney and W. Light. A Course in Approximation Theory. Brook/Cole Publishing Company, 2000.
[19] R. Cheng, J. Chen, and X. Xie. Cleaning uncertain data with quality guarantees. Proceedings of the
VLDB Endowment, 1(1):722–735, 2008.
[20] C. Daskalakis, A. De, I. Diakonikolas, A. Moitra, and R. A. Servedio. A polynomial-time approximation scheme for fault-tolerant distributed storage. In SODA, pages 628–644. SIAM, 2014.

27

[21] B. Dean, M. Goemans, and J. Vondrák. Adaptivity and approximation for stochastic packing problems.
In ACM-SIAM symposium on Discrete algorithms, pages 395–404, 2005.
[22] B. Dean, M. Goemans, and J. Vondrak. Approximating the Stochastic Knapsack Problem: The Benefit
of Adaptivity. Mathematics of Operations Research, 33(4):945, 2008.
[23] M. Fazel and M. Chiang. Network utility maximization with nonconcave utilities using sum-of-squares
method. In Decision and Control, 2005 and 2005 European Control Conference. CDC-ECC’05. 44th
IEEE Conference on, pages 1867–1874. IEEE, 2005.
[24] P. Fishburn. Utility Theory and Decision Making. John Wiley & Sons, Inc, 1970.
[25] M. Garey and D. Johnson. “Computers and Intractability: A Guide to the Theory of NPCompleteness”. W.H. Freeman, 1979.
[26] S. Geetha and K. Nair. On stochastic spanning tree problem. Networks, 23(8):675–679, 1993.
[27] A. Goel and P. Indyk. Stochastic load balancing and related problems. In Annual Symposium on
Foundations of Computer Science, page 579, 1999.
[28] D. Gottlieb and C. Shu. On the Gibbs phenomenon and its resolution. SIAM review, 39(4):644–668,
1997.
[29] V. Goyal and R. Ravi. Chance constrained knapsack problem with random item sizes. To appear in
Operation Research Letter, 2009.
[30] S. Guha and K. Munagala. Adaptive Uncertainty Resolution in Bayesian Combinatorial Optimization
Problems. To appear in ACM Transactions on Algorithms, 2008.
[31] A. Gupta, M. Pál, R. Ravi, and A. Sinha. Boosted sampling: approximation algorithms for stochastic
optimization. In ACM Symposium on Theory of Computing, pages 417–426. ACM, 2004.
[32] M. Henig. Risk criteria in a stochastic knapsack problem. Operations Research, 38(5):820–825, 1990.
[33] L. Huang and J. Li. Approximating the expected values for combinatorial optimization problems over
stochastic points. In Automata, Languages, and Programming, pages 910–921. Springer, 2015.
[34] L. Huang, J. Li, J. M. Phillips, and H. Wang. -kernel coresets for stochastic points. arXiv preprint
arXiv:1411.0194, 2014.
[35] H. Ishii, S. Shiode, and T. Nishida Yoshikazu. Stochastic spanning tree problem. Discrete Applied
Mathematics, 3(4):263–273, 1981.
[36] D. Kahneman and A. Tversky. Prospect theory: An analysis of decision under risk. Econometrica:
Journal of the Econometric Society, pages 263–291, 1979.
[37] O. Kariv and S. Hakimi. An algorithmic approach to network location problems. II: The p-medians.
SIAM Journal on Applied Mathematics, 37(3):539–560, 1979.
[38] J. Kleinberg, Y. Rabani, and É. Tardos. Allocating bandwidth for bursty connections. In ACM Symposium on Theory of Computing, page 673, 1997.

28

[39] J. Li and A. Deshpande. Consensus answers for queries over probabilistic databases. In ACM
SIGMOD-SIGACT-SIGART Symposium on Principles of Database Systems, 2009.
[40] J. Li and A. Deshpande. Ranking continuous probabilistic datasets. Proceedings of the VLDB Endowment, 3(1), 2010.
[41] J. Li, B. Saha, and A. Deshpande. A unified approach to ranking in probabilistic databases. In Proceedings of the VLDB Endowment, 2009.
[42] J. Li and T. Shi. A fully polynomial-time approximation scheme for approximating a sum of random
variables. Operations Research Letters, 42(3):197–202, 2014.
[43] J. Li and W. Yuan. Stochastic combinatorial optimization via poisson approximation. In ACM Symposium on Theory of Computing, 2013.
[44] R. Loui. Optimal paths in graphs with stochastic or multidimensional weights. Communications of the
ACM, 26(9):670–676, 1983.
[45] R. Martin. The St. Petersburg Paradox. The Stanford Encyclopedia of Philosophy, 2004. http:
//plato.stanford.edu/archives/fall2004/entries/paradox-stpetersburg.
[46] S. Mittal and A. Schulz. A general framework for designing approximation schemes for combinatorial
optimization problems with many objectives combined into one. Approximation, Randomization and
Combinatorial Optimization. Algorithms and Techniques, pages 179–192, 2008.
[47] A. Munteanu, C. Sohler, and D. Feldman. Smallest enclosing ball for probabilistic data. In Proceedings
of the thirtieth annual symposium on Computational geometry, page 214. ACM, 2014.
[48] I. Murthy and S. Sarkar. Exact algorithms for the stochastic shortest path problem with a decreasing
deadline utility function. European Journal of Operational Research, 103(1):209–229, 1997.
[49] I. Murthy and S. Sarkar. Stochastic shortest path problems with piecewise-linear concave utility functions. Management Science, 44(11):125–136, 1998.
[50] E. Nikolova. Approximation Algorithms for Reliable Stochastic Combinatorial Optimization. International Workshop on Approximation Algorithms for Combinatorial Optimization Problems, pages
338–351, 2010.
[51] E. Nikolova, M. Brand, and D. Karger. Optimal route planning under uncertainty. In Proceedings of
International Conference on Automated Planning and Scheduling, 2006.
[52] E. Nikolova, J. Kelner, M. Brand, and M. Mitzenmacher. Stochastic shortest paths via quasi-convex
maximization. In European Symposium on Algorithms, pages 552–563, 2006.
[53] F. Oberhettinger. Fourier transforms of distributions and their inverses: a collection of tables. Academic press, 1973.
[54] M. R. Osborne and G. K. Smyth. A modified prony algorithm for fitting sums of exponential functions.
SIAM Journal of Scientific Computing, 1995.
[55] C. Papadimitriou and M. Yannakakis. On the approximability of trade-offs and optimal access of web
sources. In Annual Symposium on Foundations of Computer Science, 2000.
29

[56] M. J. D. Powell. Approximation theory and methods. Cambridge University Press, 1981.
[57] A. Ralston and R. Rabinowitz. A First Course in Numerical Analysis. 2001.
[58] H. Safer, J. B. Orlin, and M. Dror. Fully polynomial approximation in multi-criteria combinatorial
optimization, 2004. MIT Working Paper.
[59] P. A. Samuelson. St. petersburg paradoxes: Defanged, dissected, and historically described. Journal
of Economic Literature, 15(1):24–55, 1977.
[60] D. Shmoys and C. Swamy. An approximation scheme for stochastic linear programming and its application to stochastic integer programs. J. ACM, 53(6):1012, 2006.
[61] C. Sigal, A. Pritsker, and J. Solberg. The stochastic shortest route problem. Operations Research,
28(5):1122–1129, 1980.
[62] E. Stein and R. Shakarchi. Fourier analysis: an introduction. Princeton University Press, 2003.
[63] C. Swamy. Risk-Averse Stochastic Optimization: Probabilistically-Constrained Models and Algorithms for Black-Box Distributions. ACM-SIAM symposium on Discrete algorithms, 2010.
[64] C. Swamy and D. Shmoys. Approximation algorithms for 2-stage stochastic optimization problems.
ACM SIGACT News, 37(1):46, 2006.
[65] J. von Neumann and O. Morgenstern. Theory of Games and Economic Behavior. Princeton Univ.
Press, 2nd edition, 1947.

A

Computing E[φwe ]

If X is a random variable, then the characteristic function of X is defined as
G(z) = E[eizX ].
We can see E[φwe ] is nothing but the value of the characteristic function of we evaluated at −i ln φ (here ln
is the complex logarithm function). For many important distributions, including negative binomial, Poisson,
exponential, Gaussian, Chi-square and Gamma, a closed-form characteristic function is known. See [53] for
a more comprehensive list.
Example 3 Consider the Poisson distributed we with mean λ, i.e., Pr(we = k) = λk e−λ /k! . Its characiz
teristic function is known to be G(z) = eλ(e −1) . Therefore,
E[φwe ] = G(−i ln φ) = eλ(φ−1) .
1

Example 4 For Gaussian distribution N (µ, σ 2 ), we know its characteristic function is G(z) = eizµ− 2 σ
Therefore,
1 2
E[φwe ] = G(−i ln φ) = φu+ 2 σ ln φ .

30

2 z2

.

For some continuous distributions, no closed-form characteristic function is known and we need proper
numerical approximation method.
If the support of the distribution is bounded, we can use for example Gauss-Legendre quadrature [57]. If
the support is infinite, we can truncate the distribution and approximate
finite
R b the integral over the remaining
Pk
interval; Generally speaking a quadrature method approximates a f (x)dx by a linear sum i=1 ci f (xi )
where ci and xi are some constants independent of the function f . A typical practice is to use composite rule,
that is to partition [a, b] into N subintervals and approximate the integral using some quadrature formula over
each subinterval. For the example of Gauss-Laguerre quadrature, assuming continuity of the 2kth derivative
of f (x) for some constant k, if we partition [a, b] into M subintervals and apply Gauss-Legendre quadrature
of degree k to each subinterval, the approximation error is
Error =

(b − a)2k+1
(k!)4
f (2k) (ξ)
(2k + 1)[(2k)!]3
M 2k

where ξ is some point in (a, b) [57, pp.116]. Let ∆ = b−a
M . If we treat k as a constant, the behavior of the
error (in terms of ∆) is Error(∆) = O(∆2k maxξ f (2k) (ξ)). Therefore, if the support and maxξ f (2k) (ξ)
are bounded by a polynomial, we can approximate the integral, in polynomial time, such that the error is
O(1/nβ ) for any fixed integer β.
The next lemma shows that we do not lose too much even though we can only get an approximation of
w
E[φ e ].
e
Lemma 13 Suppose in Theorem 5, we can only compute an approximate value of E[φw
i ], denoted by Ee,i ,
we
for each e and i, such that |E[φi ] − Ee,i | ≤ O(n−β ) for some positive integer β. Denote E(S) =
Q
PL
k=1 ck
e∈S Ee,i . For any solution S, we have that

|E[e
µ(w(S))] − E(S)| ≤ O(n1−β ).
Proof: We need the following simple result (see [40] for a proof): a1 , . . . , an and e1 , . . . , en are complex
numbers such that |ai | ≤ 1 and |ei | ≤ n−β for all i and some β > 1. Then, we have
n
n
Y
Y
(ai + ei ) −
Ei ≤ O(n1−β ).
i=1

i=1

Since |φi | ≤ 1, we can see that
e
|E[φw
i ]|

Z

φxi pe (x)dx| ≤ 1.

=|
x≥0

The lemma simply follows by applying the above result and noticing that L and all ck s are constants.

2

We can show that Theorem 1 still holds even though we only have the approximations of the E[φwe ]
values. The proof is straightforward and omitted.

31

