The Empirical Likelihood Approach to Quantifying Uncertainty in Sample Average Approximation\footnote{A preliminary conference version of this work has appeared in \cite{lam2015quantifying
Henry Lam\thanks{Department of Industrial and Operations Engineering
University of Michigan
Ann Arbor
MI. Email: {\tt khlam@umich.edu
We study the empirical likelihood approach to construct confidence intervals for the optimal value and the optimality gap of a given solution, henceforth quantify the statistical uncertainty of sample average approximation, for optimization problems with expected value objectives and constraints where the underlying probability distributions are observed via limited data. This approach 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. %Numerical experiments show that our approach compares favorably with central limit asymptotics in terms of finite-sample performance. \\ \noindentKeywords: empirical likelihood; sample average approximation; confidence interval; constrained optimization; stochastic program; statistical uncertainty %We consider stochastic optimization with expected value objectives and constraints, where the underlying probability distributions are observed via limited data. Sample average approximation (SAA) is a common approach to obtain approximate solutions in such situations. In this paper, we study a new approach to construct confidence intervals for the optimal values, henceforth quantify the statistical uncertainty of SAA, by using the empirical likelihood method in statistics. This approach casts the task of confidence interval construction as a pair of distributionally robust optimizations posited over the uncertain distribution, with a divergence-based uncertainty set that is suitably calibrated to provide asymptotic statistical guarantee. %Numerical experiments show that our approach compares favorably with central limit asymptotics in terms of finite-sample performance.
ELuncertainty2_Oct17.tex

\documentclass[11pt]{article}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{amsmath,amsfonts,amsthm,amssymb,graphicx,bm,color}
%\usepackage{subfigure}
%\usepackage[pdftex]{graphicx}
%\usepackage{fullpage}

\setcounter{MaxMatrixCols}{10}
%TCIDATA{OutputFilter=LATEX.DLL}
%TCIDATA{Version=5.50.0.2960}
%TCIDATA{<META NAME="SaveForMode" CONTENT="1">}
%TCIDATA{BibliographyScheme=Manual}
%TCIDATA{LastRevised=Monday, June 28, 2010 19:04:16}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}

\theoremstyle{plain}
\newtheorem{thm}{Theorem}
\newtheorem{proposition}{Proposition}
\newtheorem{lemma}{Lemma}
\newtheorem{corollary}{Corollary}
\newtheorem{assumption}{Assumption}
\newtheorem{remark}{Remark}
\newtheorem{example}{Example}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\cP}{{\mathcal P}}
\newcommand{\p}{\partial}
\newcommand{\vc}[1]{\mathbf{#1}}
\newcommand{\ept}[1]{\left\langle #1 \right\rangle}
\newcommand{\abs}[1]{\left|#1\right|}

\newcommand{\blue}{\textcolor{blue}}
%\addtolength{\textheight}{.4in}
%\renewcommand{\baselinestretch}{1.35}
\textwidth 6.52in \oddsidemargin -.01in \evensidemargin -.01in
\textheight 8.8in \topmargin -.2in \pagestyle{myheadings}

%\input{tcilatex}
\begin{document}

\title{The Empirical Likelihood Approach to Quantifying Uncertainty in Sample Average Approximation\footnote{A preliminary conference version of this work has appeared in \cite{lam2015quantifying}. 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.}}
\author{Henry Lam\thanks{Department of Industrial and Operations Engineering, University of Michigan, Ann Arbor, MI. Email: {\tt khlam@umich.edu}}
\and
Enlu Zhou\thanks{H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA. Email: {\tt enlu.zhou@isye.gatech.edu}}}
\date{}
\maketitle

\begin{abstract}
We study the empirical likelihood approach to construct confidence intervals for the optimal value and the optimality gap of a given solution, henceforth quantify the statistical uncertainty of sample average approximation, for optimization problems with expected value objectives and constraints where the underlying probability distributions are observed via limited data. This approach 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. %Numerical experiments show that our approach compares favorably with central limit asymptotics in terms of finite-sample performance.
\\

\noindent\emph{Keywords:} empirical likelihood; sample average approximation; confidence interval; constrained optimization; stochastic program; statistical uncertainty

%We consider stochastic optimization with expected value objectives and constraints, where the underlying probability distributions are observed via limited data. Sample average approximation (SAA) is a common approach to obtain approximate solutions in such situations. In this paper, we study a new approach to construct confidence intervals for the optimal values, henceforth quantify the statistical uncertainty of SAA, by using the empirical likelihood method in statistics. This approach casts the task of confidence interval construction as a pair of distributionally robust optimizations posited over the uncertain distribution, with a divergence-based uncertainty set that is suitably calibrated to provide asymptotic statistical guarantee. %Numerical experiments show that our approach compares favorably with central limit asymptotics in terms of finite-sample performance.
\end{abstract}

\section{Introduction}
We consider a stochastic optimization problem in the form
\begin{equation} \label{obj}
\min_{x\in\Theta}\{h(x):=E[H(x;\xi)]\},
\end{equation}
where $x=(x_1,\ldots,x_p)$ is a continuous decision variable in the deterministic feasible region $\Theta\subseteq\mathbb{R}^p$, and $\xi$ is a random vector on $\mathbb{R}^d$. We are interested in situations where the underlying probability distribution that controls the expectation $E[\cdot]$ is not fully known and can only be accessed via limited data $\xi_1,\ldots,\xi_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., \cite{shapiro2014lectures}):
\begin{equation} \label{obj empirical}
\min_{x\in\Theta}\frac{1}{n}\sum_{i=1}^nH(x;\xi_i).
\end{equation}

We further consider problems with expected value constraints, in the form
\begin{equation} \label{obj constrained}
\begin{array}{ll}\min&h(x)=E[H(x;\xi)]\\
\text{subject to}&f_k(x)=E[F_k(x;\xi)]\leq0,\ k=1,\ldots,m\\
&g_k(x)\leq0,\ k=1,\ldots,s
\end{array}
\end{equation}
where $g_k(\cdot)$'s are deterministic functions. Thus \eqref{obj constrained} can include both stochastic and deterministic constraints. Again, under limited data $\xi_1,\ldots,\xi_n$, an SAA version of \eqref{obj constrained} is (e.g., \cite{wang2008sample})
\begin{equation} \label{obj constrained empirical}
\begin{array}{ll}\min&\frac{1}{n}\sum_{i=1}^nH(x;\xi_i)\\
\text{subject to}&\frac{1}{n}\sum_{i=1}^nF_k(x;\xi_i)\leq0,\ k=1,\ldots,m\\
&g_k(x)\leq0,\ k=1,\ldots,s
\end{array}
\end{equation}

Our premise is that beyond the $n$ observations, new samples are not easily accessible because of either a lack of data or limited computational capacity in running further Monte Carlo simulation. The optimal value and solution obtained from \eqref{obj empirical} and \eqref{obj constrained empirical} thus deviate from those under the genuine distribution in \eqref{obj} and \eqref{obj constrained}. Moreover, the error of the solution implies a non-zero optimality gap with the true optimal value, resulting in suboptimal decisions. Estimating these errors is important and has been studied over the years (e.g., \cite{kleywegt2002sample}, \cite{mak1999}, Chapter 5 in \cite{shapiro2014lectures}).

%The focus of this paper is to quantify these uncertainties and estimate these errors.

%The focus of this paper is to quantify the uncertainty in SAA by constructing confidence intervals (CIs) for the true optimal values. Moreover, we will show how the technique we develop can be similarly applied to estimate the optimality gap for a given solution.

Our main contribution is to bring in a new approach to rigorously quantify the uncertainty in \eqref{obj empirical} and \eqref{obj constrained empirical} 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 numerical examples 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 \ref{EL section} and \ref{optimality gap section} study the theory of our approach applied to the optimal value and the optimality gap, and Section \ref{sec:numerics} shows some numerical results and comparison with previous methods.

%as an alternative to directly using the CLT, that has asymptotic statistical guarantees and could possess better finite-sample performance.
%(REWRITE) To demonstrate our proposed approaches, we consider conditional value-at-risk (CVaR) as a specific example of \eqref{obj}. CVaR is widely used in portfolio risk management in finance and insurance, among other applications \cite{rockafellar2000optimization}. The distributions of these portfolios are unknown and can only be represented via data in any real-life contexts. We shall provide some numerics on quantifying the uncertainty of the empirical optimization of CVaR.

%In the rest of this paper, Section \ref{EL section} will


%on a new approach to quantify these deviations in the form of confidence bounds, via the so-called empirical likelihood method in statistics.

  %studies a new method for quantifying the uncertainty of the optimal value from sample average approximation in the form of confidence bounds. Obviously, \eqref{obj empirical} is subject to approximation errors with respect to \eqref{obj}, and this error is the main focus of this paper.

%In this setup, previous results on quantifying the approximation errors have fallen in two directions. The first is central limit convergence of the optimality gap, i.e. difference between \eqref{obj} and \eqref{obj empirical}, typically to some Gaussian random variable under suitable scaling (\cite{shapiro2014lectures}; \cite{kim2015guide}). Such results can be used to construct confidence interval (CI) that has asymptotically correct coverage probability as sample size increases. %, and typically require the condition of Lipschitz continuity of $h(X;\cdot)$ (REF).
%The second approach consists of sampling complexity results, in the form of concentration inequalities on the optimality gap.
%%Previous work on this topic mostly surrounds sampling complexity results, in the form of concentration inequalities on the optimality gap, i.e. difference between \eqref{obj} and \eqref{obj empirical}, or central limit convergence.
%For instance, \cite{shapiro2005complexity} present large deviations inequalities that depend on certain size measures of $\Theta$ and the Lipschitz constant of $h$. In classification-type learning problems, where $h$ is often taken as an indicator function, deviation bounds have been obtained based on the entropy or the so-called Vapnik-Chervonenkis dimension of $h$ (e.g. \cite{mendelson2003few}).
%
%%\cite{kim2015guide} reviews some convergence results on the empirical optimal value, including convergence in distribution of the optimality gap to some function of Gaussian process. All these results serve as valuable guidelines for choosing the sample size $n$; however, they are either of worst-case nature (hence can be loose), or depend on parameters that are hard to estimate.
%
%The above two approaches are both valuable in evaluating the reliability of an adopted solution or choosing a good sample size $n$ in advance. However, they have antagonistic limitations: Central limit convergence typically only works well for the large-sample case and has no guarantee for smaller sample size, while sampling complexity results are typically loose because of their worst-case nature, even though they hold for any sample size. Driven by these limitations, we are interested in exploring techniques that can potentially yield improvement on both ends, by giving tight quantification of approximation errors even in situations of limited samples. % and have the potential to work for $h$ that have minimal regularity requirement. This is motivated by the limitations of the two existing approaches that either only have asymptotic guarantee, or are
%
%%Our focus is finite-sample performance, take a different, post-hoc, perspective in quantifying the errors made by \eqref{obj empirical}. In particular, we investigate numerical schemes that can generate accurate statistical estimates of the optimality gap or the true optimal value, given the samples $X_i$'s have already been collected. For modelers, this sort of estimates would be useful to assess the quality and the confidence in using an adopted solution.
%%We consider three distinct approaches: the empirical likelihood (EL) method, the Dirichlet method, and nonparametric bootstrapping. These approaches differ in  statistical philosophy and have distint benefits and disadvantages in terms of both the uncertainty assessment scope and computational load. To give some initial highlight, the EL method can be viewed as a nonparametric analog of classical likelihood inference, and can generate CI for the true optimal value or the optimality gap in the frequentist sense. %It is known to be competitive against plain use of normal asymptotic theory in terms of the CI's coverage error under certain finite-sample situations.
%%On the other hand, the Dirichlet method is utilized from a Bayesian perspective to approximate a corresponding posterior distribution, whereas nonparametric bootstrapping can be viewed as a limit of the Dirichlet method by taking suitable asymptotic on the prior distribution.%, and also possesses certain frequentist properties.



%Our paper is organized as follows. In Section \ref{methods}, we first present the three approaches for quantifying uncertainty of the optimal value output by \eqref{obj empirical}. We then adapt these approaches for optimality gap in Section \ref{methods gap}. We discuss some comparisons among the approaches in Section \ref{sec:comparisons}, followed by some numerical results on CVaR in Section \ref{numerics}.

%\section{METHODS FOR QUANTIFYING UNCERTAINTY OF THE OPTIMAL VALUE}\label{methods}

%We will introduce three methods (empirical likelihood method, Dirichlet method, and bootstrap sampling method) and discuss their relations.
\section{The Empirical Likelihood Method for Constructing Confidence Bounds for Optimal Values}\label{EL section}
%\subsubsection{An Optimization Formulation}
This section studies in detail the EL method in constructing CIs for the optimal values. Section \ref{sec:unconstrained} focuses on \eqref{obj} that only has deterministic constraints, and Section \ref{sec:constrained} generalizes to the stochastically constrained case \eqref{obj constrained}.

\subsection{Deterministically Constrained Optimization}\label{sec:unconstrained}
Let us first fix some notations. Given the set of i.i.d. data $\xi_1,\xi_2,\ldots,\xi_n$, we denote a probability vector over $\{\xi_1,\ldots,\xi_n\}$ as $w=(w_1,\ldots,w_n)\in\mathbb R^n$, where $\sum_{i=1}^nw_i=1$ and $w_i\geq0$ for all $i=1,\ldots,n$. We denote $\chi^2_{q,\beta}$ as the $1-\beta$ quantile of a $\chi^2$ distribution with degree of freedom $q$. We use ``$\Rightarrow$" to denote convergence in distribution, and ``a.s." to denote ``almost surely".

Our method utilizes the optimization problems
\begin{equation}
\begin{array}{ll}
\max/\min_w&\min_{x\in\Theta}\sum_{i=1}^nw_iH(x;\xi_i)\\
\text{subject to}&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{p+1,\beta}\\
&\sum_{i=1}^nw_i=1\\
&w_i\geq0\text{\ for all\ }i=1,\ldots,n
\end{array}\label{EL opt}
\end{equation}
where ``$\max/\min$" denotes a pair of maximization and minimization. Note that the optimal value of the SAA problem \eqref{obj empirical} lies between those of \eqref{EL opt}.%, contain the true value of $\max_{x\in\Theta}H(x)$ with probability at least $1-\beta$ asymptotically.

The quantity $-(1/n)\sum_{i=1}^n\log(nw_i)$ can be interpreted as the Burg-entropy divergence (\cite{pardo2005statistical}, \cite{ben2013robust}) between the probability distributions represented by the weights $w$ and by the uniform weights $(1/n)_{i=1,\ldots,n}$ on the support $\{\xi_1,\ldots,\xi_n\}$. %In fact, this distance belongs to the class of so-called $\phi$-divergence, where the $\phi$ function can be identified as $\phi(x)=-\log x$ here \cite{pardo2005statistical}.
Thus, the first constraint in \eqref{EL opt} is a Burg-entropy divergence ball centered at the uniform weights, with radius $\chi^2_{p+1,\beta}/(2n)$. From the viewpoint of DRO (e.g., \cite{delage2010distributionally,ben2013robust,wiesemann2014distributionally}), the optimization problems in \eqref{EL opt} output the worst-case estimates of $\min_{x\in\Theta}\{h(x)=E[H(x;\xi)]\}$ when $E[\cdot]$ is uncertain and its underlying distribution is believed to lie inside the divergence ball. We should point out, however, that this DRO interpretation differs from those in the existing literature (e.g., \cite{bertsimas2016}), as our divergence ball (i.e. the ``uncertainty set" in the terminology of robust optimization) may have low coverage of the true distribution $P$. This can be seen particularly when $P$ is a continuous distribution, in which case the coverage of the divergence ball is zero because of the violation of the absolute continuity requirement needed in properly defining the divergence.

%when the underlying probability distribution of a stochastic problem is not fully known, one can impose an uncertainty set in which the probability distribution is believed to lie, and calculate the worst-case scenario among all distributions within the set. In our case, \eqref{EL opt} is computing the worst-case scenarios as the max or min for the quantity $\max_{x\in\Theta}E[h(x;\xi)]$, where the uncertainty set can be thought as the neighborhood around the uniform distribution on $\{\xi_1,\ldots,\xi_n\}$ measured by a $\phi$-divergence with $\phi(x)=-\log x$.


%\subsection{Statistical Guarantees}
The EL method is a mechanism to endow statistical meaning to \eqref{EL opt}. In particular, it asserts that using the ball size $\chi^2_{p+1,\beta}/(2n)$ in \eqref{EL opt} gives rise to statistically valid $1-\beta$ confidence bounds for the optimal value of \eqref{obj} (despite that the ball may under-cover the true distribution). This method originates as a nonparametric analog of maximum likelihood estimation first proposed by \cite{owen1988empirical}. On the data set $\{\xi_1,\ldots,\xi_n\}$, we first define a ``nonparametric likelihood" $\prod_{i=1}^nw_i$, where $w_i$ is a probability weight applied to each datum. It is straightforward to see that the maximum value of $\prod_{i=1}^nw_i$, among all $w$ in the probability simplex, is $\prod_{i=1}^n(1/n)$. In fact, the same conclusion holds even if one allows putting weights outside the support of the data, which could only make the likelihood $\prod_{i=1}^nw_i$ smaller. In this sense, $\prod_{i=1}^n(1/n)$ can be viewed as a maximum likelihood in the nonparametric space. Correspondingly, we define the nonparametric likelihood ratio between the weights $w$ and the maximum likelihood weights as $\prod_{i=1}^nw_i/\prod_{i=1}^n(1/n)=\prod_{i=1}^n(nw_i)$. %While Wilks' Theorem stipulates that such ratio

%The statistical guarantee for \eqref{EL opt} and the choice of the ``neighborhood size" $\chi^2_{p+1,\beta}$ can be made rigorous through the theory of the EL method. The key of the method is an analog of the celebrated Wilks' Theorem \cite{cox1979theoretical} in parametric likelihood inference, namely that the ratio between the maximum likelihood and the true likelihood, i.e. the so-called likelihood ratio, converges to a chi-square distribution in a suitable logarithmic scale.


The key of the EL method is a nonparametric counterpart of the celebrated Wilks' Theorem \cite{wilks1938large} 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 $\chi^2$-distribution in a suitable logarithmic scale. To develop this analog, we first incorporate a target parameter of interest, i.e. the quantity whose statistical uncertainty is to be assessed (or to be ``estimated"). Say this parameter is $\theta\in\mathbb R^p$. %In our case this is $\max_{x\in\Theta}E[h(x;\xi)]$. The EL method bridges to Wilks' Theorem via a quantity known as the nonparametric profile likelihood ratio, defined as the maximum ratio between the empirical likelihood and the nonparametric maximum likelihood, optimized over all probability weights with support over the data set and satisfying some given moment conditions incorporating the target parameter of interest. The crux is that nonparametric profile likelihood ratio satisfies similar asymptotic properties as the likelihood ratio in the parametric context. %We shall illustrate these in detail in the context where the target parameter is $\max_{x\in\Theta}E[h(x;\xi)]$.
%We state the above mathematically.
Suppose the true parameter is known to satisfy the set of equations $E[t(\theta;\xi)]=0$ where $E[\cdot]$ is the expectation for the random object $\xi\in\mathbb R^d$, and $t(\theta;\xi),0\in\mathbb R^b$. We define the nonparametric profile likelihood ratio as
\begin{equation}
\mathcal R(\theta)=\max\left\{\prod_{i=1}^nnw_i:\sum_{i=1}^nw_it(\theta;\xi_i)=0,\ \sum_{i=1}^nw_i=1,\ w_i\geq0\text{\ for all\ }i=1,\ldots,n\right\}\label{profile1}
\end{equation}
where profiling refers to the categorization of all weights that respect the set of equations $E[t(\theta;\xi)]=0$.

With the above definitions, the crux is the empirical likelihood theorem (ELT):
\begin{thm}[Theorem 3.4 in \cite{owen2001empirical}]
Let $\xi_1,\ldots,\xi_n\in\mathbb R^d$ be i.i.d. data. Let $\theta_0\in\mathbb R^p$ be a value of the parameter that satisfies $E[t(\theta;\xi)]=0$, where $t(\theta;\xi),0\in\mathbb R^b$. Assume the covariance matrix $Var(t(\theta_0;\xi))$ is finite and has rank $q>0$. Then $-2\log\mathcal R(\theta_0)\Rightarrow\chi^2_q$, where $\mathcal R(\theta)$ is defined in \eqref{profile1}.\label{ELT estimating equations}
\end{thm}

The quantity $-2\log\mathcal R(\theta)$ is defined as $\infty$ if the optimization in \eqref{profile1} is infeasible.

We now explain how \eqref{EL opt} provides confidence bounds for optimization problem \eqref{obj}. We make the following assumptions:
\begin{assumption}
\begin{enumerate}
%\item $\xi\in\mathcal \xi$ where $\mathcal \xi$ a continuous space.
\item $h(x)$ is differentiable in $x$ with $\nabla_xh(x)=E[\nabla_xH(x;\xi)]$ for all $x\in\Theta$.\label{differentiability}
\item $x^*\in\text{argmin}_{x\in\Theta}h(x)$ if and only if $\nabla_xh(x^*)=0$. Moreover, this relation is \emph{distributionally stable}, meaning that $\tilde x^*\in\text{argmin}_{x\in\Theta}\tilde h(x)$ if and only if $\nabla_x\tilde h(\tilde x^*)=0$ for any $\tilde h(x)=\tilde E[H(\tilde x;\xi)]$ that has the expectation $\tilde E[\cdot]$ generated under an arbitrary distribution $\tilde P$ such that
$$\sup_{x\in\Theta}|\tilde h(x)-h(x)|<\epsilon$$
for small enough $\epsilon>0$.
\label{first order}
\item There exists an $x^*\in\text{argmin}_{x\in\Theta}h(x)$ such that the covariance matrix of the random vector $(\nabla_x H(x^*;\xi),H(x^*;\xi))\in\mathbb R^{p+1}$ is finite and has positive rank.\label{rank}
\item
$\frac{1}{n}\sum_{i=1}^nH(x;\xi_i)\to h(x)$ uniformly over $x\in\Theta$ a.s..\label{uniform convergence}
\item $E\left[\sup_{x\in\Theta}H(x;\xi)^2\right]<\infty$\label{second moment}
\end{enumerate}\label{EL assumptions}
\end{assumption}

Assuming the existence of $\nabla_xH(x;\xi)$ a.s., the interchangeability of derivative and expectation in Assumption \ref{EL assumptions}.\ref{differentiability} can generally be justified by the pathwise Lipschitz continuity condition
$$|H((x_1,\ldots,x_{j-1},u,x_{j+1},\ldots,x_p);\xi)-H((x_1,\ldots,x_{j-1},v,x_{j+1},\ldots,x_p);\xi)|\leq M_j|u-v|\text{\ \ a.s.}$$
for any $u,v$ in a nonrandom neighborhood around the point $x_j$ to be differentiated and $M_j$ measurable with $EM_j<\infty$ (e.g., \cite{asmussen2007stochastic}). Another sufficient condition is that $H(x;\xi)$ is a.s. continuous and piecewise differentiable in $x_j$ and $\sup_{u\in D}|(\partial/\partial x_j)H((x_1,\ldots,x_{j-1},u,x_{j+1},\ldots,x_p);\xi)|$ is integrable where $D$ is a neighborhood around $x_j$ \cite{Glasserman:1988:PCD:318123.318245}. Assumption \ref{EL assumptions}.\ref{first order} states that the first order condition for optimality is both sufficient and necessary. Assumptions \ref{EL assumptions}.\ref{first order} and \ref{EL assumptions}.\ref{uniform convergence} together ensure that this first order condition is unchanged when the true distribution is replaced by a (weighted) empirical version as the sample size gets large. Assumption \ref{EL assumptions}.\ref{second moment} is a technical condition required to bound the error between the empirical distribution and its weighted version within the divergence ball. Assumption \ref{EL assumptions}.\ref{rank} is used to invoke Theorem \ref{ELT estimating equations}. Note that $x^*$ is not necessarily unique.

As our subsequent development will reveal, both the necessity and the sufficiency of the first order condition in Assumption \ref{EL assumptions}.\ref{first order} are required; in particular, we need the necessity of $\nabla_xh(x^*)=0$ for $x^*\in\text{argmin}_{x\in\Theta}h(x)$ and the sufficiency of $\nabla_x\tilde h(\tilde x^*)=0$ for $\tilde x^*\in\text{argmin}_{x\in\Theta}\tilde h(x)$ in order for our argument on statistical guarantee to go through. %We have chosen to strengthen the assumption for simplicity.
Assumptions \ref{EL assumptions}.\ref{first order}, \ref{EL assumptions}.\ref{uniform convergence} and \ref{EL assumptions}.\ref{second moment} can be replaced by a single condition
\begin{assumption}
$x^*\in\text{argmin}_{x\in\Theta}h(x)$  if and only if $\nabla_xh(x^*)=0$, and $\tilde x^*\in\text{argmin}_{x\in\Theta}\sum_{i=1}^nw_iH(x;\xi_i)$ if and only if $\sum_{i=1}^nw_i\nabla_xH(\tilde x^*;\xi_i)=0$ for any support set $\{\xi_1,\ldots,\xi_n\}\subset\Theta$ and arbitrary probability vector $w$.\label{first order alternate}
\end{assumption}
Assumption \ref{first order alternate} is satisfied by, for instance, $H(\cdot;\xi)$ that is coersive and convex for any $\xi$ and $\Theta=\mathbb R^p$.

 
%Condition \ref{rank} states that all partial derivatives of $H(x^*;\xi)$ together with $H(x^*;\xi)$ itself are neither perfectly correlated or nondegenerate.

We have the following statistical guarantee:
\begin{thm}
Suppose $\xi_1,\ldots,\xi_n\in\mathbb R^d$ are i.i.d. data. Let $z^*$ be the optimal value of \eqref{obj}, and $\overline z$ and $\underline z$ be the maximum and minimum values of \eqref{EL opt} respectively. Then, under Assumption \ref{EL assumptions}, we have
$$\liminf_{n\to\infty}P\left(z^*\in[\underline z,\overline z]\right)\geq1-\beta.$$\label{EL thm}
\end{thm}
%\vspace{-2mm}
\begin{proof}
By Assumption \ref{EL assumptions}.\ref{rank}, there exists an $x^*\in\text{argmin}_{x\in\Theta}h(x)$ such that the covariance matrix of the random vector $(\nabla_x H(x^*;\xi),H(x^*;\xi))$ has a positive rank, call it $r$. Also, by Assumptions \ref{EL assumptions}.\ref{differentiability} and \ref{EL assumptions}.\ref{first order}, $x^*$ satisfies $E[\nabla_x H(x^*;\xi)]=0$. %, where $\mathbf 0$ is a length $p$ vector of zeros.

We define the nonparametric profile likelihood ratio as
\begin{equation}
\mathcal R(x,z)=\max\left\{\prod_{i=1}^nnw_i:\sum_{i=1}^nw_i\nabla_x H(x;\xi_i)=0,\ \sum_{i=1}^nw_iH(x;\xi_i)=z,\ \sum_{i=1}^nw_i=1,\ w_i\geq0\text{\ for all\ }i=1,\ldots,n\right\}.\label{profile}
\end{equation}
%This is the maximum ratio between the empirical likelihood and the nonparametric maximum likelihood, among all probability weights that respect the condition $E[\nabla_x h(x_0;\xi)]=\mathbf 0$ and $E[h(x_0;\xi)]=z$. The parameter $z$ is our target. We
Let $z^*=\min_{x\in\Theta}h(x)=h(x^*)$. From Theorem \ref{ELT estimating equations}, the nonparametric profile likelihood ratio \eqref{profile} satisfies $-2\log\mathcal R(x^*,z^*)\Rightarrow\chi^2_r$ as $n\to\infty$. %, where $\chi^2_q$ denotes the $\chi^2$-distribution with degree of freedom $q$.
This implies $P(-2\log\mathcal R(x^*,z^*)\leq\chi^2_{r,\beta})\to1-\beta$.

The rest of the proof focuses on the event $-2\log\mathcal R(x^*,z^*)\leq\chi^2_{r,\beta}$. Write
\begin{eqnarray}
&&-2\log\mathcal R(x^*,z^*)\notag\\
&=&\min\Bigg\{-2\sum_{i=1}^n\log(nw_i):\sum_{i=1}^nw_i\nabla_x H(x^*;\xi_i)=0,\ \sum_{i=1}^nw_iH(x^*;\xi_i)=z^*,\ \sum_{i=1}^nw_i=1,\ w_i\geq0{}\notag\\
&&{}\text{\ for all\ }i=1,\ldots,n\Bigg\}\label{profile updated}
\end{eqnarray}
We argue that $-2\log\mathcal R(x^*,z^*)\leq\chi^2_{r,\beta}$ implies the existence of a probability vector $w$ such that \begin{equation}
\sum_{i=1}^nw_i\nabla_x H(x^*;\xi_i)=0,\ \sum_{i=1}^nw_iH(x^*;\xi_i)=z^*,\ -2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{r,\beta}\label{properties}
\end{equation}
Notice that $-2\sum_{i=1}^n\log(nw_i)=\infty$ if $w_i=0$ for any $i$. Hence it suffices to replace, in \eqref{profile updated}, $w_i\geq0$ with $w_i\geq\epsilon$ for all $i$, for some small enough $\epsilon>0$. In this modified, compact, feasible set, $-2\sum_{i=1}^n\log(nw_i)$ is bounded and hence must possess an optimal solution $w$, which is a probability vector that satisfies \eqref{properties}.

This further implies that $z^*$ is bounded from above and below by the optimization problems
%So the set $\mathcal C=\{(x,z):-2\log\mathcal R(x,z)\leq\chi^2_{p+1,\beta}\}$ forms an asymptotically valid $1-\beta$ level confidence region for $(x_0,z_0)$, i.e. $P((x_0,z_0)\in\mathcal C)\to1-\beta$. This implies that with probability $1-\beta$ asymptotically,
%This further implies that with probability $1-\beta$ asymptotically, $z_0$ is contained by the pair of optimizations
\begin{equation}
\begin{array}{ll}
\max/\min_{w}&\sum_{i=1}^nw_iH(x^*;\xi_i)\\
\text{subject to}&\sum_{i=1}^nw_i\nabla_x H(x^*;\xi_i)=0\\
&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{r,\beta}\\
&\sum_{i=1}^nw_i=1\\
&w_i\geq0\text{\ for all\ }i=1,\ldots,n
\end{array}\label{opt interim}
\end{equation}
We argue that as $n\to\infty$, \eqref{opt interim} is equivalent to
\begin{equation}
\begin{array}{ll}
\max/\min_{w}&\sum_{i=1}^nw_iH(x^*;\xi_i)\\
\text{subject to}&w\in\{(w_1,\ldots,w_n):x^*\in\text{argmin}_{x\in\Theta}\sum_{i=1}^nw_iH(x;\xi_i)\}\\
&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{r,\beta}\\
&\sum_{i=1}^nw_i=1\\
&w_i\geq0\text{\ for all\ }i=1,\ldots,n
\end{array}\label{opt interim1}
\end{equation}
eventually (i.e. with probability 1). Note that the first constraint in \eqref{opt interim1} states that the probability vector $w$ must be chosen such that $x^*$, a minimizer of $h(x)$ picked at the beginning of this proof, also minimizes $\sum_{i=1}^nw_iH(x;\xi_i)$.

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 $\{\xi_1,\ldots,\xi_n\}$. Denote $E^w[\cdot]$ as the associated expectation and $h^w(x)=E^w[H(x;\xi)]$. We will show that
\begin{equation}
\sup_{x\in\Theta,w\in\mathcal W_r}|h^w(x)-h(x)|\to0\text{\ \ a.s.}\label{opt interim2}
\end{equation}
where
\begin{equation}
\mathcal W_r=\left\{(w_1,\ldots,w_n)\in\mathbb R^n:-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{r,\beta},\ \sum_{i=1}^nw_i=1,\ w_i\geq0\text{\ for all\ }i=1,\ldots,n\right\}\label{constraint set}
\end{equation}
Assumption \ref{EL assumptions}.\ref{first order} then implies that with probability 1, for sufficiently large $n$, $\sum_{i=1}^nw_i\nabla_x H(x^*;\xi_i)=0$ if and only if $x^*\in\text{argmin}_{x\in\Theta}\sum_{i=1}^nw_iH(x;\xi_i)$ for any $w\in\mathcal W_r$, leading to the equivalence.

We now show \eqref{opt interim2}. Consider
\begin{eqnarray}
\sup_{x\in\Theta,w\in\mathcal W_r}|h^w(x)-h(x)|&=&\sup_{x\in\Theta,w\in\mathcal W_r}\left|\int H(x;\xi)d(P^w-P)(\xi)\right|\notag\\
&\leq&\sup_{x\in\Theta,w\in\mathcal W_r}\left|\int H(x;\xi)d(P^w-\hat P)(\xi)\right|+\sup_{x\in\Theta,w\in\mathcal W_r}\left|\int H(x;\xi)d(\hat P-P)(\xi)\right|{}\notag\\
&&{}\text{\ \ where $\hat P$ denotes the empirical distribution generated from $\{\xi_1,\ldots,\xi_n\}$}\notag\\
&\leq&\sup_{x\in\Lambda,1\leq i\leq n,w\in\mathcal W_r}|H(x;\xi_i)|d_{TV}(P^w,\hat P)+\sup_{x\in\Theta}\left|\int H(x;\xi)d(\hat P-P)(\xi)\right|{}\label{interim new3}\\
&&{}\text{\ \ where $d_{TV}$ denotes the total variation distance}\notag
\end{eqnarray}
Now by Lemma 11.5 in \cite{owen2001empirical} (restated in the Appendix) and Assumption \ref{EL assumptions}.\ref{second moment}, we have
\begin{equation}
\sup_{x\in\Theta,1\leq i\leq n}|H(x;\xi_i)|=\max_{1\leq i\leq n}\sup_{x\in\Theta}|H(x;\xi_i)|=o(n^{1/2})\text{\ \ a.s.}\label{interim new1}
\end{equation}
On the other hand, by Pinsker's inequality, for any $w\in\mathcal W_r$,
\begin{equation}
d_{TV}(P^w,\hat P)\leq\sqrt{\frac{d_{KL}(P^w,\hat P)}{2}}=\sqrt{\frac{-\sum_{i=1}^n\log(nw_i)}{2n}}\leq\sqrt{\frac{\chi^2_r}{4n}}\label{interim new2}
\end{equation}
where $d_{KL}$ denotes the Kullback-Leibler divergence. Combining \eqref{interim new1} and \eqref{interim new2}, the first term in \eqref{interim new3} goes to 0 a.s.. The second term in \eqref{interim new3} converges to 0 a.s. by Assumption \ref{EL assumptions}.\ref{uniform convergence}. Hence $\sup_{x\in\Theta}|h^w(x)-h(x)|\to0$ a.s.. Therefore \eqref{opt interim} is equivalent to \eqref{opt interim1} eventually as $n\to\infty$.

Consider \eqref{opt interim1}. With the first constraint, the objective function must be equal to $\min_{x\in\Theta}\sum_{i=1}^nw_iH(x;\xi_i)$. Thus \eqref{opt interim1} is equivalent to
\begin{equation}
\begin{array}{ll}
\max/\min_w&\min_{x\in\Theta}\sum_{i=1}^nw_iH(x;\xi_i)\\
\text{subject to}&w\in\{(w_1,\ldots,w_n):x^*\in\text{argmin}_{x\in\Theta}\sum_{i=1}^nw_iH(x;\xi_i)\}\\
&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{r,\beta}\\
&\sum_{i=1}^nw_i=1\\
&w_i\geq0\text{\ for all\ }i=1,\ldots,n
\end{array}\label{opt}
\end{equation}
Let $\overline v$ and $\underline v$ be the maximum and minimum values of \eqref{opt}. Note that $r\leq p+1$ since $r$ is the rank of a $\mathbb R^{(p+1)\times(p+1)}$ matrix. This implies $\chi^2_{r,\beta}\leq\chi^2_{p+1,\beta}$. Together with a relaxation by removing the first constraint in \eqref{opt}, we have $\underline v\geq\underline z$ and $\overline v\leq\overline z$ where $\overline z$ and $\underline z$ are the maximum and minimum values of \eqref{EL opt}. From this we conclude that
$$\liminf_{n\to\infty}P(\underline z\leq z^*\leq\overline z)\geq\liminf_{n\to\infty}P(\underline v\leq z^*\leq\overline v)\geq\liminf_{n\to\infty}P(-2\log\mathcal R(x^*,z^*)\leq\chi^2_{r,\beta})=1-\beta$$


%$-2\log\mathcal R(x^*,z^*)\leq\chi^2_{r,\beta}$ implies that $z^*$ is bounded by the optimization pair \eqref{EL opt}. Therefore, we have
%$$\liminf_{n\to\infty}P\left(\underline z\leq z^*\leq\overline z\right)\geq\lim_{n\to\infty}P(-2\log\mathcal R(x^*,z^*)\leq\chi^2_{p+1,\beta})=1-\beta$$
%
%with probability at least $1-\beta$ asymptotically, $z_0$ is contained in \eqref{EL opt}, which proves Theorem \ref{EL thm}.
\end{proof}

Note that we have obtained bounds by relaxing the constraints in \eqref{opt}, and the degree of freedom in the $\chi^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 exists techniques (e.g., bootstrap calibration or Bartlett correction; \cite{owen1988empirical,diciccio1991empirical}) that can improve the coverage of the EL method in estimation problems. Investigation of these techniques in the optimization context is delegated to future work.

Note that the equivalence of \eqref{opt interim} and \eqref{opt interim1} holds for all $n$ if Assumption \ref{first order alternate} replaces Assumptions \ref{EL assumptions}.\ref{first order}, \ref{EL assumptions}.\ref{uniform convergence} and \ref{EL assumptions}.\ref{second moment}. %  proof Theorem \ref{EL thm} is an asymptotic statement. For small samples, accuracy (measured in terms of the empirical coverage probability being correctly larger than $1-\beta$) can be enhanced by other calibration methods than using $\chi^2_{p+1,\beta}$ in the optimization \eqref{EL opt}. Chapters 2 and 13 in \cite{owen2001empirical} provide some discussion.

\subsection{Stochastically Constrained Optimization}\label{sec:constrained}
We generalize the EL method to the stochastically constrained problem \eqref{obj constrained}. In this setting, we construct CI via the following optimization problems
\begin{equation}
\begin{array}{ll}
\max/\min_w&\left\{\begin{array}{ll}\min_x&\sum_{i=1}^nw_iH(x;\xi_i)\\
\text{subject to}&\sum_{i=1}^nw_iF_k(x;\xi_i)\leq0,\ k=1,\ldots,m\\
&g_k(x)\leq0,\ k=1,\ldots,s
\end{array}\right\}\\
\text{subject to}&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{p+m+1,\beta}\\
&\sum_{i=1}^nw_i=1\\
&w_i\geq0\text{\ for all\ }i=1,\ldots,n
\end{array}\label{opt constrained}
\end{equation}
While resembling \eqref{EL opt}, we note that the degree of freedom in the $\chi^2$-distribution is now $p+m+1$, which includes the number of stochastic constraints compared to \eqref{EL opt}.

For convenience, we denote
$$\Lambda=\left\{x\in\mathbb R^p:g_k(x)\leq0,\ k=1,\ldots,s\right\}$$
as the set of $x$ satisfying the deterministic constraints in \eqref{obj constrained}.

%To describe our main results on statistical guarantees for stochastically constrained problems, we first
We make the following assumptions in parallel to Assumption \ref{EL assumptions}:
\begin{assumption}
We assume:
\begin{enumerate}
%\item $h(x;\xi)$ and $f_k(x;\xi),\ k=1,\ldots,m$ are uniformly bounded over $x\in\Theta$.
\item $h(x)=E[H(x;\xi)]$, $f_k(x)=E[F_k(x;\xi)],\ k=1,\ldots,m$ and $g_k(x),\ k=1,\ldots,s$ are all differentiable in $x\in\Lambda$, and
$$\nabla_xh(x)=E[\nabla_xH(x;\xi)],\ \nabla_xf_k(x)=E[\nabla_xF_k(x;\xi)]$$
%Moreover, $\nabla_xh(\xi;x)$, $\nabla_xf_k(x;\xi)$ are all uniformly bounded over $x\in\Theta$.
\item Let $S^*$ be the set of all optimal solutions for \eqref{obj constrained}. $x^*\in S^*$ if and only if $x^*$ satisfies the KKT condition, where the active set of the KKT condition (i.e. equalities) is unique among all $x^*\in S^*$ and is sufficient for determining $S^*$. This relation is \emph{distributionally stable}, meaning that $\tilde x^*\in\tilde S^*$, where $\tilde S^*$ is the set of optimal solutions for
\begin{equation} \label{obj constrained perturbed}
\begin{array}{ll}\min&\tilde h(x)\\
\text{subject to}&\tilde f_k(x)\leq0,\ k=1,\ldots,m\\
&\tilde g_k(x)\leq0,\ k=1,\ldots,s
\end{array}
\end{equation}
if and only if $\tilde x^*$ satisfies the corresponding KKT condition, where $\tilde h(x)=\tilde E[H(x;\xi)]$ and $\tilde f_k(x)=\tilde E[F_k(x;\xi)]$, with $\tilde E$ denoting the expectation under an arbitrary distribution $\tilde P$ such that
$$\sup_{x\in\Lambda}|\tilde h(x)-h(x)|<\epsilon$$
$$\sup_{x\in\Lambda}|\tilde f_k(x)-f_k(x)|<\epsilon\text{\ \ for all\ }k=1,\ldots,m$$
for small enough $\epsilon>0$. Moreover, for any such $\epsilon>0$, the active set of the KKT condition at any $\tilde x^*\in\tilde S^*$ for \eqref{obj constrained perturbed} is the same as that at any $x^*\in S^*$ for \eqref{obj constrained} and is sufficient for determining $\tilde S^*$.%Moreover, the set of active constraints in the KKT condition remains the same as that of \eqref{obj constrained}.
%\item The
%
%     $\bm\theta^*$ be the $E[h(X;\bm\theta)]$, $E[f_k(X;\bm\theta)],\ k=1,\ldots,m$ and $E[g_k(X;\bm\theta)],\ k=1,\ldots,m'$ are all linearly independent in $\bm\theta$.
\item
  There exists an optimal solution $x^*$ for \eqref{obj constrained}, with associated Lagrange multipliers for the stochastic constraints in \eqref{obj constrained} given by $\lambda^*=(\lambda_1^*,\ldots,\lambda_m^*)$, such that the covariance matrix of the variables $H(x^*;\xi)$, $\frac{\partial}{\partial x_j}H(x^*;\xi)+\sum_{k=1}^m\lambda_k^*\frac{\partial}{\partial x_j}F_k(x^*;\xi)$, and $F_k(x^*;\xi)$, for all indices $j$ and $k$ corresponding to the active set of the KKT condition, is finite and has positive rank.
\item
$$\frac{1}{n}\sum_{i=1}^nH(x;\xi_i)\to h(x)\text{\ \ and\ \ }\frac{1}{n}\sum_{i=1}^nF_k(x;\xi_i)\to f_k(x),\ \ k=1,\ldots,m$$
uniformly over $x\in\Lambda$ a.s..
\item $E\left[\sup_{x\in\Theta}H(x;\xi)^2\right]<\infty$ and $E\left[\sup_{x\in\Theta}F_k(x;\xi)^2\right]<\infty$ for $k=1,\ldots,m$.
\end{enumerate}\label{assumptions constrained}
\end{assumption}

Denote $\nu^*=(\nu_1^*,\ldots,\nu_s^*)$ as the Lagrange multiplier for the deterministic constraints in \eqref{obj constrained}. In Assumptions \ref{assumptions constrained}.2 and \ref{assumptions constrained}.3 above, the active set of the KKT condition satisfied by $(x^*,\lambda^*,\nu^*)$ is in the form
$$\frac{\partial}{\partial x_j}h(x^*)+\sum_{k=1}^m\lambda_k^*\frac{\partial}{\partial x_j}f_k(x^*)+\sum_{k=1}^s\nu_k^*\frac{\partial}{\partial x_j}g_k(x^*)=0,\ j\in\mathcal A_1^*\equiv\{1,\ldots,p\}$$
  $$f_k(x^*)=0,\ k\in\mathcal A_2^*\subset\{1,\ldots,m\}$$
  $$g_k(x^*)=0,\ k\in\mathcal A_3^*\subset\{1,\ldots,s\}$$
  %$$x_j^*=0,\ j\in\{1,\ldots,p\}\setminus\mathcal A_1^*$$
  $$\lambda_k^*=0\ k\in\{1,\ldots,m\}\setminus\mathcal A_2^*$$
  $$\nu_k^*=0,\ \ k\in\{1,\ldots,s\}\setminus\mathcal A_3^*$$
%  The optimization \eqref{obj constrained} has an optimal solution $x^*$ as the unique solution that satisfies all the active constraints in the KKT condition.
where $\mathcal A_1^*$, $\mathcal A_2^*$ and $\mathcal A_3^*$ denote the sets of indices that correspond to the equalities in the condition, which are unique among any optimal solutions of \eqref{obj constrained} by Assumption \ref{assumptions constrained}.2. The $j$ and $k$ described in Assumption \ref{assumptions constrained}.3 refer to the indices in $\mathcal A_1^*$ and $\mathcal A_2^*$. Assumption \ref{assumptions constrained}.2 further enforces the sets $\mathcal A_1^*$, $\mathcal A_2^*$ and $\mathcal A_3^*$ to remain as the active sets under a perturbation to $\tilde P$ described therein, and the equalities indexed via these sets are enough to determine $S^*$ and $\tilde S^*$. Assumptions \ref{assumptions constrained}.2 and \ref{assumptions constrained}.3 generalize Assumptions \ref{EL assumptions}.\ref{first order} and \ref{EL assumptions}.\ref{rank} from the first order condition to the KKT condition. Similar to Section \ref{sec:unconstrained}, we require the necessity of the KKT and the active set conditions regarding \eqref{obj constrained} and the sufficiency regarding \eqref{obj constrained perturbed} for our development to go through. Constraint qualification for the validity of the KKT condition is implicitly assumed in Assumption \ref{assumptions constrained}.2.

%Note that Conditions 1 and 3 are similar to those of Assumption \ref{EL assumptions}. Conditions 2 and 4 together ensure that the active set of constraints is unchanged when the true distribution is replaced by its (weighted) empirical version. Condition 5 is a technical condition required to bound the error between the empirical distribution and its weighted version within the divergence ball.


We have the following result:
\begin{thm}
Suppose $\xi_1,\ldots,\xi_n$ are i.i.d. data. Under Assumption \ref{assumptions constrained}, we have
$$\liminf_{n\to\infty}P(z^*\in[\underline z,\overline z])\geq1-\beta$$
where $z^*$ is the optimal value of \eqref{obj constrained}, and $\underline z$ and $\overline z$ are the minimum and maximum values of \eqref{opt constrained}.
\end{thm}
%
%
%where $q$ is the sum of the total number of decision variables $p+p'$ and the total number of stochastic constraints in \eqref{opt constrained} (which is at most $m+m'$) contain the true optimal value of \eqref{opt constrained} with probability at least $\beta$.

\begin{proof}
%Let $\lambda_1,\ldots,\lambda_m,\nu_1,\ldots,\nu_s$ be the Lagrange multipliers for the functional constraints in \eqref{obj constrained}. The active set of the KKT condition satisfied by $x^*$ is in the form
%$$E\left[\frac{\partial}{\partial x_j}h(x;\xi)-\sum_{k=1}^m\lambda_k\frac{\partial}{\partial x_j}f_k(x;\xi)\right]-\sum_{k=1}^s\nu_k\frac{\partial}{\partial x_j}g_k(x)=0,\ j\in\mathcal A_1\subset\{1,\ldots,d\}$$
%  $$E[f_k(x;\lambda)]=0,\ k\in\mathcal A_2\subset\{1,\ldots,m\}$$
%  $$g_k(x)=0,\ k\in\mathcal A_3\subset\{1,\ldots,s\}$$
%  $$x_j=0,\ j\in\{1,\ldots,d\}\setminus\mathcal A_1$$
%  $$\lambda_k=0\ k\in\{1,\ldots,m\}\setminus\mathcal A_2$$
%  $$\nu_k=0,\ \ k\in\{1,\ldots,s\}\setminus\mathcal A_3$$
%with the rest of the KKT conditions being inequalities.

%    $$E\left[\frac{\partial}{\partial\theta_j'}h(X;\bm\theta)-\sum_{k=1}^m\lambda_k\frac{\partial}{\partial\theta_j'}f_k(X;\bm\theta)-\sum_{k=1}^{m'}\lambda_k'\frac{\partial}{\partial\theta_j'}g_k(X;\bm\theta)\right]=0,\ j=1,\ldots,p'$$
%
%    $$E[g_k(X;\bm\theta)]=0,\ k=1,\ldots,m'$$
%    $$\theta_j\geq0,\ j=1,\ldots,p$$
%    $$\lambda_k\geq0,\ k=1,\ldots,m$$
%    The number of active expressions (i.e. equalities being satisfied) from above at the optimal primal-dual pair $(\bm\theta^*,\bm\lambda^*)$ is $p+p'+m+m'$. Let $E[\mathbf a(X_i;\bm\theta,\bm\lambda)]=\textbf 0$ be this active set of equalities.


%To see the theorem above, note that KKT condition and Assumption 2 imply that the optimal $(\bm\theta^*,\bm\lambda^*)$ must satisfy the equalities
%$$\theta_j=0\text{\ \ or\ \ }E\left[\frac{\partial}{\partial\theta_j}h(X;\bm\theta)-\sum_{k=1}^m\lambda_k\frac{\partial}{\partial\theta_j}f_k(X;\bm\theta)-\sum_{k=1}^{m'}\lambda_k'\frac{\partial}{\partial\theta_j}g_k(X;\bm\theta)\right]=0,\text{\ \ for each\ \ }j=1,\ldots,p$$
%$$E\left[\frac{\partial}{\partial\theta_j'}h(X;\bm\theta)-\sum_{k=1}^m\lambda_k\frac{\partial}{\partial\theta_j'}f_k(X;\bm\theta)-\sum_{k=1}^{m'}\lambda_k'\frac{\partial}{\partial\theta_j'}g_k(X;\bm\theta)\right]=0,\ j=1,\ldots,p'$$
%$$\lambda_k=0\text{\ \ or\ \ }E[f_k(X;\bm\theta)]=0,\text{\ \ for each\ \ }k=1,\ldots,m$$
%$$E[g_k(X;\bm\theta)]=0,\ k=1,\ldots,m'$$
%and the inequalities
%$$E\left[\frac{\partial}{\partial\theta_j}h(X;\bm\theta)-\sum_{k=1}^m\lambda_k\frac{\partial}{\partial\theta_j}f_k(X;\bm\theta)-\sum_{k=1}^{m'}\lambda_k'\frac{\partial}{\partial\theta_j}g_k(X;\bm\theta)\right]\leq0,\ j=1,\ldots,p$$
%$$E[f_k(X;\bm\theta)]\leq0,\ k=1,\ldots,m$$
%    $$\theta_j\geq0,\ j=1,\ldots,p$$
%    $$\lambda_k\geq0,\ k=1,\ldots,m$$

%By Assumption 3, $(\bm\theta^*,\bm\lambda^*)$ must satisfy $p+p'+m+m'$ of the equalities written above. Let the profile likelihood
Consider the nonparametric profile likelihood ratio
\begin{equation}
\mathcal R(x,\lambda,\nu,z)=\max\left\{\prod_{i=1}^nnw_i:\begin{array}{l}\sum_{i=1}^nw_iH(x;\xi_i)=z\\
\sum_{i=1}^nw_i\left(\frac{\partial}{\partial x_j}H(x;\xi_i)+\sum_{k=1}^m\lambda_k\frac{\partial}{\partial x_j}F_k(x;\xi_i)\right)+\sum_{k=1}^s\nu_k\frac{\partial}{\partial x_j}g_k(x)=0,\ j\in\mathcal A_1^*\\
\sum_{i=1}^nw_iF_k(x;\xi_i)=0,\ k\in\mathcal A_2^*\\
g_k(x)=0,\ k\in\mathcal A_3^*\\
%x_j=0,\ j\in\{1,\ldots,p\}\setminus\mathcal A_1^*\\
\lambda_k=0\ k\in\{1,\ldots,m\}\setminus\mathcal A_2^*\\
\nu_k=0,\ \ k\in\{1,\ldots,s\}\setminus\mathcal A_3^*\\
\sum_{i=1}^nw_i=1\\
w_i\geq0\text{\ for all\ }i=1,\ldots,n
\end{array}\right\}\label{profile constrained}
\end{equation}
Let $x^*$ be an optimal solution for \eqref{obj constrained} satisfying Assumption \ref{assumptions constrained}.3, and $\lambda^*=(\lambda_1^*,\ldots,\lambda_m^*),\nu^*=(\nu_1^*,\ldots,\nu_s^*)$ be its associated Lagrange multipliers. %We have $(x^*,\lambda^*,\nu^*)$ satisfies the equations
%$$E\left[\frac{\partial}{\partial x_j}H(x;\xi)+\sum_{k=1}^m\lambda_k\frac{\partial}{\partial x_j}F_k(x;\xi)\right]+\sum_{k=1}^s\nu_k\frac{\partial}{\partial x_j}g_k(x)=0,\ j\in\mathcal A_1$$
%  $$E[F_k(x;\xi)]=0,\ k\in\mathcal A_2$$
  By Assumption \ref{assumptions constrained}.3, the covariance of the random vector concatenated by
  \begin{eqnarray*}
  &&H(x^*;\xi)\\
  &&\frac{\partial}{\partial x_j}H(x^*;\xi)+\sum_{k=1}^m\lambda_k^*\frac{\partial}{\partial x_j}F_k(x^*;\xi)+\sum_{k=1}^s\nu_k^*\frac{\partial}{\partial x_j}g_k(x^*)\text{\ \ for\ }j\in\mathcal A_1^*\\
  &&F_k(x^*;\xi)\text{\ \ for\ }k\in\mathcal A_2^*\\
  \end{eqnarray*}
  has rank $r$ for some $r>0$. Let $z^*$ be the optimal value of \eqref{obj constrained} equal to $h(x^*)$. Since the other active KKT conditions are deterministic, Theorem \ref{ELT estimating equations} implies that $-2\log\mathcal R(x^*,\lambda^*,\nu^*,z^*)\Rightarrow\chi^2_r$, %Note that, by scrutinizing the set of possible equalities in the KKT condition, we must have $r\leq q+1$ where $q$ is as defined in the theorem.
which further implies $P(-2\log\mathcal R(x^*,\lambda^*,\nu^*,z^*)\leq\chi^2_{r,\beta})\to1-\beta$.
    %&&g_k(x)\text{\ \ for\ }k\in\mathcal A_3\\
%  &&x_j\text{\ \ for\ }j\in\{1,\ldots,p\}\setminus\mathcal A_1\\
%  &&\lambda_k\text{\ \ for\ }k\in\{1,\ldots,m\}\setminus\mathcal A_2\\
%  &&\nu_k\text{\ \ for\ }k\in\{1,\ldots,s\}\setminus\mathcal A_3


%So the set $\mathcal C=\{(x,\lambda,\nu,\mu):-2\log\mathcal R(x,\lambda,\nu,\mu)\leq\chi^2_{r,\beta}\}$ forms an asymptotically valid $1-\beta$ level confidence region for $(x^*,\lambda^*,\nu^*,\mu^*)$, i.e. $P((x,\lambda,\nu,\mu)\in\mathcal C)\to1-\beta$. %This implies that with probability $1-\beta$ asymptotically, there exists a probability vector $\mathbf w$ such that $\sum_{i=1}^nw_i\nabla_x h(x_0;\xi_i)=\mathbf 0$, $\sum_{i=1}^nw_ih(x_0;\xi_i)=z_0$ and $-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{p+1,\beta}$. This further implies that with probability $1-\beta$ asymptotically, $z_0$ is contained by the pair of optimizations

%Note that $(x,\lambda,\nu,\mu)\in\mathcal C$ implies that $\underline V\leq\mu\leq\overline V$ where $\underline V$ and $\overline V$ are the min and max values of the programs
Similar to the proof of Theorem \ref{EL thm}, $-2\log\mathcal R(x^*,\lambda^*,\nu^*,z^*)\leq\chi^2_{r,\beta}$ implies the existence of a $w$ that satisfies $-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{r,\beta}$ and all constraints in \eqref{profile constrained} evaluated at $x^*,\lambda^*,\nu^*,z^*$. This in turn implies that $z^*$ is bounded by
\begin{equation}
\begin{array}{ll}
\max/\min_w&\sum_{i=1}^nw_iH(x^*;\xi_i)\\
\text{subject to}&\sum_{i=1}^nw_i\left(\frac{\partial}{\partial x_j}H(x^*;\xi_i)+\sum_{k=1}^m\lambda_k^*\frac{\partial}{\partial x_j}F_k(x^*;\xi_i)\right)+\sum_{k=1}^s\nu_k^*\frac{\partial}{\partial x_j}g_k(x^*)=0,\ j\in\mathcal A_1^*\\
&\sum_{i=1}^nw_iF_k(x^*;\xi_i)=0,\ k\in\mathcal A_2^*\\
&g_k(x^*)=0,\ k\in\mathcal A_3^*\\
%&x_j=0,\ j\in\{1,\ldots,p\}\setminus\mathcal A_1^*\\
&\lambda_k^*=0\ k\in\{1,\ldots,m\}\setminus\mathcal A_2^*\\
&\nu_k^*=0,\ \ k\in\{1,\ldots,s\}\setminus\mathcal A_3^*\\
&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{r,\beta}\\
&\sum_{i=1}^nw_i=1\\
&w_i\geq0\text{\ for all\ }i=1,\ldots,n
\end{array}\label{interim opt constrained}
\end{equation}
Using the same argument as in the proof of Theorem \ref{EL thm}, we obtain from Assumptions \ref{assumptions constrained}.4 and \ref{assumptions constrained}.5 that
$$\sup_{x\in\Lambda,w\in\mathcal W_r}|h^w(x)-h(x)|\to0\text{\ \ a.s.}$$
$$\sup_{x\in\Lambda,w\in\mathcal W_r}|f_k^w(x)-f_k(x)|\to0\text{\ \ a.s. for all\ }k=1,\ldots,m$$
where $\mathcal W_r$ is defined in \eqref{constraint set}, and $h^w(x)=E^w[H(x;\xi)]$, $f_k^w(x)=E^w[F_k(x;\xi)]$ with $E^w$ denoting the expectation with respect to $P^w$, the probability distribution represented by the weights $w$ on the support $\{\xi_1,\ldots,\xi_n\}$. %the active constraints in the KKT condition remains the same as that of %Thus, denoting $\mathcal A_j^w$ as the set of indices for the active constraints in the KKT condition under the sample counterpart weighted by $w$, we have
Thus, by Assumption \ref{assumptions constrained}.2, the set of active KKT conditions for an optimal solution of the weighted sample problem
$$\begin{array}{ll}\min_x&\sum_{i=1}^nw_iH(x;\xi_i)\\
\text{subject to}&\sum_{i=1}^nw_iF_k(x;\xi_i)\leq0,\ k=1,\ldots,m\\
&g_k(x)\leq0,\ k=1,\ldots,s
\end{array}$$
for any $w\in\mathcal W_r$ is identical to that for $x^*$ for \eqref{obj constrained} eventually as $n\to\infty$, and Assumption \ref{assumptions constrained}.2 further implies that \eqref{interim opt constrained} is equivalent to
\begin{equation}
\begin{array}{ll}
\max/\min_w&\sum_{i=1}^nw_iH(x^*;\xi_i)\\
\text{subject to}&w\in\left\{(w_1,\ldots,w_n):x^*\in\left\{\begin{array}{ll}\text{argmin}_x&\sum_{i=1}^nw_iH(x;\xi_i)\\
\text{subject to}&\sum_{i=1}^nw_iF_k(x;\xi_i)\leq0,\ k=1,\ldots,m\\
&g_k(x)\leq0,\ k=1,\ldots,s
\end{array}\right\}\right\}\\
&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{r,\beta}\\
&\sum_{i=1}^nw_i=1\\
&w_i\geq0\text{\ for all\ }i=1,\ldots,n
\end{array}\label{interim opt new1 constrained}
\end{equation}
eventually as $n\to\infty$. With the first constraint, the objective function in \eqref{interim opt new1 constrained} must be equal to $\min_x\{\sum_{i=1}^nw_iH(x;\xi_i):\sum_{i=1}^nw_iF_k(x;\xi_i)\leq0,k=1,\ldots,m,\ g_k(x)\leq0,k=1,\ldots,s\}$. Note that
\begin{equation}
r\leq1+|\mathcal A_1^*|+|\mathcal A_2^*|\leq1+p+m\label{df}
\end{equation}
where $|\cdot|$ denotes cardinality. This implies that $\chi^2_{r,\beta}\leq\chi^2_{p+m+1,\beta}$. Thus, together with a relaxation of the first constraint in \eqref{interim opt new1 constrained}, the same argument as in the proof of Theorem \ref{EL thm} stipulates that the maximum and minimum values of \eqref{interim opt new1 constrained} are bounded from above and below respectively by those of \eqref{opt constrained} and concludes the theorem.
%\begin{equation*}
%\begin{array}{ll}
%\max/\min_w&\left\{\begin{array}{ll}\min_x&\sum_{i=1}^nw_iH(x;\xi_i)\\
%\text{subject to}&\sum_{i=1}^nw_iF_k(x;\xi_i)\leq0,\ k=1,\ldots,m\\
%&g_k(x)\leq0,\ k=1,\ldots,s
%\end{array}\right\}\\
%\text{subject to}&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{p+m+1,\beta}\\
%&\sum_{i=1}^nw_i=1\\
%&w_i\geq0\text{\ for all\ }i=1,\ldots,n
%\end{array}
%\end{equation*}



%$\underline v\geq\underline z$ and $\overline v\leq\overline z$, where $\underline v$ and $\overline v$ denote the minimum and maximum values of \eqref{opt constrained}. This concludes that
%$$\liminf_{n\to\infty}P(\underline z\leq z^*\leq\overline z)\geq\liminf_{n\to\infty}P(\underline v\leq z^*\leq\overline v)\geq\lim_{n\to\infty}P(-2\log\mathcal R(x^*,\lambda^*,\nu^*,z^*)\leq\chi^2_{r,\beta})=1-\beta$$



% it suffices to show that the distribution represented by $\mathbf w$ on the data support, which we denote as $P^{\mathbf w}$, with $\mathbf w$ such that $-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{r,\beta}$, satisfies

%Note that
%\begin{eqnarray}
%\sup_{x\in\Lambda}|h^{\mathbf w}(x)-h(x)|&=&\sup_{x\in\Lambda}\left|\int H(x;\xi)d(P^{\mathbf w}-P)(\xi)\right|\notag\\
%&\leq&\sup_{x\in\Lambda}\left|\int H(x;\xi)d(P^{\mathbf w}-\hat P)(\xi)\right|+\sup_{x\in\Lambda}\left|\int H(x;\xi)d(\hat P-P)(\xi)\right|{}\notag\\
%&&{}\text{\ \ where $\hat P$ denotes the empirical distribution}\notag\\
%&\leq&\sup_{x\in\Lambda,\ 1\leq i\leq n}|H(x;\xi_i)|d_{TV}(P^{\mathbf w},\hat P)+\sup_{x\in\Lambda}\left|\int H(x;\xi)d(\hat P-P)(\xi)\right|{}\notag\\
%&&{}\text{\ \ where $d_{TV}$ denotes the total variation distance}\label{interim new3}
%\end{eqnarray}
%Now by Lemma 11.5 in \cite{owen2001empirical} (restated in the Appendix) and Assumption \ref{assumptions constrained} Condition 5, we have
%\begin{equation}
%\sup_{x\in\Lambda,\ 1\leq i\leq n}|H(x;\xi_i)|=\max_{1\leq i\leq n}\sup_{x\in\Lambda}|H(x;\xi)|=o(n^{1/2})\text{\ \ a.s.}\label{interim new1}
%\end{equation}
%On the other hand, by Pinsker's inequality,
%\begin{equation}
%d_{TV}(P^{\mathbf w},\hat P)\leq\sqrt{\frac{d_{KL}(P^{\mathbf w},\hat P)}{2}}=\sqrt{\frac{-\sum_{i=1}^n\log(nw_i)}{2n}}\leq\sqrt{\frac{\chi^2_r}{4n}}\label{interim new2}
%\end{equation}
%where $d_{KL}$ denotes the Kullback-Leibler divergence. Combining \eqref{interim new1} and \eqref{interim new2}, the first term in \eqref{interim new3} goes to 0 a.s.. The second term in \eqref{interim new3} converges to 0 a.s. by Assumption \ref{assumptions constrained} Condition 4. Hence $\sup_{x\in\Lambda}|h^{\mathbf w}(x)-h(x)|\to0$ a.s.. The same argument shows that $\sup_{x\in\Lambda}|f_k^{\mathbf w}(x)-f_k(x)|\to0$ a.s. for all $k=1,\ldots,m$. Therefore, \eqref{interim opt} is equivalent to \eqref{opt constrained}.
%
%\end{align}
%\leq Cd_W(P^{\mathbf w},P)\label{Wasserstein}
%\end{equation}
%where $d_W(P^{\mathbf w},P)$ denotes the Wasserstein metric of order 1, given by
%$$d_W(Q,P)=\sup\left\{\int f(\xi)d(Q-P)(\xi):f\in\text{Lip}(1)\right\}$$
%with $\text{Lip}(1)$ denotes the set of all functions having Lipschitz constant 1, and $C$ is the uniform bound on the derivative of $h$ as guaranteed in Assumption \ref{assumptions constrained} Condition 1.
%
%Using triangle inequality, \eqref{Wasserstein} is bounded from above by
%$$d_W(P^{\mathbf w},\hat P)+d_W(\hat P,P)$$
%where $\hat P$ denotes the empirical distribution. Since $\mathbf w$ satisfies $-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{r,\beta}$

%Let $\textbf w=(w_1,\ldots,w_n)$. The two numbers solved by the programs
%\begin{equation}
%\begin{array}{ll}
%\max/\min_{\textbf w,\bm\theta,\bm\lambda,\mu}&\mu\\
%\text{subject to}&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{q+1,\beta}\\
%&\sum_{i=1}^nw_ih(X_i;\bm\theta)=\mu\\
%&\sum_{i=1}^nw_i\mathbf a(X_i;\bm\theta,\bm\lambda)=\textbf 0\\
%&\sum_{i=1}^nw_i=1\\
%&w_i\geq0\text{\ for all\ }i=1,\ldots,n
%\end{array}\label{opt1}
%\end{equation}
%contain the true optimal value of \eqref{opt constrained} with probability at least $\beta$ asymptotically. Because of Assumption 3, solving

\end{proof}

Note that, much like the proof of Theorem \ref{EL thm}, we have relaxed constraints and placed a conservative bound on the degree of freedom of the $\chi^2$-distribution in \eqref{df}, which could potentially be improved with more refined analysis.

%\section{Computational Load}
%We shall discuss how to compute \eqref{EL opt}, which consists of a min-max problem and a max-max problem. Note that the outer objective function $\max_{x\in\Theta}\sum_{i=1}^nw_ih(x;\xi_i)$ is a convex function in $\mathbf w$. Suppose that $h(\xi;\cdot)$ is concave in $x$. Then evaluating $\max_{x\in\Theta}\sum_{i=1}^nw_ih(x;\xi_i)$ for fixed $\mathbf w$ boils down to a convex optimization, and hence both the outer and inner optimizations in the min-max problem are convex and efficiently solvable. Moreover, to speed up computation, by observing that the feasible region for $\mathbf w$ is compact one can consider an exchange of min-max to max-min via Sion's minimax theorem \cite{sion1958general}, which gives
%\begin{equation}
%\min_{\mathbf w\in\mathcal A}\max_{x\in\Theta}\sum_{i=1}^nw_ih(x;\xi_i)=\max_{x\in\Theta}\min_{\mathbf w\in\mathcal A}\sum_{i=1}^nw_ih(x;\xi_i)\label{Sion}
%\end{equation}
%where $\mathcal A$ denotes the feasible region in \eqref{EL opt}. Fixing $x$, the inner minimization $\min_{\mathbf w\in\mathcal A}\sum_{i=1}^nw_ih(x;\xi_i)$ in the maximin formulation can be solved by Lagrangian relaxation
%\begin{eqnarray}
%&&\max_{\nu\geq0,\lambda\in\mathbb R}\min_{\mathbf w\geq\mathbf 0}\sum_{i=1}^nw_ih(x;\xi_i)+\nu\left(-2\sum_{i=1}^n\log(nw_i)-\chi^2_{p+1,\beta}\right)-\lambda\left(\sum_{i=1}^nw_i-1\right)\notag\\
%&=&\max_{\nu\geq0,\lambda\in\mathbb R}\sum_{i=1}^n\min_{w_i\geq0}\{w_i(h(x;\xi_i)-\lambda)-2\nu\log(nw_i)\}-2\nu n\log n-\nu\chi^2_{p+1,\beta}+\lambda.\label{interim Lagrangian}
%\end{eqnarray}
%Straightforward calculation reveals that
%$$\min_{w_i\geq0}\{w_i(h(x;\xi_i)-\lambda)-2\nu\log(nw_i)\}=\left\{\begin{array}{ll}2\nu\left(1-\log\left(\frac{2\nu}{h(x;\xi_i)-\lambda}\right)\right)&\text{\ if\ }h(x;\xi_i)>\lambda\\
%-\infty&\text{\ otherwise}\end{array}\right.$$
%with the optimal solution of $w_i$ being $2\nu/(h(x;\xi_i)-\lambda)$ or $\infty$ correspondingly. This implies that \eqref{interim Lagrangian} is equal to
%$$\max_{\nu\geq0,\ \lambda<\min_ih(x;\xi_i)}-\sum_{i=1}^n2\nu\log\left(\frac{2\nu}{h(x;\xi_i)-\lambda}\right)+2\nu n(1-\log n)-\nu\chi^2_{p+1,\beta}+\lambda$$
%where the optimal solution of $w_i$ is $2\nu/(h(x;\xi_i)-\lambda)$. Hence the inner minimization of the right hand side of \eqref{Sion} becomes a convex maximization problem with two variables (instead of having $n$ variables in the original formulation).
%
%On the other hand, the max-max problem in \eqref{EL opt} can be more challenging, because the outer optimization involves maximizing the convex function $\max_{x\in\Theta}\sum_{i=1}^nw_ih(x;\xi_i)$ over $\mathbf w$. A quick heuristic for obtaining a solution for the max-max problem is to do alternating maximization, namely repeatedly fixing $\mathbf w$ and maximizing over $x$, and fixing $x$ and maximizing over $\mathbf w$, until no improvement is observed. \cite{ct1984} contains global convergence results for this sort of procedure under some improvement bound assumptions uniformly on the alternating steps. Similar as above, when $x$ is fixed, one can write $\max_{\mathbf w\in\mathcal A}\sum_{i=1}^nw_ih(x;\xi_i)$ via Lagrangian relaxation as
%$$\min_{\nu\geq0,\ \lambda>\max_ih(x;\xi_i)}2\nu\sum_{i=1}^n\log\left(\frac{2\nu}{\lambda-h(x;\xi_i)}\right)+2\nu n(\log n-1)+\nu\chi^2_{p+1,\beta}+\lambda$$
%with the optimal solution of $w_i$ being $2\nu/(\lambda-h(x;\xi_i))$.
%
%The computation of (\ref{opt constrained}) is similar, except that the constraints $\sum_{i=1}^n w_iF_k(x;\xi_i)\leq 0, k=1,\ldots,m'$ and $\sum_{i=1}^n w_iF_k(x;\xi_i) = 0, k=m'+1,\ldots,m$ involve both $x$ and $w$. Hence, we put these constraints into the objective function by introducing Lagarangian multipliers $\lambda \in\Omega = \{\lambda_i\geq 0, i=1,\ldots,m'; \lambda_i\in\mathbb{R}, i=m'+1,\ldots,m\}$.
%\begin{equation}
%\min_{\mathbf w\in\mathcal{A}} \min_{\lambda\in\Omega}\max_{x\in\Lambda}~~\sum_{i=1}^nw_i\left(H(x;\xi_i)+\sum_{k=1}^{m}\lambda_kF_k(x;\xi_i)\right)
%\end{equation}
%By exchanging the minimizations, we have
%\begin{equation}
%\min_{\lambda\in\Omega}\min_{\mathbf w\in\mathcal{A}}\max_{x\in\Lambda}~~\sum_{i=1}^nw_i \tilde{H}(x;\xi_i),
%\end{equation}
%where $\tilde{H}(x;\xi_i) = H(x;\xi_i)+\sum_{k=1}^m\lambda_kF_k(x;\xi_i)$.
%For a fixed $\lambda$, note the objective function for $w$, i.e., $\max_{x\in\Lambda}~~\sum_{i=1}^nw_i\tilde{H}(x;\xi_i)$ is a convex function for $w$. Hence, the same simplification described above for the unconstrained problem can be used, leading to an optimization problem
%$$
%\min_{\lambda\in\Omega}\max_{x\in\Lambda} \max_{\nu\geq0,\mu<\min_i \tilde{H}(x;\xi_i)}~~-\sum_{i=1}^n2\nu\log\left(\frac{2\nu}{\tilde{H}(x;\xi_i)-\mu}\right)+2\nu n(1-\log n)-\nu\chi^2_{p+1,\beta}+\mu.
%$$

\section{The Empirical Likelihood Method for Constructing Confidence Bounds for Optimality Gaps}\label{optimality gap section}
We study the construction of CI for the optimality gap of a given solution using the EL method. Suppose $\hat x$ is obtained from some procedure independently of the data $\xi_1,\ldots,\xi_n$. The optimality gap of $\hat x$ is given by $\mathcal G(\hat x)=h(\hat x)-z^*$ where $z^*$ is the optimal value of either \eqref{obj} or \eqref{obj constrained}. We will show how we can apply the results in Section \ref{EL section} to find the CI for $\mathcal G(\hat x)$.

Consider the optimization problems
\begin{equation}
\begin{array}{ll}
\max/\min_w&\max_{x\in\Theta}\sum_{i=1}^nw_i[H(\hat x;\xi_i)-H(x;\xi_i)]\\
\text{subject to}&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{p+1,\beta}\\
&\sum_{i=1}^nw_i=1\\
&w_i\geq0\text{\ for all\ }i=1,\ldots,n
\end{array}\label{EL optimality gap}
\end{equation}

We have the following guarantee in using \eqref{EL optimality gap} to construct the CI for $\mathcal G(\hat x)$ for \eqref{obj}:
\begin{thm}
Suppose $\xi_1,\ldots,\xi_n\in\mathbb R^d$ are i.i.d. data independent of a given solution $\hat x$. Let $\mathcal G(\hat x)$ be the optimality gap of $\hat x$ for \eqref{obj}, and $\overline z$ and $\underline z$ be the maximum and minimum values of the programs in \eqref{EL optimality gap} respectively. Suppose Assumption \ref{EL assumptions} holds except that in Condition \ref{rank}, we consider the covariance matrix of $(\nabla_x H(x^*;\xi),H(x^*;\xi)-H(\hat x;\xi))$ instead. We have
\begin{equation}
\liminf_{n\to\infty}P\left(\mathcal G(\hat x)\in[\underline z,\overline z]\right)\geq1-\beta.\label{guarantee optimality gap}
\end{equation}
\label{EL thm optimality gap}
\end{thm}

\begin{proof}
Let $\bar H(x;\xi)=H(x;\xi)-H(\hat x;\xi)$, and $\bar h(x)=E[\bar H(x;\xi)]=h(x)-h(\hat x)$. We verify that Assumption \ref{EL assumptions}, with the change that the covariance matrix of $(\nabla_x H(x^*;\xi),H(x^*;\xi)-H(\hat x;\xi))$ is considered instead in Condition \ref{rank}, implies that $\bar h$ and $\bar H$ satisfies Assumption \ref{EL assumptions} too with $h$ and $H$ replaced by $\bar h$ and $\bar H$ and $\epsilon$ replaced by $2\epsilon$.

\noindent\underline{Condition 1:} We have $\nabla_x\bar h(x)=\nabla_x(h(x)-h(\hat x))=\nabla_xh(x)=E[\nabla_xH(x;\xi)]=E[\nabla_x(H(x;\xi)-H(\hat x;\xi))]=E[\nabla_x\bar H(x;\xi)]$.

\noindent\underline{Condition 2:} We have $x^*\in\text{argmin}_{x\in\Theta}\bar h(x)\Leftrightarrow x^*\in\text{argmin}_{x\in\Theta}h(x)\Leftrightarrow\nabla_xh(x^*)=0\Leftrightarrow\nabla_x\bar h(x^*)=0$. Similarly, $x^*\in\text{argmin}_{x\in\Theta}\bar{\tilde h}(x)\Leftrightarrow x^*\in\text{argmin}_{x\in\Theta}\tilde h(x)\Leftrightarrow\nabla_x\tilde h(x^*)=0\Leftrightarrow\nabla_x\bar{\tilde h}(x^*)=0$ for any $\bar{\tilde h}(x)=\tilde h(x)-\tilde h(\hat x)$ that satisfies
$$\sup_{x\in\Theta}|\bar{\tilde h}(x)-\bar h(x)|\leq\sup_{x\in\Theta}|\tilde h(x)-h(x)|+|\tilde h(\hat x)-h(\hat x)|<2\epsilon$$

\noindent\underline{Condition 3:} By our modification of this condition we have the covariance of $(\nabla_x\bar H(x^*;\xi),\bar H(x^*;\xi))=(\nabla_xH(x^*;\xi),H(x^*;\xi)-H(\hat x;\xi))$ finite and having a positive rank.

\noindent\underline{Condition 4:} It is straightforward to show that $\frac{1}{n}\sum_{i=1}^n\bar H(x;\xi_i)=\frac{1}{n}\sum_{i=1}^nH(x;\xi_i)-\frac{1}{n}\sum_{i=1}^nH(\hat x;\xi_i)\to0$ a.s. uniformly over $x\in\Theta$.

\noindent\underline{Condition 5:} We have $E[\sup_{x\in\Theta}\bar H(x;\xi)^2]=E[\sup_{x\in\Theta}(H(x;\xi)-H(\hat x;\xi))^2]\leq 4(E[\sup_{x\in\Theta}H(x;\xi)^2]+E[H(\hat x;\xi)^2])<\infty$.

We have therefore verified our claim. Using Theorem \ref{EL thm}, we get that
$$\liminf_{n\to\infty}P(\underline v\leq\bar h(x^*)\leq\overline v)\geq1-\beta$$
where $\overline v$ and $\underline v$ are the maximum and minimum values of
\begin{equation}
\begin{array}{ll}
\max/\min_w&\min_{x\in\Theta}\sum_{i=1}^nw_i\bar H(x;\xi_i)\\
\text{subject to}&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{p+1,\beta}\\
&\sum_{i=1}^nw_i=1\\
&w_i\geq0\text{\ for all\ }i=1,\ldots,n
\end{array}\label{EL optimality gap1}
\end{equation}
Noting that $\mathcal G(x^*)=-\bar h(x^*)$, we get \eqref{guarantee optimality gap} immediately.
\end{proof}

For the optimality gap of the stochastically constrained problem \eqref{obj constrained}, we use the following optimization problems
\begin{equation}
\begin{array}{ll}
\max/\min_w&\left\{\begin{array}{ll}\max_x&\sum_{i=1}^nw_i[H(\hat x;\xi_i)-H(x;\xi_i)]\\
\text{subject to}&\sum_{i=1}^nw_iF_k(x;\xi_i)\leq0,\ k=1,\ldots,m\\
&g_k(x)\leq0,\ k=1,\ldots,s
\end{array}\right\}\\
\text{subject to}&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{p+m+1,\beta}\\
&\sum_{i=1}^nw_i=1\\
&w_i\geq0\text{\ for all\ }i=1,\ldots,n
\end{array}\label{optimality gap constrained}
\end{equation}

We have the following guarantee in parallel to the deterministically constrained case:
\begin{thm}
Suppose $\xi_1,\ldots,\xi_n\in\mathbb R^d$ are i.i.d. data independent of a given solution $\hat x$. Let $\mathcal G(\hat x)$ be the optimality gap of $\hat x$ for \eqref{obj constrained}, and $\overline z$ and $\underline z$ be the maximum and minimum values of the programs \eqref{optimality gap constrained} respectively. Suppose Assumption \ref{assumptions constrained} holds except that in Condition 3, we consider the covariance matrix of $H(x^*;\xi)-H(\hat x;\xi)$, $\frac{\partial}{\partial x_j}H(x^*;\xi)+\sum_{k=1}^m\lambda_k^*\frac{\partial}{\partial x_j}F_k(x^*;\xi)$, and $F_k(x^*;\xi)$, for all indices $j$ and $k$ corresponding to the active set of the KKT condition for \eqref{obj constrained}. We have
$$\liminf_{n\to\infty}P\left(\mathcal G(\hat x)\in[\underline z,\overline z]\right)\geq1-\beta.$$\label{EL thm optimality gap constrained}
\end{thm}

\begin{proof}
The proof follows verbatim from that of Theorem \ref{EL thm optimality gap}, and noting that the operations involving $f_k$, $F_k$ and $g_k$ are unaffected by the substitution of $h$ and $H$ with $\bar h$ and $\bar H$.
\end{proof}

\section{Numerical Examples}\label{sec:numerics}
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 a stochastically constrained portfolio optimization problem. The latter examples strictly speaking do not satisfy our assumptions, since the first order condition or the KKT condition does not hold for their sample counterparts. 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 (\cite{shapiro2014lectures}, Theorem 5.7). For deterministically constrained problems in the form \eqref{obj}, the ($1-\beta$) CI on the optimal value is given by
\begin{equation}\label{CLT CI}
\left[\hat z_n^*\pm z_{1-\beta/2}\frac{\hat\sigma(\hat x_n^*)}{\sqrt n}\right]
\end{equation}
where $z_{1-\beta/2}$ is the critical value of the standard normal distribution with confidence $1-\beta$, $\hat x^*_n$ is the empirical optimal solution obtained from \eqref{obj empirical}, $\hat z^*_n = (1/n)\sum_{i=1}^nH(\hat x^*;\xi_i)$ is the empirical optimal value, and $\hat\sigma(\hat x^*_n) = \sqrt{(1/(n-1))\sum_{i=1}^n(H(\hat x^*_n; \xi_i)-\hat z^*_n)^2}$ is the empirical standard deviation of $H(\hat x^*_n;\xi)$. Since $\hat z^*_n$ is a low biased estimator of $z^*$, the CI (\ref{CLT CI}) suffers from the under coverage issue. So we also compare with a 2-sample CLT (CLT2) method, as suggested by \cite{mak1999}, which uses first half 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
\begin{equation} \label{CLT2 CI}
\left[\hat z^*_{n/2} - z_{1-\beta/2}\frac{\hat\sigma(\hat x^*_{n/2})}{\sqrt {n/2}}, \bar{z}^*_{n/2} + z_{1-\beta/2}\frac{\bar\sigma(\hat x^*_{n/2})}{\sqrt {n/2}}\right]
\end{equation}
where $z^*_{n/2}, \hat x^*_{n/2}, \sigma(\hat x^*_{n/2})$ are computed as before using first half of the data $\{\xi_1,\ldots,\xi_{n/2}\}$, $\bar{z}^*_{n/2}=(2/n)\sum_{i=\frac{n}{2}+1}^nH(\hat x^*_{n/2};\xi_i)$ is the evaluation of $\hat x^*_{n/2}$ using the remaining half of the data, and $\bar\sigma(\hat x^*_{n/2})=\sqrt{(1/(n/2-1))\sum_{i=\frac{n}{2}+1}^n(H(\hat x^*_{n/2}; \xi_i)-\bar z^*_{n/2})^2}$ is the empirical standard deviation at $\hat{x}^*_{n/2}$. Note that $ \bar{z}^*_{n/2}$ is a high biased estimator of $z^*$, and thus the CI (\ref{CLT2 CI}) 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 \cite{Bayraksan2006} to estimate CIs on the optimality gap. For a given solution $\hat{x}$ that is independent of the data, the SRP outputs a one-sided ($1-\beta$) CI on the optimality gap given by
$$
\left[0, \hat{\mathcal G}_n(\hat{x}) + z_{1-\beta}\frac{\tilde\sigma(\hat x_n^*)}{\sqrt n}\right],
$$
where as before $\hat x^*_n$ is the empirical optimal solution, $\hat{\mathcal G}_n(\hat{x})= (1/n)\sum_{i=1}^n\left(H(\hat x, \xi_i)-H(\hat x^*_n,\xi_i)\right)^2$, and $$\tilde \sigma^2(\hat x_n^*) = \frac{1}{n-1}\sum_{i=1}^n\left[(H(\hat x,\xi_i)-H(\hat x_n^*,\xi_i))-(\hat h(\hat x) - \hat z^*_n) \right]^2,$$ where $\hat z^*_n$ is the empirical optimal value and $\hat h(\hat x) = (1/n)\sum_{i=1}^nH(\hat x,\xi_i)$. In all the examples considered below, we set $\beta = 0.05$. Note that all the above discussion holds for deterministically constrained problems. Nonetheless, we also apply these methods in a stochastically constrained problem as a benchmark (which is heuristic since there has been no formal proof of their validity in this case).

%Before going into the numerical details, we discuss the tractability of the optimization pairs \eqref{EL opt} and \eqref{opt constrained}, each of which consists of a max-min and a min-min problem. For the unconstrained case \eqref{EL opt}, if $H$ is a convex function and $\Theta$ is a convex set, then
Note that the EL method consists of solving a max-min and a min-min problem. Supposing that the original problem \eqref{obj} or \eqref{obj constrained} 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 challenging because the outer optimization involves minimizing the concave function $\min_{x\in\Theta}\sum_{i=1}^nw_iH(x;\xi_i)$ over $w$. This is not a convex problem in general. However, fixing either $w$ or $x$, optimizing over the other 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. Such type of schemes has appeared in chance-constrained programming (e.g., \cite{chen2010cvar,zymler2013distributionally,jiang2012data}), and it appears to work well in our examples despite a lack of global convergence guarantee. %All the discussion here holds for the constrained case \eqref{opt constrained} if the constraints are all convex.


\subsection{Quadratic Optimization}\label{sec:quad opt}
We consider a simple unconstrained problem of minimizing a quadratic function
\begin{equation}\label{quad}
\min_{x} E[(x-\xi)^2],
\end{equation}
where $\xi$ follows an unknown distribution $F^c$. It is easy to see that the optimal solution is $x^* = E[\xi]$ and the optimal value is $z^* = Var(\xi)$.  We set $F^c$ as a standard normal distribution, and thus $x^*=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 \eqref{quad}. We test on three cases where we randomly generate $n= 10, 50, 100$ data points from $F^c$. For each case, we repeat the experiment $100$ times, and note down 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 \ref{comparisons simple}.

\begin{table}[h]
\centering
 % \begin{subtable}[b]{0.45\textwidth}\centering
 % {
 {\small

 \begin{tabular}{cc|ccccc}
 & &Coverage & Mean lower & Mean upper & Mean interval & Standard deviation
\\
& &probability & bound & bound & width & of interval width\\
\hline
$n=10$ & EL& $0.79$ &$0.32$ &$1.89$ &$1.57$ &$0.95$\\
       & CLT&$0.72$ &$0.19$ &$1.69$ &$1.50$ &$0.93$ \\
       & CLT2&$0.80$ &$0.01$ &$2.42$ &$2.41$ &$1.75$ \\ \hline
$n=50$ & EL&$0.97$ &$0.29$ &$1.96$ &$1.67$ &$1.07$\\
       %& CLT&$0.86$&$1.22$&$2.23$&$1.01$&$0.42$ \\
       & CLT&$0.84$ &$0.60$ &$1.31$ &$0.71$ &$0.19$ \\
       &CLT2&$0.87$ &$0.45$ &$1.53$ &$1.08$ &$0.51$ \\ \hline
$n=100$& EL&$0.99$ &$0.66$ &$1.43$ &$0.77$ &$0.27$\\
      % &CLT&$0.89$&$1.36$&$2.07$&$0.71$&$0.21$ \\
       & CLT&$0.84$ &$0.70$ &$1.22$ &$0.52$ &$0.10$ \\
       &CLT2&$0.88$ &$0.60$ &$1.33$ &$0.74$ &$0.29$ \\
\end{tabular}
}
%}
%\end{subtable}
\caption{Confidence intervals on optimal values of the quadratic optimization problem}
\label{comparisons simple}
\end{table}

To compare EL and SRP on optimality gap, we first generate a solution $\hat{x}$ and evaluate its true optimality gap using a large sample size ($10^8$). Then for each of the three cases $n=10, 50, 100$, we repeat the experiment  $100$ times for each method to obtain $95\%$ confidence bounds and estimate their empirical coverage probabilities. The results are summarized in Table \ref{tab:quad_optgap}, where the suboptimal solution $\hat x = 0.62$ and its corresponding optimality gap is $0.39$.

\begin{table}[h]
\centering
 % \begin{subtable}[b]{0.45\textwidth}\centering
 % {
 {\small

 \begin{tabular}{cc|ccccc}
 & &Coverage & Mean lower & Mean upper & Mean interval & Standard deviation
\\
& &probability & bound & bound & width & of interval width\\
\hline
$n=10$ & EL& $0.95$ &$0.06$ &$1.95$ &$1.89$ &$0.94$\\
       & CLT-SRP&$0.86$ &$0$ &$1.16$ &$1.16$ &$0.72$ \\ \hline
%       & CLT2&$$ &$$ &$$ &$$ &$$ \\ \hline
$n=50$ & EL&$0.99$ &$0.09$ &$1.42$ &$1.32$ &$1.54$\\
       %& CLT&$0.86$&$1.22$&$2.23$&$1.01$&$0.42$ \\
       & CLT-SRP&$0.93$ &$0$ &$0.72$ &$0.72$ &$0.26$ \\ \hline
%       &CLT2&$$ &$$ &$$ &$$ &$$ \\ \hline
$n=100$& EL&$0.97$ &$0.13$ &$0.83$ &$0.70$ &$0.50$\\
      % &CLT&$0.89$&$1.36$&$2.07$&$0.71$&$0.21$ \\
       & CLT-SRP&$0.92$ &$0$ &$0.57$ &$0.57$ &$0.16$ \\
%       &CLT2&$$ &$$ &$$ &$$ &$$ \\
\end{tabular}
}
%}
%\end{subtable}
\caption{Confidence intervals on optimality gaps of the quadratic optimization problem}
\label{tab:quad_optgap}
\end{table}

%In this example,  EL in general has the highest coverage probability on optimal values, and its interval width is wider than the direct CLT method but comparable to the 2-sample CLT method. On the optimality gap, EL has higher coverage probability than SRP, while its interval width is also wider. %This is not surprising since both methods relax the bounds.

\subsection{CVaR Estimation}
In this example, we consider estimating CVaR$_{\alpha, F^c}(\xi)$, the $\alpha$-level conditional-value-at-risk of a random variable $\xi$, which we assume follows an unknown distribution $F^c$. This can be rewritten as a stochastic optimization problem:
\begin{equation}
\min_{x\in\mathbb{R}} \left\{x + \frac{1}{1-\alpha}E[(\xi-x)^+]\right\},\label{CVaR opt}
\end{equation}
where $(\cdot)^+$ is short for $\max(\cdot,0)$. We set $F^c$ as a standard normal distribution and $\alpha=0.9$. As the previous example in Section \ref{sec:quad opt}, we run the experiment $100$ times for each method and each case of $n=10, 50, 100$. The results are summarized in Table \ref{tab:CVaR} and \ref{tab:CVaR_optgap}. Note that the true optimal value can be accurately calculated and is equal to $1.755$; the suboptimal solution in this experiment is $0.71$ with optimality gap $0.36$.

%Moreover, we also compare with the CIs obtained from the CLT and the delta method (\cite{shapiro2014lectures}, Theorem 5.7), given by
%\begin{equation}\label{CLT CI}
%\left[\hat H_n^*\pm z_{1-\beta/2}\frac{\hat\sigma(\hat x_n^*)}{\sqrt n}\right]
%\end{equation}
%where $z_{1-\beta/2}$ is the critical value of the standard normal distribution with confidence $1-\beta$, $\hat x^*_n$ is the empirical optimal solution, $\hat H^*_n$ is the empirical optimal value given by $(1/n)\sum_{i=1}^nH(\hat x^*;\xi_i)$, and $\hat\sigma^2(\hat x^*_n)$ is the empirical standard deviation of the objective value at the optimal solution given by $\sqrt{(1/(n-1))\sum_{i=1}^n(h(\hat x^*_n; \xi_i)-\hat H^*_n)^2}$. Since $\hat H^*_n$ is a low biased estimator of $H^*$, the confidence interval (\ref{CLT CI}) suffers from the under coverage issue. So we also compare with a 2-sample CLT method, as suggested by \cite{mak1999}, which uses first half of the data to compute the empirical optimal value and the corresponding solution and then evaluates the solution using the remaining half of the data to generate an upper bound. The 2-sample CLT CI is given by
%\begin{equation} \label{CLT2 CI}
%\left[\hat H^*_{n/2} - z_{1-\beta/2}\frac{\hat\sigma(\hat x^*_{n/2})}{\sqrt {n/2}}, \bar{H}^*_{n/2} + z_{1-\beta/2}\frac{\bar\sigma(\hat x^*_{n/2})}{\sqrt {n/2}}\right]
%\end{equation}
%where $H^*_{n/2}, \hat x^*_{n/2}, \sigma(\hat x^*_{n/2})$ are computed as before using first half of the data $\{\xi_1,\ldots,\xi_{n/2}\}$, $\bar{H}^*_{n/2}$ is the evaluation of $\hat x^*_{n/2}$ using the remaining half of the data, i.e., $(2/n)\sum_{i=\frac{n}{2}+1}^nH(\hat x^*_{n/2};\xi_i)$, and $\bar\sigma^2(\hat x^*_{n/2})$ is the empirical standard deviation at $\hat{x}^*_{n/2}$ given by $\sqrt{(1/(n/2-1))\sum_{i=\frac{n}{2}+1}^n(h(\hat x^*_{n/2}; \xi_i)-\bar H^*_{n/2})^2}$. Note that $ \bar{H}^*_{n/2}$ is a high biased estimator of $H^*$, and thus the confidence interval (\ref{CLT2 CI}) aims to alleviate 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 confidence interval.

%We test on three cases when the data size $n$ is $10$, $50$ and $100$, respectively. For each case, we repeat the experiment $100$ times, and note down the empirical coverage probability, mean upper and lower bounds, and the mean and standard deviation of the interval width for each method.

\begin{table}[h]
\centering
 % \begin{subtable}[b]{0.45\textwidth}\centering
 % {
 {\small

 \begin{tabular}{cc|ccccc}
 & &Coverage & Mean lower & Mean upper & Mean interval & Standard deviation
\\
& &probability & bound & bound & width & of interval width\\
\hline
$n=10$ & EL&$0.39$ &$0.80$&$1.65$&$0.85$&$0.56$\\
       & CLT&$0.50$&$0.95$&$2.16$&$1.21$&$1.02$ \\
       & CLT2&$0.47$&$1.14$&$4.03$&$2.89$&$4.17$ \\ \hline
$n=50$ & EL&$0.90$ &$1.21$&$2.31$&$1.10$&$0.40$\\
       %& CLT&$0.86$&$1.22$&$2.23$&$1.01$&$0.42$ \\
       & CLT&$0.81$&$1.23$&$2.16$&$0.93$&$0.35$ \\
       &CLT2&$0.78$&$0.99$&$2.55$&$1.56$&$1.22$ \\ \hline
$n=100$& EL&$0.98$ &$1.34$&$2.28$&$0.94$&$0.27$\\
      % &CLT&$0.89$&$1.36$&$2.07$&$0.71$&$0.21$ \\
       & CLT&$0.86$&$1.35$&$2.06$&$0.71$&$0.20$ \\
       &CLT2&$0.88$&$1.19$&$2.35$&$1.15$&$0.63$ \\
\end{tabular}
}
%}
%\end{subtable}
\caption{Confidence intervals on optimal values of the CVaR estimation problem}
\label{tab:CVaR}
\end{table}

\begin{table}[h]
\centering
 % \begin{subtable}[b]{0.45\textwidth}\centering
 % {
 {\small

 \begin{tabular}{cc|ccccc}
 & &Coverage & Mean lower & Mean upper & Mean interval & Standard deviation
\\
& &probability & bound & bound & width & of interval width\\
\hline
$n=10$ & EL& $0.96$ &$0.02$ &$3.49$ &$3.47$ &$2.08$\\
       & CLT-SRP&$0.83$ &$0$ &$1.70$ &$1.70$ &$1.38$ \\ \hline
%       & CLT2&$$ &$$ &$$ &$$ &$$ \\ \hline
$n=50$ & EL&$1.00$ &$0.03$ &$1.46$ &$1.43$ &$0.50$\\
       %& CLT&$0.86$&$1.22$&$2.23$&$1.01$&$0.42$ \\
       & CLT-SRP&$0.85$ &$0$ &$0.81$ &$0.81$ &$0.43$\\ \hline
%       &CLT2&$$ &$$ &$$ &$$ &$$ \\ \hline
$n=100$& EL&$0.99$ &$0.07$ &$1.05$ &$0.98$ &$0.24$\\
      % &CLT&$0.89$&$1.36$&$2.07$&$0.71$&$0.21$ \\
       & CLT-SRP&$0.91$ &$0$ &$0.71$ &$0.71$ &$0.26$ \\
%       &CLT2&$$ &$$ &$$ &$$ &$$ \\
\end{tabular}
}
%}
%\end{subtable}
\caption{Confidence intervals on optimality gaps of the CVaR estimation problem}
\label{tab:CVaR_optgap}
\end{table}

%We note that in this example EL, CLT, and CLT2 are relatively comparable. When there are relatively few data ($n=10$), the CLT and 2-sample CLT methods have higher empirical coverage probability, but their CIs are also much wider on average with a much larger standard deviation in interval widths. When the data size increases, the coverage probability of the EL method improves significantly, while the mean interval width is slightly larger than the CLT but smaller than the 2-sample CLT. As in the previous example, EL has higher coverage probability with a wider interval width than SRP in all three cases.

\subsection{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 denote by $x = [x^1,\ldots,x^d]'$ the vector of holding proportions in $d$ assets, $\xi=[\xi^1,\ldots,\xi^d]'$ the random vector of asset returns, and $r_b$ the threshold for expected return.  We assume short selling is not allowed. The problem can be written as
\begin{equation} \label{CVaRmean}
\begin{array}{ll}
  \min_{x} &  CVaR_{\alpha}(-\xi'x)  \\
  \text{subject to} & E[\xi'x] \geq r_b \\
   & \sum_{i=1}^dx_i = 1 \\
   & x_i\geq 0, i=1,\ldots,d \\
\end{array}
\end{equation}
We can rewrite the problem in the form of (\ref{obj constrained}) as
\begin{equation} \label{CVaRmean}
\begin{array}{ll}
  \min_{x,c} & c+ \frac{1}{1-\alpha}E[(-\xi'x-c)^+]  \\
  \text{subject to} & E[\xi'x] \geq r_b \\
   & \sum_{i=1}^dx_i = 1 \\
   & x_i\geq 0, i=1,\ldots,d \\
\end{array}
\end{equation}

%the objective function can be rewritten as $\min_{x, c}~~~ c + \frac{1}{1-\alpha}E[(-\xi'x-c)^{+}].$ When the distribution of $\xi$ is unknown and we only have data $\xi_1,\ldots,\xi_n$,  the SAA version of (\ref{CVaRmean}) is
%\begin{equation}
%\begin{array}{ll}
%  \min_{x,c} &  c + \frac{1}{(1-\alpha)n}\sum_{i=1}^n (-\xi_i'x-c)^{+}  \\
% \text{subject to}  &\frac{1}{n}\sum_{i=1}^n \xi_i'x \geq r_b \\
%   & \sum_{i=1}^dx_i = 1 \\
%   & x_i\geq 0, i=1,\ldots,d
%\end{array}
%\end{equation}
%By introducing new variables $b_i\geq (-\xi_i x - c)^{+}, i=1,\ldots,n$, and denoting $b = [b_1,\ldots,b_n]'$, the above problem can be converted to a linear program:
%\begin{equation}
%\begin{array}{ll}
%  \min_{x,c,b} &  c + \frac{1}{1-\alpha}\sum_{i=1}^n b_i  \\
%  \text{subject to} & b_i \geq -\xi_i x - c,~ i=1,\ldots,n \\
%   &(\frac{1}{n}\sum_{i=1}^n \xi_i)'x \geq r_b \\
%   & \sum_{i=1}^dx_i = 1 \\
%   & x_i\geq 0, i=1,\ldots,d \\
%   & b_i \geq 0, i=1,\ldots,n
%\end{array}
%\end{equation}
%To quantify the uncertainty in the function value associated with the SAA solution, we apply the EL method, i.e., solve
%\begin{equation}
%\begin{array}{ll}
%\max/\min_{\mathbf w}&\left\{\begin{array}{ll}\min_{x,c} &  c + \frac{1}{1-\alpha}\sum_{i=1}^n w_i(-\xi_i'x-c)^{+} \\
%\text{subject to}& \sum_{i=1}^n w_i\xi_i'x \geq r_b \\
%   & \sum_{i=1}^dx_i = 1 \\
%   & x_i\geq 0,  i=1,\ldots,d \\
%\end{array}\right\}\\
%\text{subject to}&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{r,\beta}\\
%&\sum_{i=1}^nw_i=1\\
%&w_i\geq0, i=1,\ldots,n
%\end{array}\label{portfolio_EL}
%\end{equation}

%To quantify the uncertainty in the function value associated with the SAA solution, we apply the EL method, i.e., solve
%\begin{equation}
%\begin{array}{ll}
%\max/\min_{\mathbf w}&\left\{\begin{array}{ll}\min_{x,c,b} &  c + \frac{1}{1-\alpha}\sum_{i=1}^n w_ib_i \\
%\text{subject to}& b_i \geq -\xi_i x - c,~ i=1,\ldots,n\\
% & (\sum_{i=1}^n w_i\xi_i)'x \geq r_b \\
%   & \sum_{i=1}^dx_i = 1 \\
%   & x_i\geq 0,  i=1,\ldots,d \\
%   & b_i \geq 0, i=1,\ldots,n
%\end{array}\right\}\\
%\text{subject to}&-2\sum_{i=1}^n\log(nw_i)\leq\chi^2_{r,\beta}\\
%&\sum_{i=1}^nw_i=1\\
%&w_i\geq0, i=1,\ldots,n
%\end{array}\label{opt constrained}
%\end{equation}


The parameter setting is as follows: $\xi$ follows a normal distribution with mean $\mu = [0.8, 1.2]'$ and covariance $\Sigma = [1~ 0; 0~ 4]$; the minimum expected return is $r_b = 1$; the CVaR level is $\alpha = 0.9$, and the confidence level is $1-\beta = 0.95$. It is easy to verify that the optimal solution to (\ref{CVaRmean}) is $x^* = [0.5, 0.5]'$, and the associate optimal value can be evaluated by Monte Carlo simulation with a large number ($10^8$) of samples, which yields $z^*\approx 0.96$. For comparison, we also implement the CLT and 2-sample CLT methods by computing the CIs according to (\ref{CLT CI}) or (\ref{CLT2 CI}); though the validity of these schemes has not been proved, we use them as heuristic to provide a benchmark. For each case of $n=10, 50, 100$, we repeat the experiment 100 times, and summarize the numerical results in Table \ref{tab:portfolio} and \ref{tab:portfolio_optgap}. In this experiment, the suboptimal solution is $[0.21 0.79]$ with optimality gap $0.73$.
\begin{table}[h]
\centering
 % \begin{subtable}[b]{0.45\textwidth}\centering
 % {
 {\small

% \begin{tabular}{cc|ccccc}
% & &Coverage & Mean lower & Mean upper & Mean interval & Standard deviation
%\\
%& &probability & bound & bound & width & of interval width\\
%\hline
%$n=10$ & EL & $0.55$ & $-0.10$ & $4.11$ & $4.21$ & $11.51$\\
%       & CLT& $0.25$ & $0.29$ & $1.29$ & $1.00$ & $1.72$ \\ \hline
%$n=50$ & EL &$0.71$  & $0.08$  & $2.40$ &$ 2.32$ & $3.03$ \\
%       & CLT& $0.40$    & $0.57$ & $1.72$ & $1.16$ & $0.67$ \\ \hline
%$n=100$& EL & $0.78$ &$0.30$   &$2.20$  & $1.90$ & $1.87$\\
%       &CLT & $0.54$    &$0.70$  & $1.57$ & $0.87$ & $0.41$ \\ \hline
%$n=200$& EL & $0.91$ &$0.48$   &$2.31$  & $1.83$ & $2.33$\\
%       &CLT & $0.60$    &$0.86$  & $1.52$ & $0.66$ & $0.23$ \\
%\end{tabular}
%}
%%}
%%\end{subtable}
%\caption{Comparison of EL and CLT methods on the portfolio optimization problem}
%\label{tab:portfolio}
%\end{table}

 \begin{tabular}{cc|ccccc}
 & &Coverage & Mean lower & Mean upper & Mean interval & Standard deviation
\\
& &probability & bound & bound & width & of interval width\\
\hline
$n=10$ & EL & $0.26$ & $-0.01$ & $1.04$ & $1.05$ & $4.66$\\
       & CLT& $0.21$ & $0.33$ & $1.21$ & $0.88$ & $1.61$ \\
       &CLT2& $0.52$ & $0.32$ & $5.46$ & $5.14$ & $6.90$ \\ \hline
$n=50$ & EL &$0.69$  & $0.14$ & $1.20$ & $1.06$ & $0.49$ \\
       & CLT& $0.58$ & $0.60$ & $1.85$ & $1.24$ & $0.64$ \\
       &CLT2& $0.67$ & $0.34$ & $2.45$ & $2.12$ & $1.59$\\ \hline
$n=100$& EL & $0.74$ & $0.28$ & $1.38$ & $1.09$ & $0.60$\\
       &CLT & $0.59$ & $0.73$ & $1.63$ & $0.90$ & $0.34$ \\
       &CLT2& $0.64$ & $0.63$ & $2.02$ & $1.39$ & $0.80$ \\
%$n=200$& EL & $0.91$ &$0.48$   &$2.31$  & $1.83$ & $2.33$\\
%       &CLT & $0.60$    &$0.86$  & $1.52$ & $0.66$ & $0.23$ \\
\end{tabular}
}
%}
%\end{subtable}
\caption{Confidence intervals on optimal values of the portfolio optimization problem}
\label{tab:portfolio}
\end{table}

\begin{table}[h]
\centering
 % \begin{subtable}[b]{0.45\textwidth}\centering
 % {
 {\small

 \begin{tabular}{cc|ccccc}
 & &Coverage & Mean lower & Mean upper & Mean interval & Standard deviation
\\
& &probability & bound & bound & width & of interval width\\
\hline
$n=10$ & EL& $0.56$ &$0.57$ &$2.44$ &$1.87$ &$1.58$\\
       & CLT-SRP& $0.79$ &$0$ &$1.92$ &$1.92$ &$1.56$ \\ \hline
%       & CLT2&$$ &$$ &$$ &$$ &$$ \\ \hline
$n=50$ & EL&$0.90$ &$0.33$ &$1.91$ &$1.59$ &$1.60$\\
       %& CLT&$0.86$&$1.22$&$2.23$&$1.01$&$0.42$ \\
       & CLT-SRP&$0.57$ &$0$ &$0.92$ &$0.92$ &$0.82$\\ \hline
%       &CLT2&$$ &$$ &$$ &$$ &$$ \\ \hline
$n=100$& EL&$0.92$ &$0.48$ &$1.50$ &$1.02$ &$0.50$\\
      % &CLT&$0.89$&$1.36$&$2.07$&$0.71$&$0.21$ \\
       & CLT-SRP&$0.59$ &$0$ &$0.87$ &$0.87$ &$0.61$ \\
%       &CLT2&$$ &$$ &$$ &$$ &$$ \\
\end{tabular}
}
%}
%\end{subtable}
\caption{Confidence intervals on optimality gaps of the portfolio optimization problem}
\label{tab:portfolio_optgap}
\end{table}

\subsection{Summary of Numerical Results}

We note in all three examples EL in general has the highest coverage probability on optimal values except when the data size is very small ($n=10$). Although EL has wider intervals than the direct CLT method, its interval widths are 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 optimal 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 converge probability in the last example is in general smaller compared to the previous two examples. A potential reason (in addition to the assumptions or the validity of the compared methods not being rigorously justified) is the higher dimensionality that naturally requires more data to achieve a similar level of accuracy in the SAA solution. %The constrained problem is also more challenging to solve than the unconstrained problem, which may lead to more numerical imprecisions when applying these methods.
Nevertheless, we can see that the EL method still produces CIs with a higher coverage probability than CLT methods when the data size is not too small, and the coverage probability improves as the data size increases.

% We have similar observation for the results on optimal gaps as the previous examples.


\section{Conclusion}
We have studied the EL method to construct statistically valid CIs for the optimal value and the 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 ball constraints, with the ball size suitably calibrated by a $\chi^2$-quantile with a suitable degree of freedom. We have studied the theory leading to the statistical guarantees and numerically compared our method to approaches suggested by the CLT. Built on a rigorous foundation, our method provides a competitive method for evaluating the statistical uncertainty for stochastic optimization problems under limited data. In future work, we plan to further refine the accuracy of our method.

\section*{Appendix}
\begin{lemma}[Lemma 11.2 in \cite{owen2001empirical}]
Let $Y_i$ be i.i.d. random variables in $\mathbb R$ with $EY_i^2<\infty$. We have $\max_{1\leq i\leq n}|Y_i|=o(n^{1/2})$ a.s..%\label{auxiliary}
\end{lemma}

%\section*{Acknowledgements}
%The research of the first author was partially supported by the National Science Foundation under grants CMMI-1400391/1542020 and CMMI-1436247/1523453. ... A preliminary conference version of this work has appeared in \cite{lam2015quantifying}.

\bibliographystyle{abbrv}
\bibliography{bibliography2}

\end{document} 
ELuncertainty2_Oct17.bbl
\begin{thebibliography}{10}

\bibitem{asmussen2007stochastic}
S.~Asmussen and P.~W. Glynn.
\newblock {\em Stochastic Simulation: Algorithms and Analysis}, volume~57.
\newblock Springer Science \& Business Media, 2007.

\bibitem{Bayraksan2006}
G.~Bayraksan and D.~Morton.
\newblock Assessing solution quality in stochastic programs.
\newblock {\em Mathematical Programming}, 108:495--514, 2006.

\bibitem{ben2013robust}
A.~Ben-Tal, D.~den Hertog, A.~De~Waegenaere, B.~Melenberg, and G.~Rennen.
\newblock Robust solutions of optimization problems affected by uncertain
  probabilities.
\newblock {\em Management Science}, 59(2):341--357, 2013.

\bibitem{bertsimas2016}
D.~Bertsimas, V.~Gupta, and N.~Kallus.
\newblock Robust sample average approximation.
\newblock {\em available at arXiv:1408.4445}, 2016.

\bibitem{chen2010cvar}
W.~Chen, M.~Sim, J.~Sun, and C.-P. Teo.
\newblock From {CV}a{R} to uncertainty set: Implications in joint
  chance-constrained optimization.
\newblock {\em Operations Research}, 58(2):470--485, 2010.

\bibitem{delage2010distributionally}
E.~Delage and Y.~Ye.
\newblock Distributionally robust optimization under moment uncertainty with
  application to data-driven problems.
\newblock {\em Operations Research}, 58(3):595--612, 2010.

\bibitem{diciccio1991empirical}
T.~DiCiccio, P.~Hall, and J.~Romano.
\newblock Empirical likelihood is {B}artlett-correctable.
\newblock {\em The Annals of Statistics}, 19(2):1053--1061, 1991.

\bibitem{Glasserman:1988:PCD:318123.318245}
P.~Glasserman.
\newblock Performance continuity and differentiability in {M}onte {C}arlo
  optimization.
\newblock In {\em Proceedings of the 20th Winter Simulation Conference}, pages
  518--524, New York, NY, USA, 1988. ACM.

\bibitem{jiang2012data}
R.~Jiang and Y.~Guan.
\newblock Data-driven chance constrained stochastic program.
\newblock {\em Mathematical Programming, Series A}, 158(1):291–--327, 2012.

\bibitem{kleywegt2002sample}
A.~J. Kleywegt, A.~Shapiro, and T.~Homem-de{-}Mello.
\newblock The sample average approximation method for stochastic discrete
  optimization.
\newblock {\em SIAM Journal on Optimization}, 12(2):479--502, 2002.

\bibitem{lam2015quantifying}
H.~Lam and E.~Zhou.
\newblock Quantifying uncertainty in sample average approximation.
\newblock In {\em Proceedings of the 2015 Winter Simulation Conference}, pages
  3846--3857. IEEE Press, 2015.

\bibitem{mak1999}
W.-K. Mak, D.~Morton, and R.~Wood.
\newblock Monte {C}arlo bounding techniques for determining solution quality.
\newblock {\em Operations Research Letters}, 24:47--56, 1999.

\bibitem{owen1988empirical}
A.~B. Owen.
\newblock Empirical likelihood ratio confidence intervals for a single
  functional.
\newblock {\em Biometrika}, 75(2):237--249, 1988.

\bibitem{owen2001empirical}
A.~B. Owen.
\newblock {\em Empirical Likelihood}.
\newblock CRC press, 2001.

\bibitem{pardo2005statistical}
L.~Pardo.
\newblock {\em Statistical Inference Based on Divergence Measures}.
\newblock CRC Press, 2005.

\bibitem{shapiro2014lectures}
A.~Shapiro, D.~Dentcheva, and A.~Ruszczy{\'n}ski.
\newblock {\em Lectures on Stochastic Programming: Modeling and Theory},
  volume~16.
\newblock SIAM, 2014.

\bibitem{wang2008sample}
W.~Wang and S.~Ahmed.
\newblock Sample average approximation of expected value constrained stochastic
  programs.
\newblock {\em Operations Research Letters}, 36(5):515--519, 2008.

\bibitem{wiesemann2014distributionally}
W.~Wiesemann, D.~Kuhn, and M.~Sim.
\newblock Distributionally robust convex optimization.
\newblock {\em Operations Research}, 62(6):1358--1376, 2014.

\bibitem{wilks1938large}
S.~S. Wilks.
\newblock The large-sample distribution of the likelihood ratio for testing
  composite hypotheses.
\newblock {\em The Annals of Mathematical Statistics}, 9(1):60--62, 1938.

\bibitem{zymler2013distributionally}
S.~Zymler, D.~Kuhn, and B.~Rustem.
\newblock Distributionally robust joint chance constraints with second-order
  moment information.
\newblock {\em Mathematical Programming, Series A}, 137(1):167--198, 2013.

\end{thebibliography}
