Proceedings of the 2019 Winter Simulation Conference
N. Mustafee, K.-H.G. Bae, S. Lazarova-Molnar, M. Rabe, C. Szabo, P. Haas, and Y.-J. Son, eds.

VALIDATING OPTIMIZATION WITH UNCERTAIN CONSTRAINTS
Henry Lam
Huajie Qian
Department of Industrial Engineering and Operations Research
Columbia University
500 W. 120th Street
New York, NY 10027, USA

ABSTRACT
We consider optimization with uncertain or probabilistic constraints under the availability of limited data orMonte Carlo samples. In this situation, the obtained solutions are subject to statistical noises that affect boththe feasibility and the objective performance. To guarantee feasibility, common approaches in data-driven
optimization impose constraint reformulations that are “safe” enough to ensure solution feasibility with highconfidence. Often times, selecting this safety margin relies on loose statistical estimates, in turn leading
to overly conservative and suboptimal solutions. We propose a validation-based framework to balance
the feasibility-optimality tradeoff more efficiently, by leveraging the typical low-dimensional structure ofsolution paths in these data-driven reformulations instead of estimates based on the whole decision space
utilized by past approaches. We demonstrate how our approach can lead to a feasible solution with lessconservative safety adjustment and confidence guarantees.
1

INTRODUCTION

We consider stochastically constrained optimization problems in the form
min

x∈X

f (x)
(1)

subject to H(x) := EF [h(x, ξ )] ≥ γ
where x ∈ Rd is the decision variable, X ⊂ Rd is the deterministic decision space, ξ ∈ Rm is a random vector
following an unknown distribution F, and EF [·] and correspondingly PF (·) denotes the expectation andprobability under F. h(·, ·) : Rd × Rm → R is a known function, which represents some notion of gain whose
expected value is constrained above a given threshold γ. We assume f (·) to be a deterministic function, as
our focus is on handling the stochasticity in the constraint (this assumption can be relaxed). Formulation(1) appears as chance constrained programs (CCPs) (Prékopa 2003) when h(x, ξ ) := 1((x, ξ ) ∈ A), in whichcase we want the solution to satisfy the event with high probability. It also appears in convex expectedvalue constrained programs (Atlason et al. 2004; Krokhmal et al. 2002) when h(x, ξ ) is concave in x forevery ξ .
We focus on the situation where the governing distribution F is not fully known but only observed
through a set of i.i.d. data or Monte Carlo sample {ξ1 , . . . , ξn }. The goal is to compute a solution, basedon the available data, that is both feasible and possesses good objective performance. Note that, due tothe statistical noise from the data, feasibility can be guaranteed at best with a high confidence, which inturn also imposes a price on the optimality of the obtained solution. Finding a good solution thus requires
striking a balance between feasibility and optimality: If a procedure is overly protective (i.e., removing
much of the feasible region), then the obtained solution will be surely feasible, but optimality will suffer,Lam and Qianand vice versa. The efficacy of this balancing depends on the efficiency in assimilating the data into
the optimization and in the corresponding estimation procedure. The aim of this paper is to construct atheoretically sound and implementable framework that makes such balancing highly efficient.
2

EXISTING APPROACHES AND CHALLENGES

Before introducing our framework, let us first discuss the established approaches and point out the major
challenges for this problem, which has been studied prominently in the data-driven optimization literature.
The common practice is to reformulate the unknown constraint in (1) into a data-driven constraint that
depends only on the data, so that this reformulated constraint is “safe”. To be more precise, let F be thefeasible region of (1). We would like to construct a region Fˆ that depends only on {ξ1 , . . . , ξn }, such thatan optimal solution x̂∗ (or a best obtainable feasible solution to (2) if exact computation is hard) of the
problem
min
f (x)
x∈X
(2)
subject to x ∈ Fˆsatisfies

Pdata (x̂∗ ∈ F ) ≥ 1 − β(3)where Pdata refers to the probability generated from the data. In other words, the data-driven solution x̂∗is indeed feasible for the original problem, with a confidence level 1 − β .
Note that simply replacing H(·) with a naive point estimate (e.g., the sample mean from the data) is
typically inadequate to guarantee (3). To see this, suppose the true optimal solution x∗ is at the boundaryof F , i.e., H(x∗ ) = γ. If we consider Fˆ = {x : (1/n) ∑ni=1 h(x, ξi ) ≥ γ}. Then, roughly speaking, with half
probability the obtained solution x̂∗ will have H(x̂∗ ) below γ , which leads to infeasibility for the original
problem. This issue may not arise if x∗ or x̂∗ is in the interior of the feasible region. However, a priori
we do not know our decision. Thus, in general, to guarantee (3), one would impose a safety margin on
an estimation of H(·), such that any solution obtained from (2) is also feasible for (1), with the required
confidence. That is, we want
Pdata (Fˆ ⊂ F ) ≥ 1 − β .
(4)We contend that most approaches in data-driven optimization rely on the above reasoning and are basedon (4). In particular, (4) provides a convenient way to certify feasibility, by requiring that all solutions
feasible for (2) are also feasible for (1) with high confidence. This set-level guarantee generally leads toa simultaneous estimation problem across all x in the decision space X , for which a proper control ofthe statistical error can lead to a substantial shrinkage of the size of Fˆ that exacerbates with problem
dimension (either of the decision space or the probability space).
We provide several examples to illustrate the phenomenon above. Our discussion considers the (single)
CCP, where H(x) is in the form PF (G(x, ξ ) ≤ b) with G(x, ξ ) : Rd × Rm → R.
Sample average approximation (SAA): The SAA approach sets
(
Fˆ =

1 n
x ∈ X : ∑ 1(G(x, ξi ) + ε ≤ b) ≥ γ + δ
n i=1

)

where εpand δ are suitably tuned parameters. For example, when G is Lipschitz continuous in x, selectingδ = Ω( (d/n) log(1/ε)) can guarantee (4) (Luedtke and Ahmed 2008), and similar relations also holdin discrete decision space (Luedtke and Ahmed 2008) and expected value constraints (Wang and Ahmed
2008). These estimates come from concentration inequalities in which union bounds are needed and give√rise to the dependence on the dimension d. Note that the resulting margin δ scales in the order of d,and also that to get any reasonably small δ n must be of higher order than d.

Lam and Qian
Robust optimization (RO) and safe convex approximation (SCA): RO setsFˆ = {x ∈ X : G(x, ξ ) ≤ b, for all ξ ∈ U }(5)where U is known as the uncertainty set, and ξ in (5) is viewed as a deterministic unknown (Bertsimas et al.
2011; Ben-Tal et al. 2009). A common example of U is an ellipsoidal set {ξ : (ξ − µ̂)T Σ̂−1 (ξ − µ̂) ≤ ρ}
where µ̂ ∈ Rd , Σ̂ ∈ Rd×d a positive semidefinite matrix, and ρ ∈ R. Here the center µ̂ and shape Σ̂ typicallycorrespond to the mean and covariance of the data, and ρ controls the set size. A duality argument showsthat, in the case of linear chance constraint in the form G(x, ξ ) = xT ξ , (5) is equivalent to the quadratic√
constraint µ̂ T x + ρkΣ̂1/2 xk2 ≤ b. Using such type of convex constraints as inner approximations for
intractable chance constraints is also known as SCA (e.g., Nemirovski 2003; Nemirovski and Shapiro
2006).
If the random variable ξ is known to be bounded, the above approach guarantees the obtained solution
has a satisfaction probability of order 1 − e−ρ/2 via Hoeffding’s inequality, and ρ is chosen by matching
this expression with the tolerance level γ. Although ρ calibrated this way may not explicitly depend on
the problem dimension, the level of conservativeness measured by the shrinkage of F to Fˆ is governedby concentration bounds whose tightness can be challenging to quantify. Another viewpoint that has beenutilized recently in data-driven RO (Bertsimas et al. 2018; Tulabandhula and Rudin 2014; Goldfarb andIyengar 2003; Hong et al. 2017) is to take U to be a set that contains γ-content of the distribution of ξ , i.e.,
PF (ξ ∈ U ) ≥ γ, with a confidence level 1 − β . In this case, a solution x̂∗ obtained from solving (5) would
satisfy PF (G(x̂∗ , ξ ) ≤ b) ≥ γ with at least 1 − β confidence, thus achieving (4) as well. Such generated
uncertainty set however typically has a size that scales with the dimension of the probability space. For
example, consider G(x, ξ ) = xT ξ with ξ ∈ Rm being standard multivariate Gaussian and the uncertainty
set U is an ellipsoid with µ̂ and Σ̂ being the true mean and covariance, i.e., U = {ξ ∈ Rd : kξ k22 ≤ ρ}.
Then, in order to make U a γ-content set the radius ρ√has to be at least of order m since kξ k22 has a mean√m, resulting in the robust counterpart ρ kxk2 = Θ( m) kxk2 ≤ b. However, the exact chance constraintzγ kxk2 ≤ b, where zγ is the γ-quantile of the univariate standard normal, is independent of the dimension.
Distributionally robust optimization (DRO): DRO setsFˆ = x ∈ X : inf EQ [h(x, ξ )] ≥ γ
Q∈Uwhere U is a set in the space of probability measures that is constructed from data, and is often known as
the ambiguity set or uncertainty set. The rationale here is similar to RO, but views the uncertainty in termsof the distribution. If U is constructed such that it contains the true distribution F with high confidence,i.e., Pdata (F ∈ U ) ≥ 1 − β , then a solution x̂∗ obtained from the DRO will satisfy PF (G(x̂∗ , ξ ) ≤ b) ≥ γwith at least 1 − β confidence so that (4) holds.
Popular choices of U include moment sets, i.e., specifying the moments of Q (to be within a range for
instance) (Delage and Ye 2010; Wiesemann et al. 2014; Goh and Sim 2010), and distance-based sets, i.e.,
specifying Q in the neighborhood ball surrounding a baseline distribution, where the ball size is measured
by a statistical distance such as φ -divergence (Petersen et al. 2000; Ben-Tal et al. 2013; Glasserman and
Xu 2014; Lam 2016; Hu and Hong 2013; Jiang and Guan 2016) or Wasserstein distance (Esfahani andKuhn 2018; Blanchet and Murthy 2019; Gao and Kleywegt 2016).
Ensuring Pdata (F ∈ U ) ≥ 1 − β means that U is a confidence region for F. In the moment set case, thisboils down to finding confidence regions for the moments whose sizes in general scale with the probability
space dimension. To explain, when only the mean EF [ξ ] is estimated, the confidence region constructed12from the central limit theorem (CLT) takes the form {µ̂ + Σ̂ 2 v : v ∈ Rm , kvk22 ≤ χm,1−β}, where µ̂ and Σ̂ are2
the sample mean and covariance and χm,1−β (which is of order m) is the 1 − β quantile of the χ 2 distribution√
with degree of freedom m, therefore the diameter of the confidence region scales as m. When the mean

Lam and Qianand covariance are jointly estimated, the dimension dependence scales up further. In the distance-based
set case, one needs to estimate statistical distances. If the Wasserstein distance is used to construct theball surrounding the empirical distribution, results from measure concentration (Fournier and Guillin 2015)
1indicate that the ball size needs to be of order n− m to ensure Pdata (F ∈ U ) ≥ 1 − β . Alternatively, if U isconstructed as a φ -divergence ball surrounding some nonparametric kernel-type density estimate, resultsfrom kernel density estimation (see Section 4.3 in Wand and Jones 1994) suggests that the estimation error
4
is of order n− m+4 . In either case, the required size of the uncertainty set exhibits exponential dependence
on the dimension. Recently, the empirical likelihood method has also been proposed to calibrate the ballsize such that U can be (much) smaller than what is needed in being a confidence region for F, while atthe same time (4) still holds (Lam and Zhou 2017;Duchi et al. 2016; Lam 2019; Blanchet and Kang 2016).
However, the ball size in this approach scales as the supremum of a so-called χ 2 -process over the decision
space (e.g., Lam 2019). An analysis using metric entropy (e.g., Example 2 in Section 14 in Lifshits (1995))
shows that the χ 2 -process supremum can scale linearly in the decision space dimension d, a much betterbut still considerable dependence on the dimension.
Finally, we mention that the only exceptional paradigm to our knowledge that provides an alternate
guarantee for (3) in the case of CCP, without using (4), is scenario optimization (SO) (e.g., Calafiore andCampi 2005; Campi and Garatti 2008). In its basic form, this approach setsFˆ = {x ∈ X : G(x, ξi ) ≤ b for all i = 1, . . . , n}i.e., using sampled constraints formed from the data. As the number of constraints increases, Fˆ is postulated
to populate the decision space in some sense and ensure the obtained solution x̂∗ lies in F . While the sample
size required in the basic SO is linear in the decision dimension d, recent works reduce this dependence
by using regularization (Campi and Carè 2013), tighter support rank estimates (Schildbach et al. 2013;
Campi and Garatti 2018) and validation-type schemes (Carè et al. 2014; Calafiore 2017). The approachthat we propose next is closest to some of the validation-type schemes suggested for SO, but substantiallymore general as it applies beyond CCP and to all the exemplified methods mentioned above.
In the following sections, we will overview our main idea, procedures and guarantees, leaving the fulldemonstration of our results and mathematical analyses to a journal version of this work. Section 3 explains
the rationale of our framework and thereby invokes our two-phase framework. Sections 4 and 5 present
two methods to be used in the framework. Section 6 shows a numerical example. Section 7 concludes thepaper.
3

FRAMEWORK AND RATIONALE

Our primary goal is to create a framework to obtain solutions from data-driven optimization reformulations
that are less conservative than past proposals, especially in terms of the dimension dependence that someof them face. Our framework works for general stochastically constrained problems in (1) including CCP.
Our key observation is the following. In all the described approaches, the data-driven reformulation
requires a parameter that controls the level of conservativeness:
1.
2.
3.
4.

SAA: safety margin δRO and SCA: uncertainty set size ρ
DRO: neighborhood ball size or moment set size
SO: number of constraints

These parameters have the properties that setting it to one extreme (e.g., 0) would signal no uncertainty inthe formulation, leading to a solution very likely infeasible, while setting it to another extreme (e.g., ∞)would cover the entire decision space, leading to a very conservative solution. Besides SO, the rationale
undertaken in Section 2 is to select these parameters so that the resulting feasible region Fˆ is containedLam and Qian
entirely in F , which involves a vast simultaneous estimation problem on the whole decision space andpotentially makes the choice of this parameter overly protective.
On the other hand, given a specific data-driven reformulation, it is easy to see that no matter how wechoose this “conservativeness” parameter, the solution must lie in a low-dimensional manifold. To be more
precise, suppose the given data-driven reformulation is
min

x∈X

f (x)

subject to x ∈ Fˆ (s)(6)

where s ∈ R denotes the conservativeness parameter that determines the size of feasible region Fˆ . We
denote the obtained solution from (6) as x∗ (s). The solution path {x∗ (s) : s ∈ R} contains all possible
obtainable solutions from the data-driven reformulation (6). Intuitively, it is only necessary to focus on
this solution path, instead of a simultaneous estimation for the whole decision space.
Nonetheless, besides the conservativeness parameter, a data-driven reformulation could have otherparameters playing various roles (e.g., center and shape of the set in ellipsoidal RO, baseline distribution
in distance-based DRO etc.). The flexibility of these parameter values can enlarge the obtainable solution
space and elevate dimension dependence. Suppose we want to contain this enlargement, and at the sametime be able to select the optimal candidate among the low-dimensional manifold {x∗ (s) : s ∈ R}. We
propose the following two-phase framework to achieve this.
Algorithm 1 The Two-Phase Framework
Input: data ξ 1:n = {ξ1 , . . . , ξn }; numbers of data n1 , n2 allocated to each phase (n1 + n2 = n); a confidence
level 1 − β ; a given method to construct data-driven reformulation with a (possibly multi-dimensional)
parameter s ∈ S; a set of candidate parameter values {s1 , s2 , . . . , s p } ⊆ S.
Phase one:
1. Use n1 observations, which we index as {ξn2 +1 , . . . , ξn } for convenience, to construct the data-driven
reformulation OPT (s) in the form (6) parameterized by s ∈ S.
2. For each j = 1, . . . , p, compute the optimal solution x∗ (s j ) of OPT (s j ).
Phase two:
Use a validator V to select (ŝ∗ , x∗ (ŝ∗ )) = V ({ξ1 , . . . , ξn2 }, {x∗ (s1 ), . . . , x∗ (s p )}, 1 − β ), where x∗ (ŝ∗ ) is asolution and ŝ∗ is the associated parameter value.
Output: x∗ (ŝ∗ ).
Our procedure (Algorithm 1) splits the data into two groups. With the first group of data, we construct
a given data-driven reformulation parametrized by a conservativeness parameter s that varies over a space
S, which we call OPT (s). We obtain the optimal solution x∗ (s) for a range of values s = s j , j = 1, . . . , p.
This step assumes the availability of an efficient solver for OPT (s). Next, the second group of data is
fed into a validator V that aims to identify the best feasible solution x∗ (ŝ∗ ) among {x∗ (s j ) : j = 1, . . . , p}.
The number of points p required to validate depends crucially on the size of S, which is constructed to
be low-dimensional. Here for simplicity we assign the last n1 observations to Phase two and all other toPhase one, and note that it is equivalent to uniformly choose n1 observations for Phase one because of thei.i.d. nature. There are multiple ways to set up the validator V , each with its own benefits which we will
describe precisely in the next two sections. Because of space limits, we solely focus on chance constraintsof the form
min
f (x)
x∈X
(7)
subject to P(x) := PF ((x, ξ ) ∈ A) ≥ 1 − αwhere A is a deterministic subset of Rd × Rm , and 1 − α is the tolerance level (i.e., 1 − α = γ in (1)).

Lam and Qian
4

VALIDATION VIA MULTIVARIATE GAUSSIAN SUPREMUM

Our first validator uses a simultaneous estimation of the satisfaction probability P(x) over the discretizedsolution path of x∗ (s). More precisely, given the solution set {x∗ (s j ) : j = 1, . . . , p}, we use a sample
2
1((x, ξi ) ∈ A) − ε, to replace theaverage with an appropriately calibrated safety margin, i.e., (1/n2 ) ∑ni=1
unknown P(·) in (7) and output the best solution among the set. The margin ε is calibrated via the limiting
2
1((x∗ (s j ), ξi ) ∈ A)) j=1,...,p which is multivariate Gaussian. It contains a critical
distribution of ((1/n2 ) ∑ni=1
value q1−β that is the quantile of a Gaussian supremum. Algorithms 2 and 3 describe two variants of thisvalidator, one unnormalized while another one normalized by the standard deviation at each s j . In thefollowing, we denote Np (0, Σ) as a p-dimensional Gaussian vector with mean zero and covariance Σ.
Algorithm 2 V : Unnormalized Gaussian Supremum
Input: {ξ1 , . . . , ξn2 }, {x∗ (s1 ), . . . , x∗ (s p )}, 1 − β
2
1. For each j = 1, . . . , p compute the sample mean P̂j = (1/n2 ) ∑ni=1
1((x∗ (s j ), ξi ) ∈ A) and sample
n2
covariance matrix Σ̂ with Σ̂( j1 , j2 ) = (1/n2 ) ∑i=1 (1((x∗ (s j1 ), ξi ) ∈ A) − P̂j1 )(1((x∗ (s j2 ), ξi ) ∈ A) − P̂j2 ).
2. Compute q1−β , the (1 − β )-quantile of max{Z1 , . . . , Z p } where (Z1 , . . . , Z p ) ∼ Np (0, Σ̂), and let

n
o
q1−β
ŝ∗ = argmin f (x∗ (s j )) : P̂j ≥ 1 − α + √ , 1 ≤ j ≤ p .
n2(8)Output: ŝ∗ , x∗ (ŝ∗ ).
Algorithm 3 V : Normalized Gaussian Supremum
Input: {ξ1 , . . . , ξn2 }, {x∗ (s1 ), . . . , x∗ (s p )}, 1 − β
1. Same as in Algorithm 2.
2. Denote σ̂ 2j = Σ̂( j, j). Compute q1−β , the (1 − β )-quantile of max{Z j /σ̂ j : σ̂ 2j > 0, 1 ≤ j ≤ p} where(Z1 , . . . , Z p ) ∼ Np (0, Σ̂), and letnoq1−β σ̂ jŝ∗ = argmin f (x∗ (s j )) : P̂j ≥ 1 − α + √,1 ≤ j ≤ p .
n2(9)Output: ŝ∗ , x∗ (ŝ∗ ).
The first Gaussian supremum validator (Algorithm 2) is reasoned from a joint CLT that governs the√convergence of n2 (P̂1 − P(x∗ (s1 )), . . . , P̂p − P(x∗ (s p ))) to Np (0, Σ), where Σ( j1 , j2 ) = Cov(1((x∗ (s j1 ), ξ ) ∈A), 1((x∗ (s j2 ), ξ ) ∈ A)). Using the sample covariance Σ̂ from Step 1 of Algorithm 2 as an approximationof Σ, we have, by the continuous mapping theorem,√
max n2 (P̂j − P(x∗ (s j ))) ≈ max Z j in distribution
1≤ j≤p1≤ j≤pwhere (Z1 , . . . , Z p ) ∼ Np (0, Σ̂). Therefore using the 1 − β quantile q1−β of the Gaussian supremum in themargin leads toq1−β
P(x∗ (s j )) ≥ P̂j − √ for all j = 1, . . . , p, with probability ≈ 1 − β .
n2The second validator (Algorithm 3) uses an alternate version of the CLT that is normalized by the√componentwise standard deviation σ j , i.e., n2 ((P̂1 − P(x∗ (s1 )))/σ1 , . . . , (P̂p − P(x∗ (s p )))/σ p ) converges

Lam and Qianto Np (0, DΣD), where D is a diagonal matrix of 1/σ j , j = 1, . . . , p. Note that the quantile q1−β in both
validators can be computed to high accuracy via Monte Carlo.
The following two theorems describe the feasibility guarantees of the obtained solutions output by thevalidators in Algorithms 2 and 3:
Theorem 1 Let ᾱ = 1 − max1≤ j≤p P(x∗ (s j )). For every solution set {x∗ (s j ) : 1 ≤ j ≤ p}, every n2 , andβ ∈ (0, 12 ), the solution output by Algorithm 2 satisfies
Pξ 1:n (x∗ (ŝ∗ ) is feasible for (7)) ≥ 1 − β −C log7 (pn )  16
2

n2 α

2

with

s

ε = α − ᾱ −C+ exp − cn2 min{ε,log(p/β ) n2+ε 2 
}ᾱ(10)where C and c are universal constants, and Pξ 1:n denotes the probability with respect to Phase two data2{ξ1 , . . . , ξn2 } conditioned on Phase one data.
Theorem 2 Under the same conditions of Theorem 1, the solution output by Algorithm 3 satisfies
Pξ 1:n (x∗ (ŝ∗ ) is feasible for (7)) ≥ 1 − β −C
2

with

s

ε = α − ᾱ −C

 log7 (pn )  162n2 α

log2 (pn2 )ε 2 
+ √+ exp − cn2 min{ε, }
n2 α
ᾱ


ᾱ + log(n2 α)/n2 log(p/β ) 
n2
+

(11)

where C and c are universal constants.
In both Theorems 1 and 2, the finite-sample coverage probability consists of two sources of errors.
The first source comes from the CLT approximation that decays polynomially in the Phase 2 sample sizen2 . The second error arises from the possibility that none of the solutions {x∗ (s1 ), . . . , x∗ (s p )} satisfies thecriterion in (8) or (9), which vanishes exponentially fast. When ε in (10) and (11) are of constant order,
the CLT error dominates. In this case the finite-sample error depends logarithmically on p, the numberof candidate parameter values, and the bounds dictate a coverage tending to 1 − β when p is as large as
1/7
exp(o(n2 )). The derivation of this logarithmic dependence on p utilizes a high-dimensional CLT recently
developed in (Chernozhukov et al. 2017).
We explain the implication on the dimensionality of the problem. Note that to sufficiently cover thewhole solution path, p is typically exponential in the dimension of S, denoted dim(S) (this happens when
we uniformly discretize the parameter space S). The discussion above thus implies a requirement that n2 isof higher order than dim(S)7 . Here the low dimensionality of S is crucial; for instance, a one-dimensional
conservativeness parameter s would mean dim(S) = 1, so that a reasonably small n2 can already ensure
adequate feasibility coverage. Moreover, the margin adjustments in Algorithms 2 and 3 both depend only
on dim(S). Thus, the choice of ŝ∗ also depends only on dim(S), but not the dimension of the whole decision
space. This indicates a substantial reduction in conservativeness compared to the approaches described in
Section 2.
Comparing between the two validators, we also see that the normalized one (Algorithm 3) is statisticallymore efficient than the unnormalized one (Algorithm 2) when the tolerance level 1 − α is large (i.e., thecommon case). More specifically, in order to make the exponential error non-trivial, one needs at least
ε > 0. In the case of Algorithm 2 expression (10) suggests that, after ignoring the logarithmic factor
log(p/β ), this requires an n2 of order (α − ᾱ)−2 , whereas in (11) it can be seen to need only an n2 oforder α(α − ᾱ)−2 , a much smaller size when 1 − α is close to 1.
Theorems 1 and 2 also give immediately the following asymptotic feasibility guarantee:

Lam and Qian
Corollary 1 Let ᾱ = 1 − max1≤ j≤p P(x∗ (s j )). For every solution set {x∗ (s j ) : 1 ≤ j ≤ p} such that ᾱ < αand every β ∈ (0, 21 ), the solution output by Algorithm 2 or 3 satisfies
lim inf

1/7
n2 →∞ and p exp(−n2 )→0

5

Pξ 1:n (x∗ (ŝ∗ ) is feasible for (7)) ≥ 1 − β .
2

VALIDATION VIA UNIVARIATE GAUSSIAN MARGIN

We offer an alternate validator that can perform more efficiently than Algorithms 2 and 3, provided thatsome further assumptions are in place. This is a scheme that simply uses a standard univariate Gaussian
critical value to calibrate the margin (Algorithm 4).
Algorithm 4 outputs a solution with objective value no worse than Algorithms 2 and 3. Comparing thecriteria to choose ŝ∗ , we see that, thanks to the stochastic dominance of the maximum among a multivariate
Gaussian vector over each of its individual components, the margin in (8) satisfies q1−β ≥ z1−β σ̂ j for allj, and similarly the margin in (9) satisfies q1−β σ̂ j ≥ z1−β σ̂ j , so that both are bounded from below by themargin of (12). Consequently the solution from (12) achieves an objective value no worse than the othertwo.
Algorithm 4 V : Univariate Gaussian Validator
Input: {ξ1 , . . . , ξn2 }, {x∗ (s1 ), . . . , x∗ (s p )}, 1 − β
2
1. For each j = 1, . . . , p compute the sample mean P̂j = (1/n2 ) ∑ni=1
1((x∗ (s j ), ξi ) ∈ A) and sample varianceσ̂ 2j = P̂j (1 − P̂j ).
2. Compute()
z1−β σ̂ j∗∗ŝ = argmin f (x (s j )) P̂j ≥ 1 − α + √,1 ≤ j ≤ p(12)n2

where z1−β is the 1 − β quantile of the standard Gaussian distribution.
Output: ŝ∗ , x∗ (ŝ∗ ).
The validity of the univariate Gaussian critical value is based on the statistical consistency of theobtained solution x∗ (ŝ∗ ) to some limiting solution (correspondingly ŝ∗ to some limiting optimal parametervalue) as n2 increases. Intuitively, this implies that with sufficient sample size one can focus feasibility
validation on a small neighborhood of ŝ∗ , which further suggests that we need to control only the statisticalerror at effectively one solution parametrized at ŝ∗ . For this argument to hold, however, we would needadditional assumptions on the CCP (7):Assumption 1 (Continuous objective) The objective function f (x) is continuous on X .
Assumption 2 (Linear constraint) 1((x, ξ ) ∈ A) = 1(aTk x ≤ bk for k = 1, . . . , K), where each ak has a density
on Rd and each bk is a non-zero constant.
We comment that Assumption 2 is a sufficient but not necessary condition, and in fact some of theresults presented below holds for more general constraints. For instance, results in Theorem 3 remain valid
for constraints of the form 1(aTk Ak (x) ≤ bk for k = 1, . . . , K) where each Ak (·) : Rd → Rm is a continuousfunction. We further assume a smoothness condition for the data-driven reformulation OPT (s), s ∈ S:
Assumption 3 (Piecewise continuous solution curve) The parameter space S is a finite interval [sl , su ].
The optimal solution x∗ (s) of OPT (s) exists and is unique except for a finite number of parameter valuess̃i , i = 1, . . . , M − 1 such that sl = s̃0 < s̃1 < · · · < s̃M−1 < s̃M = su , and the parameter-solution mapping x∗ (s)is uniformly continuous on [s̃0 , s̃1 ), (s̃M−1 , s̃M ], and (s̃i−1 , s̃i ) for i = 2, . . . , M − 1.
Continuity of the solution path allows approximating the whole solution curve by discretizing theparameter space S. Also note that under Assumption 3 the solution x∗ (s) exists and is unique for almost
surely every s ∈ S with respect to the Lebesgue measure. Therefore, if one discretizes the parameter

Lam and Qian
space by randomizing via a continuous distribution over S, then with probability one the solution x∗ (s) is
unique at all sampled parameter values. This provides an easy way to ensure the assumption that noneof the parameter values {s1 , . . . , sm } used in Phase one of Algorithm 1 belongs to the discontinuity set{s̃1 , . . . , s̃M−1 }.
To explain the superior performance of Algorithm 4, we introduce a notion of optimality within thesolution path {x∗ (s) : s ∈ S}. First we fill in the holes of of the solution curve. Under piecewise uniform
continuity, the parameter-solution mapping x∗ (s) on each piece (s̃i−1 , s̃i ) can be continuously extended tothe closure [s̃i−1 , s̃i ]. Specifically, we define the extended parameter-solution mapping at the endpoints of
the i-th piecexi∗ (s̃i−1 ) := lim x∗ (s), xi∗ (s̃i ) := lim x∗ (s)s→s̃i−1 +s→s̃i −∗ (s̃ ) and x∗ (s̃ ) take different values if the i-th and (i + 1)-th pieces
and xi∗ (s) = x∗ (s) for s ∈ (s̃i−1 , s̃i ). xi+1i
i i
are disconnected. With the extended parameter-solution mappings xi∗ (·)’s, we define the following as the
optimal solution set among the solution set parameterized by s

XS∗ := argmin{ f (x) : P(x) ≥ 1 − α, x = xi∗ (s) for some s ∈ [s̃i−1 , s̃i ] and i = 1, . . . , M}
and the define the set of optimal parameter as
S∗ := {s : s ∈ [s̃i−1 , s̃i ] and xi∗ (s) ∈ XS∗ for some i = 1, . . . , M}.
We need several additional technical assumptions. The first is that the chance constraint is not bindingat the endpoints of each piece of the solution path:∗ (s̃ )) 6= 1 − α for all i = 1, . . . , M − 1, P(x∗ (s )) 6= 1 − α,Assumption 4 P(xi∗ (s̃i )) 6= 1 − α and P(xi+1il∗∗P(x (su )) 6= 1 − α, and sups∈{
P(x
(s))
>
1 − α.
/ s̃1 ,...,s̃M−1 }Since the solution path {x∗ (s) : s ∈ S} depends on the Phase one data ξ n2 +1:n , the path and hence the∗ (s̃ ) are random objects, and so the first part of Assumption 4 is expected to holdendpoints xi∗ (s̃i ), xi+1ialmost surely provided that the set {x ∈ X : P(x) = 1 − α} is a null set under the Lebesgue measure. The
second part states that the solution path contains a strictly feasible solution which in turn ensures that the
optimal solution set XS∗ is non-empty. Note that this can be achieved by simply including very conservative
parameter values in S.
Another property we assume regards the monotonicity of the feasible set size with respect to theparameter s in the reformulation OPT (s):Assumption 5 Denote by Sol(s) ⊆ X the feasible set of OPT (s). Assume Sol(s) is a closed set for alls ∈ S and Sol(s2 ) ⊆ Sol(s1 ) for all s1 , s2 ∈ S such that s1 < s2 .
Assmption 5 holds for almost all common reformulations because of a monotonic relation between theparameter s and the conservativeness level. For instance, in RO with ellipsoidal uncertainty set, the RO
feasible region shrinks with the radius of the ellipsoid, and similar relations hold for DRO, SAA, and SO.
Lastly, we assume the following technical assumption for the set of optima:
Assumption 6 For any ε > 0 there exists an s ∈/ {s̃1 , . . . , s̃M−1 } such that P(x∗ (s)) > 1−α and d(x∗ (s), XS∗ ) <∗∗∗ε, where d(x (s), XS ) := infx∈XS∗ kx (s) − xk2 .
This assumption trivially holds if we have maxx∈XS∗ P(x) > 1 − α. Otherwise, if P(x) = 1 − α for allx ∈ XS∗ , it rules out the case that the solution path x∗ (s) passes through the optima without entering theinterior of the feasible set of (7). This extreme case typically happens with zero probability, again in viewof the fact that the solution path is itself random with respect to Phase one data.
With these assumptions, we have the following asymptotic performance guarantee for Algorithm 4:
Theorem 3 Under Assumptions 1-6 hold, the optimal solution set is a singleton, i.e., XS∗ = {xS∗ }. Denoteby εs = sups∈S inf1≤ j≤p s − s j the mesh size. We also have that, with respect to {ξ1 , . . . , ξn2 }, the solution

Lam and Qianand parameter output by Algorithm 4 satisfy
lim

n2 →∞,εs →0x∗ (ŝ∗ ) = xS∗ andlim

n2 →∞,εs →0d(ŝ∗ , S∗ ) = 0

almost surely, where d(ŝ∗ , S∗ ) = infs∈S∗ |ŝ∗ − s|. Moreover(
lim infn2 →∞,εs →0 Pξ 1:n (x∗ (ŝ∗ ) is feasible for (7)) ≥ 1 − β2limn2 →∞,εs →0 Pξ 1:n (x∗ (ŝ∗ ) is feasible for (7)) = 1
2

if P(xS∗ ) = 1 − α
if P(xS∗ ) > 1 − α.

Theorem 3 states that as the mesh {s1 , . . . , s p } gets increasingly fine and and the data size grows, thesolution given by Algorithm 4 enjoys performance guarantees concerning both feasibility and optimality.
In particular, the estimated solution and the parameter converge to the unique optimal solution xS∗ and the
optimal parameter set S∗ respectively, and simultaneously the obtained solution x∗ (ŝ∗ ) is feasible with thedesired confidence level 1 − β .
6

A NUMERICAL EXAMPLE

We present a numerical example to demonstrate the performances of our framework. We consider thefollowing linear CCP
min cT x subject to PF (ξ T x ≤ b) ≥ 1 − αwhere c ∈ Rd , b ∈ R are deterministic, the distribution F of the randomness ξ ∈ Rd is multivariate Gaussianwith mean µ and covariance Σ, and the tolerance level 1 − α is set to 90%.
We test the proposed
framework on RO with ellipsoid uncertainty set that leads to a robust counterpart√
in the form µ̂ T x + skΣ̂1/2 xk2 ≤ b where µ̂ and Σ̂ are the sample mean and covariance for ξ computedfrom Phase one data. The benchmark (“SCA” in the table)pis set to an SCA (equation 2.4.11 of Ben-Tal
et al. 2009), which in our case can be expressed as µ 0 x + 2 log(1/α)kΣ1/2 xk2 ≤ b. Here, we give this
SCA or RO the advantage of knowing the true mean µ and covariance Σ of the randomness. Note that thisis equivalent to setting s = 2 log α1 . To implement our framework, we take the (1 − α)n1 -th order statisticŝ1−α of {(ξn2 +i − µ̂)0 Σ̂−1 (ξn2 +i − µ̂) : i = 1, . . . , n1 }, where ξn2 +i , i = 1, . . . , n1 are the Phase one data, sothat {ξ : (ξ − µ̂)0 Σ̂−1 (ξ − µ̂) ≤ ŝ1−α } is roughly a (1 − α)-content set for ξ (such type of quantile-based
selection has been used in Hong et al. 2017). We then set the values s j = (ŝ1−α + 20) j/50 for j = 1, . . . , 50.
Table 1 summarizes the results under dimension d = 10 and data sizes n = 200, 500. “unnorm. GS”,
“norm. GS” and “uni. Gaussian” denote Algorithms 2, 3 and 4 respectively. For each setting we repeat theexperiments 1000 times each with an independently generated data set and a data-driven solution output,
and then take down the average objective value (“mean obj. val.”) achieved by these solutions and the
proportion (“feasibility level”) of feasible solutions as the empirical feasibility coverage. Therefore, thesmaller the “mean obj. val.” is, the better is the solution in terms of optimality, and “feasibility level”≥ 95% indicates that the desired feasibility confidence level 95% is achieved and otherwise not.
Table 1: RO with ellipsoidal uncertainty set. d = 10. Data are split to n1 = n2 = n/2.
SCAmean obj. val.
feasibility level

−3.57
100%

unnorm. GS
n = 200 n = 500
−3.68
−4.42
99.9%
99.8%

norm. GSn = 200 n = 500
−4.20
−4.58
98.5%
99.6%

uni. Gaussiann = 200 n = 500
−4.43
−4.80
97.5%
98.8%

We highlight two observations. First, our framework with the three proposed validators outperforms
the SCA benchmark. In terms of the objective performance, all our validators achieve lower objective
value than SCA (with a difference ≥ 0.6), while at the same time retain the feasibility confidence to above
95% in all the three tables. Second, among the three proposed validators, the univariate Gaussian validator

Lam and Qian
appears less conservative than the Gaussian supremum counterparts in achieving better objective values,
and relatedly tighter feasibility confidence levels (i.e., closer to 95%).
7

CONCLUSION

We have studied a validation-based framework to reduce the conservativeness in data-driven optimization
with uncertain constraints, faced by the several conventional approaches that rely on implicit estimationof the entire feasible region. Our framework leverages the low dimensionality of the possible solutions
output from a data-driven reformulation class, by extracting the associated parametrized solution path andselecting the best candidate within the path. We have proposed several validators for solution selection, andhave presented results on the feasibility guarantees on the obtained solutions under our validators that scalefavorably with the problem dimension. We have illustrated the benefits of our approach with a numericalexample.
ACKNOWLEDGMENTS
We gratefully acknowledge support from the National Science Foundation under grants CAREER CMMI1653339/1834710 and IIS-1849280.
REFERENCES
Atlason, J., M. A. Epelman, and S. G. Henderson. 2004. “Call Center Staffing with Simulation and Cutting Plane Methods”.
Annals of Operations Research 127(1-4):333–358.
Ben-Tal, A., D. Den Hertog, A. De Waegenaere, B. Melenberg, and G. Rennen. 2013. “Robust Solutions of Optimization
Problems Affected by Uncertain Probabilities”. Management Science 59(2):341–357.
Ben-Tal, A., L. El Ghaoui, and A. Nemirovski. 2009. Robust Optimization. Princeton, New Jersey: Princeton University Press.
Bertsimas, D., D. B. Brown, and C. Caramanis. 2011. “Theory and Applications of Robust Optimization”. SIAM Review 53(3):464–
501.
Bertsimas, D., V. Gupta, and N. Kallus. 2018. “Data-Driven Robust Optimization”. Mathematical Programming 167(2):235–292.
Blanchet, J., and Y. Kang. 2016. “Sample Out-Of-Sample Inference Based on Wasserstein Distance”. arXiv preprint
arXiv:1605.01340.
Blanchet, J., and K. Murthy. 2019. “Quantifying Distributional Model Risk via Optimal Transport”. Mathematics of Operations
Research 44(2):565–600.
Calafiore, G., and M. C. Campi. 2005. “Uncertain Convex Programs: Randomized Solutions and Confidence Levels”. Mathematical
Programming 102(1):25–46.
Calafiore, G. C. 2017. “Repetitive Scenario Design”. IEEE Transactions on Automatic Control 62(3):1125–1137.
Campi, M. C., and A. Carè. 2013. “Random Convex Programs with L1 -Regularization: Sparsity and Generalization”. SIAM
Journal on Control and Optimization 51(5):3532–3557.
Campi, M. C., and S. Garatti. 2008. “The Exact Feasibility of Randomized Solutions of Uncertain Convex Programs”. SIAM
Journal on Optimization 19(3):1211–1230.
Campi, M. C., and S. Garatti. 2018. “Wait-And-Judge Scenario Optimization”. Mathematical Programming 167(1):155–189.
Carè, A., S. Garatti, and M. C. Campi. 2014. “Fast—Fast Algorithm for the Scenario Technique”. Operations Research 62(3):662–
671.
Chernozhukov, V., D. Chetverikov, and K. Kato. 2017. “Central Limit Theorems and Bootstrap in High Dimensions”. The
Annals of Probability 45(4):2309–2352.
Delage, E., and Y. Ye. 2010. “Distributionally Robust Optimization Under Moment Uncertainty with Application to Data-Driven
Problems”. Operations Research 58(3):595–612.
Duchi, J., P. Glynn, and H. Namkoong. 2016. “Statistics of Robust Optimization: A Generalized Empirical Likelihood Approach”.
arXiv preprint arXiv:1610.03425.
Esfahani, P. M., and D. Kuhn. 2018. “Data-Driven Distributionally Robust Optimization Using the Wasserstein Metric:
Performance Guarantees and Tractable Reformulations”. Mathematical Programming 171(1-2):115–166.
Fournier, N., and A. Guillin. 2015. “On the Rate of Convergence in Wasserstein Distance of the Empirical Measure”. Probability
Theory and Related Fields 162(3-4):707–738.
Gao, R., and A. J. Kleywegt. 2016. “Distributionally Robust Stochastic Optimization with Wasserstein Distance”. arXiv preprint
arXiv:1604.02199.
Glasserman, P., and X. Xu. 2014. “Robust Risk Measurement and Model Risk”. Quantitative Finance 14(1):29–58.

Lam and Qian
Goh, J., and M. Sim. 2010. “Distributionally Robust Optimization and Its Tractable Approximations”. Operations Research 58(4part-1):902–917.
Goldfarb, D., and G. Iyengar. 2003. “Robust Portfolio Selection Problems”. Mathematics of Operations Research 28(1):1–38.
Hong, L. J., Z. Huang, and H. Lam. 2017. “Learning-Based Robust Optimization: Procedures and Statistical Guarantees”.
arXiv preprint arXiv:1704.04342.
Hu, Z., and L. J. Hong. 2013. “Kullback-Leibler Divergence Constrained Distributionally Robust Optimization”. Available at
Optimization Online.
Jiang, R., and Y. Guan. 2016. “Data-Driven Chance Constrained Stochastic Program”. Mathematical Programming 158(12):291–327.
Krokhmal, P., J. Palmquist, and S. Uryasev. 2002. “Portfolio Optimization with Conditional Value-At-Risk Objective andConstraints”. Journal of Risk 4(2):43–68.
Lam, H. 2016. “Robust Sensitivity Analysis for Stochastic Systems”. Mathematics of Operations Research 41(4):1248–1275.
Lam, H. 2019. “Recovering Best Statistical Guarantees via the Empirical Divergence-Based Distributionally Robust Optimization”.
Operations Research 67(4):1090–1105.
Lam, H., and E. Zhou. 2017. “The Empirical Likelihood Approach to Quantifying Uncertainty in Sample Average Approximation”.
Operations Research Letters 45(4):301–307.
Lifshits, M. A. 1995. Gaussian Random Functions, Volume 322 of Mathematics and Its Applications. Dordrecht, Netherlands:
Kluwer Academic Publisher.
Luedtke, J., and S. Ahmed. 2008. “A Sample Approximation Approach for Optimization with Probabilistic Constraints”. SIAM
Journal on Optimization 19(2):674–699.
Nemirovski, A. 2003. “On Tractable Approximations of Randomly Perturbed Convex Constraints”. In 42nd IEEE International
Conference on Decision and Control, Volume 3, 2419–2422. Piscataway, New Jersey: Institute of Electrical and Electronics
Engineers, Inc.
Nemirovski, A., and A. Shapiro. 2006. “Convex Approximations of Chance Constrained Programs”. SIAM Journal on
Optimization 17(4):969–996.
Petersen, I. R., M. R. James, and P. Dupuis. 2000. “Minimax Optimal Control of Stochastic Uncertain Systems with Relative
Entropy Constraints”. IEEE Transactions on Automatic Control 45(3):398–412.
Prékopa, A. 2003. “Probabilistic Programming”. In Stochastic Programming, Volume 10 of Handbooks in Operations Researchand Management Science, 267–351. New York: Elsevier.
Schildbach, G., L. Fagiano, and M. Morari. 2013. “Randomized Solutions to Convex Programs with Multiple Chance Constraints”.
SIAM Journal on Optimization 23(4):2479–2501.
Tulabandhula, T., and C. Rudin. 2014. “Robust Optimization Using Machine Learning for Uncertainty Sets”. arXiv preprint
arXiv:1407.1097.
Wand, M. P., and M. C. Jones. 1994. Kernel Smoothing. London, UK: Chapman and Hall/CRC.
Wang, W., and S. Ahmed. 2008. “Sample Average Approximation of Expected Value Constrained Stochastic Programs”.
Operations Research Letters 36(5):515–519.
Wiesemann, W., D. Kuhn, and M. Sim. 2014. “Distributionally Robust Convex Optimization”. Operations Research 62(6):1358–
1376.

AUTHOR BIOGRAPHIESHENRY LAM is an Associate Professor in the Department of Industrial Engineering and Operations Research at Columbia
University. His research focuses on Monte Carlo simulation, uncertainty quantification, risk analysis, and stochastic and robust
optimization. His email address is khl2114@columbia.edu.
HUAJIE QIAN is a Ph.D. student in the Department of Industrial Engineering and Operations Research at Columbia University.
His research interest lies in simulation uncertainty quantification, data-driven simulation analysis, and stochastic optimization.
His email address is hq2157@columbia.edu.