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.

OPTIMALLY TUNING FINITE-DIFFERENCE ESTIMATORS
Haidong Li

Henry Lam

Department of Industrial Engineering andManagement
Peking University
5 Yiheyuan Road
Beijing 100871, P. R. CHINA

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

ABSTRACT
We consider stochastic gradient estimation when only noisy function evaluations are available. Central
finite-difference scheme is a common method in this setting, which involves generating samples under
perturbed inputs. Though it is widely known how to select the perturbation size to achieve the optimalorder of the error, exactly achieving the optimal first-order error, which we call asymptotic optimality,
is considered much more challenging and not attempted in practice. In this paper, we provide evidence
that designing asymptotically optimal estimator is practically possible. In particular, we propose a new
two-stage scheme that first estimates the required parameter in the perturbation size, followed by running
finite-difference based on the estimated parameter in the first stage. Both theory and numerical experiments
demonstrate the optimality of the proposed estimator and the robustness over conventional finite-difference
schemes based on ad hoc tuning.
1

INTRODUCTION

In this paper, we consider finite-difference stochastic gradient estimation (e.g., Glasserman 2013, Asmussen
and Glynn 2007, Fu 2006, L’Ecuyer 1991), commonly used when only noisy simulation observations areavailable to evaluate the function value or model output. In stochastic optimization, such setting is known
as black-box or zeroth-order (Ghadimi and Lan 2013, Nesterov and Spokoiny 2017). Finite-difference
estimators are in contrast to unbiased derivative estimators, which include the infinitesimal perturbation
analysis (Ho et al. 1983, Heidelberger et al. 1988), the likelihood ratio or the score function method (Glynn
1990, Rubinstein 1986, Reiman and Weiss 1989), measure-valued or weak differentiation (Heidergott andVázquez-Abad 2008, Heidergott et al. 2010), and other variants such as the generalized likelihood ratio
method (Peng et al. 2018).
Finite-difference estimators consist of generating samples under perturbed input parameters. This
perturbation size is chosen with consideration of controlling the bias and variance that contributes tothe overall mean squared error (MSE). As the perturbation size increases, bias increases while variancedecreases (and vice versa). To minimize the MSE, the perturbation size is determined by balancing themagnitudes of the two error sources. It is widely known (Zazanis and Suri 1993, Fox and Glynn 1989) that,for twice continuously differentiable functions, the optimal perturbation size of the central finite-difference(CFD) scheme turns out to be of order n−1/6 which leads to an optimal MSE order n−2/3 , where nrefers to the number of differencing pairs in the simulation. On the other hand, in forward or backward
finite-difference scheme, the optimal perturbation size is of order n−1/4 and the corresponding optimal
MSE order deteriorates to n−1/2 . Note that although the optimal order of error, in terms of n, is widely
known and achievable, it is considered much more challenging to obtain an estimator that achieves theexact optimal first-order MSE (i.e., the “constant” in front of the n). This is because it requires additionalLi and Lam
model information such as higher-order derivatives and the noise variance, which are unknown a priori.
Indeed, as far as we know, there have been no mainstream estimators that attempt to optimize the MSE at
this level of accuracy.
In this paper, we provide some evidence that it is, in fact, practically possible to obtain finite-difference
estimators that guarantee to exactly optimize the first-order MSE. We call such estimator asymptotically
optimal. More precisely, we propose a new, conceptually simple, two-stage scheme to obtain such an
estimator. This scheme first estimates the parameters required to tune the perturbation size that optimizes
the first-order MSE (the “Estimation” phase). Then, at the second stage, we use the estimated parameters
to obtain a nearly optimal perturbation size (the “Minimization” phase) and run standard finite-difference(we will focus on CFD due to its superiority over other variants, though our framework can be readily
generalized). The key reason for the implementability of this approach is that the parameter estimationthat allows a nearly optimal CFD only requires a loose accuracy, so that few samples need to be allocatedto the Estimation phase. We materialize the above intuition and show that this two-stage scheme, which
we call Estimation-Minimization Central Finite-Difference (EM-CFD), is nearly asymptotically optimal(where “nearly” means that asymptotic optimality can be achieved when the allocation to the Estimation
phase is relatively negligible). Compared to conventional CFD using the right order of perturbation size,
but with the underlying constant chosen in an ad hoc fashion, our estimator is more robust as it performs
consistently close to an “oracle” CFD that assumes knowledge of unknown model parameters. We support
our theory with empirical results.
The rest of the paper is organized as follows. In Section 2, we introduce the setting of CFD andthe challenge in achieving asymptotic optimality. Section 3 proposes a two-stage scheme to optimize theperformance of CFD. Sections 4 and 5 discuss the asymptotic properties of parameter estimation and theoverall MSE, respectively. Section 6 presents numerical results, and Section 7 concludes the paper andoutlines future directions.
2

SETTING AND MOTIVATIONIn this paper, we focus our discussions on the single-dimensional case. Let f (·) : R → R be a performancemeasure of interest, where we have access to an unbiased estimate fˆ(x) for any chosen x ∈ R. In otherwords, fˆ(·) is a family of random variables indexed by x such that E[fˆ(x)] = f (x) and V ar(fˆ(x)) = σ 2 (x)for any x ∈ R. Suppose we do not apply common random numbers (CRNs) in generating fˆ(x), and thusfˆ(x)’s are assumed to be independent across different points x. We would like to estimate the first-orderderivative f (1) (x0 ) where x0 ∈ R is the point of interest.
In estimating f (1) (x0 ), the CFD scheme elicits the output
Y (δ) =fˆ(x0 + δ) − fˆ(x0 − δ),2δ

where δ > 0 is the perturbation size. Suppose that f (x) is thrice continuously differentiable with non-zerothird-order derivative f (3) (x0 ), we have as δ → 0
Y (δ) = f (1) (x0 ) + (Bδ 2 + o(δ 2 )) +

(δ)
,
δ

where B = f (3) (x0 )/6 and (δ) ∈ R is a random variable such that E[(δ)] = 0 and V ar((δ)) = η 2 (δ).
Suppose we do not apply common random numbers (CRNs) in generating fˆ(x0 + δ) and fˆ(x0 − δ), andthat V ar(fˆ(x0 ± δ)) → V ar(fˆ(x0 )) as δ → 0. Then η 2 (δ) → σ 2 (x0 )/2 as δ → 0. Given the capability
to output independent runs of Yi (δ), i = 1, . . . , n, the CFD scheme obtains an estimate of f (1) (x0 ) byP
taking the sample average, i.e., θb = (1/n) ni=1 Yi (δ).

Li and Lam
The MSE of θb can be expressed as
MSE = E[(θb − f (1) (x0 ))2 ]= (E[θb − f (1) (x0 )])2 + V ar(θb − f (1) (x0 ))= (B + o(1))2 δ 4 + (σ 2 (x0 )/2 + o(1))/(nδ 2 ),where o(1) means a term that goes to zero as δ goes to zero. Consider choosing the optimal δ, in relation ton, in order to minimize the MSE. First, it is well-known and easy to see that, in order to obtain the optimalorder of the MSE, one should balance the squared bias term B 2 δ 4 and the variance term (σ 2 (x0 )/2)/(nδ 2 )to the same order in terms of n. This amounts to setting δ to be order n−1/6 (since otherwise by perturbing
the order of δ either one of the two terms would increase), which leads to an optimal MSE order n−2/3 .
Next, to get the optimal δ more precisely, note that the MSE depends on the priori unknown parameters
B and σ 2 (x0 ). Note that by Holder’s inequality
B 2 δ 4 + (σ 2 (x0 )/2)/(nδ 2 ) ≥ 3 B 2 σ 4 (x0 )/(16n2 )

1/3

1/6
and equality holds if δ = σ 2 (x0 )/(4nB 2 ). That is, the MSE with an exactly optimal first-order termis
1/3
MSEopt = 3 B 2 σ 4 (x0 )/(16n2 )(1 + o(1))(1)
with the dominating order n−2/3 and the frontal constant 3(B 2 σ 4 (x0 )/16)1/3 . We call a scheme asymptotically optimal if it achieves (1).
Thus, though it is easy to find the optimal order of δ by setting δ = αn−1/6 for some α > 0 as in
the conventional CFD scheme, exactly achieving the optimal first-order MSE, or asymptotic optimality,1/6
requires α = σ 2 (x0 )/(4B 2 )and hence knowing the model information B and σ 2 (x0 ). In the literature,
extracting this information is often viewed as challenging, in the sense that doing so is not easier thanestimating the gradient directly. Our main contribution is to propose a scheme that extracts and utilizes thisinformation efficiently. Our key insight is that, in order to tune α close to optimal, it suffices to obtain only
moderately accurate estimates of B and σ 2 (x0 ). Thus, one do not need to devote too much computational
effort into this task, in order to achieve (1). This makes our procedure readily implementable.
3

OVERVIEW OF THE TWO-STAGE PROCEDURE

We propose a new two-stage scheme, which we call Estimation-Minimization Central Finite-Difference(EM-CFD), to obtain an asymptotically optimal estimator. The first stage is the Estimation stage (E-stage).
At the E-stage, we allocate n1 = bλnc, 0 < λ < 1 samples, each with a possibly different perturbation
size, to estimate B and σ 2 (x0 ). The detailed estimation method will be discussed in Section 4. The secondb and σstage is the Minimization stage (M-stage). Here, we plug the estimated parameters B
b2 (x0 ) into

1/6
the optimal perturbation size δ = σ 2 (x0 )/(4nB 2 ). Then allocate n2 = d(1 − λ)ne samples and run
standard CFD to estimate f (1) (x0 ). The full implementation of EM-CFD is shown in Algorithm 1.
4

ESTIMATING MODEL PARAMETERS

We present in detail the E-stage, namely on the estimation of B and σ 2 (x0 ). We use linear regression.
Since B is related to the third order derivative of f (x), we take the Taylor expansion of f (x) to justify theestimation accuracy of B. Suppose that f (x) is five-times continuously differentiable with non-zero finite
f (5) (x0 ) (the non-zero assumption can be dropped, but we put this to streamline our results), we have as
δ→0fˆ(x0 + δ) − fˆ(x0 − δ) = f (1) (x0 )2δ + B2δ 3 + (D2δ 5 + o(δ 5 )) + (δ)

Li and Lam
Algorithm 1: Estimation-Minimization Central Finite-Difference
The total number of samples is n, the allocation ratio for the E-stage is λ > 0;E-stage: Allocate n1 = bλnc samples whose perturbation sizes δi , i = 1, . . . , n1 , are i.i.d.
b and σgenerated from a distribution P. Then estimate B and σ 2 (x0 ) by B
b2 (x0 ).

1/6
b2)
.
M-stage: Allocate n2 = d(1 − λ)ne samples with perturbation size δ = σ
b2 (x0 )/(4n2 B
Then estimate f (1) (x0 ) by1
θbEM =n2n2 ˆ
Xfi (x0 + δ) − fˆi (x0 − δ)
i=1

2δ

.

return the gradient estimate θbEM .
where B = f (3) (x0 )/6, D = f (5) (x0 )/120, and (δ) ∈ R is a random variable such that E[(δ)] = 0 andV ar((δ)) → 2σ 2 (x0 ) as δ → 0. Lethi0
y = fˆ(x0 + δ1 ) − fˆ(x0 − δ1 ), . . . , fˆ(x0 + δn1 ) − fˆ(x0 − δn1 ) ,02δ1 · · · 2δn1
X=,
2δ13 · · · 2δn3 1
hi0
β = f (1) (x0 ), B ,0
E = D2δ15 + (δ1 ), . . . , D2δn5 1 + (δn1 ) ,where δ1 , . . . , δn1 are the perturbation sizes of samples at the E-stage. Then, y = Xβ + E, and we canestimate β given by βb = (X 0 X)−1 X 0 y. Specifically, we havePn1 2 Pn1 3
P 1 4 Pn1
δi i=1 δi yi − ni=1δ
δ i yii=1Pb
P i 1 i=1
P
(2)
B=
1
1
δi4 )2 )δi6 − ( ni=1
δi2 ni=12( ni=1
where yi represents the i-th component in y.
Correspondingly, we havePn1 2 Pn1 3
P 1 4 Pn1
δi i=1 δi Ei − ni=1δ
δ i Eii=1Pb
Pn1 6
Pni 1 i=1
B−B =
.
n1
2
4
2
2( i=1 δi i=1 δi − ( i=1 δi ) )Without loss of generality, suppose δi is i.i.d. generated from a distribution with E[δ0 ] = 0 and V ar[δ0 ] = Σ0 ,and Σ0 → 0 as n1 → +∞. Additionally, any order moment of δi is assumed to be finite, which allows us
to apply the Central Limit Theorem. Then we have
P 1
 

(1/n1 ) P n
2DE[δ06 ]i=1 δi Ein1(1/n1 ) Pi=1
δi3 Ei  2DE[δ08 ]

√
d

 

n1
δ 2  −  E[δ 2 ] ) −−−−−−→n1 ( (1/n1 ) Pi=1 (1/n ) n1 δi4   E[δ04 ]  n1 →+∞1 Pi=1 i
0
6
1(1/n1 ) n
E[δ06 ]
i=1 δi   2
4D Cov(δ06 , δ06 ) + 2σ 2 (x0 )E[δ02 ]0
0 4D2 Cov(δ08 , δ06 ) + 2σ 2 (x0 )E[δ04 ]  
2DCov(δ02 , δ06 )N (0 , 0 2DCov(δ04 , δ06 )0
2DCov(δ06 , δ06 )4D2 Cov(δ06 , δ08 ) + 2σ 2 (x0 )E[δ04 ]4D2 Cov(δ08 , δ08 ) + 2σ 2 (x0 )E[δ06 ]2DCov(δ02 , δ08 )2DCov(δ04 , δ08 )2DCov(δ06 , δ08 )2DCov(δ06 , δ02 )2DCov(δ08 , δ02 )Cov(δ02 , δ02 )Cov(δ04 , δ02 )Cov(δ06 , δ02 )2DCov(δ06 , δ04 )2DCov(δ08 , δ04 )Cov(δ02 , δ04 )Cov(δ04 , δ04 )Cov(δ06 , δ04 )
2DCov(δ06 , δ06 )8
6
2DCov(δ0 , δ0 )

Cov(δ02 , δ06 ) ).
4
6
Cov(δ0 , δ0 ) 
Cov(δ06 , δ06 )Li and Lam
With the multivariate delta method, we obtain√
D(E[δ02 ]E[δ08 ] − E[δ04 ]E[δ06 ])σ 2 (x0 )E[δ02 ]dbn1 (B − B) −
−
−
−
−
−
→
N
(0,
)
n1 →+∞
E[δ02 ]E[δ06 ] − E[δ04 ]E[δ04 ]2(E[δ02 ]E[δ06 ] − E[δ04 ]E[δ04 ])where we use Σ0 → 0 to argue the negligibility of higher-order terms.
b when using δi that is randomly drawn
The following theorem gives the expression for the MSE of B
from a normal distribution in an i.i.d. fashion. It shows that the MSE vanishes to zero as n1 → +∞. In
b is consistent.
other words, the estimate B
Theorem 1 Suppose f (·) is five-times continuously differentiable, and fˆ(·) has a finite second moment at
any point with V ar(fˆ(x0 ± δ)) → V ar(fˆ(x0 )) as δ → 0. If we randomly draw δ1 , . . . , δn1 from N (0, Σ0 ),b in (2), we haveand consider B
2 (x )
σ022b − B) ] = (10DΣ0 ) +E[(B(1 + o(1)),
12n1 Σ30
where o(1) means a term that goes to zero as n1 → +∞. Moreover, for any Σ0 = Θ(nk1 ), −1/3 < k < 0,b − B)2 ] = 0.
lim E[(Bn1 →+∞

Proof.

For δi ∼ N (0, Σ0 ), i = 1, . . . , n1 ,D(1 + o(1))(E[δ02 ]E[δ08 ] − E[δ04 ]E[δ06 ])E[δ02 ]E[δ06 ] − E[δ04 ]E[δ04 ]σ 2 (x0 )2= (10DΣ0 ) +(1 + o(1)).
12n1 Σ30b − B)2 ] =E[(BWhen Σ0 = Θ(nk1 ), −1/3 < k < 0, we have2+

σ 2 (x0 )(1 + o(1))E[δ02 ]2n1 (E[δ02 ]E[δ06 ] − E[δ04 ]E[δ04 ])lim (10DΣ0 )2 = 0 and

n1 →+∞σ 2 (x0 )= 0. Therefore,
n1 →+∞ 12n1 Σ30lim

b − B)2 ] = 0.
lim E[(Bn1 →+∞

Note that we can generalize the above to non-normal mean-zero δi , though normal distribution is a
−1/5
natural choice. Moreover, we can show further that choosing Σ0 = Θ(n1 ) leads to an optimal order−2/5
b − B)2 ]. Therefore, at the E-stage, we randomly take perturbation size δi = αi n−1/10 wheren1of E[(B1
αi is i.i.d. generated from a fixed normal distribution.

As for σ 2 (x0 ), we use the sample variance of outputs fˆ(x0 + δi ) − fˆ(x0 − δi ) to estimate it, i.e.,σb2 (x0 ) =1 ky − ȳ 1n1 ×1 k22
2
n1 − 1(3)

P 1
where ȳ = ni=1
yi /n1 and 1n1 ×1 is a column vector with all elements equal to one. The following
theorem provides the consistency of σb2 (x0 ). This holds as long as Σ0 → 0 as n1 → ∞, and is not limited−1/10
to the case that δi = αi n1.
Theorem 2 Suppose fˆ(·) has a finite second moment at any point with V ar(fˆ(x0 ± δ)) → V ar(fˆ(x0 ))as δ → 0. Suppose we randomly draw δ1 , . . . , δn1 from N (0, Σ0 ), with Σ0 → 0 as n1 → ∞. Considerσb2 (x0 ) in (3). We havelim σ
b2 (x0 ) = σ 2 (x0 ), a.s.
n1 →+∞

Li and Lam
Proof.
This follows quite immediately from the consistency of sample variance. To give more details,
since Σ0 → 0 as n1 → +∞, we have lim δi = 0, a.s. for any i. For each n1 ∈ N+ ,n1 →+∞

n

Combining

n

1
1
1 X
1 X
yi2 − (
yi )2
n1
n1

n1
σ
b (x0 ) =2(n1 − 1)
2

i=1

!
.

i=1

lim δi = 0, a.s. and the strong Law of Large Numbers, we haven1 →+∞

n

1
1 X
lim
yi2 = E[2 (0)], a.s.
n1 →+∞ n1

i=1

and

n

1
1 X
lim
yi = E[(0)], a.s.
n1 →+∞ n1

i=1

Therefore,lim σ
b2 (x0 ) =n1 →+∞5

 1
1
E[2 (0)] − (E[(0)])2 = V ar((0)) = σ 2 (x0 ), a.s.
2
2

OPTIMIZING FINITE-DIFFERENCE PERFORMANCEAt the E-stage, we estimate unknown model parameters B and σ 2 (x0 ) by (2) and (3). Then, at the M-stage,
we use these parameter estimates to obtain a nearly optimal perturbation size and estimate f (1) (x0 ) bystandard CFD. The error in parameter estimation will propagate to the MSE of CFD. The following theorem
provides the expression for the MSE of the resulting gradient estimator from EM-CFD. It also shows that thisMSE closely matches the optimal in the first order, in turn implying that EM-CFD is nearly asymptotically
optimal, when the allocation assigned to the E-stage is relatively small.
Theorem 3 Suppose f (·) is five-times continuously differentiable, and fˆ(·) has a finite second moment at
any point with V ar(fˆ(x0 ± δ)) → V ar(fˆ(x0 )) as δ → 0. Suppose we randomly draw δ1 , . . . , δn1 fromN (0, Σ0 ), with Σ0 = Θ(nk ), −1/3 < k < 0. Consider EM-CFD in Algorithm 1. The MSE of θbEM is
1

MSEEM = E[(θbEM − f (1) (x0 ))2 ]
=


  
2
1/3

−1/3
b2)
E B σ
b2 (x0 )/(4B
n2
+ E


σ 2 (x0 )

2/3
2n2




1/3 b2)σb2 (x0 )/(4B

 
!
1/3
−1/3
2
2
b )
+V ar B σb (x0 )/(4Bn2(1 + o(1)),
which is a function of B and σ 2 (x0 ). Moreover,lim

lim n2/3 MSEEM − MSEopt = 0,λ→0 n→+∞

where MSEopt is defined in (1).

(4)

Li and Lam
The M-stage uses the standard CFD scheme, and thus conditional on any given δ we have



√ b
σ 2 (x0 )d(1)2n2 θEM − (f (x0 ) + Bδ ) −−−−−→ N 0,.
n2 →+∞
2δ 2

Proof.

1/6−1/6
b and σb2)n2where both B
b2 (x0 ) are random variables estimated from theWe use δ = σ
b2 (x0 )/(4B
E-stage. Therefore, for the overall EM-CFD scheme, we have
 

h
i
1/3
−1/3
2
(1)
2
2
2
b
b
b
b
E[θEM ] = E E[θEM |B, σb (x0 )] = f (x0 ) + E B σb (x0 )/(4B )n2+ o(δ )
andV ar(θbEM )
h
i


b σ
b σ
= E V ar(θbEM |B,b2 (x0 )) + V ar E[θbEM |B,b2 (x0 )]


 

1/3
2
σ (x0 )

−1/3
4 
2
2
2
b
= E
b (x0 )/(4B )n2+ o(δ ) .

1/3 + o(δ ) + V ar B σ
2/3
2
2
b
2n2
σ
b (x0 )/(4B )Here |o(δ 2 )/δ 2 | → 0 a.s. and |o(δ 4 )/δ 4 | → 0 a.s. when δ → 0 a.s. Then we obtain
MSEEM = E[(θbEM − f (1) (x0 ))2 ]

2


=
E[θbEM − f (1) (x0 )] + V ar θbEM − f (1) (x0 )

  

2
1/3

−1/3
b2)
=
E B σ
b2 (x0 )/(4Bn2+ o(δ 2 )+ Eσ 2 (x0 )
2/3

2n2



b2)
σ
b2 (x0 )/(4B4 
1/3 + o(δ )

 

1/3
−1/3
2
2
2
b
+V ar B σb (x0 )/(4B )n2+ o(δ ) .
We havelim E(o(δ 2 )/δ 2 ) = 0,n2 →+∞

lim E(o(δ 4 )/δ 4 ) = 0, andn2 →+∞

lim V ar(o(δ 2 )/δ 2 ) = 0. With then2 →+∞

results in Theorem 1 and 2, we have, asymptotically, MSEEM satisfies
lim n2/3 MSEEM − 3

n→+∞



σ 2 (x0 )B
4n

2/3


=3

σ 2 (x0 )B
4(1 − λ)

2/3


−3

Further,
lim

lim n

λ→0 n→+∞

Recall that MSEopt = 3

 2

σ (x0 )B
4n

2/3

2/3


MSEEM − 3

σ 2 (x0 )B
4n

2/3
= 0.

(1 + o(1)), which concludes the theorem.

σ 2 (x0 )B
4

2/3
.

Li and Lam
6

NUMERICAL RESULTS

In this section, we conduct numerical experiments to test the performance of the proposed EM-CFD scheme.
We consider three examples. Example 1 aims to estimate the first-order derivative of f (x) = k sin(x) at
the point x0 = 0, where k is an a priori unknown parameter. Observations at δ follow a normal distribution,i.e., fˆ(δ) ∼ N (f (δ), 1). Therefore in this example, we have B = −k/6 and σ 2 (x0 ) = 1 but they are
unknown for EM-CFD scheme. Example 2 concerns the derivative estimation of a polynomial functionfˆ(x) = 1 − 6x + 36x2 − 53x3 + 22x5 + N (0, 0.05) at different points x. We have B = 220x2 − 53 andσ 2 (x) = 0.05 accordingly. Example 3 considers a generic M/M/1 queueing system. Observations are theaveraged system time of the first 10 customers, and the derivative of their expectation with respect to thearrival rate is of interest to us. The true derivative is 0.0946 when both arrival and service rates are 4.
The proposed EM-CFD scheme is compared with an “oracle” CFD that uses the optimal perturbation
size assuming knowledge of B and σ 2 (x0 ) (CFD-opt), and the conventional CFD scheme that uses anarbitrary α:••

CFD-opt: Both B and σ 2 (x0 ) are known for this scheme. Then it outputs n independent runs
1/6 −1/6
of Yi (δ) = (fˆ(x0 + δ) − fˆ(x0 − δ))/(2δ) where δ = σ 2 (x0 )/(4B 2 )n, and estimatesPn(1)bf (x0 ) by θ = (1/n) i=1 Yi (δ).
CFD: This scheme uses αn−1/6 for some arbitrarily chosen α. We do it by arbitrarily choosing
e and σsome value (Be2 (x0 )) for each unknown parameter, and set the perturbation size δ =

1/6




e2)
e σ
e σ
σ
e2 (x0 )/(4Bn−1/6 . We take B,
e2 (x0 ) = (−100/6, 1) in example 1, B,
e2 (x0 ) =e σ(5, 0.05) in example 2 while B,e2 (x0 ) = (1/2, 1) in example 3. For each example, we output
P
n independent runs of Yi (δ), and estimate f (1) (x0 ) by θb = (1/n) n Yi (δ).
i=1

Additionally, in order to examine the performance of the proposed EM-CFD scheme under different
allocation ratio λ, three settings are tested, namely λ = 0.1, λ = 0.2, and λ = 0.3. In all numericalexperiments, the statistical efficiency of the estimation schemes is measured by the MSE estimated by
10,000 independent experimental replications. The MSE is reported as a function of the (unknown) model
parameter k, the value of x, or the number of samples n in each experiment.
In Figure 1, we can see that EM-CFD performs better than the conventional CFD when the modelb andparameter k is less than 60, which could be attributed to the reason that parameter estimation (B2σb (x0 )) in EM-CFD is much closer to optimal than the arbitrary choice in CFD. As the model parameterk increases, the arbitrary choice in CFD approaches the optimal one and is equal to the true value whenk = 100. Of course, a priori we do not know this k, so that the arbitrary CFD can perform well or
can perform poorly. In other words, EM-CFD is more robust. Moreover, the performance of EM-CFD
with λ = 0.1 is still comparable to that of CFD-opt, showing that EM-CFD is close to having the best
possible performance using any CFD for the whole considered range of k. Another observation is that theperformance of EM-CFD gets better as the allocation ratio λ decreases. Although allocating more samples
in the E-stage can lead to a better parameter estimation, it would limit the sample size to estimate thegradient, ultimately resulting in a worse MSE. Such observation is in accord with the asymptotic propertyin (4).
Figure 2 further examines the performance of EM-CFD scheme with lower allocation ratio λ. Similar
to results in Figure 1, EM-CFD remains more efficient than the conventional CFD, and CFD-opt is slightly
better than EM-CFD. However, given a finite number of samples, the performance of EM-CFD becomesworse when decreasing the allocation ratio λ from 0.01 to 0.005. It could be because the number of samples
allocated to estimate the unknown parameters is too few. In other words, we need enough samples toguarantee the consistency of parameter estimation. The best balance in obtaining best overall performanceseems to be obtained at λ = 0.01.

Li and Lam

Figure 1: MSE of all tested CFD schemes in example 1.

Figure 2: MSE of all tested CFD schemes in example 1.

Li and Lam
Figure 3 also demonstrates the robustness of EM-CFD. The MSE of EM-CFD at x ∈ [0, 1] is generally
below 0.5, and EM-CFD performs better than the conventional CFD when x < 0.3 or x > 0.6. In contrast,
the conventional CFD can perform well when its arbitrary choice of B is close to the true value, i.e., x
is around 0.5. Unfortunately, the conventional CFD performs poorly at other times. At point x = 1, theMSE of the conventional CFD is greater than 2.5, which is five times as large as that of EM-CFD.

Figure 3: MSE of all tested CFD schemes in example 2.
Figure 4 illustrates the performance of CFD and EM-CFD in example 3. When n = 1000, we cansee that the proposed EM-CFD achieves MSE = 3.8 × 10−4 whereas the conventional CFD scheme leads
to MSE = 7.7 × 10−4 . That is to say EM-CFD enhances the estimation efficiency by more than 50%
compared to this arbitrarily tuned CFD. In addition, note that unlike examples 1 and 2, simulation varianceσ 2 (x) varies with different points x in this example, which can be viewed as a more challenging setting.
The favorable results in Figure 4 thus demonstrate potential of EM-CFD to apply in more sophisticatedand practical settings.
7

CONCLUSION

While it is well-known how to find perturbation sizes with the optimal order for finite-difference estimatorsin gradient estimation when only noisy function evaluations are available, finding the perturbation size thatmatches the exact first-order MSE, or what we called asymptotic optimality, is generally considered much
more challenging as it relies on unknown model information. In this paper we derived and implemented an
asymptotically optimal CFD, henceforth providing evidence that designing such estimators is practicallypossible. This estimator comprises a two-stage scheme, called EM-CFD, that first estimates the needed
but unknown model parameters and then, based on these estimates, chooses a nearly optimal perturbation
size. Its implementability hinges on the main insight that these model parameter estimates only needmoderate accuracy in order to achieve optimality, thus allowing us to allocate few samples in the estimation
stage. Besides theoretically proving near asymptotic optimality, we conducted some numerical studies toLi and Lam

Figure 4: MSE of CFD and EM in example 3.
demonstrate that the proposed estimator works closely compared to an oracle benchmark and more robustly
compared to CFD with arbitrarily chosen constant in the perturbation size despite the right order. In future
work, we will investigate the multi-dimensional generalizations, investigate the higher-order convergence
rate of the MSE that allows us to further enhance the choice of perturbation size, and conduct more extensive
numerical studies.
ACKNOWLEDGMENTS
This work was supported in part by the National Science Foundation under CAREER CMMI-1834710
and IIS-1849280, and by the China Scholarship Council (CSC) that supports Haidong’s visit at Columbia
University.
REFERENCES
Asmussen, S., and P. W. Glynn. 2007. Stochastic simulation: algorithms and analysis, Volume 57. Springer
Science & Business Media.
Fox, B. L., and P. W. Glynn. 1989. “Replication schemes for limiting expectations”. Probability in theEngineering and Informational Sciences 3(3):299–318.
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.
Glasserman, P. 2013. Monte Carlo methods in financial engineering, Volume 53. Springer Science &Business Media.
Glynn, P. W. 1990. “Likelihood ratio gradient estimation for stochastic systems”. Communications of the
ACM 33(10):75–84.

Li and Lam
Heidelberger, P., X.-R. Cao, M. A. Zazanis, and R. Suri. 1988. “Convergence properties of infinitesimal
perturbation analysis estimates”. Management Science 34(11):1281–1302.
Heidergott, B., and F. Vázquez-Abad. 2008. “Measure-valued differentiation for Markov chains”. Journalof Optimization Theory and Applications 136(2):187–209.
Heidergott, B., F. J. Vázquez-Abad, G. Pflug, and T. Farenhorst-Yuan. 2010. “Gradient estimation for discreteevent systems by measure-valued differentiation”. ACM Transactions on Modeling and Computer
Simulation (TOMACS) 20(1):1–28.
Ho, Y.-C., X. Cao, and C. Cassandras. 1983. “Infinitesimal and finite perturbation analysis for queueing
networks”. Automatica 19(4):439–445.
L’Ecuyer, P. 1991. “An overview of derivative estimation”. In Winter Simulation Conference, 207–217.
Nesterov, Y., and V. Spokoiny. 2017. “Random gradient-free minimization of convex functions”. Foundationsof Computational Mathematics 17(2):527–566.
Peng, Y., M. C. Fu, J.-Q. Hu, and B. Heidergott. 2018. “A new unbiased stochastic derivative estimator for
discontinuous sample performances with structural parameters”. Operations Research 66(2):487–499.
Reiman, M. I., and A. Weiss. 1989. “Sensitivity analysis for simulations via likelihood ratios”. Operations
Research 37(5):830–844.
Rubinstein, R. Y. 1986. “The score function approach for sensitivity analysis of computer simulation
models”. Mathematics and Computers in Simulation 28(5):351–379.
Zazanis, M. A., and R. Suri. 1993. “Convergence rates of finite-difference sensitivity estimates for stochastic
systems”. Operations research 41(4):694–703.
AUTHOR BIOGRAPHIESHAIDONG LI is a Ph.D. candidate in the Department of Industrial Engineering and Management, Peking
University, Beijing, China. He received his B.S. Degree from the Department of Engineering Mechanics
at Peking University. His research interests include simulation optimization and network analysis. He
is a visiting scholar under the supervision of Prof. Henry Lam during 2019-2020. His email address ishaidong.li@pku.edu.cn.
HENRY 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, andstochastic and robust optimization. His email address is henry.lam@columbia.edu.