Proceedings of the 2020 Winter Simulation Conference
K.-H. Bae, B. Feng, S. Kim, S. Lazarova-Molnar, Z. Zheng, T. Roeder, and R. Thiesing, eds.

DISTRIBUTIONALLY CONSTRAINED STOCHASTIC GRADIENT ESTIMATION USINGNOISY FUNCTION EVALUATIONSHenry Lam

Junhui Zhang

Department of Industrial Engineering &Operations Research
Columbia University
500 W. 120th Street
New York, NY 10027 USA

Department of Applied Physics &Applied Mathematics
Columbia University
500 W. 120th Street
New York, NY 10027 USA

ABSTRACT
We consider gradient estimation with only noisy function evaluation, where the function can only beevaluated at values lying within a probability simplex. We are interested in obtaining gradient estimators
where each (pair of) data collection or simulation run applies simultaneously to all directions at once. Our
problem is motivated from the use of stochastic approximation in distributionally robust simulation analysis,
which involves solving for worst-case input distributions in a black-box simulation model. In this context,
conventional gradient schemes such as simultaneous perturbation face challenges as the required moment
conditions that allow the “canceling” of higher-order error terms cannot be satisfied without violating thesimplex constraints. We investigate a new set of required conditions on the probability distribution thatgoverns the perturbation, which leads us to a class of implementable gradient estimators using Dirichlet
mixtures. We study the statistical properties of these estimators and demonstrate their effectiveness withnumerical results.
1

INTRODUCTION

We study stochastic gradient estimation with only noisy function evaluation, where the function can only
be evaluated at values lying within a probability simplex. More precisely, given a real-valued function
Z(·) that defines on an n-dimensional set of probability weights, we want to compute its (properly defined)
gradient ∇Z, by using only noisy observations of Z at suitably chosen points. Our particular focus is on
schemes that can simultaneously estimate all components of the gradient with only a single or a pair ofobservations on Z.
Gradient estimation with only noisy function evaluation arises ubiquitously in simulation optimization,
when one uses stochastic approximation (SA) or gradient descent and the simulation model is only
accessible as a “black-box” (e.g., Glasserman 2013; Asmussen and Glynn 2007; Fu 2006; L’Ecuyer 1991).
In derivative-free optimization, this is also known as the zeroth-order gradient oracle (Ghadimi and Lan2013). The classical Kiefer-Wolfowitz (KW) scheme (Kushner and Yin 2003) is applicable precisely inthe above optimization setting. At each iteration, KW uses a pair of noisy evaluation to compute a finite
difference that approximates the partial derivative at each dimension, and convergence is provably attained bysuitably controlling the perturbation size and iteration step size. The number of evaluations at each iteration
required by this scheme, however, grows with the problem dimension. To address this, schemes that allow
a simultaneous estimation of derivatives at all dimension have been investigated. These schemes randomly
generate a perturbation vector according to some well-chosen distribution, and evaluate the function at
the realized perturbation. Then, through multiplying by a factor related to the generating distribution (andpossibly subtracting another evaluation at the original point), an estimate of the entire gradient is obtained.

Lam and Zhang
Prominent examples of such approaches include the simultaneous perturbation stochastic approximation(SPSA) (Spall 1992; Spall 1997) that uses distributions with enough masses at zero, Gaussian smoothing(Nesterov and Spokoiny 2017), and uniform sampling (Flaxman et al. 2004) to generate the perturbation
vector.
Our investigation is motivated from the use of randomized gradient estimation schemes described abovein distributionally robust simulation analysis, a recent line of work to address input uncertainty problems.
The latter refers to the estimation of confidence bounds or output variability for the performance measureof interest, when the input distributions that feed into a simulation model are corrupted by uncertainty or
statistical noises. In the presence of input data, such problems can be handled by using various statisticalmethods including bootstrap resampling (Barton and Schruben 2001; Barton et al. 2014), delta method(Cheng and Holland 2004), empirical likelihood (Lam and Qian 2016) and Bayesian approaches (Chick
2001; Zouaoui and Wilson 2003). Distributionally robust simulation analysis targets especially at highuncertainty, nonparametric situations where the modeler may only have minimal information on the input
distributions, such as moment ranges based on expert opinions or with a small amount of data (Glasserman
and Xu 2014; Ghosh and Lam 2019; Ghosh and Lam 2015; Lam 2016). The method has also been proposed
to address uncertainty in directions beyond typical assumptions, such as dependency in the input processes(Lam 2018). In these situations, the method aims to find worst-case bounds on the performance measure
under the imposed information. It does so by solving optimization problems with objective being the target
performance measure and decision variables being the uncertain input distributions. The feasible region
is represented by the imposed information as constraints on the input distributions, often known as an
uncertainty set or ambiguity set. If this set contains the true input distribution with high confidence, thenthe resulting optimal value would bound the true value of the performance measure with at least the samelevel of confidence guarantee.
Our studied gradient estimator is used to solve the aforementioned worst-case optimization problems.
The objective function of these problems, which is the performance measure of interest, can be a blackbox simulation model in general where we can only observe noisy outputs given the input distribution.
Computing the worst-case performance measure requires running stochastic approximation or gradientdescent, and in this case thus falls into the setting of the zeroth-order gradient oracle where a randomized
perturbation is required at each iteration. Note that, since the decision variable is the input distribution, it
lies in the space of probability distributions, and the objective is typically undefined or incomputable when
the decision variable is outside this region. Thus, in order to run the zeroth-order gradient-based iteration,
we need to obtain a gradient estimator on the constrained space, which is precisely our focus in this paper.
Though not primarily motivated, other potential applications of our investigation include SA approaches
used in solving bilinear games and stochastic utility problems (Nemirovski et al. 2009).
A challenge in obtaining gradient estimator via random perturbation, compared to previously studied
unconstrained settings, lies in the satisfaction of the distributional and moment conditions required forthe random perturbation vector. These conditions ensure that the higher-order error terms in the function
evaluation or finite difference under the perturbation can be sufficiently “canceled out”, so that the resulting
gradient estimator enjoys good bias properties. For instance, the standard version of SPSA requires theperturbation vector to have independent, mean-zero, symmetric components with finite reciprocal moments(Spall 1992). Because of the constraints in our problem, however, the perturbation vector we generatemust have correlated and unsymmetric components, and with restrictive moment behaviors, thus making
these previous conditions difficult to satisfy. Moreover, challenges also arise in using the idea of central
finite-differencing in these schemes, as the probability constraint may, depending on the evaluation pointof interest, only allow movement in either the forward or backward direction. To overcome all thesechallenges, in this paper we present a new set of conditions to ensure that the perturbation vector within
the probability simplex leads to higher-order cancellation. We then propose a specific set of estimators
using the Dirichlet mixture that satisfies these new conditions. We show how these estimators are readily

Lam and Zhang
implementable by offering precise guidelines. We also demonstrate their effectiveness in both estimating
gradients and in applying to SA.
2

PROBLEM SETTING AND MOTIVATION

Consider the function Z(·) : S → R, whereS = {p = (p1 , . . . , pn ) ∈ Rn : p1 + · · · + pn = 1, pi ≥ 0 ∀i = 1, . . . , n}
is an n-dimensional probability simplex. We are interested in estimating ∇Z(p), the gradient of Z. To be
more precise, note that, because of the probability-simplex constraint, an arbitrary perturbation of p needed
to define a gradient may shoot outside the domain of Z, and thus, when speaking of gradient, we mean adirectional gradient ∇Z(p) ∈ Rn defined such that, for any perturbation of p ∈ S to q ∈ S , we haveZ(q) − Z(p) = ∇Z(p)0 (p − q) + o(kp − qk)(1)

Equivalently, we can define each component of ∇Z(p) = (Z1 (p), . . . , Zn (p))0 as
Zi (p) = lim

ε→0+Z((1 − ε)p + εei ) − Z(p)ε(2)

where ei is a vector that takes value 1 at component i and 0 otherwise. That is, Zi (p) is defined by mixing
p with a point mass at component i and evaluating the increment as the mixture parameter shrink to 0. It is
straightforward to see that (1) and (2) coincide, up to a translation of c1 for any constant c and 1 being thevector of 1’s (note that (1) as defined is non-unique up to this translation, i.e., for any ψ(·) satisfying (1),
ψ(·) + c1 also satisfies (1) for any c). In robust statistics, definition (2) represents the so-called influencefunction (Hampel et al. 2011) that measures the sensitivity with respect to changes in the data distribution,
when viewing Z(p) as a statistical functional on the data distribution p.
Our goal is to estimate ∇Z(p) when given only the capability to output noisy function evaluations.
That is, given any q, we can output Ẑ(q) such that E[Ẑ(q)] = Z(q). Here we take the liberal view that
∇Z(p) can refer to any candidate that satisfies (1). We are particularly interested in estimating ∇Z(p) byrandomly generating a probability vector, say δ , and mix it with the original probability vector p. This canbe thought as using a perturbation vector δ − p. Our scheme then follows by evaluating Z at the perturbed
probability vector and “projecting” to get derivative estimates simultaneously for all dimensions.
More precisely, we consider several variants of estimator. The most basic scheme is to useẐ((1 − c)p + cδ )S(p, δ )c(3)

where c > 0 and S(p, δ ) ∈ Rn are chosen by user. When R simulation budget is available, we would
randomly draw δ j , j = 1, . . . , R, evaluate (3) each time with an independent evaluation of Ẑ, and at the end
output their average, i.e.,
1 R Ẑ j ((1 − c)p + cδ j )S(p, δ j )
(4)
∑
R j=1
c
where Ẑ j (·), j = 1, . . . , R denotes R independent function evaluations. We call (3) and (4) the singlefunction-evaluation (SFE) estimator.
In (3) and (4), c > 0 is a perturbation size that could depend on a given R. The vector S(p, δ ) represents
a factor or “score function” that aims to cancel out higher-order terms in the associated Taylor series inorder to control the bias.
Next, consider the use of two function evaluations, one at (1 − c)p + cδ and one at p itself. We output
Ẑ 1 ((1 − c)p + cδ ) − Ẑ 2 (p)S(p, δ )c(5)

Lam and Zhang
where Ẑ 1 (·) and Ẑ 2 (·) represent two independent evaluations, and c > 0 and S(p, δ ) ∈ Rn are chosen byuser like before. Analogously, when 2R simulation budget is available, we would use
1 R Ẑ 2 j−1 ((1 − c)p + cδ j ) − Ẑ 2 j (p)S(p, δ j )
∑
R j=1
c

(6)

where Ẑ 2 j−1 (·), Ẑ 2 j (·), j = 1, . . . , R are 2R independent function evaluations. We call (5) and (6) the
forward-function-evaluation (FFE) estimator.
Lastly, we also consider a variant similar to central finite-differencing. Consider using two function
evaluations, one at (1 − c)p + cδ and one at (1 + c)p − cδ . We output
Ẑ 1 ((1 − c)p + cδ ) − Ẑ 2 ((1 + c)p − cδ )S(p, δ )2c(7)

where Ẑ 1 (·) and Ẑ 2 (·) again represent two independent evaluations. Like FFE, when 2R simulation budget
is available, we would use
1 R Ẑ 2 j−1 ((1 − c)p + cδ j ) − Ẑ 2 j ((1 + c)p − cδ j )S(p, δ j )
∑
R j=1
2c

(8)

where Ẑ 2 j−1 (·), Ẑ 2 j (·), j = 1, . . . , R are 2R independent function evaluations. We call (7) and (8) the
central-function-evaluation (CFE) estimator.
Let us investigate and compare the statistical errors of the above estimators in estimating the gradient.
In doing so we will also highlight the involved challenges. Typically, all the estimators above have variancesthat scale in order 1/(Rc2 ), as the c appears in the denominator of these estimators. To understand thebias, supposing Z is sufficiently smooth, consider the Taylor series
Z((1 − c)p + cδ ) = Z(p) + c∇Z(p)0 (δ − p) +c2
(δ − p)0 ∇2 Z(p)(δ − p) + O(c3 )2(9)where ∇2 Z(p) is a properly defined Hessian. Now, since δ and the function evaluations are independent,
we can write the expectation of SFE as
Z((1 − c)p + cδ )E
S(p, δ )
c
 


c2
1
0
0 2
3
Z(p) + c∇Z(p) (δ − p) + (δ − p) ∇ Z(p)(δ − p) + O(c ) S(p, δ )= E
c2E[S(p, δ )]c
= Z(p)+ ∇Z(p)0 E[(δ − p)S(p, δ )] + tr(E[∇2 Z(p)(δ − p)(δ − p)0 S(p, δ )]) + O(c2 )(10)c2Thus, if we can ensure that
E[S(p, δ )] = 0
E[(δ − p)S(p, δ )] = I
2tr(E[∇ Z(p)(δ − p)(δ − p)0 S(p, δ )]) = 0

(11)
(12)
(13)

where I is the identity matrix and 0 are zero matrices each in the suitable dimension, then SFE gives rise
to a gradient estimator that has bias O(c2 ). If, however, (13) cannot be satisfied, then the bias could beO(c), which is less desirable.

Lam and Zhang
It is easy to see that FFE follows a similar, though slightly different, behavior. Its expectation isZ((1 − c)p + cδ ) − Z(p)E
S(p, δ )cc= ∇Z(p)0 E[(δ − p)S(p, δ )] + tr(E[∇2 Z(p)(δ − p)(δ − p)0 S(p, δ )]) + O(c2 )2Compared with (10), the term Z(p)E[S(p, δ )]/c disappears. Thus, if we can ensure that
E[(δ − p)S(p, δ )] = I
tr(E[∇ Z(p)(δ − p)(δ − p)0 S(p, δ )]) = 0
2

(14)
(15)

then we have bias O(c2 ), and similarly if only (14) is enforced, then bias is O(c). As we will see, (11) canbe easily satisfied, and thus it may appear that SFE is preferable as it only requires one function evaluation
to get one “gradient sample”. However, the first term Z(p)E[S(p, δ )]/c has a c in the denominator, so witha finite simulation run this term could add significantly to the variance. This is intuitive as SFE does notresemble a finite-difference scheme, the latter capturing the differencing needed in evaluating a derivative.
Next we discuss CFE. Thanks to the central finite-differencing, the first and third-order terms in the
Taylor series of Z((1 − c)p + cδ ) and Z((1 + c)p − cδ ) cancel out, so that
Z((1 − c)p + cδ ) − Z((1 + c)p − cδ )S(p, δ ) = ∇Z(p)0 E[(δ − p)S(p, δ )] + O(c2 )E
2cThus, to ensure a bias O(c2 ), we only requireE[(δ − p)S(p, δ )] = I
This appears to be desirable. However, note that (1 + c)p − cδ is not a proper mixture and may notcorrespond to a probability distribution.
Thus, taking into account the variance and the implementability, FFE appears the best one among thethree candidates. We will see that (11) and (12) are quite easy to enforce. In this sense, CFE can result inO(c2 ) bias, but it is not implementable generally due to the probability simplex constraint. SFE and FFE,on the other hand, can also ensure bias O(c) easily, and FFE has significantly better variance with a priceof needing twice as many function evaluations.
The main challenge, however, is that without the central finite-differencing trick, it is not obvious how
to lift FFE and SFE to a bias O(c2 ). This would be the main investigation in the next section.
To close this section, we describe how the above estimators are applied. They are primarily motivatedfrom distributionally robust simulation analysis. Here, think of Z(P) as a performance measure that depends
on some input distribution P. For instance, P could denote the interarrival or service times in a queueing
system, and Z(P) represents some output performance measure such as the expected wait time. When P
is uncertain but some information is available, we can consider computing the worst case bounds for the
performance measure given by
min Z(P) and max Z(P)
(16)
P∈U

P∈U

where U is the uncertainty set that contains information on P. For instance, U can denote momentand support constraints, so that U = {P : EP [ fl (X)] ≤ µl , l = 1, 2, ..., s, supp P = A} for some moment
functions fl (·)’s, µl ∈ R, X denotes a generic random variable under distribution P, and A is a subset
in the domain of P. As another example, one could also consider a neighborhood of a baseline model
measured by statistical distances such as the φ -divergence, in which case U = {P : dφ (P, Pb ) ≤ η} wheredφ (P, Pb ) denotes the φ -divergence from some baseline distribution Pb . Depending on the information,these constraints or combinations of them can be used in a given context.

Lam and Zhang
One motivation of our investigated gradient estimator is to iteratively solve (16) using SA when Z isaccessible only via noisy function evaluation. Note that, in many simulation analysis, P can be a continuous
distribution. However, one can suitably discretize the distribution by using support points generated froma heavy-tailed distribution and solving the problem on the sampled support points (see Ghosh and Lam2019). This reduces to solving
min Z(p) and max Z(p)
(17)
p∈U

p∈U

Then, to solve say the minimization problem in (17), we can use for instance the Frank-Wolfe (FW) method
0 (q − p). Here, in view of the discussion above, weˆwhich requires at each iteration solves minq∈U ∇Z(p)can equivalently replace the gradient ∇Z(p) with any ψ(p) = c1 ∇Z(p) + c2 1 for some c1 > 0 and c2 ∈ R.
Taking into account the error in the gradient estimation, Frank-Wolfe stochastic approximation (FWSA)(Ghosh and Lam 2019) would solve minq∈U ψ̂(p)0 (q − p) where ψ̂(p) is an estimator of ψ(p).
3IMPLEMENTABLE GRADIENT ESTIMATORS WITH DIRICHLET MIXTURES

In view of the discussion in Section 2, we aim to get an estimator for ψ(p), an equivalent to ∇Z(p), by
using random perturbation that satisfies (11) to (12). We will focus mostly on SFE and FFE as they are
more readily implementable, but we would also discuss briefly how CFE can also be used if one impose
extra requirement on the initial solution in the associated SA scheme. In particular, we will base our choiceon the use of a mixture of Dirichlet distributions, i.e.,Kδ = ∑ θi ∆i ,θi > 0, ∆i ∼ Dir(α i ), α i = (αi,1 , αi,2 , ..., αi,n )

i=1

where (θ1 , . . . , θK ) ∈ RK are probability weights. Moreover, we set the score function S(p, δ ) = γ(δ − p)for some γ ∈ R.
For our development, the following facts on Dirichlet distributions will be useful. Let ∆i = (∆i,1 , ∆i,2 , ..., ∆i,n ) ∼Dir(αi,1 , αi,2 , ..., αi,n ), and αi,0 = ∑nj=1 αi, j . We haveE(∆i, j ) =Var(∆i, j ) =

αi, j
αi,0

αi, j /αi,0 (αi, j /αi,0 )2−
αi,0 + 1
αi,0 + 1

Cov(∆i, j , ∆i,k ) = −Thus,
E((∆i − E(∆i ))(∆i − E(∆i ))0 ) =(αi, j /αi,0 )(αi,k /αi,0 )
αi,0 + 1

diag(α i /αi,0 ) (α i /αi,0 )(α i /αi,0 )0−αi,0 + 1
αi,0 + 1

where α i = (αi,1 , αi,2 , ..., αi,n ).
3.1 Explicit Moment Conditions
We begin with a more explicit set of conditions to replace (11) to (13). To streamline presentation, fromnow on we skip the argument p in ∇Z(p) and ∇2 Z(p) and simply use ∇Z and ∇2 Z. We have the following:
Theorem 1 If δ satisfies
E[δ ] = pγE[(δ − p)0 (δ − p)] = I −

(18a)

11

0

n
E[(δ − p)i (δ − p) j (δ − p)k ] = µ ∀ (i, j, k) ∈ [n] × [n] × [n]

(18b)
(18c)

Lam and Zhang
for some µ, γ ∈ R, assume Z is trice differentiable, and δ has bounded fourth moments, and we let our
gradient estimator beγ
ψ̂(p) = Ẑ((1 − c)p + cδ )(δ − p)cThencµγ 0 2
10 ∇Z
E[ψ̂(p)] = ∇Z + (1 ∇ Z1 −)1 + O(c2 )2nProof.

Using the Taylor expansion of Z((1 − c)p + cδ ) in (9), we getγγcE[ψ̂(p)] = Z(p)E[(δ − p)] + γE[(δ − p)(δ − p)0 ]∇Z + E[(δ − p)0 ∇2 Z(δ − p)(δ − p)] + O(c2 )c2110cµγ 0 2= (I −)∇Z + (
1 ∇ Z 1)1 + O(c2 )n2cµγ 0 210 ∇Z
= ∇Z + (1 ∇ Z1 −)1 + O(c2 )2nThus, Theorem 1 implies that if we can enforce the third-order constraint in (18) with µ = 0, then we
can have an estimator with bias O(c2 ) in estimating ψ(p).
3.2 An Eligible Dirichlet Mixture
We will find a specific implementable δ that satisfies the conditions in (18), by using a mixture of
K = n(n−1)/2+1 of Dirichlet distributions. Without loss of generality, we assume p = (p1 , p2 , ..., pn ), p1 ≤p2 ≤ ... ≤ pn . Moreover, each Dirichlet only perturbs two variables, and putting all of them together will
have enough degree of freedom to control the covariance matrix and the third-order constraints in (18).
Now considern−1n∑ θ l, j ∆l, j + θ n ∆nδ=∑

l=1 j=l+1

where ∆l,i ∼ Dir(α l,i ) and ∆n ∼ Dir(α n ), with
θ 1 = θ 1,i = 2p1 /(n − 1), i = 2, 3, ..., nθ 2 = θ 2,i = (2p2 − θ 1 )/(n − 2), i = 3, 4, ..., n...
l−1

θ l = θ l,i = (2pl − ∑ θ j )/(n − l), i = l + 1, ..., n
j=1

...
n−2

θ n−1 = θ n−1,n = (2pn−1 − ∑ θ j )
j=1

n−1

θ n = pn − ∑ θ j /2
j=1

(19)

Lam and Zhang
and0
α l,i /αl,i= (ei + el )/20αl,i= 10 α l,i = C(θ l,i )2 − 1

α n = en
where C is large constant such that C(θ l,i )2 > 1, and γ = 4C/n.
Theorem 2 The distribution δ chosen as (19) with the depicted parameter choices satisfies (18) with µ = 0.
Theorem 2 stipulates that the choice of δ in (19) can be used in SFE and FFE to obtain a bias O(c2 ).
Its proofs is skipped due to space limit.
The Dirichlet mixture δ for SFE and FFE works for CFE, assuming the “backward” perturbed distribution(1 + c)p − cδ is well-defined. Moreover, in such situations, we can obtain a workable δ that mixes only nDirichlet distributions instead of order n2 , which is depicted as follows. Let δ = ∑nl=1 θ l ∆l , ∆l ∼ Dir(α l )where(
np1 (l = 1)θl =pl − p1 (l = 2, 3, ..., n)(
1/n (l = 1)αl =el (l = 2, 3, ..., n)and γ = np2 2 . We then have
1

n

n

l=2

l=1

E[δ ] = p1 1 + ∑ (pl − p1 )el = ∑ pl el = pandnp2
γCov(δ ) = γ 1
2



110
110
diag(1) −=I−nn

We briefly indicate how one may avoid the possibility of the “backward” perturbation lying outside theprobability simplex in the CFE. If we start with an initial solution (i.e., distribution) that has positive massat each support point. Then, by properly tuning the step size and perturbation size at each iteration, theprobability mass at each support point can still remain sufficiently positive, so that moving down a massby the next perturbation size would still be allowable. These choices can be obtained in such a way thatstill allows convergence under standard assumptions, via an analysis on the constants needed in the scalingof these sizes.
4

NUMERICAL RESULTS

In this section, we conduct numerical experiments to compare the performance of the three gradient estimators.
We also experiment the FWSA and its variant, a mirror descent (MD) SA algorithm (Nemirovski et al.
2009; Ghosh and Lam 2015). We consider a deterministic example: Z(p) = p0 A0 Ap + B0 p + C wherep ∈ R10 . (Here we generate A ∈ R10×10 , B ∈ R10 , C ∈ R randomly such that each entry of them follows
Uni f orm(0, 1) independently). Thus, we can evaluate Z exactly and so var(Ẑ) = σ 2 = 0. Note that, though
the example is deterministic, the bias and variance discussed in previous sections still hold due to therandomization in the perturbation direction.
We consider solving the minimization problem in (17) with moment constraints
i−1∑ pi xi ∈ [0.3, 0.7], ∑ pi xi2 ∈ [0.2815, 0.4222], xi = 9 , i = 1, ..., 10iiLam and Zhang
We choose these bounds because setting p1 = ... = p1 = 1/10, we have ∑i pi xi = 0.5, ∑i pi xi2 ≈ 0.3519, so
this set of moment constraints describes a set that is close to the uniform distribution.
The experiment has 2 parts: first, we show that the proposed Dirichlet mixture distribution indeed
satisfies the constraints in (18) ((c) for SFE, FFE) and we compare the mean squared error (MSE) of
the three gradient estimators under different R and c. Second, we test the proposed FWSA and MDSA
algorithm for the convex quadratic function Z.
4.1 Performances of Gradient Estimators
First, we test the distributions for the proposed Dirichlet mixture distributions for Estimator I (SFE), II(FFE) and III (CFE) (For I and II we use the same distribution of δ , so we only show the plot for II,
and for III we use the simpler version that only uses n Dirichlets in the mixture). To see this, we show
E[δ ], γE[(δ − p)(δ − p)0 ], E[(δ − p)i (δ − p) j (δ − p)k ] as estimated in the 200th iteration during the SPSA
scheme minimizing Z. Here, we choose R = 500 for each iteration, and the perturbation size at k = 200 is−4
ck = 0.1
k = 5 × 10 .

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: Distribution of δ . (abc): Estimator I, II; (def): Estimator III.
From Figure 1 we see that for all three estimators, E[δ ] = p and γE[(δ − p)(δ − p)0 ] = I − 110 /n aresatisfied. The upper and lower rows in (a) and (d) are almost the same, and the diagonals of (b) and (e) areapproximately 0.9, while the off-diagonals are approximately -0.1. By design, however, only Estimator I
and II satisfy E[(δ − p)i (δ − p) j (δ − p)k ] = 0. Indeed, from the histogram over all combination of i, j, k,(c) shows that (δ − p)i (δ − p) j (δ − p)k ∼ 10−7 are close to 0, while in (f) the values are much larger, andthere are outliers. Thus, in Estimator I and II, we need to generate more Dirichlet for each δ , but then thegenerated δ appears well behaved such that it “cancels out” the O(c2 ) terms in the Taylor expansion.
Second, we use Estimators I, II, and III to estimate ∇Z(p0 ) where p0 is randomly generated in theprobability simplex. We vary c (the perturbation size) and R (the number of δ used in each estimation).
We run 100 experiments for each set of parameters and use MSE to measure the performance.

Lam and Zhang
Notice that we know the true gradient at p0 : ∇Z(p0 ) = 2A0 Ap0 + B. And from analysis we know thatthe bias for each estimator is
10 (2A0 Ap0 + B)
10 ∇Z(p0 )1=−1bias ≈ −
n
n

(a)

(b)

(c)

(d)

Figure 2: MSE against c.
Figure 2 tells us several things:
1. From Figures (b) and (d), we see that the MSE of Estimator I is much larger than Estimator II and
III. In addition, after choosing R = O(c−2 ) as a function of c rather than a constant, the MSE stays at thesame level, indicating that its MSE is of order R−1 c−2 .
2. The difference between the MSE and the squared bias can be seen as an estimate of the variance ofthe estimator. So from figure (a) we see that simply changing c has no effect on variance reduction; while
figure (c) shows that to reduce the variance, we can decrease c and increase R together.
3. Based on Figure 2, when choosing which estimator to use in the optimization problem, we should
never choose Estimator I due to its large variance. Estimators II and III are comparable in terms of theMSE, but II requires more Dirichlets in the distribution for δ , and III puts constraints on the perturbation
size c such that the perturbed p is still in the simplex. So which one to use might depend on the actual
function Z we are interested in. Generally speaking, II is preferable as the overhead in using more Dirichlets
is relatively negligible, and its being free of the constraints on c makes it more flexible against III.
4.2 Performances of FWSA and MDSAIn this section, we use FWSA and MDSA to minimize Z. For FWSA, we choose the update stepsize at
iteration k to be γk = 0.05/k, the perturbation stepsize ck = 0.1/k1/4 . For MDSA, we choose γk = 1/k1.5 in
the proximal operator, and the same perturbation size ck as FWSA. (The MDSA scheme with estimator I
gives “nan” during optimization, so we only show the result for FW+Estimator I,II,III and MD+Estimator
II,III.) We choose R = 500 in each iteration.

Lam and Zhang
We evaluate the performance based on the Frank-Wolfe (FW) gap (for FWSA)
g(p) = − min ψ(p)0 (q − p)q∈U

and on the function value Z.

Figure 3: FWSA vs. MDSA.
From Figure 3 we see that all 4 experiments using Estimator II/III are converging, but Estimator I either
blows up the optimization (as in MDSA) or does not decrease the function at all (as in FWSA). This might
be explained by the poor gradient estimation it provides. However, it is interesting that MDSA converges
much faster than FWSA. After 300 iterations, FWSA is still well above MDSA (and is still decreasing).
This might be because we did not choose the best set of parameters for the optimization, or might be dueto the superiority of MDSA itself. In future work, we plan to test on various functions with different setsof parameters for more extensive investigations.
5

CONCLUSION

In this paper, we investigated gradient estimators using only noisy function evaluations when the input
is constrained to be a probability distribution, a problem motivated from running gradient descent for
distributionally robust simulation analysis. Traditional random perturbation approaches encounter difficulties
in satisfying the conditions needed to cancel out higher-order bias due to the probability simplex constraint.
We proposed a new set of conditions for distributionally constrained random perturbation, and studied
three types of gradient estimators, SFE, FFE and CFE, by properly mixing the input distribution witha perturbation vector. We showed that while CFE, like in the classical setting, is the easiest to cancelout higher-order bias, the backward perturbation involved may not satisfy the simplex constraint. SFE
and FFE can cancel out higher-order bias with more sophisticated perturbation designs, and FFE tends tohave a lower variance than SFE due to first-order term canceling. Based on these, we suggested FFE as
the recommended choice. We provided two explicit implementations of these estimators using carefully
tuned Dirichlet mixtures. Experimental results supported our theoretical findings, and also showed that the
proposed estimators are effective for use in FWSA and MDSA. In particular, FFE or CFE together withMDSA appeared to have the best performance in terms of convergence speed.
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CAREER CMMI1834710 and IIS-1849280.

Lam and Zhang
REFERENCES
Asmussen, S., and P. W. Glynn. 2007. Stochastic simulation: algorithms and analysis, Volume 57. Springer Science & BusinessMedia.
Barton, R. R., B. L. Nelson, and W. Xie. 2014. “Quantifying input uncertainty via simulation confidence intervals”. INFORMS
journal on computing 26(1):74–87.
Barton, R. R., and L. W. Schruben. 2001. “Resampling methods for input modeling”. In Proceeding of the 2001 Winter
Simulation Conference (Cat. No. 01CH37304), Volume 1, 372–378. IEEE.
Cheng, R. C., and W. Holland. 2004. “Calculation of confidence intervals for simulation output”. ACM Transactions on Modelingand Computer Simulation (TOMACS) 14(4):344–362.
Chick, S. E. 2001. “Input distribution selection for simulation experiments: accounting for input uncertainty”. Operations
Research 49(5):744–758.
Flaxman, A. D., A. T. Kalai, and H. B. McMahan. 2004. “Online convex optimization in the bandit setting: gradient descent
without a gradient”. arXiv preprint cs/0408007.
Fu, M. C. 2006. “Gradient estimation”. Handbooks in operations research and management science 13:575–616.
Ghadimi, S., and G. Lan. 2013. “Stochastic first-and zeroth-order methods for nonconvex stochastic programming”. SIAM
Journal on Optimization 23(4):2341–2368.
Ghosh, S., and H. Lam. 2015. “Mirror descent stochastic approximation for computing worst-case stochastic input models”.
In 2015 Winter Simulation Conference (WSC), 425–436. IEEE.
Ghosh, S., and H. Lam. 2019. “Robust analysis in stochastic simulation: Computation and performance guarantees”. Operations
Research 67(1):232–249.
Glasserman, P. 2013. Monte Carlo methods in financial engineering, Volume 53. Springer Science & Business Media.
Glasserman, P., and X. Xu. 2014. “Robust risk measurement and model risk”. Quantitative Finance 14(1):29–58.
Hampel, F. R., E. M. Ronchetti, P. J. Rousseeuw, and W. A. Stahel. 2011. Robust statistics: the approach based on influencefunctions, Volume 196. John Wiley & Sons.
Kushner, H., and G. G. Yin. 2003. Stochastic approximation and recursive algorithms and applications, Volume 35. Springer
Science & Business Media.
Lam, H. 2016. “Robust sensitivity analysis for stochastic systems”. Mathematics of Operations Research 41(4):1248–1275.
Lam, H. 2018. “Sensitivity to serial dependency of input processes: A robust approach”. Management Science 64(3):1311–1327.
Lam, H., and H. Qian. 2016. “The empirical likelihood approach to simulation input uncertainty”. In 2016 Winter Simulation
Conference (WSC), 791–802. IEEE.
L’Ecuyer, P. 1991. “An overview of derivative estimation”. In Winter Simulation Conference, 207–217.
Nemirovski, A., A. Juditsky, G. Lan, and A. Shapiro. 2009. “Robust stochastic approximation approach to stochastic programming”.
SIAM Journal on optimization 19(4):1574–1609.
Nesterov, Y., and V. Spokoiny. 2017. “Random gradient-free minimization of convex functions”. Foundations of Computational
Mathematics 17(2):527–566.
Spall, J. C. 1992. “Multivariate stochastic approximation using a simultaneous perturbation gradient approximation”. IEEE
transactions on automatic control 37(3):332–341.
Spall, J. C. 1997. “A one-measurement form of simultaneous perturbation stochastic approximation”. Automatica 33(1):109–112.
Zouaoui, F., and J. R. Wilson. 2003. “Accounting for parameter uncertainty in simulation input modeling”. Iie Transactions 35(9):781–792.

AUTHOR BIOGRAPHIESHENRY LAM is an Associate Professor in the Department of Industrial Engineering and Operations Research at Columbia
University. He received his Ph.D. degree in statistics from Harvard University in 2011, and was on the faculty of Boston
University and the University of Michigan before joining Columbia in 2017. His research focuses on Monte Carlo simulation,
uncertainty quantification, risk analysis, and stochastic and robust optimization. His email address is henry.lam@columbia.edu.
JUNHUI ZHANG is an undergraduate junior student at Columbia University, majoring in Applied Mathematics. Her email
address is jz2903@columbia.edu.