Proceedings of the 2014 Winter Simulation Conference
A. Tolk, S. D. Diallo, I. O. Ryzhov, L. Yilmaz, S. Buckley, and J. A. Miller, eds.

ROBUST RARE-EVENT PERFORMANCE ANALYSIS WITH NATURAL NON-CONVEXCONSTRAINTSJose Blanchet

Christopher Dolan

IEOR Department
Columbia University
340 W. Mudd Building
500 W 120th Street
New York, NY 10027, USA

Department of Statistics
Columbia University
1021 School of Social Work1255 Amsterdam Avenue
New York, NY 10027, USA

Henry Lam
Department of Mathematics and Statistics
Boston University
111 Cummington Mall
Boston, MA 02215, USA

ABSTRACT
We consider a common type of robust performance analysis that is formulated as maximizing an expectation
among all probability models that are within some tolerance of a baseline model in the Kullback-Leibler
sense. The solution of such concave program is tractable and provides an upper bound which is robust
to model misspecification. However, this robust formulation fails to preserve some natural stochastic
structures, such as i.i.d. model assumptions, and as a consequence, the upper bounds might be pessimistic.
Unfortunately, the introduction of i.i.d. assumptions as constraints renders the underlying optimization
problem very challenging to solve. We illustrate these phenomena in the rare event setting, and propose alarge-deviations based approach for solving this challenging problem in an asymptotic sense for a naturalclass of random walk problems.
1

INTRODUCTION

Robust performance analysis is concerned with the problem of evaluating the worst case performancemeasure of interest (typically described as an expectation) among all plausible probability models, such as
those within certain tolerance of a baseline model which is believed to be reflective of reality. Taken literally,
this problem formulation can be challenging because it gives rise to an infinite dimensional optimization
problem (note that we mentioned “all models” within certain tolerance). When the tolerance region is
described in terms of Kullback-Leibler divergence (and other related notions), this apparently daunting
optimization problem is often tractable, and this tractability feature has been exploited in a range of literaturein recent years, for example in control theory (Iyengar (2005), Nilim and El Ghaoui (2005), Petersen et al.
(2000)), distributionally robust optimization (Ben-Tal et al. 2013), finance (Glasserman and Xu 2014),
economics (Hansen and Sargent 2008) and queueing (Jain et al. 2010).
Tolerance regions based on the Kullback-Leibler divergence, however, fail to incorporate informationthat is often quite natural to assume in common stochastic settings, and that should be added in terms ofBlanchet, Dolan and Lam
constraints in the underlying robust performance analysis formulation. One such natural and important
constraint is the i.i.d. property, often arising in models involving random walk input. Failing to inform thei.i.d. property even in simple situations involving random walk models can have important consequences interms of the accurate assessment of worst case performance measures of interest. Unfortunately, however,
the i.i.d. property as an additional constraint on top of the Kullback-Leibler imposed tolerance gives riseto an optimization problem that is no longer easy to handle.
The main contribution of this paper is to show that in the context of performance analysis for performancemeasures associated to large deviations events, such robust formulation yields to a problem for which
asymptotically optimal solutions can be constructed. We illustrate our ideas in the setting of i.i.d. random
walks.
The rest of the paper is organized as follows. In Section 2 we provide a precise mathematical formulationof the robust performance analysis problem with i.i.d. constraints and explain why the problem is very
challenging. In Section 3 we provide a strategy that allows to solve this challenging problem asymptoticallyin a large deviations regime. Finally, we provide numerical examples which illustrate the performance ofour proposed solution and the impact of adding i.i.d. constraints in the robust formulation.
2

PROBLEM FORMULATION

Let {Xk : k ≥ 0} be a sequence of zero mean i.i.d. random variables. Define S0 = 0 and put Sn = X1 + ... + Xn .
Let us use F (·) to denote the CDF (Cumulative Distribution Function) of Xi , that is, P (Xi ≤ x) = F (x)and we use PF (·) to denote the product measure generated by F (·). We use PFn (·) to denote the projectionof PF (·) into its n first coordinates. Simply put, PFn describes the joint distribution of the random variables(X1 , ..., Xn ). The expectation operator associated to PF (·) and PFn (·) is denoted by EF (·) and EFn (·),
respectively. We define ψF (θ ) = log EF1 exp (θ X1 ) and assume that ψF (θ ) < ∞ for θ in a neighborhoodof the origin.
Now, define An = {Sn /n ∈ A} for a closed set A which does not contain the origin. We are concerned withthe problem of estimating PF (An ). Observe that PF (An ) → 0 as n → ∞ because of the law of large numbers.
Moreover, because ψF (·) is finite in a neighborhood of the origin we have that PF (An ) ≤ exp (−δ n) forsome δ > 0 for all n sufficiently large.
In contrast to standard rare event estimation problems, however, here we assume that F (·) is unknown.
Nevertheless, based on some evidence (for example based on data or expert knowledge) let us assume
that we have obtained a CDF G (·), which approximates F (·) in a suitable sense, for example in the
Kullback-Leibler divergence sense which we shall review momentarily. Let us write PG (·) to denote theproduct measure associated to G (·) and we use EG (·) for the expectation operator corresponding to PG (·).
Similarly as before, PGn (·) is the projection of PG (·) onto its n first coordinates and we use EGn (·) to denotethe expectation operator associated to PGn (·).
We assume that the likelihood ratio dPFn /dPGn is well defined and therefore the Kullback-Leibler of PFn
with respect to PGn is defined via
R (PFn ||PGn ) = EFn log



dPFn
dPGn



= nEF1 log






Z
dPF1
dF
(X1 ) = n log(x) dF (x) .
dG
dPG1

If dPFn /dPGn fails to exist (i.e. PFn is not absolutely continuous with respect to PGn ), then the Kullback-Leibler
divergence is defined as infinity. It is elementary to verify that R(· ||PGn ) is convex (actually R (·||·) is convex
in both of its arguments (Dupuis and Ellis 1997)). The associated robust performance analysis problem
with Kullback-Leibler constraint consists in solvingmax{Qn (An ) : R (Qn ||PGn ) ≤ ηn },
n
Q

(1)

Blanchet, Dolan and Lam
where ηn should be chosen to satisfy
ηn ≈ nZdF(x) dF (x) .
log
dG

R

One might select ηn by estimating log (dF/dG (x)) dF (x) using available data.
The optimization problem (1) is a concave mathematical program; the objective function to maximize
is linear (in particular concave) in the variable Qn and, as mentioned earlier, the constraint is convex.
Moreover, as we shall see in the body of the paper, the optimal solution to (1), Qn∗ (·), can be characterized
as a suitable mixture between PGn (·|An ) and PGn (·|Acn ). As the next result shows, it turns out that Qn∗ (An )might differ substantially from PFn (An ) even if F is close to G in the Kullback-Leibler sense, that is, even incases in which ηn = o (n) as n → ∞. In more detail, typically we will have PFn (An ) = exp (−nI¯F (A) + o (n)),
for some positive constant I¯F (A), whereas the next result indicates that typically Qn∗ (An ) ≥ δ ηn /n for someδ > 0 and large enough n. So, for example, if one builds an approximation G to F from data, one would
need an exponentially large sample size (in n) in order to obtain an accurate estimate of the probability ofinterest using only the relative entropy constraint without recognizing that the data might have come froman i.i.d. model.
Theorem 1 Suppose that ηn = o(n) and that ηn > δ > 0 for some δ > 0 uniformly over n. Assume alsothat PG (An ) ∈ (exp (−n/δ ) , exp (−δ n)) for some δ > 0 and all n sufficiently large. Then the optimal valueof (1), Qn∗ (An ), satisfies
ηn(1 + o(1))
Qn∗ (An ) =− log PGn (An )as n → ∞.
One of the main reasons for such a disparity, as we shall establish in the next section, is that the feasible
region (i.e. {Qn : R Qn ||PGn ≤ ηn }) fails to recognize that we are interested only in models for which thei.i.d. property of the Xi ’s is preserved. So, introducing the i.i.d. constraint transforms problem (1) into thealternative form





Z
Z
Z 
dH
x1 + ... + xnn
∈ A dH (x1 ) ...dH (xn ) : n log(x) dH (x) ≤ ηn . (2)max PH (An ) = ... I
HndG
Observe that the previous problem is not a concave program because the objective function to maximize
is no longer concave. Unfortunately, in general (2) is very challenging to solve. In the next secion we
explain how to use large deviations theory to solve problem (2) in an asymptotic sense. We finish thissection with a proof of our first theorem.
2.1 Proof of Theorem 1
Proof.
asTo solve (1), we rewrite it in terms of the likelihood ratio between Qn and Gn , namely L = dQn /dPGn ,max
EGn [L; An ]subject to EGn [L log L] ≤ ηn ,(3)

where the maximization is over L ∈ L = {L ≥ 0 : EGn L = 1} and consider the Lagrangian relaxation
max EGn [L; An ] − α(EGn [L log L] − ηn ).
L∈L

(4)

Our goal is to find α ∗ ≥ 0 such that there is an L∗ that solves (1) and moreover that EGn [L∗ log L∗ ] = ηn .
Then this L∗ will be optimal for (1) (c.f. p.220, Theorem 1 in Luenberger (1997)). Note that when α = 0,
the optimal value to the maximization in (4) is 1, by choosing L to concentrate only on the event A. So theBlanchet, Dolan and Lamcase α ∗ = 0 is easily ruled out. Now, given any fixed α > 0, it can be verified by a convexity argumentthat the solution to the maximization (1) is given byL∗ ∝ eI(An )/α = e1/α I(An ),(5)

where I (An ) denotes the indicator function of the set An ; see, for example, Hansen and Sargent (2008).
Now we write αn for α to highlight the role of n, and introduce βn = 1/αn for convenience. We also write
pn = PG (An ) and put qn = PG Ān = 1 − pn . Then (5) can be written as(eβnβn +q on An
∗
p
e
n
n
L =
.
(6)
1
on Acnp eβn +qnn

We now proceed to find αn∗ > 0, or βn∗ = 1/αn∗ , such thatEGn [L∗ log L∗ ] = ηn .

(7)

pn eβn− log(pn eβn + qn ) = ηn .
pn eβn + qn

(8)

Using the form of (6), (7) becomes
βn

Since ηn > δ > 0 and pn → 0 as n → ∞, we must have that all βn satisfying (8) must also satisfy βn → ∞as n → ∞. Otherwise the left hand side converges to zero on some subsequence while the right hand sidestays positively bounded away from zero. Now, we claim that lim sup pn eβn = 0. Let us proceed assuming
this claim for the moment and come back to this issue at the end of our proof. Then, by a Taylor series
expansion applied to the left hand side of (8), we have that
βn pn eβn (1 + o(1)) = ηn ,(9)which gives
log βn + log pn + βn + o(1) = log ηn .
Heuristically, we must have
βn = log ηn − log pn − log βn + o(1) = (log ηn − log pn )(1 + o(1)).

(10)

To verify (10) rigorously, note that when we choose βn = log (ηn /pn ), the left hand side of (9) becomesηn (log(ηn /pn ))(1 + o(1)) which is much more larger than ηn for n large enough. On the other hand, setting
βn = 0 gives the left hand side o(1). Therefore, by continuity there must be a solution to (9) in the range[0, log(ηn /pn )]. Consequently, we have thatβn = log ηn − log pn − log βn + o(1) = log ηn − log pn + rn

(11)where the remainder term rn satisfies |rn | ≤ log(log (ηn /pn )) + o(1), or equivalently we obtain that rn =o(log (ηn /pn )), and hence (10).
Iterating the first equality in (10) using (11), we get further thatβn = log ηn − log pn − log(log ηn − log pn + rn ) + o(1) = log ηn − log pn − log(log ηn − log pn ) + o(1).
Finally, the optimal value is
EGn [L; An ] =ηnηnpn eβn∼ pn eβn =(1 + o(1)) ∼.
βnlog ηn − log pn− log pn
pn e + qn

Now we must verify that indeed lim sup pn eβn = 0. Suppose that lim sup pn eβn > 0, this implies thatβnk ≥ δ nk along a subsequence nk → ∞ for δ > 0. But then we must have from (8) that ηnk ≥ δ 0 nk forsome δ 0 > 0, contradicting our assumption that ηn = o (n). We therefore conclude the statement of ourtheorem.

Blanchet, Dolan and Lam
3

OUR MAIN RESULT

3.1 A Large Deviations Rate Characterization
In order to prove our main result we shall impose additional technical conditions. We assume that ψG (·) is
steep in the sense that for all a ∈ (−∞, ∞) there is θa such that ψG0 (θa ) = a. Under this assumption we
have that
PGn (An ) = exp (−nI¯G (A) + o (n)) ,where
I¯G (A) = inf IG (x) = inf sup (θ x − ψG (θ )) .
x∈A

x∈A θ

Similarly, for each H, we define ψH (θ ) := log EH1 exp (θ X) and define IH (x) = supθ (θ x − ψH (θ )). Similarly
as for the definition of I¯G (A), we write I¯H (A) = infx∈A IH (x). Now we are ready to state our main result.
Theorem 2 Let int (A) denote the interior of A. Suppose that I¯H (A) = I¯H (int (A)) ∈ (0, ∞) for any PH ∈ P
for some feasible set P. We have1max
PH (A) = − min I¯H (A)lim logn→∞ n
i.i.d.
PH ∈P
PH ∈P,Xi ∼ H
Proof.

The proof follows from a LD argument. First,
11
lim inf log
maxlog PH (An )PH (An ) = lim inf
max
n→∞ n
n→∞
i.i.d.
i.i.d.
PH ∈P,Xi ∼ PH
PH ∈P,Xi ∼ PH n≥

1
lim inf log PH (An )n→∞ ni.i.d.
PH ∈P,Xi ∼ PHmax

≥ − min I¯H (int (A)) = − min I¯H (A).
PH ∈P

PH ∈P

Next, since A is closed, using Chebycheff inequality (as in the proof of Cramer’s theorem; see p.27,
Remark (c) in Dembo and Zeitouni (1998)) gives
PH (An ) ≤ 2 exp (−nI¯H (A))
and hence

1max
PH (An ) ≤ − min I¯H (A).
lim sup log
i.i.d.
PH ∈Pn→∞ n
PH ∈P,Xi ∼ PH

Combining the upper and lower bounds we get our conclusion.
The significance of the previous result is that the optimization problem that must be solved can now becast as a concave program and thus it is more tractable. In order to have a concrete class of examples, let
us focus on the case in which A = [a, ∞). Really the key property that holds using this specific selection isthat we can identify a specific element a ∈ A such that IH (a) = I¯H (A). We consider the general problemof finding the minimum rate function over a class of distributions, namely
inf IH (a) = inf sup{θ a − ψH (θ )}.

PH ∈P

PH ∈P θ(12)

Lemma 1 If P is a convex set, then the optimization program (12) is convex.
Proof.
The inner objective function in (12) is concave as a function of θ , since ψH (·) is convex. For
the outer objective function, note that ψH (θ ) is concave in PH , and so θ a − ψH (θ ) is convex in PH . As the
maximum over a set of convex functions (indexed by θ ), the outer objective function supθ {θ a − ψH (θ )}is also convex as a function of PH . Therefore both the inner and outer optimizations in (12) are convex
programs.

Blanchet, Dolan and Lam
3.2 Numerical Procedure
For numerical procedure, in order to avoid the issue of optimization over infinite-dimensional variables,
we concentrate on the case of discrete PH . For convenience, we write p = (p1 , . . . , pm ) as the weights oversupport points {x1 , . . . , xm }. Moreover, we write()mZ(p) = max θ a − log ∑ pi eθ xi
θ

(13)

i=1

as the outer objective in (12). Note that for discrete distribution p, the optimal solution of θ in (13) canbe solved simply by finding the root ofθ xi∑mi=1 pi xi e
=a
θ xi
∑m
i=1 pi e

We now focus on minp∈P Z(p). Suppose that P is a feasible region dictated by KL-divergence constraint,
m
0
0
0
0
i.e. P = {p ≥ 0 : ∑m
i=1 pi log(pi /pi ) ≤ η, ∑i=1 pi = 1} for some baseline distribution p = (p1 , . . . , pm ).
A particularly convenient procedure to approximate the optimal solution is to use conditional gradient (orFrank-Wolfe) method (Frank and Wolfe 1956). This lies on the stepwise optimization, given the current
solution pk = (pk1 , . . . , pkm ) at step k,(14)
min ∇Z(pk )(p − pk )p∈PThis subroutine can be easily solved. In fact, we have

!
m
k
d
log ∑ p j eθ x j∇Z(pk ) = −
d pi
j=1



k

 =
pi =pki

i

eθ xi
− m
k
∑ j=1 p j eθ x j!
i

by simple arithmetic or by the use of the envelope theorem, where θ k is the solution to
k
θ xi
∑m
i=1 pi xi e=ak θ xi∑m
i=1 pi e

For convenience, we letkeθ xiξi (p ) = − m k θ k x∑ j=1 p j e jkbe the i-th coordinate of ∇Z(pk ). The solution to (14) is given by qk+1 = (q1k+1 , . . . , qk+1
m ), wherekqik+1 =

p0i eβ ξi (p )k∑mj=1 p0j eβ ξ j (p )and β < 0 satisfies the equationkm0
k β ξi (p )k∑mi=1 β pi ξi (p )e−log
p0j eβ ξ j (p ) = η
∑
k)
m
0
β
ξ
(p
∑ j=1 p j e j
j=1

If there is no negative root to this equation, then qk+1 is plainly a degenerate mass on argmin{∇Z(pk )}.
Therefore, the iterative procedure is the following:
Iterative Procedure: Start from the baseline distribution p0 (or any other distribution). At each iteration k,
given pk , do the following:

Blanchet, Dolan and Lam
1. Compute the root θ that solves
k
θ xi
∑m
i=1 pi xi e=ak θ xi∑m
i=1 pi e

2. Computeξi = −eθ xi∑mj=1 pkj eθ x j

for i = 1, . . . , m.
3. Compute, if any, the negative root of
m
0
β ξi
∑m
i=1 β pi ξi e−
log
∑ p0j eβ ξ j = η∑mj=1 p0j eβ ξ j
j=1

4. If there is a negative root β , thenqik+1 =

p0i eβ ξi∑mj=1 p0j eβ ξ j

for i = 1, . . . , m. Otherwise qk+1
= 1 if i = argmin{ξi }, and 0 for all other i’s.
i
k+1
k+1
k
k+1
5. Update p= (1 − ε )p + ε qk+1 for some step size ε k+1 (by using limited minimization
rule or Armijo Rule; see p.217, Bertsekas (1999)).
4

NUMERICAL EXPERIMENTS

We will apply our algorithm in the case of two standard discrete distributions with finite support, namely
the binomial distribution and a discrete distribution with random weights. We compare the performanceof the upper bound with i.i.d. constraints with the bound without i.i.d constraints. Figure ?? shows the
log-probabilities associated with a binomial model with parameters m = 10, p = 0.5 as n increases. The truemodel is assumed to be binomial with m = 10, p = 0.55. Figure ?? shows the log-probabilities associatedwith a discrete distribution with random weights on the integer support from 1 to 10, using MLE froma random sample as the base model. As can be seen, in both cases the upper bound with i.i.d constraint
provides a much tighter bound to the real model than otherwise.
Binomial Model
−1
−2
−3

log(probability)

−4
−5
−6
−7
−8
Upper Bound i.i.d constraints
Real Model
Base Model
Upper Bound no constraints

−9
−10
−11

3

4

5

6

7
n

Figure 1: Binomial

8

9

10

Blanchet, Dolan and Lam
Random Weight Model
0−5log(probability)

−10

−15

−20
Upper bound i.i.d constraints
Real Model
Base Model
Upper bound no constraints

−25

−30

5

10

15

20

25

30

35

40

n

Figure 2: Random Weights
REFERENCES
Ben-Tal, A., D. Den Hertog, A. De Waegenaere, B. Melenberg, and G. Rennen. 2013. “Robust solutionsof optimization problems affected by uncertain probabilities”. Management Science 59 (2): 341–357.
Bertsekas, D. P. 1999. “Nonlinear programming”.
Dembo, A., and O. Zeitouni. 1998. Large deviations techniques and applications, Volume 2. Springer.
Dupuis, P., and R. Ellis. 1997. A Weak Convergence Approach to the Theory of Large Deviations. New
York: John Wiley and Sons.
Frank, M., and P. Wolfe. 1956. “An algorithm for quadratic programming”. Naval Research Logistics
Quarterly 3:95–110.
Glasserman, P., and X. Xu. 2014. “Robust risk measurement and model risk”. Quantitative Finance 14 (1):29–58.
Hansen, L., and T. Sargent. 2008. Robustness. Princeton, NJ: Princeton University Press.
Iyengar, G. N. 2005. “Robust dynamic programming”. Mathematics of Operations Research 30 (2): 257–280.
Jain, A., A. E. Lim, and J. G. Shanthikumar. 2010. “On the optimality of threshold control in queues withmodel uncertainty”. Queueing Systems 65 (2): 157–174.
Luenberger, D. 1997. Optimization by Vector Space Methods. New York: John Wiley and Sons.
Nilim, A., and L. El Ghaoui. 2005. “Robust control of Markov decision processes with uncertain transition
matrices”. Operations Research 53 (5): 780–798.
Petersen, I. R., M. R. James, and P. Dupuis. 2000. “Minimax optimal control of stochastic uncertain systems
with relative entropy constraints”. Automatic Control, IEEE Transactions on 45 (3): 398–412.
AUTHOR BIOGRAPHIESJOSE BLANCHET is an Associate Professor in the Departments of IEOR and Statistics at Columbia
University. Jose holds a Ph.D. in Management Science and Engineering from Stanford University. Prior
to joining Columbia he was a faculty member in the Statistics Department at Harvard University. Jose isa recipient of the 2009 Best Publication Award given by the INFORMS Applied Probability Society andof the 2010 Erlang Prize. He also received a PECASE award given by NSF in 2010. He worked as ananalyst in Protego Financial Advisors, a leading investment bank in Mexico. He has research interests inapplied probability and Monte Carlo methods. He serves in the editorial board of Advances in Applied
Probability, Journal of Applied Probability, Mathematics of Operations Research and QUESTA. His email is .

Blanchet, Dolan and Lam
CHRISTOPHER DOLAN is a PhD student in the Department of Statistics at Columbia University. He
received his BA and MS in Mathematics from the Courant Institute at NYU. His research interests are inapplied probability and rare event simulation. His email address is dolan@stat.columbia.edu.
HENRY LAM is an Assistant Professor in the Department of Mathematics and Statistics at Boston University.
He graduated from Harvard University with a Ph.D. degree in statistics in 2011. His research focuses onlarge-scale stochastic simulation, rare-event analysis, and simulation optimization, with application interests
in service systems and risk management. His email is khlam@bu.edu.