Proceedings of the 2015 Winter Simulation Conference
L. Yilmaz, W. K. V. Chan, I. Moon, T. M. K. Roeder, C. Macal, and M. D. Rossetti, eds.

MIRROR DESCENT STOCHASTIC APPROXIMATION FOR COMPUTING WORST-CASE
STOCHASTIC INPUT MODELS
Soumyadip Ghosh
Math Sciences Department
T.J. Watson IBM Research Center
Yorktown Heights, NY 10598, USA

Henry Lam
Department of Industrial and Operations Engineering
University of Michigan
1205 Beal Ave.
Ann Arbor, MI 48109, USA

ABSTRACT
Performance analysis via stochastic simulation is often subject to input model uncertainty, meaning that the
input model is unknown and needs to be inferred from data. Motivated especially from situations with limited
data, we consider a worst-case analysis to handle input uncertainty by representing the partially available
input information as constraints and solving a worst-case optimization problem to obtain a conservative
bound for the output. In the context of i.i.d. input processes, such approach involves simulation-based
nonlinear optimizations with decision variables being probability distributions. We explore the use ofa specialized class of mirror descent stochastic approximation (MDSA) known as the entropic descentalgorithm, particularly effective for handling probability simplex constraints, to iteratively solve for the local
optima. We show how the mathematical program associated with each iteration of the MDSA algorithm canbe efficiently computed, and carry out numerical experiments to illustrate the performance of the algorithm.
1

INTRODUCTION

We consider stochastic simulation with uncertain input specifications. This arises in many contexts when
only finite data is available for calibrating the input model. Conventional approaches to handle such
uncertainty include goodness-of-fit tests (e.g. Banks et al. 2000), Bayesian model selection or averaging(e.g. Chick 2001), and bootstrapping (e.g. Barton et al. 2013). In this paper, we consider an alternative,
worst-case approach, which can be particularly relevant when data is limited or unavailable, or when themodeler wants to assess risk beyond what data could tell. This approach comprises solving optimization
programs on the space of input models, with constraints representing the partially known information onthese inputs. The optimal values of these programs give the worst- or the best-case performance measures
under the uncertainty.
More concretely, consider a stochastic simulation problem of computing Ep [h(X)], where X =(X1 , . . . , XT ) are the i.i.d. replications from a common input model represented by the probability disi.i.d.

tribution p. In other words, Xt ∼ p. h(·) is some function that captures the system logic applied to X,
and Ep [·] denotes the expectation with respect to the i.i.d. replications of Xt ’s each with distribution p. As
an example, one can think of X as a sequence of interarrival or service times, and p is the corresponding
interarrival or service time distribution. h(X) can be the indicator function of long waiting time for customer
T , so that Ep [h(X)] is the corresponding long waiting time probability.
Our premise is that p is unknown, but some information is available. These information can be
represented by constraints on p, denoted A . Our goal is to perform a worst-case analysis to calculate
min Ep [h(X)] or max Ep [h(X)].

p∈A

p∈A

(1)

Ghosh and Lam
For instance, A can specify that the first two moments of p are within certain ranges. In this case, we
can write
A = {p : µ ≤ Ep X1 ≤ µ, σ ≤ Ep X12 ≤ σ }for some numbers µ ≤ µ and σ ≤ σ .
As another example, suppose that the modeler adopts a baseline model b, e.g. exponential distribution,and believes that it is close to the truth. Then we can use the constraint
A = {p : d(p, b) ≤ η}
where d(·, ·) is some statistical distance defined between two distributions, i.e. a notion of distance thatis not restricted to any specific parametric models, and η > 0 is a threshold parameter that describes thelevel of uncertainty.
The key is that in both examples, the modeler represents the information on the input model nonparametrically, without committing to a specific parametric class that can be potentially misspecified. Later inthis paper, we shall use the above two examples of A to illustrate our optimization procedure.
The worst-case formulation like (1) is inspired from the literature of distributionally robust optimization(DRO). This literature considers decision-making under ambiguity of the underlying probabilistic models(e.g. Delage and Ye 2010, Goh and Sim 2010, Ben-Tal et al. 2013), often cast in the form of minimax
problems. Glasserman and Xu (2014) further studies the solutions to these problems via simulation.
However, unlike the standard DRO literature that possess convex formulations and linear (or stage-wiselinear) objectives, the objective function in (1) is not linear, or even convex in p in general, since computing
the expectation Ep [·] actually involves a T -fold convolution of p due to the i.i.d. structure. The objectivefunction in (1) is also accessible only via simulation in our context of interest. These challenges obstruct
the direct usage of previous techniques. Moreover, because of non-convexity, we would set our goal as toobtain a local optimum for (1).
Formulation (1) also faces two additional challenges. The first is that p can be high-dimensional. In
fact, for a continuous input distribution, p is infinite-dimensional. In this paper we shall restrict ourselvesto discrete finite-support distributions, with the understanding that discretization is used a priori to handle
continuous distributions. Second, (1) consists of constraints that can be computationally costly to various
extents.
Ghosh and Lam (2015) discusses some of the above challenges and proposes a constrained stochastic
approximation (SA) algorithm based on the Frank-Wolfe (FW) method in nonlinear programming (which
they call FWSA), by successively solving linear-objective subprograms that give the best feasible direction
at each iteration. They show almost sure convergence and analyze the convergence rate of the algorithm.
This paper describes an alternative to FWSA by using the idea of mirror descent (MD) (e.g. Nemirovski
and Yudin 1983, Nemirovski et al. 2009). In fact, we will focus on a specialized class of mirror descentstochastic approximation (MDSA) that is based on a so-called entropic distance-generating function. The
deterministic version of this algorithm is known as the entropic descent algorithm (EDA) (e.g. Beck andTeboulle 2003). The motivation of this alternative is two-fold. The first is that MD method, and inparticular EDA, is known to scale favorably with the decision space’s dimension when the space is a
probability simplex. Since we consider potentially discrete distributions that have many support points (orvery fine discretization of continuous distributions), this property of MD is attractive for our problem. Whiletheoretical guarantees of this sort have been shown only for convex programs, our preliminary experimental
results show the algorithm works well for some potentially non-convex settings. Second, unlike FWSA,MDSA does not require a growing sample size for gradient estimation. This can potentially save huge
computational effort as iteration increases.
The rest of the paper is organized as follows. Section 2 introduces MDSA. Section 3 discusses gradientestimation. Section 4 analyzes the solution to the mathematical program required in each iteration of
MDSA. Finally, Section 5 shows some numerical results for an M/GI/1 example.

Ghosh and Lam
2

MIRROR DESCENT STOCHASTIC APPROXIMATIONThe MD procedure was first proposed by Nemirovski and Yudin (1983) for convex deterministic optimization.
One of its earlier motivations is to define a proper gradient descent type algorithm for decision variable thatlies on function spaces. In general, the gradient of the objective lies in the dual of the decision space, which
is different from the primal, and so the usual iteration routine by moving the current solution along thegradient direction is not well-defined (the standard Euclidean space is an exception because it is isometricto its dual space). MD resolves this issue by first mapping the solution into the dual space, followed byapplying the gradient descent in the dual space and at the end mapping back into the primal. This method
can as well be used for decision variables in the Euclidean space: By appropriately defining the space andits norm that are pertained to the geometry of the constraints, it can achieve a rate of convergence that isfavorably scaled with the dimension of the problem.
We first give a general description of MD, focusing on a generic minimization problem minx∈X Z(x).
Adopting Nemirovski et al. (2009), we first define a distance-generating function ω(·) and a norm k · k on
the space X . The function ω(·) needs to be strongly convex with respect to the norm k · k, meaning that
ω(y) − ω(x) ≥ ∇ω(x)0 (y − x) +αky − xk2
2

for some constant α > 0 and any x, y ∈ X , which is equivalent to(y − x)0 (∇ω(y) − ∇ω(x)) ≥ αky − xk2 .
Next we define the prox-function V (·, ·) as
V (x, z) = ω(z) − ω(x) − ∇ω(x)0 (z − x).

(2)

In other words, V (·, ·) is the remainder of the first order approximation of ω(·). We then define theprox-mapping proxx (y) as
proxx (y) = argminz∈X y0 (z − x) +V (x, z).
For an optimization minx∈X Z(x), the MD method iteratively updates xk+1 from xk by using xk+1 =proxxk (γk ∇Z(xk )) where γk is some step size. If ∇Z(xk ) can only be assessed through unbiased samples,c k ) is used, in which case it becomes MDSA.
then an estimate ∇Z(xWe put the above discussion into our context. Our objective function is Z(p) = Ep [h(X)]. We consider
p in the probability simplex on n fixed support points {u1 , . . . , un }, denoted by P = {(p1 , . . . , pn ) : ∑ni=1 pi =1, pi ≥ 0 for all i = 1, . . . , n}. The feasible region A ⊂ P is our space in consideration. We shall use theentropic distance-generating functionnω(p) = − ∑ pi log pi .

(3)

i=1

The corresponding prox-function isn

V (p, q) = ∑ qi log
i=1

qi
pi

(4)

where p = (p1 , . . . , pn ), q = (q1 , . . . , qn ) ∈ A , which can be seen by plugging (3) into (2). The MD iteration
using the prox-mapping is
pk+1 = proxpk (γk ∇Z(pk )) = argminp∈A γk ∇Z(pk )0 (p − pk ) +V (pk , p).
This procedure is known as the entropic descent algorithm (EDA).

Ghosh and Lam
Since ∇Z(pk ) is not known exactly in the stochastic simulation context, one needs to estimate it bysampling. In the next sections, we shall discuss two issues with implementing a stochastic version of EDA:
The first is how to obtain an estimate of the gradient, and the second is on efficient methods for solving
the subprogram that defines the prox-mapping.
In the convex optimization setting where unbiased samples for ∇Z is available,√it is known that theconvergence rate of the optimality gap, i.e. E[Z(xk ) − Z(x∗ )] for MDSA is of order 1/ k. Compared to theprojective gradient method, i.e. carrying out an Euclidean projection on the feasible region at each iteration,2
which is in fact a special case of MDSA
by taking ω(x)p= (1/2)kxk2 where k · k2 denotes L2 -norm, the√convergence constant in front of 1/ k for EDA can be n/ log n times smaller (Nemirovski et al. 2009).
This makes EDA favorable for large n, and motivates our use of EDA in this paper.
Another motivation for using MDSA is that, in contrast to FWSA proposed in Ghosh and Lam (2015),
we do not increase the sample size with the iteration number. To explain the difference in the required
sample sizes on an intuitive level, we note that FWSA consists of solving a subprogram at each iteration
with a linearized objective in order to obtain the best feasible direction. This subprogram introduces bias
in estimating the best feasible direction, even when the gradient itself can be estimated unbiasedly. The
increase in sample size with respect to iteration is therefore designed to offset this bias. On the other
hand, prox-mapping possesses contraction properties similar to a projection, which facilitates the analysis
for MDSA using a notion of distance derived from the prox-function rather than by looking at the bias
separately.
3

GRADIENT ESTIMATION

Note that since p ∈ P is a discrete distribution, we can writeZ(p) = Ep [h(X)] = ∑ · · · ∑ h(xi1 , . . . , xiT )pi1 · · · piT
i1

(5)

iT

so that Z(p) is a high-dimensional polynomial in the variables p. Although explicit form (5) is available,
there are often numerous summands in (5) that renders exact computation infeasible, and so simulation
has to be used. In carrying out gradient estimation for (5), one issue is that naive differentiation of (5)with respect to each coordinate of p does not lead to any implementable procedure, because it results in ahigh-dimensional multivariate polynomial that has no probabilistic meaning.
To get around the above issue, Ghosh and Lam (2015) proposes using the Gateaux derivative on afunctional of probability distribution for simulating the gradient of (5), and we shall use it here. Given ap, consider a mixture of p with the point mass at each coordinate, i.e. ei = (0, . . . , 1, . . . , 0)0 where the 1 isat the i-th coordinate. This mixture is written as (1 − ε)p + εei for i = 1, . . . , n. Then consider
ψi (p) =dZ((1 − ε)p + εei )
dε

(6)
ε=0

and define ψ (p) = (ψ1 (p), . . . , ψn (p))0 . This vector ψ (p) captures all the necessary gradient informationof Z(p) within the space P, in the following sense:
Theorem 1 (Adapted from Ghosh and Lam (2015)) The vector ψ (p) defined in (6) satisfies∇Z(p)0 (q − p) = ψ (p)0 (q − p)(7)

for any p, q ∈ P, and can be evaluated by writing as
ψi (p) = Ep [h(X)si (X)]
where

(8)

T

I(xt = ui )
− T.
pi
t=1

si (x1 , . . . , xT ) = ∑(9)Ghosh and Lam
Here I(·) is the indicator function, and (u1 , . . . , un ) are the support points for the probability simplex P.
Since the operations in MDSA only involve ∇Z(p)0 (q − p) but not the gradient itself directly, (7)dictates that one can use ψ (p) to replace ∇Z(p). The functions si in (9) can be regarded as a nonparametricversion of the score function (Ghosh and Lam 2015), and with (8) one can use h(X)si (X) as an unbiased
estimator for the gradient information vector ψ (p). In summary, given pk , our MDSA computes at iterationkψ (pk ))
pk+1 = proxpk (γk ψ̂whereψ (pk ) =ψ̂1 M
∑ (h(X)si (X)) j
M j=1

and (h(X)si (X)) j , j = 1, . . . , M denote the samples using M independent sample paths of X.
4

COMPUTING PROX-MAPPING

We focus on two examples in this section.
4.1 Kullback-Leibler (KL) Divergence Constraint
Let

(
A =

n

pi
p ∈ P : ∑ pi log ≤ η
bi
i=1

)
(10)

where b = (b1 , . . . , bn ) is some baseline distribution. This type of constraints has been used commonly inrobust control (e.g. Hansen and Sargent 2008, Petersen et al. 2000) and mathematical finance (e.g. Glasserman
and Xu 2014, Glasserman and Xu 2013), among other disciplines. The prox-mapping proxp (ξξ ) forp = (p1 , . . . , pn ) ∈ P and ξ = (ξ1 , . . . , ξn ) ∈ Rn is described as follows:
Proposition 1 For the feasible region defined in (10), the prox-mapping for the minimization problem in(1) using prox-function defined in (4) is given by proxp (ξξ ) = ((proxp (ξξ ))1 , . . . , (proxp (ξξ ))n ) as below:
First definen−ξi log pi − ξnpe∑i=1 iibi
ϕ(p, ξ ) =−log∑ pi e−ξi∑ni=1 pi e−ξi
i=1
and

1
1+β

κ(β ; p, ξ ) =∑ni=1 pi


β
ξ 
− i
bi1+β e 1+β log bpii − ξi
1

β

−

ξi

(1 + β ) ∑ni=1 pi1+β bi1+β e 1+β

n

1

β

−

ξi

− log ∑ pi1+β bi1+β e 1+β .
i=1If ϕ(p, ξ ) ≤ η, then(proxp (ξξ ))i =pi e−ξi∑nj=1 p j e−ξ jelse if ϕ(p, ξ ) > η, then1

(proxp (ξξ ))i =

β

−

ξi

pi1+β bi1+β e 1+β1
1+β∑nj=1 p j

β

−

ξj

1+β
b 1+β
j e

where β is the positive root of the equation
κ(β ; p, ξ ) = η.

(11)

Ghosh and Lam
Proposition 1 reveals that the prox-mapping under KL constraint is either an exponential twisting with
1β

exponent −ξi , or an exponential twisting together with a “geometric mixture” of p and b given by pi1+β bi1+β .
Note that the main computation in the prox-mapping is merely a one-dimensional root-finding problem forβ.
Proof of Proposition 1.

Consider0

proxp (ξξ ) = argminq∈A ξ (q − p) +V (p, q)
where V (p, q) is defined in (4), which can be rewritten as the solution ofξ 0 (q − p) + ∑ni=1 qi log qpii
minsubject to ∑ni=1 qi log bqii ≤ η∑ni=1 qi = 1
qi ≥ 0 for all i = 1, . . . , n.

(12)Relaxing the first constraint in (12), we have the Lagrangiannqiξ (q − p) + ∑ qi log + β
pi
i=1
0

n

qi
∑ qi log bi − η
i=1

!
(13)

where β is the Lagrange multiplier. By Theorem 1, P.220 in Luenberger (1968), since the program (12)is convex, if we can find β ≥ 0 and q ∈ A such that q minimizes
 (13) over P for the fixed β ≥ 0, andqinthat β and q satisfy the complementary slackness condition β ∑i=1 bi log bi − η = 0, then this q is an
optimal solution for (12).
Suppose that β = 0. Then (13) reduces ton0ξ (q − p) + ∑ qi log
i=1

qi
.
pi

To minimize (14) over q ∈ P, we consider a further relaxation of the constraint ∑ni=1 qi = 1 to get
!
n
n
qi
0
ξ (q − p) + ∑ qi log + λ ∑ qi − 1 .
pi
i=1
i=1

(14)

(15)

Differentiate (15) with respect to each qi and set it to zero to get
ξi + log qi + 1 − log pi + λ = 0or

qi ∝ pi e−ξi .

The constraint ∑ni=1 qi = 1 implies thatqi =pi e−ξi.
∑nj=1 p j e−ξ j(16)

It is easy to check that q defined by (16) satisfies the KKT condition as the minimizer of the convex
minimization program of (14) over q ∈ P.
Now, we go back to the original problem (12). From the discussion above, q defined by (16) is an
optimal solution for (12) if it is feasible in A . In other words, we have to check that ∑ni=1 qi log bqii ≤ η.

Ghosh and Lam
Substituting (16) into ∑ni=1 qi log qbii gives ϕ(p, ξ ), and so the condition becomes ϕ(p, ξ ) ≤ η. This gives
rise to the first case in Proposition 1.
Suppose that ϕ(p, ξ ) > η. Then the q considered in (16) is not feasible, and we must choose β > 0.
Relaxing the constraint ∑ni=1 qi = 1 on the Lagrangian (13), we have
!
!
n
n
n
q
q
i
i
0
ξ (q − p) + ∑ qi log + β ∑ qi log − η + λ ∑ qi − 1 .
(17)
p
b
i
i
i=1
i=1
i=1
Differentiate with respect to each qi and set it to zero gives
ξi + log qi + 1 − log pi + β (log qi + 1 − log bi ) + λ = 0
or

1

β

−

ξi

−

ξi

qi ∝ pi1+β bi1+β e 1+β .
The constraint ∑ni=1 qi = 1 gives1β

pi1+β bi1+β e 1+β

qi =

1

β

ξj

.

(18)

1+β
1+β −
∑nj=1 p j b j e 1+β

q defined by (18) satisfies the KKT condition as the minimizer of (17) over q ∈ P. Going back to theoriginal problem (12), this q is an optimal solution if β is chosen such that ∑ni=1 qi log qbii = η. Substituting(18) into ∑ni=1 qi log bqii gives κ(β ; p, ξ ) and so ∑ni=1 qi log qbii = η is equivalent to (11). Note that we canalways find a positive root for (11) in this case because when β = 0,κ(0; p, ξ ) = ϕ(p, ξ ) > η
whileκ(β ; p, ξ ) → 0
as β → ∞. Hence there must be a positive root for β since κ(β ; p, ξ ) is continuous in β .
4.2 Moment Constraint
Let

(
A =

)

n

p ∈ P : ∑ pi rl (xi ) ≤ µl , l = 1, . . . , m

(19)

i=1

where rl (·)’s are some basis functions. For instance, r1 (x) = x and r2 (x) = x2 correspond to the first two
moments. The prox-mapping proxp (ξξ ) is described as follows:
Proposition 2 For the feasible region defined in (19), the prox-mapping for the minimization problem in(1) using the prox-function defined in (4) is given by proxp (ξξ ) = ((proxp (ξξ ))1 , . . . , (proxp (ξξ ))n ), wherempi e−ξi −∑l=1 βl rl (xi )(proxp (ξξ ))i = nm∑ j=1 p j e−ξ j −∑l=1 βl rl (x j )and {βl , l = 1, . . . , m} is the minimizer ofnminβl ≥0,l=1,...,mmmlog ∑ pi e−ξi −∑l=1 βl rl (xi ) + ∑ βl µl .
i=1

l=1

(20)

Ghosh and Lam
The prox-mapping applies an exponential twisting on the probability distribution, with exponent−ξi − ∑ml=1 βl rl (xi ). Each application of prox-mapping requires solving a convex optimization with only
non-negativity constraints in a dimension equal to the number of moment constraints. If the inequalities are
replaced by equalities, the procedure is exactly the same except that the corresponding Lagrange multipliers
βl in (20) become unconstrained.
Proof of Proposition 2.

The prox-mapping proxp (ξξ ) is found by solving0minξ (q − p) + ∑ni=1 qi log qpii
subject to ∑ni=1 qi rl (xi ) ≤ µl , l = 1, . . . , m∑ni=1 qi = 1
qi ≥ 0 for all i = 1, . . . , n.

(21)

Relaxing the first constraint, (21) has the same optimal value as the Lagrangian formulationn0

min ξ (q − p) + ∑ qi log

max

βl ≥0,l=1,...,m q∈P

i=1

m
qi
+ ∑ βl (rl (xi )qi − µl )pi l=1(22)

since (21) is convex. Moreover, the saddle point solution of (22) is also optimal for (21). Consider the
inner minimization in (22). Relaxing further the constraint ∑ni=1 qi = 1, we have
!
m
n
n
qi
0
ξ (q − p) + ∑ qi log + ∑ βl (rl (xi )qi − µl ) + λ ∑ qi − 1 .
pi l=1i=1
i=1Differentiate with respect to each qi and set it to zero, we havem

ξi + log qi + 1 − log pi + ∑ βl rl (xi ) + λ = 0
l=1

or

m

qi ∝ pi e−ξi −∑l=1 βl rl (xi ) .

The constraint ∑ni=1 qi = 1 implies thatmpi e−ξi −∑l=1 βl rl (xi )qi = n.
m∑ j=1 p j e−ξ j −∑l=1 βl rl (x j )(23)

q defined by (23) satisfies the KKT condition as the minimizer in the inner optimization in (22). Substituting(23) into the inner objective function in (22) gives
n

m

m

0

− log ∑ pi e−ξi −∑l=1 βl rl (xi ) − ∑ βl µl − ξ pi=1

so that (22) becomes

l=1

n

max

βl ≥0,l=1,...,m

m

m

0

− log ∑ pi e−ξi −∑l=1 βl rl (xi ) − ∑ βl µl − ξ pi=1

whose optimal solution is the same as that of (20).

l=1Ghosh and Lam
5

NUMERICAL EXPERIMENTS

This section describes numerical experiments conducted on the canonical simulation model of the M/GI/1
single-server queue to characterize the performance of the MDSA algorithm. This queueing model serves
as a good testing ground because a good approximation to the solution of the optimization problem (1)where the set A is of form (10) is available via a deterministic re-formulation, which lets us observe therate of convergence in the objective value for our algorithm directly.
Consider an M/GI/1 queue where the arrival process is Poisson known with high accuracy to haverate λ = 1. On the other hand, the time Xt taken to provide service to the t−th customer is uncertain but
is assumed i.i.d.. A simulation model estimates Z(p) = Ep [h(XT )], the long-run average of the waitingtimes, i.e. the time between the customer’s arrival and for service to commence, of the first T customers:#"
1 T
where Wt = max{0,Wt−1 + Xt − At }.
Wt ,Z(p) = E
T∑
1
Note that the t-th customer’s waiting time Wt is calculated by iterating over the celebrated Lindley’srecursion, and At ∼ exp(1) is the inter-arrival time of the (t − 1)-th customer.
We carry out optimization (1) that captures the worst-case performance for Z(p) amongst service time
distributions p that is η distance away from b in terms of KL divergence (focusing on the minimization
formulation):
 
pi
≤η
(24)
min Z(p) s.t. ∑ pi logp∈Pbi
iwhere P denotes the probability simplex on {u1 , . . . , un }. Here we choose ui ’s as the uniform discretizationof the interval [0, 1], i.e. ui = (i + 1)/n, and the nominal distribution b is the corresponding discretization ofa mixture distribution 0.3 × Beta(2, 6) + 0.7 × Beta(6, 2) over the interval [0, 1], i.e. bi is the probability
mass of this continuous distribution assigned to the interval [ui−1 , ui ].
As T grows, the average waiting time converges to the corresponding steady-state value, which is
given in closed form for M/GI/1 queues when the traffic intensity ρp = Ep [Xt ] is less than 1, by the
Pollaczek-Khinchine formula (Khintchine 1932) as:
Z∞ (p) =ρp Ep [X1 ] +Varp (X1 ).
2(1 − ρp )So, when T is large, an approximation Z∞∗ to the worst-case performance estimate can be obtained byreplacing Z(p) in the problem (24) with Z∞ (p). (In experiments, a choice of T = 2000 seems to show close
agreement.) With Ep [X1 ] = ∑ pi xi and Ep [X12 ] = ∑ pi xi2 , an approximation for the minimization problem (24)is given by the steady-state formulation (SS) (shown below) which is in the form of a so-called convexconcave-fractional program. As Ghosh and Lam (2015) has noted (also see Boyd and Vandenberghe 2009),(SS) can be reformulated as an equivalent convex optimization problem (SS0 ) using the variable substitutions
t = 1/(2(1 − ∑ pi xi )) and yi = pit. An approximation for the worst-case maximization problem can besimilarly derived.
min
p

s.t.

∑ pi xi2(SS)2(1 − ∑ pi xi ) pi
∑ pi log bi ≤ ηi∑ pi = 1
0 ≤ pi ≤ 1,

∀i = 1, . . . , nminp(SS0 )∑ yi xi2
yis.t. ∑ yi log≤ ηttbi
2t − 2 ∑ yi xi = 1∑ yi = t
0 ≤ yi ≤ t∀i = 1, . . . , nGhosh and Lam
The key parameters in the MDSA algorithm are the fixed sample-size M and the step-size sequence
γk . We take the latter to be of the form γk = k−α . We shall experiment on the convergence behavior of thealgorithm for varied M and α. Note that Nemirovski et al. (2009) study the case α = 0.5 in particular.
All the experimental results were terminated for the MDSA algorithm at iteration k when at least oneof the following criteria are met:
1.
2.
3.
4.

The relative difference between objective value Z(pk ) and the average of the observed values in 30−6previous iterations, (∑30v=1 Z(pk−v ))/30, is below 5 × 10 , or
The gradient estimate ψ̂(pk ) has l2 -norm smaller than 1 × 10−3 , orThe change in iterates pk and pk+1 has l1 -norm smaller than 1 × 10−6 , or
The iteration count k reaches 200.

(a) α = 1.5, M varied as shown(b) M = 76, 800, α varied as shown

Figure 1: Approximated Optimality gap of the MDSA algorithm iterates as function of iteration count(both in log-scale). On the left, the step-size exponent α is held constant at 1.5 while the sample size M isvaried. Conversely, the plot on the right holds the sample-size M at 76, 800 while the step-size exponent
α is varied.
Fig. 1 captures the performance of the MDSA algorithm as a function of the sample size M and step-size
exponent α parameters. It plots the (approximate) optimality gap as a function of the iteration count of
the algorithm, where the parameters are set to be η = 0.025, n = 100 and T = 2000. As can be expected,
Fig. 1a shows that the performance of the MDSA algorithm generally improves as the the sample size M
used to evaluate the gradient ψ(p) is increased. Note that Fig. 1a is plotted in log-scale, and the relatively
constant slopes among different M’s show that M seems to have primarily an effect on the constant factor
that accompanies the first order rate of convergence of the MDSA method as a function of the iteration
count. Note that for very large M values, the numerical accuracy issues intervene, resulting in the breaks
observed in the data series.
Fig. 1b plots the performance as the step-size exponent α is varied. Compared to M, α affects a bit
more the slopes of the trends (although still not very significantly), which demonstrates some effects on
the rate of convergence of MDSA in the power order of the iteration count. The result shows that the best
performance is achieved when α ≈ 1.5. When α is too aggressive (say α ≥ 2.0) the step-sizes become too
small which limits the movement of the solution iterates. On the other hand, too small an α (say α ≤ 1.25)

Ghosh and Lam
may lead to overshooting the best next iterate; this can be inferred from the noisier sample paths for small
α. (The noise observed in the sample paths for α = 1.5 or 1.75 at the lower-right end of the plot is more
likely due to numerical accuracy issues in deriving the approximation from (SS’).)(a) α = 1.5, M varied as shown(b) M = 76, 800, α varied as shown

Figure 2: Approximated Optimality gap of the rolling simple average of the previous thirty MDSA algorithm
iterates as function of iteration count (both in log-scale). On the left, the step-size exponent α is held
constant at 1.5 while the sample size M is varied. Conversely, the plot on the right holds the sample-size
M at 76, 800 while the step-size exponent α is varied.
Next, we also experiment some averaging on the MDSA scheme. This sort of averaging has been
considered in Nemirovski et al. (2009), where the weights are defined as a specific function of the step-size,
and a uniform version has also been widely used in the classical SA literature (e.g. Polyak 1990 and Ruppert
1991). Conceptually, averaging helps tamp down the effect of the noise induced by the fixed sample size at
every iteration. Fig. 2 provides the performance of a simple rolling average of the previous thirty iterates.
We see that the trends noticed in Fig. 1 emerge more clearly in Fig. 2, with the noise of the trends being
smaller across the board.
ACKNOWLEDGMENTS
The second author is grateful to the support by the National Science Foundation under Grants CMMI-1400391
and CMMI-1436247.
REFERENCES
Banks, J., J. Carson, B. Nelson, and D. Nicol. 2000. Discrete-Event System Simulation. Prentice hall
Englewood Cliffs, NJ, USA.
Barton, R. R., B. L. Nelson, and W. Xie. 2013. “Quantifying Input Uncertainty via Simulation ConfidenceIntervals”. INFORMS Journal on Computing 26 (1): 74–87.
Beck, A., and M. Teboulle. 2003. “Mirror Descent and Nonlinear Projected Subgradient Methods forConvex Optimization”. Operations Research Letters 31 (3): 167–175.
Ben-Tal, A., D. Den Hertog, A. De Waegenaere, B. Melenberg, and G. Rennen. 2013. “Robust Solutions of
Optimization Problems Affected by Uncertain Probabilities”. Management Science 59 (2): 341–357.

Ghosh and Lam
Boyd, S., and L. Vandenberghe. 2009. Convex Optimization. Cambridge university press.
Chick, S. E. 2001. “Input Distribution Selection for Simulation Experiments: Accounting for Input Uncertainty”. Operations Research 49 (5): 744–758.
Delage, E., and Y. Ye. 2010. “Distributionally Robust Optimization under Moment Uncertainty with
Application to Data-Driven Problems”. Operations Research 58 (3): 595–612.
Ghosh, S., and H. Lam. 2015. “Computing Worst-Case Input Models in Stochastic Simulation”. submitted,
available at: http://www-personal.umich.edu/ khlam/files/robustSAOR3.pdf .
Glasserman, P., and X. Xu. 2013. “Robust Portfolio Control with Stochastic Factor Dynamics”. Operations
Research 61 (4): 874–893.
Glasserman, P., and X. Xu. 2014. “Robust Risk Measurement and Model Risk”. Quantitative Finance 14(1): 29–58.
Goh, J., and M. Sim. 2010. “Distributionally Robust Optimization and its Tractable Approximations”.
Operations Research 58 (4-Part-1): 902–917.
Hansen, L. P., and T. J. Sargent. 2008. Robustness. Princeton university press.
Khintchine, A. Y. 1932. “Mathematical Theory of a Stationary Queue”. Matematicheskii Sbornik 39:7384.
Luenberger, D. G. 1968. Optimization by Vector Space Methods. John Wiley & Sons.
Nemirovski, A., A. Juditsky, G. Lan, and A. Shapiro. 2009. “Robust Stochastic Approximation Approachto Stochastic Programming”. SIAM Journal on Optimization 19 (4): 1574–1609.
Nemirovski, A., and D. Yudin. 1983. Problem Complexity and Method Efficiency in Optimization. Wiley,
New York.
Petersen, I. R., M. R. James, and P. Dupuis. 2000. “Minimax Optimal Control of Stochastic Uncertain
Systems with Relative Entropy Constraints”. IEEE Transactions on Automatic Control 45 (3): 398–412.
Polyak, B. T. 1990. “New Stochastic Approximation Type Procedures”. Automat. Remote Control 7:98107.
Ruppert, D. 1991. “Stochastic Approximation”. In Handbook in Sequential Analysis, edited by B. K. Ghoshand P. K. Sen. Marcel Dekker.
AUTHOR BIOGRAPHIESSOUMYADIP GHOSH is a Research Staff Member in the Business Analytics and Mathematical Sciences
Department at the IBM T.J. Watson Research Center. His current research interests lie in simulation based
optimization techniques for stochastic optimization problems, with a focus on applications in Energy andPower systems and supply chain management. His email is ghoshs@us.ibm.com and his web page is athttps://researcher.ibm.com/researcher/view.php?person=us-ghoshs.
HENRY LAM is an Assistant Professor in the Department of Industrial and Operations Engineering at
the University of Michigan, Ann Arbor. He graduated from Harvard University with a Ph.D. degree instatistics in 2011, and has been an Assistant Professor in the Department of Mathematics and Statistics at
Boston University until 2014. His research focuses on stochastic simulation, risk analysis, and simulation
optimization. His email address is khlam@umich.edu.