The Empirical Likelihood Approach to Quantifying Uncertainty in
Sample Average Approximation∗
Henry Lam†

Enlu Zhou‡

Abstract
We study the empirical likelihood approach to construct confidence intervals for the optimalvalue and the optimality gap of a given solution, henceforth quantify the statistical uncertainty
of sample average approximation, for optimization problems with expected value objectives andconstraints where the underlying probability distributions are observed via limited data. Thisapproach relies on two distributionally robust optimization problems posited over the uncertain distribution, with a divergence-based uncertainty set that is suitably calibrated to provide
asymptotic statistical guarantees.
Keywords: empirical likelihood; sample average approximation; confidence interval; constrained
optimization; stochastic program; statistical uncertainty

1

Introduction

We consider a stochastic optimization problem in the form
min{h(x) := E[H(x; ξ)]},
x∈Θ

(1)

where x = (x1 , . . . , xp ) is a continuous decision variable in the deterministic feasible region Θ ⊆ Rp ,and ξ is a random vector on Rd . We are interested in situations where the underlying probability
distribution that controls the expectation E[·] is not fully known and can only be accessed via
limited data ξ1 , . . . , ξn . It is customary in this setting to work on an empirical counterpart of the
problem, namely by solving the sample average approximation (SAA) (e.g., [16]):n1X
H(x; ξi ).
x∈Θ n

min

(2)

i=1

We further consider problems with expected value constraints, in the formminh(x) = E[H(x; ξ)]subject to fk (x) = E[Fk (x; ξ)] ≤ 0, k = 1, . . . , m
gk (x) ≤ 0, k = 1, . . . , s∗(3)

A preliminary conference version of this work has appeared in [11]. Research of the first author was partially
supported by the National Science Foundation under Grants CMMI-1400391/1542020 and CMMI-1436247/1523453.
Research of the second author was partially supported by the National Science Foundation under Grant CAREER
CMMI-1453934, and Air Force Office of Scientific Research under Grant YIP FA- 9550-14-1-0059.
†
Department of Industrial and Operations Engineering, University of Michigan, Ann Arbor, MI. Email:khlam@umich.edu‡
H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA.
Email: enlu.zhou@isye.gatech.edu

12
where gk (·)’s are deterministic functions. Thus (3) can include both stochastic and deterministic
constraints. Again, under limited data ξ1 , . . . , ξn , an SAA version of (3) is (e.g., [17])1 Pn
minn Pi=1 H(x; ξi )subject to n1 ni=1 Fk (x; ξi ) ≤ 0, k = 1, . . . , m(4)
gk (x) ≤ 0, k = 1, . . . , s
Our premise is that beyond the n observations, new samples are not easily accessible because ofeither a lack of data or limited computational capacity in running further Monte Carlo simulation.
The optimal value and solution obtained from (2) or (4) thus deviate from those under the genuine
distribution in (1) or (3). Moreover, the error of the solution implies a non-zero optimality gap withthe true optimal value, resulting in suboptimal decisions. Estimating these errors is important andhas been studied over the years (e.g., [10], [12], Chapter 5 in [16]).
Our main contribution is to bring in a new approach to rigorously quantify the uncertainty in (2)and (4) through constructing confidence intervals (CIs) for the true optimal value and the optimality
gap for a given solution. The machinery underlying our framework uses the so-called empirical
likelihood (EL) method in statistics, and culminates at a reformulation of the problem of finding
the upper and lower bounds of a CI into solving two optimization problems that closely resemble
distributionally robust optimization (DRO). The uncertainty set in the DRO is a divergence-based
ball cast over an uncertain probability distribution, where the size of the ball is suitably calibrated
so that it provides asymptotic guarantees for the coverage probability of the resulting CI.
We study the theory giving rise to such guarantees. We demonstrate through several numericalexamples that our method compares favorably with some existing methods, such as bounds using
the central limit theorem (CLT) and the delta method, in terms of finite-sample performance. In
the remainder of this paper, Sections 2 and 3 study the theory of our approach applied to the
optimal value and the optimality gap, and our online Supplemental Material shows the numericalresults and comparison with previous methods.

2

The Empirical Likelihood Method for Constructing Confidence
Bounds for Optimal Values

This section studies in detail the EL method in constructing CIs for the optimal values. Section 2.1 focuses on (1) that only has deterministic constraints, and Section 2.2 generalizes to thestochastically constrained case (3).

2.1

Deterministically Constrained Optimization

Let us first fix some notations. Given the set of i.i.d. dataP ξ1 , ξ2 , . . . , ξn , we denote a probability
vector over {ξ1 , . . . , ξn } as w = (w1 , . . . , wn ) ∈ Rn , where ni=1 wi = 1 and wi ≥ 0 for all i = 1, . . . , n.
We denote χ2q,β as the 1 − β quantile of a χ2 distribution with degree of freedom q. We use “⇒” todenote convergence in distribution, and “a.s.” to denote “almost surely”.
Our method utilizes the optimization problems
Pnmax / minw minPx∈Θ
i=1 wi H(x; ξi )nsubject to−2 i=1 log(nwi ) ≤ χ2p+1,βPn(5)
i=1 wi = 1
wi ≥ 0 for all i = 1, . . . , n

3where “max / min” denotes a pair of maximization and minimization. Note that the optimal valueof the SAA problem (2) lies
Pnbetween those of (5).
The quantity −(1/n) i=1 log(nwi ) can be interpreted as the Burg-entropy divergence ([15],[3]) between the probability distributions represented by the weights w and by the uniform weights(1/n)i=1,...,n on the support {ξ1 , . . . , ξn }. Thus, the first constraint in (5) is a Burg-entropy divergence ball centered at the uniform weights, with radius χ2p+1,β /(2n). From the viewpointof DRO (e.g., [6, 3, 18]), the optimization problems in (5) output the worst-case estimates ofminx∈Θ {h(x) = E[H(x; ξ)]} when E[·] is uncertain and its underlying distribution is believed tolie inside the divergence ball. We should point out, however, that this DRO interpretation differsfrom those in the existing literature (e.g., [4]), as our divergence ball (i.e. the “uncertainty set” inthe terminology of robust optimization) may have low coverage of the true distribution P . Thiscan be seen particularly when P is a continuous distribution, in which case the coverage of thedivergence ball is zero because of the violation of the absolute continuity requirement needed inproperly defining the divergence.
The EL method is a mechanism to endow statistical meaning to (5). In particular, it assertsthat using the ball size χ2p+1,β /(2n) in (5) gives rise to statistically valid 1 − β confidence bounds
for the optimal value of (1) (despite that the ball may under-cover the true distribution). Thismethod originates as a nonparametric analog of maximum likelihood estimationQfirst proposed by[13]. On the data set {ξ1 , . . . , ξn }, we first define a “nonparametric likelihood” ni=1 wi , where wi
is aQprobability weight applied to each datum. It is straightforward
to see that the maximum valueQnnof i=1 wi , among all w in the probability simplex, is i=1 (1/n). In fact, the same conclusion
holds even if one
Q the support of the data, which could only makeQ allows putting weights outside
the likelihood ni=1 wi smaller. In this sense, ni=1 (1/n) can be viewed as a maximum likelihood in
the nonparametric space. Correspondingly, we define thelikelihood
ratio between
Q nonparametricQQ
the weights w and the maximum likelihood weights as ni=1 wi / ni=1 (1/n) = ni=1 (nwi ).
The key of the EL method is a nonparametric counterpart of the celebrated Wilks’ Theorem[19] in parametric likelihood inference. The latter states that the ratio between the maximum
likelihood and the true likelihood (the parametric likelihood ratio) converges to a χ2 -distribution
in a suitable logarithmic scale. To develop this analog, we first incorporate a target parameter ofinterest, i.e. the quantity whose statistical uncertainty is to be assessed (or to be “estimated”).
Say this parameter is θ ∈ Rp . Suppose the true parameter is known to satisfy the set of equations
E[t(θ; ξ)] = 0 where E[·] is the expectation for the random object ξ ∈ Rd , and t(θ; ξ), 0 ∈ Rb . We
define the nonparametric profile likelihood ratio as
( n
)
n
n
Y
X
X
R(θ) = maxnwi :wi t(θ; ξi ) = 0,
wi = 1, wi ≥ 0 for all i = 1, . . . , n
(6)
i=1

i=1

i=1

where profiling refers to the categorization of all weights that respect the set of equations E[t(θ; ξ)] =0.
With the above definitions, the crux is the empirical likelihood theorem (ELT):
Theorem 1 (Theorem 3.4 in [14]). Let ξ1 , . . . , ξn ∈ Rd be i.i.d. data. Let θ0 ∈ Rp be a valueof the parameter that satisfies E[t(θ; ξ)] = 0, where t(θ; ξ), 0 ∈ Rb . Assume the covariance matrix
V ar(t(θ0 ; ξ)) is finite and has rank q > 0. Then −2 log R(θ0 ) ⇒ χ2q , where R(θ) is defined in (6).
The quantity −2 log R(θ) is defined as ∞ if the optimization in (6) is infeasible.
We now explain how (5) provides confidence bounds for optimization problem (1). We makethe following assumptions:4
Assumption 1.

1. h(x) is differentiable in x with ∇x h(x) = E[∇x H(x; ξ)] for all x ∈ Θ.

2. x∗ ∈ argminx∈Θ h(x) if and only if ∇x h(x∗ ) = 0. Moreover, this relation is distributionally stable, meaning that x̃∗ ∈ argminx∈Θ h̃(x) if and only if ∇x h̃(x̃∗ ) = 0 for any h̃(x) = Ẽ[H(x; ξ)]that has the expectation Ẽ[·] generated under an arbitrary distribution P̃ such thatsup |h̃(x) − h(x)| < x∈Θ

for small enough  > 0.
3. There exists an x∗ ∈ argminx∈Θ h(x) such that the covariance matrix of the random vector(∇x H(x∗ ; ξ), H(x∗ ; ξ)) ∈ Rp+1 is finite and has positive rank.
P
4. n1 ni=1 H(x; ξi ) → h(x) uniformly over x ∈ Θ a.s..

5. E supx∈Θ H(x; ξ)2 < ∞Assuming the existence of ∇x H(x; ξ) a.s., the interchangeability of derivative and expectationin Assumption 1.1 can generally be justified by the pathwise Lipschitz continuity condition|H((x1 , . . . , xj−1 , u, xj+1 , . . . , xp ); ξ) − H((x1 , . . . , xj−1 , v, xj+1 , . . . , xp ); ξ)| ≤ Mj |u − v| a.s.
for any u, v in a nonrandom neighborhood around the point xj to be differentiated and Mj measurable with EMj < ∞ (e.g., [1]). Another sufficient condition is that H(x; ξ) is a.s. continuous andpiecewise differentiable in xj and supu∈D |(∂/∂xj )H((x1 , . . . , xj−1 , u, xj+1 , . . . , xp ); ξ)| is integrable
where D is a neighborhood around xj [8]. Assumption 1.2 states that the first order condition
for optimality is both sufficient and necessary. Assumptions 1.2 and 1.4 together ensure that thisfirst order condition is unchanged when the true distribution is replaced by a (weighted) empiricalversion as the sample size gets large. Assumption 1.5 is a technical condition required to boundthe error between the empirical distribution and its weighted version within the divergence ball.
Assumption 1.3 is used to invoke Theorem 1. Note that x∗ is not necessarily unique.
As our subsequent development will reveal, both the necessity and the sufficiency of the firstorder condition in Assumption 1.2 are required; in particular, we need the necessity of ∇x h(x∗ ) = 0
for x∗ ∈ argminx∈Θ h(x) and the sufficiency of ∇x h̃(x̃∗ ) = 0 for x̃∗ ∈ argminx∈Θ h̃(x) in order for
our argument on statistical guarantee to go through. Assumptions 1.2, 1.4 and 1.5 can be replaced
by a single condition
P
AssumptionP
2. x∗ ∈ argminx∈Θ h(x) if and only if ∇x h(x∗ ) = 0, and x̃∗ ∈ argminx∈Θ ni=1 wi H(x; ξi )if and only if ni=1 wi ∇x H(x̃∗ ; ξi ) = 0 for any support set {ξ1 , . . . , ξn } ⊂ Θ and arbitrary probability
vector w.
Assumption 2 is satisfied by, for instance, H(·; ξ) that is coersive and convex for any ξ andΘ = Rp .
We have the following statistical guarantee:
Theorem 2. Suppose ξ1 , . . . , ξn ∈ Rd are i.i.d. data. Let z ∗ be the optimal value of (1), and z andz be the maximum and minimum values of (5) respectively. Then, under Assumption 1, we have
lim inf P (z ∗ ∈ [z, z]) ≥ 1 − β.
n→∞5Proof. By Assumption 1.3, there exists an x∗ ∈ argminx∈Θ h(x) such that the covariance matrix ofthe random vector (∇x H(x∗ ; ξ), H(x∗ ; ξ)) has a positive rank, call it r. Also, by Assumptions 1.1
and 1.2, x∗ satisfies E[∇x H(x∗ ; ξ)] = 0.
We define the nonparametric profile likelihood ratio as
( n
)
n
n
n
Y
X
X
X
R(x, z) = max
nwi :wi ∇x H(x; ξi ) = 0,
wi H(x; ξi ) = z,wi = 1, wi ≥ 0 for all i = 1, . . . , n .
i=1

i=1

i=1

i=1

(7)
Let z ∗ = minx∈Θ h(x) = h(x∗ ). From Theorem 1, the nonparametric profile likelihood ratio (7)satisfies −2 log R(x∗ , z ∗ ) ⇒ χ2r as n → ∞. This implies P (−2 log R(x∗ , z ∗ ) ≤ χ2r,β ) → 1 − β.
The rest of the proof focuses on the event −2 log R(x∗ , z ∗ ) ≤ χ2r,β . Write−2 log R(x∗ , z ∗ )
(
n
n
n
n
X
X
X
X
= min − 2log(nwi ) :wi ∇x H(x∗ ; ξi ) = 0,
wi H(x∗ ; ξi ) = z ∗ ,wi = 1, wi ≥ 0
i=1

i=1

i=1

i=1

)
for all i = 1, . . . , n(8)

We argue that −2 log R(x∗ , z ∗ ) ≤ χ2r,β implies the existence of a probability vector w such thatnX
i=1

wi ∇x H(x∗ ; ξi ) = 0,

n
X
i=1

wi H(x∗ ; ξi ) = z ∗ , −2nXlog(nwi ) ≤ χ2r,β

(9)

i=1

P
Notice that −2 ni=1 log(nwi ) = ∞ if wi = 0 for any i. Hence it suffices to replace, in (8), wi ≥ 0
with
Pwi ≥  for all i, for some small enough  > 0. In this modified, compact, feasible set,
−2 ni=1 log(nwi ) is bounded and hence must possess an optimal solution w, which is a probability
vector that satisfies (9).
This further implies that z ∗ is bounded from above and below by the optimization problemsPmax / minw Pni=1 wi H(x∗ ; ξi )n∗subject toi=1
Pnwi ∇x H(x ; ξi ) =2 0
−2 i=1 log(nwi ) ≤ χr,β
(10)
Pn
w
=
1
i=1 i
wi ≥ 0 for all i = 1, . . . , n
We argue that as n → ∞, (10) is equivalent toPn∗max / minw
i=1 wi H(x ; ξi )P
subject tow ∈P{(w1 , . . . , wn ) : x∗ ∈ argminx∈Θ ni=1 wi H(x; ξi )}−2 ni=1 log(nwi ) ≤ χ2r,βPn
i=1 wi = 1
wi ≥ 0 for all i = 1, . . . , n(11)eventually (i.e. with probability 1). Note that the first constraint in (11) states that the probability∗
vector w must be
Pnchosen such that x , a minimizer of h(x) picked at the beginning of this proof,
also minimizes i=1 wi H(x; ξi ).

6
To develop the argument for the asymptotic equivalence, let us denote P w as the distribution
represented by the probability weights w on the support {ξ1 , . . . , ξn }. Denote E w [·] as the associatedexpectation and hw (x) = E w [H(x; ξ)]. We will show that|hw (x) − h(x)| → 0 a.s.

sup

(12)

x∈Θ,w∈Wr

where
(
Wr =

(w1 , . . . , wn ) ∈ Rn : −2nXlog(nwi ) ≤ χ2r,β ,

i=1

n
X

)
wi = 1, wi ≥ 0 for all i = 1, . . . , n

(13)

i=1

Pn
∗
Assumption 1.2 then implies thatwithprobability
1,for
sufficiently
largen,i=1 wi ∇x H(x ; ξi ) = 0
P
if and only if x∗ ∈ argminx∈Θ ni=1 wi H(x; ξi ) for any w ∈ Wr , leading to the equivalence.
We now show (12). Consider
Zw
sup |h (x) − h(x)| =supH(x; ξ)d(P w − P )(ξ)
x∈Θ,w∈Wrx∈Θ,w∈Wr

Z
≤

sup

Z

w

H(x; ξ)d(P − P̂ )(ξ) + supx∈Θ,w∈Wr

H(x; ξ)d(P̂ − P )(ξ)

x∈Θwhere P̂ denotes the empirical distribution generated from {ξ1 , . . . , ξn }Zw
sup
|H(x; ξi )|dT V (P , P̂ ) + sup
H(x; ξ)d(P̂ − P )(ξ) (14)≤x∈Θ,1≤i≤n,w∈Wrx∈Θ

where dT V denotes the total variation distance
Now by Lemma 11.5 in [14] (restated in the Appendix) and Assumption 1.5, we havesupx∈Θ,1≤i≤n|H(x; ξi )| = max sup |H(x; ξi )| = o(n1/2 ) a.s.
1≤i≤n x∈ΘOn the other hand, by Pinsker’s inequality, for any w ∈ Wr ,sr Pnr
− i=1 log(nwi )dKL (P w , P̂ )
χ2r
w
=
≤
dT V (P , P̂ ) ≤
2
2n
4n

(15)

(16)

where dKL denotes the Kullback-Leibler divergence. Combining (15) and (16), the first term in(14) goes to 0 a.s.. The second term in (14) converges to 0 a.s. by Assumption 1.4. Hencesupx∈Θ |hw (x) − h(x)| → 0 a.s.. Therefore (10) is equivalent to (11) eventually as n → ∞.P
Consider (11). With the first constraint, the objective function must be equal to minx∈Θ ni=1 wi H(x; ξi ).
Thus (11) is equivalent toP
max / minw minx∈Θ ni=1 wi H(x; ξi )P
subject tow ∈P{(w1 , . . . , wn ) : x∗ ∈ argminx∈Θ ni=1 wi H(x; ξi )}−2 ni=1 log(nwi ) ≤ χ2r,β(17)Pn
i=1 wi = 1
wi ≥ 0 for all i = 1, . . . , n
Let v and v be the maximum and minimum values of (17). Note that r ≤ p + 1 since r is the rankof a R(p+1)×(p+1) matrix. This implies χ2r,β ≤ χ2p+1,β . Together with a relaxation by removing the

7first constraint in (17), we have v ≥ z and v ≤ z where z and z are the maximum and minimum
values of (5). From this we conclude that
lim inf P (z ≤ z ∗ ≤ z) ≥ lim inf P (v ≤ z ∗ ≤ v) ≥ lim inf P (−2 log R(x∗ , z ∗ ) ≤ χ2r,β ) = 1 − β
n→∞

n→∞

n→∞

We also have the following result based on Assumption 2 that can be more natural to verify insome cases:
Corollary 1. Suppose ξ1 , . . . , ξn ∈ Rd are i.i.d. data. Let z ∗ be the optimal value of (1), and z
and z be the maximum and minimum values of (5) respectively. Then, under Assumptions 1.1, 1.3
and 2, we have
lim inf P (z ∗ ∈ [z, z]) ≥ 1 − β.
n→∞Proof. The proof follows similarly as that of Theorem 2, with the observation that the equivalenceof (10) and (11) holds for all n if Assumption 2 replaces Assumptions 1.2, 1.4 and 1.5.
Note that we have obtained bounds by relaxing the constraints in (17), and the degree of freedom
in the χ2 -distribution may not be optimally chosen. Nevertheless, our numerical examples show
that, at least for small p, the EL method provides reasonably tight CIs. There exist techniques(e.g., bootstrap calibration or Bartlett correction; [13, 7]) that can improve the coverage of the EL
method in estimation problems. Investigation of these techniques in the optimization context isdelegated to future work.

2.2

Stochastically Constrained Optimization

We generalize the EL method to the stochastically constrained problem (3). In this setting, weconstruct CI via the following optimization problems


Pn
w
H(x;
ξ
)

 minxiii=1
Pnsubject tomax / minw
i=1 wi Fk (x; ξi ) ≤ 0, k = 1, . . . , m g(x)≤ 0, k = 1, . . . , skP(18)
subject to−2 ni=1 log(nwi ) ≤ χ2p+m+1,βPn
i=1 wi = 1
wi ≥ 0 for all i = 1, . . . , n
While resembling (5), we note that the degree of freedom in the χ2 -distribution is now p + m + 1,
which includes the number of stochastic constraints compared to (5).
For convenience, we denoteΛ = {x ∈ Rp : gk (x) ≤ 0, k = 1, . . . , s}as the set of x satisfying the deterministic constraints in (3).
We make the following assumptions in parallel to Assumption 1:
Assumption 3. We assume:8
1. h(x) = E[H(x; ξ)], fk (x) = E[Fk (x; ξ)], k = 1, . . . , m and gk (x), k = 1, . . . , s are all differentiable in x ∈ Λ, and∇x h(x) = E[∇x H(x; ξ)], ∇x fk (x) = E[∇x Fk (x; ξ)]
2. Let S ∗ be the set of all optimal solutions for (3). x∗ ∈ S ∗ if and only if x∗ satisfies the KKT
condition, where the active set of the KKT condition (i.e. equalities) is unique among allx∗ ∈ S ∗ and is sufficient for determining S ∗ . This relation is distributionally stable, meaningthat x̃∗ ∈ S̃ ∗ , where S̃ ∗ is the set of optimal solutions for
min
h̃(x)
subject to f˜k (x) ≤ 0, k = 1, . . . , mg̃k (x) ≤ 0, k = 1, . . . , s(19)

if and only if x̃∗ satisfies the corresponding KKT condition, where h̃(x) = Ẽ[H(x; ξ)] andf˜k (x) = Ẽ[Fk (x; ξ)], with Ẽ denoting the expectation under an arbitrary distribution P̃ suchthat
sup |h̃(x) − h(x)| < x∈Λ

sup |f˜k (x) − fk (x)| <  for all k = 1, . . . , mx∈Λ

for small enough  > 0. Moreover, for any such  > 0, the active set of the KKT condition
at any x̃∗ ∈ S̃ ∗ for (19) is the same as that at any x∗ ∈ S ∗ for (3) and is sufficient for
determining S̃ ∗ .
3. There exists an optimal solution x∗ for (3), with associated Lagrange multipliers for the∗ = (λ∗ , . . . , λ∗ ), such that the covariance matrix
stochastic constraints in (3) given by λP
m
1
∂
∗
∗ ∂
∗
∗
∗
of the variables H(x ; ξ), ∂xj H(x ; ξ) + m
k=1 λk ∂xj Fk (x ; ξ), and Fk (x ; ξ), for all indices jand k corresponding to the active set of the KKT condition, is finite and has positive rank.
4.

n

n

i=1

i=1

1X
1X
H(x; ξi ) → h(x) andFk (x; ξi ) → fk (x), k = 1, . . . , mnn
uniformly over x ∈ Λ a.s..

5. E supx∈Θ H(x; ξ)2 < ∞ and E supx∈Θ Fk (x; ξ)2 < ∞ for k = 1, . . . , m.
Denote ν ∗ = (ν1∗ , . . . , νs∗ ) as the Lagrange multiplier for the deterministic constraints in (3). In
Assumptions 3.2 and 3.3 above, the active set of the KKT condition satisfied by (x∗ , λ∗ , ν ∗ ) is in
the form
m

s

k=1

k=1

X
X
∂
∂
∂
h(x∗ ) +λ∗kfk (x∗ ) +νk∗gk (x∗ ) = 0, j ∈ A∗1 ≡ {1, . . . , p}
∂xj
∂xj
∂xj
fk (x∗ ) = 0, k ∈ A∗2 ⊂ {1, . . . , m}
gk (x∗ ) = 0, k ∈ A∗3 ⊂ {1, . . . , s}λ∗k = 0 k ∈ {1, . . . , m} \ A∗29νk∗ = 0, k ∈ {1, . . . , s} \ A∗3
where A∗1 , A∗2 and A∗3 denote the sets of indices that correspond to the equalities in the optimality
condition, which are unique among any optimal solutions of (3) by Assumption 3.2. The j and k
described in Assumption 3.3 refer to the indices in A∗1 and A∗2 . Assumption 3.2 further enforces
the sets A∗1 , A∗2 and A∗3 to remain as the active sets under a perturbation to P̃ described therein,
and the equalities indexed via these sets are enough to determine S ∗ and S̃ ∗ . Assumptions 3.2
and 3.3 generalize Assumptions 1.2 and 1.3 from the simple zero-derivative optimality condition to
the KKT condition. Similar to Section 2.1, we require the necessity of the KKT and the active
set conditions regarding (3) and the sufficiency regarding (19) for our development to go through.
Constraint qualification for the validity of the KKT condition is implicitly assumed in Assumption
3.2.
We have the following result:
Theorem 3. Suppose ξ1 , . . . , ξn are i.i.d. data. Under Assumption 3, we have
lim inf P (z ∗ ∈ [z, z]) ≥ 1 − βn→∞where z ∗ is the optimal value of (3), and z and z are the minimum and maximum values of (18).
Proof. Consider the nonparametric profile likelihood ratioPnξi ) = z
i=1 wi H(x;




Pm
Ps
Pn




∂
∂
∂
∗


H(x;
ξ
)
+
λ
F
(x;
ξ
)
+
ν
g
(x)
=
0,
j
∈
A
w
i
i

1 
k=1 k ∂xj k
k=1 k ∂xj k
i=1 i ∂xj




P


n
∗


w
F
(x;
ξ
)
=
0,
k
∈
A


i
2
i=1 i k
n

Y
∗
g
(x)
=
0,
k
∈
A
k
3
R(x, λ, ν, z) = max
nwi :


λk = 0 k ∈ {1, . . . , m} \ A∗2


i=1




∗


ν
=
0,
k
∈
{1,
.
.
.
,
s}
\
A


k
3


P


n




w
=
1
i


i=1


wi ≥ 0 for all i = 1, . . . , n
(20)
∗
∗
∗
∗
Let x be an optimal solution for (3) satisfying Assumption 3.3, and λ = (λ1 , . . . , λm ), ν ∗ =(ν1∗ , . . . , νs∗ ) be its associated Lagrange multipliers. By Assumption 3.3, the covariance of the random
vector concatenated byH(x∗ ; ξ)
m

s

k=1

k=1

X
X
∂
∂
∂
H(x∗ ; ξ) +λ∗kFk (x∗ ; ξ) +νk∗gk (x∗ ) for j ∈ A∗1
∂xj
∂xj
∂xj
∗

Fk (x ; ξ)for k ∈ A∗2

has rank r for some r > 0. Let z ∗ be the optimal value of (3) equal to h(x∗ ). Since the other active
KKT conditions are deterministic, Theorem 1 implies that −2 log R(x∗ , λ∗ , ν ∗ , z ∗ ) ⇒ χ2r , whichfurther implies P (−2 log R(x∗ , λ∗ , ν ∗ , z ∗ ) ≤ χ2r,β ) → 1 − β.
Similar to the proof of Theorem 2, −2 log R(x∗ , λ∗ , ν ∗ , z ∗ ) ≤ χ2r,β implies the existence of a w
Pthat satisfies −2 ni=1 log(nwi ) ≤ χ2r,β and all constraints in (20) evaluated at x∗ , λ∗ , ν ∗ , z ∗ . This in10turn implies that z ∗ is bounded byPn∗; ξ )
max / minwwi H(x
i

 P
Pi=1
Pm
n
s
∂
∗ ∂
∗; ξ ) +∗ ∂ F (x∗ ; ξ ) +∗∗subject to
w
H(x
λ
i
i
i
k
i=1
k=1 k ∂xj
k=1 νk ∂xj gk (x ) = 0, j ∈ A1
∂xj
Pn
∗
∗
i=1 wi Fk (x ; ξi ) = 0, k ∈ A2∗∗
gk (x ) = 0, k ∈ A3
λ∗k = 0 k ∈ {1, . . . , m} \ A∗2νk∗ =k ∈ {1, . . . , s} \ A∗3P0,n−2 i=1 log(nwi ) ≤ χ2r,βPn
i=1 wi = 1
wi ≥ 0 for all i = 1, . . . , n(21)Using the same argument as in the proof of Theorem 2, we obtain from Assumptions 3.4 and 3.5
that
sup |hw (x) − h(x)| → 0 a.s.
x∈Λ,w∈Wrsupx∈Λ,w∈Wr|fkw (x) − fk (x)| → 0 a.s. for all k = 1, . . . , mwhere Wr is defined in (13), and hw (x) = E w [H(x; ξ)], fkw (x) = E w [Fk (x; ξ)] with E w denoting the
expectation with respect to P w , the probability distribution represented by the weights w on the
support {ξ1 , . . . , ξn }. Thus, by Assumption 3.2, the set of active KKT conditions for an optimal
solution of the weighted sample problem
Pnminx
Pni=1 wi H(x; ξi )subject to
i=1 wi Fk (x; ξi ) ≤ 0, k = 1, . . . , m
gk (x) ≤ 0, k = 1, . . . , sfor any w ∈ Wr is identical to that for x∗ for (3) eventually as n → ∞, and Assumption 3.2 further
implies that (21) is equivalent toPn∗max / minw
i=1wi H(x ; ξi )


Pn

 argminx

i=1 wi H(x; ξi )Pnsubject to
w
F
(x;
ξ
)
≤
0,
k
=
1,
.
.
.
,
m
subject tow ∈ (w1 , . . . , wn ) : x∗ ∈i
i=1 i k



g
(x)
≤
0,
k
=
1, . . . , skPn−2 i=1 log(nwi ) ≤ χ2r,βPn
i=1 wi = 1
wi ≥ 0 for all i = 1, . . . , n(22)eventually
P as n → ∞. With
P the first constraint, the objective function in (22) must be equal tominx { ni=1 wi H(x; ξi ) : ni=1 wi Fk (x; ξi ) ≤ 0, k = 1, . . . , m, gk (x) ≤ 0, k = 1, . . . , s}. Note thatr ≤ 1 + |A∗1 | + |A∗2 | ≤ 1 + p + m(23)

where | · | denotes cardinality. This implies that χ2r,β ≤ χ2p+m+1,β . Thus, together with a relaxationof the first constraint in (22), the same argument as in the proof of Theorem 2 stipulates that the
maximum and minimum values of (22) are bounded from above and below respectively by those of(18) and concludes the theorem.

11
Note that, much like the proof of Theorem 2, we have relaxed constraints and placed a conservative bound on the degree of freedom of the χ2 -distribution in (23), which could potentially be
improved with more refined analysis.

3

The Empirical Likelihood Method for Constructing Confidence
Bounds for Optimality Gaps

We study the construction of CI for the optimality gap of a given solution using the EL method. We
focus on optimization problem (1) and suppose x̂ is a feasible solution obtained from some procedure
independently of the data ξ1 , . . . , ξn . The optimality gap of x̂ is given by G(x̂) = h(x̂) − z ∗ wherez ∗ is the optimal value of (1). We will show how we can apply the results in Section 2 to find the
CI for G(x̂).
Consider the optimization problems
Pnmax / minw maxPx∈Θ
i=1 wi [H(x̂; ξi ) − H(x; ξi )]subject to−2 ni=1 log(nwi ) ≤ χ2p+1,βPn(24)
i=1 wi = 1
wi ≥ 0 for all i = 1, . . . , n
We have the following guarantee in using (24) to construct the CI for G(x̂) for (1):Theorem 4. Suppose ξ1 , . . . , ξn ∈ Rd are i.i.d. data independent of a given solution x̂ that is feasiblefor (1). Let G(x̂) be the optimality gap of x̂ for (1), and z and z be the maximum and minimum
values of the programs in (24) respectively. Suppose Assumption 1 holds except that in Condition
2, the relation holds for any h̃(x) = Ẽ[H(x; ξ)] such that supx∈Θ |(h̃(x) − h(x)) − (h̃(x̂) − h(x̂))| < and in Condition 3, we consider the covariance matrix of (∇x H(x∗ ; ξ), H(x∗ ; ξ) − H(x̂; ξ)) instead.
We have
lim inf P (G(x̂) ∈ [z, z]) ≥ 1 − β.
(25)
n→∞

Proof. Let H̄(x; ξ) = H(x; ξ) − H(x̂; ξ), and h̄(x) = E[H̄(x; ξ)] = h(x) − h(x̂). We verify that
Assumption 1, with the change that supx∈Θ |(h̃(x) − h(x)) − (h̃(x̂) − h(x̂))| <  is used in Condition
2 and the covariance matrix of (∇x H(x∗ ; ξ), H(x∗ ; ξ) − H(x̂; ξ)) is considered instead in Condition
3, implies that h̄ and H̄ satisfies Assumption 1 too with h and H replaced by h̄ and H̄.
Condition 1: We have ∇x h̄(x) = ∇x (h(x) − h(x̂)) = ∇x h(x) = E[∇x H(x; ξ)] = E[∇x (H(x; ξ) −H(x̂; ξ))] = E[∇x H̄(x; ξ)].
Condition 2: We have x∗ ∈ argminx∈Θ h̄(x) ⇔ x∗ ∈ argminx∈Θ h(x) ⇔ ∇x h(x∗ ) = 0 ⇔ ∇x h̄(x∗ ) =¯¯
0. Similarly, x̃∗ ∈ argminx∈Θ h̃(x) ⇔ x̃∗ ∈ argminx∈Θ h̃(x) ⇔ ∇x h̃(x̃∗ ) = 0 ⇔ ∇x h̃(x̃∗ ) = 0 for any¯¯h̃(x) = h̃(x) − h̃(x̂) that satisfies supx∈Θ |h̃(x) − h̄(x)| <  by our modification of this condition.
Condition 3: By our modification of this condition we have the covariance of (∇x H̄(x∗ ; ξ), H̄(x∗ ; ξ)) =(∇x H(x∗ ; ξ), H(x∗ ; ξ) − H(x̂; ξ)) finite and having
P a positive rank.P
P
Condition 4: It is straightforward to show that n1 ni=1 H̄(x; ξi ) = n1 ni=1 H(x; ξi )− n1 ni=1 H(x̂; ξi ) →h̄(x) a.s. uniformly over x ∈ Θ.
Condition 5: We have E[supx∈Θ H̄(x; ξ)2 ] = E[supx∈Θ (H(x; ξ)−H(x̂; ξ))2 ] ≤ 4(E[supx∈Θ H(x; ξ)2 ]+E[H(x̂; ξ)2 ]) < ∞.
We have therefore verified our claim. Using Theorem 2, we get that
lim inf P (v ≤ h̄(x∗ ) ≤ v) ≥ 1 − βn→∞12where v and v are the maximum and minimum values ofPn
max / minw minPx∈Θ
i=1 wi H̄(x; ξi )subject to−2 ni=1 log(nwi ) ≤ χ2p+1,βPn
i=1 wi = 1
wi ≥ 0 for all i = 1, . . . , n(26)

Noting that G(x∗ ) = −h̄(x∗ ), we get (25) immediately.
In parallel to Corollary 1, we have the following result on optimality gap assessment based on
an alternative set of assumptions:
Corollary 2. Suppose ξ1 , . . . , ξn ∈ Rd are i.i.d. data independent of a given solution x̂ that is
feasible for (1). Let G(x̂) be the optimality gap of x̂ for (1), and z and z be the maximum andminimum values of the programs in (24) respectively. Under Assumption 1.1, Assumption 1.3 where
the covariance matrix of (∇x H(x∗ ; ξ), H(x∗ ; ξ) − H(x̂; ξ)) is considered instead, and Assumption 2,
we have
lim inf P (G(x̂) ∈ [z, z]) ≥ 1 − βn→∞Proof. Similar to the proof of Theorem 4, we can define H̄ and h̄ and verify that Assumption
1.1 and the modified Assumption 1.3 hold for H̄ and h̄. Assumption 2 is also satisfied for h̄∗ ) = 0 ⇔ ∇ h̄(x∗ ) = 0
h(x) ⇔ ∇x h(xPand H̄ because x∗ ∈P
argminx∈Θ h̄(x) ⇔ x∗ ∈ argminx∈Θ
x
P
n
n
n
∗) = 0 ⇔∗ ∈ argmin∗ ∈ argmin
w
∇
H(x̃
w
H(x)
⇔
w
H̄(x)
⇔
x̃
and
x̃
i
x
i
i
x∈Θ
x∈Θ
i=1
i=1
i=1
Pn
∗
i=1 wi ∇x H̄(x̃ ) = 0 for any support set {ξ1 , . . . , ξn } and arbitrary probability weight vector w.
The rest of the proof then follows as that of Theorem 4.
The statistical uncertainty of the optimality gap of the stochastically constrained problem (3)can in principle be analyzed in a similar fashion. However, validating the feasibility of a given
solution x̂ under limited data is not straightforward in such a scenario (see, e.g., Section 3 in [17]),and for this reason we skip the corresponding result in this work.

4

Conclusion

We have studied the EL method to construct statistically valid CIs for the optimal value andthe optimality gap of a given solution for stochastic optimization problems. The method builds on
positing two optimization problems that resemble DRO problems with Burg-entropy divergence ballconstraints, with the ball size suitably calibrated by a χ2 -quantile with a chosen degree of freedom.
We have studied the theory leading to the statistical guarantees and numerically compared our
method to approaches suggested by the CLT (in our online Supplemental Material). Built on
a rigorous foundation, our method provides a competitive method for evaluating the statisticaluncertainty for stochastic optimization problems under limited data. In future work, we plan tofurther refine the accuracy and extend the scope of our method.

Appendix
Lemma 1 (Lemma 11.2 in [14]). Let Yi be i.i.d. random variables in R with EYi2 < ∞. We havemax1≤i≤n |Yi | = o(n1/2 ) a.s..

13

References[1] S. Asmussen and P. W. Glynn. Stochastic Simulation: Algorithms and Analysis, volume 57.
Springer Science & Business Media, 2007.
[2] G. Bayraksan and D. Morton. Assessing solution quality in stochastic programs. Mathematical
Programming, 108:495–514, 2006.
[3] A. Ben-Tal, D. den Hertog, A. De Waegenaere, B. Melenberg, and G. Rennen. Robust solutionsof optimization problems affected by uncertain probabilities. Management Science, 59(2):341–
357, 2013.
[4] D. Bertsimas, V. Gupta, and N. Kallus. Robust sample average approximation. available atarXiv:1408.4445, 2016.
[5] W. Chen, M. Sim, J. Sun, and C.-P. Teo. From CVaR to uncertainty set: Implications in joint
chance-constrained optimization. Operations Research, 58(2):470–485, 2010.
[6] E. Delage and Y. Ye. Distributionally robust optimization under moment uncertainty withapplication to data-driven problems. Operations Research, 58(3):595–612, 2010.
[7] T. DiCiccio, P. Hall, and J. Romano. Empirical likelihood is Bartlett-correctable. The Annalsof Statistics, 19(2):1053–1061, 1991.
[8] P. Glasserman. Performance continuity and differentiability in Monte Carlo optimization. In
Proceedings of the 20th Winter Simulation Conference, pages 518–524, New York, NY, USA,1988. ACM.
[9] R. Jiang and Y. Guan. Data-driven chance constrained stochastic program. Mathematical
Programming, Series A, 158(1):291–327, 2012.
[10] A. J. Kleywegt, A. Shapiro, and T. Homem-de-Mello. The sample average approximation
method for stochastic discrete optimization. SIAM Journal on Optimization, 12(2):479–502,
2002.
[11] H. Lam and E. Zhou. Quantifying uncertainty in sample average approximation. In Proceedingsof the 2015 Winter Simulation Conference, pages 3846–3857. IEEE Press, 2015.
[12] W.-K. Mak, D. Morton, and R. Wood. Monte Carlo bounding techniques for determining
solution quality. Operations Research Letters, 24:47–56, 1999.
[13] A. B. Owen. Empirical likelihood ratio confidence intervals for a single functional. Biometrika,
75(2):237–249, 1988.
[14] A. B. Owen. Empirical Likelihood. CRC press, 2001.
[15] L. Pardo. Statistical Inference Based on Divergence Measures. CRC Press, 2005.
[16] A. Shapiro, D. Dentcheva, and A. Ruszczyński. Lectures on Stochastic Programming: Modelingand Theory, volume 16. SIAM, 2014.
[17] W. Wang and S. Ahmed. Sample average approximation of expected value constrained stochastic programs. Operations Research Letters, 36(5):515–519, 2008.

14[18] W. Wiesemann, D. Kuhn, and M. Sim. Distributionally robust convex optimization. Operations
Research, 62(6):1358–1376, 2014.
[19] S. S. Wilks. The large-sample distribution of the likelihood ratio for testing composite hypotheses. The Annals of Mathematical Statistics, 9(1):60–62, 1938.
[20] S. Zymler, D. Kuhn, and B. Rustem. Distributionally robust joint chance constraints withsecond-order moment information. Mathematical Programming, Series A, 137(1):167–198,
2013.

15

Supplemental Material
5

Numerical Examples

We test the presented method numerically on three examples. For proof of concept, the first example
is a simple unconstrained quadratic optimization problem. Then we apply the proposed method to
two more examples, including the problem of estimating Conditional-Value-at-Risk (CVaR) and astochastically constrained portfolio optimization problem. The latter examples strictly speaking do
not satisfy our assumptions, since the needed optimality conditions may not hold for their sample
counterparts (Assumption 1.2, 2 or 3.2). However, given that the EL method does not rely on these
conditions procedurally, we can still test its performance on these examples.
We compare EL with the CIs obtained from the CLT and the delta method ([16], Theorem 5.7).
For deterministically constrained problems in the form (1), the (1 − β) CI on the optimal value isgiven byσ̂(x̂∗ )σ̂(x̂∗ )ẑn∗ − z1−β/2 √ n , ẑn∗ + z1−β/2 √ n(27)nn
where z1−β/2 is the critical value of the standard normal P
distribution at level 1 − β/2, x̂∗n is the∗empirical optimal solution
obtainedPfrom (2), ẑn = (1/n) ni=1 H(x̂∗ ; ξi ) is the empirical optimalpvalue, and σ̂(x̂∗n ) =(1/(n − 1)) ni=1 (H(x̂∗n ; ξi ) − ẑn∗ )2 is the empirical standard deviation ofH(x̂∗n ; ξ). Since ẑn∗ is a negatively biased estimator of z ∗ , the CI (27) can suffer from under coverage.
So we also compare with a 2-sample CLT (CLT2) method, as suggested by [12], which uses firsthalf of the data to compute the empirical optimal value and solution, and then uses the remaining
half of the data to estimate the objective value fixed at the solution to generate an upper bound.
The 2-sample CLT CI is given by"#σ̂(x̂∗n/2 )σ̄(x̂∗n/2 )
∗
∗
ẑn/2
− z1−β/2 p, z̄n/2+ z1−β/2 p
(28)
n/2
n/2
∗ , x̂∗ , σ̂(x̂∗ ) are computed as before using first half of the data {ξ , . . . , ξ∗
where ẑn/2
1n/2 }, z̄n/2 =
n/2
n/2
Pn
(2/n) i= n +1 H(x̂∗n/2 ; ξi ) is the evaluation of x̂∗n/2 using the remaining half of the data, and2
qP∗ )2 is the empirical standard deviation at x̂∗ .
σ̄(x̂∗n/2 ) = (1/(n/2 − 1)) ni= n +1 (H(x̂∗n/2 ; ξi ) − z̄n/2
n/22∗Note that z̄n/2
is a positively biased estimator of z ∗ , and thus the CI (28) alleviates the under coverage issue; on the other hand, the effective sample size is reduced by half, and thus the estimates
are less accurate especially when the data size is small, which may in turn affect the coverage
probability of the CI.
Due to the limited data size, we use the single replication procedure (SRP) proposed in [2] toestimate CIs on the optimality gap. For a given solution x̂ that is independent of the data, the
SRP outputs a one-sided (1 − β) CI on the optimality gap given byσ̃(x̂∗n ),(29)
0, Ĝn (x̂) + z1−β √nPwhere as before x̂∗n is the empirical optimal solution, Ĝn (x̂) = (1/n) ni=1 (H(x̂, ξi ) − H(x̂∗n , ξi ))2 ,
and
n
i2
1 Xh
2 ∗
σ̃ (x̂n ) =(H(x̂, ξi ) − H(x̂∗n , ξi )) − (ĥ(x̂) − ẑn∗ ) ,
n−1
i=1

16
P
where ẑn∗ is the empirical optimal value and ĥ(x̂) = (1/n) ni=1 H(x̂, ξi ). In all the examples
considered below, we set β = 0.05; since we consider cases of small data size, we replace z1−β/2and z1−β in (27)-(29) respectively with tn−1,1−β/2 and tn−1,1−β , which are the critical values of theStudent’s t-distribution with n − 1 degree of freedom. Note that all the above discussion holds fordeterministically constrained problems. Nonetheless, we also apply these methods in a stochasticallyconstrained problem as a benchmark which, like the EL method (see the remark at the beginningof this section), are heuristic in this case without any formal validity proof.
Note that the EL method consists of solving a max-min and a min-min problem. Supposing thatthe original problem (1) or (3) is convex, then the max-min program is convex. In our examples
we use the built-in Matlab solvers. The min-min program, on the other hand, is more
Pn challenging
because the outer optimization involves minimizing the concave function minx∈Θ i=1 wi H(x; ξi )over w. This is not a convex problem in general. However, fixing either w or x, optimizing over theother variable becomes a convex problem. Thus one approach is to do alternating minimization,by iteratively minimizing w and x while fixing each others, until no improvement is observed. Suchtype of schemes has appeared in chance-constrained programming (e.g., [5, 20, 9]), and it appears
to work well in our examples despite a lack of global convergence guarantee.

5.1

Quadratic Optimization

We consider a simple unconstrained problem of minimizing a quadratic function
min E[(x − ξ)2 ],

(30)

x

where ξ follows an unknown distribution F c . It is easy to see that the optimal solution is x∗ = E[ξ]and the optimal value is z ∗ = V ar(ξ). We set F c as a standard normal distribution, and thusx∗ = 0 and z ∗ = 1.
Assuming we are given n observations from the normal distribution, we implement the different
methods to obtain 95% confidence bounds for the optimal value of (30). We test on three cases
where we randomly generate n = 10, 50, 100 data points from F c . For each case, we repeat theexperiment 1000 times, and note the empirical coverage probability, mean upper and lower bounds,
and the mean and standard deviation of the interval width for each method. The results are
summarized in Table 1.

n = 10n = 50n = 100

EL
CLT
CLT2
EL
CLT
CLT2
EL
CLT
CLT2

Coverageprobability
0.74
0.80
0.83
0.95
0.92
0.91
0.96
0.92
0.95

Mean lower
bound
0.40
0.12
−0.05
0.34
0.59
0.44
0.74
0.72
0.60

Mean upper
bound
1.57
1.68
2.54
2.11
1.36
1.59
1.39
1.26
1.41

Mean interval
width
1.17
1.56
2.59
1.78
0.77
1.15
0.65
0.54
0.81

Standard deviationof interval width
0.59
0.83
1.77
1.20
0.19
0.50
0.29
0.10
0.33

Table 1: Confidence intervals on optimal value of the quadratic optimization problem
To compare EL and SRP on optimality gap, we first generate a solution x̂ and compute itstrue optimality gap. Then for each of the three cases n = 10, 50, 100, we repeat the experiment171000 times for each method to obtain 95% confidence bounds and estimate their empirical coverage
probabilities. The results are summarized in Table 2, where the suboptimal solution x̂ = 0.624 andits corresponding optimality gap is 0.39.

n = 10n = 50n = 100

EL
CLT-SRP
EL
CLT-SRP
EL
CLT-SRP

Coverageprobability
0.95
0.89
0.99
0.95
0.98
0.94

Mean lower
bound
0.05
0
0.10
0
0.15
0

Mean upper
bound
1.89
1.36
1.45
0.77
0.79
0.64

Mean interval
width
1.84
1.36
1.35
0.77
0.64
0.64

Standard deviationof interval width
0.92
0.84
1.58
0.27
0.35
0.17

Table 2: Confidence intervals on optimality gap of the quadratic optimization problem

5.2

CVaR Estimation

In this example, we consider estimating CVaRα,F c (ξ), the α-level conditional-value-at-risk of arandom variable ξ, which we assume follows an unknown distribution F c . This can be rewritten asa stochastic optimization problem:
1+
min x +
E[(ξ − x) ] ,
(31)
x∈R
1−α
where (·)+ is short for max(·, 0). We set F c as a standard normal distribution and α = 0.9. As theprevious example in Section 5.1, we run the experiment 1000 times for each method and each caseof n = 50, 100 (the case of n = 10 is not included in this example since the data size is too small
for estimating E[(ξ − x)+ ]). To obtain confidence intervals on optimality gap, we first generatea solution x̂ and evaluate its optimality gap using a large (108 ) sample size, and then run theexperiment 1000 times for each case. The results are summarized in Table 3 and 4. Note that thetrue optimal value can be accurately calculated and is equal to 1.755; the suboptimal solution inthis experiment is 0.71 with optimality gap 0.36.

n = 50n = 100

EL
CLT
CLT2
EL
CLT
CLT2

Coverageprobability
0.87
0.85
0.76
0.95
0.90
0.85

Mean lower
bound
1.22
1.23
1.04
1.35
1.37
1.22

Mean upper
bound
2.27
2.19
2.59
2.29
2.09
2.32

Mean interval
width
1.04
0.97
1.55
0.94
0.72
1.10

Standard deviationof interval width
0.42
0.40
1.20
0.28
0.21
0.65

Table 3: Confidence intervals on optimal value of the CVaR estimation problem

5.3

Portfolio Optimization

Our last example considers minimizing the CVaR risk associated with the loss of an investment,
subject to the condition that the expected return should exceed a certain threshold. Let’s denoteby x = [x1 , . . . , xd ]0 the vector of holding proportions in d assets, ξ = [ξ 1 , . . . , ξ d ]0 the random vector

18n = 50n = 100

EL
CLT-SRP
EL
CLT-SRP

Coverageprobability
0.99
0.89
0.99
0.93

Mean lower
bound
0.04
0
0.07
0

Mean upper
bound
1.51
0.98
1.04
0.76

Mean interval
width
1.47
0.98
0.97
0.76

Standard deviationof interval width
0.51
0.50
0.26
0.29

Table 4: Confidence intervals on optimality gap of the CVaR estimation problem
of asset returns, and rb the threshold for expected return. We assume short selling is not allowed.
The problem can be written asminx
CV aRα (−ξ 0 x)
subject to E[ξ 0 x] ≥ rb
Pd
i=1 xi = 1
xi ≥ 0, i = 1, . . . , d(32)

We can rewrite the problem in the form of (3) as
1minx,cc + 1−αE[(−ξ 0 x − c)+ ]0
subject to E[ξ x] ≥ rb
Pd
i=1 xi = 1
xi ≥ 0, i = 1, . . . , d(33)

The parameter setting is as follows: ξ follows a normal distribution with mean µ = [0.8, 1.2]0
and covariance Σ = [1 0; 0 4]; the minimum expected return is rb = 1; the CVaR level is α = 0.9,
and the confidence level is 1 − β = 0.95. It is easy to verify that the optimal solution to (33) isx∗ = [0.5, 0.5]0 , and the associated optimal value can be evaluated by Monte Carlo simulation with alarge number (108 ) of samples, which yields z ∗ ≈ 0.96. For comparison, we also implement the CLT
and 2-sample CLT methods by computing the CIs according to (27) or (28); though the validity ofthese schemes has not been proved, we use them as heuristic to provide a benchmark. To compare
CIs on optimality gap, we randomly generate a suboptimal solution and evaluate its optimality gap
using a large (108 ) sample size. For each case of n = 50, 100, we repeat the experiment 1000 times,
and summarize the numerical results in Table 5 and 6. In this experiment, the suboptimal solution
is [0.09, 0.91] with optimality gap 1.09.

n = 50n = 100

EL
CLT
CLT2
EL
CLT
CLT2

Coverageprobability
0.63
0.52
0.73
0.71
0.50
0.69

Mean lower
bound
0.11
0.51
0.47
0.29
0.71
0.62

Mean upper
bound
1.22
1.71
2.57
1.36
1.61
2.01

Mean interval
width
1.11
1.20
2.09
1.07
0.90
1.39

Standard deviationof interval width
0.60
0.71
1.55
0.64
0.37
0.96

Table 5: Confidence intervals on optimal value of the portfolio optimization problem

5.4

Summary of Numerical Results

We note in all three examples EL in general has the highest coverage probability on optimal values.
Although EL in general has wider intervals than the direct CLT method, its interval widths are19n = 50n = 100

EL
CLT-SRP
EL
CLT-SRP

Coverageprobability
0.90
0.72
0.95
0.65

Mean lower
bound
0.25
0
0.63
0

Mean upper
bound
3.24
1.59
2.15
1.13

Mean interval
width
2.99
1.59
1.52
1.13

Standard deviationof interval width
2.41
1.06
0.84
2.51

Table 6: Confidence intervals on optimality gap of the portfolio optimization problem
often comparable to or smaller than the 2-sample CLT method, which usually has higher coverage
probability than the plain CLT method. EL also has higher coverage probabilities on the optimality gap than SRP, accompanied by wider intervals than SRP. Overall speaking, EL performs
competitively compared to the CLT methods.
One thing worth mentioning is that the empirical coverage probability in the last example issmaller compared to the previous two examples. A conjectured potential reason is the invalidity ofboth the CLT-based and the EL methods in this stochastically constrained problem. Here the EL
method gives roughly comparable coverage probabilities as the CLT methods on optimal values,and higher coverage probabilities on optimality gaps.