Parametric Scenario Optimization under Limited Data: A Distributionally Robust Optimization View

TL;DR

This paper proposes a distributionally robust optimization framework that leverages parametric assumptions and Monte Carlo simulation to obtain feasible solutions in small-data probabilistic constraint problems.

Contribution

It introduces a novel DRO approach that incorporates parametric information and Monte Carlo methods to improve feasibility guarantees with limited data.

Findings

Parametric assumptions can improve feasibility in small-data scenarios.

Optimal generating distributions differ from baseline in parametric DRO.

Numerical examples demonstrate the effectiveness of the proposed approach.

Abstract

We consider optimization problems with uncertain constraints that need to be satisfied probabilistically. When data are available, a common method to obtain feasible solutions for such problems is to impose sampled constraints, following the so-called scenario optimization approach. However, when the data size is small, the sampled constraints may not statistically support a feasibility guarantee on the obtained solution. This paper studies how to leverage parametric information and the power of Monte Carlo simulation to obtain feasible solutions for small-data situations. Our approach makes use of a distributionally robust optimization (DRO) formulation that translates the data size requirement into a Monte Carlo sample size requirement drawn from what we call a generating distribution. We show that, while the optimal choice of this generating distribution is the one eliciting the data or the baseline distribution in a nonparametric divergence-based DRO, it is not necessarily so in the parametric case. Correspondingly, we develop procedures to obtain generating distributions that improve upon these basic choices. We support our findings with several numerical examples.

Tables

Table 1. Table 1. Comparisons among choices of ℙ 0 subscript ℙ 0 \mathbb{P}_{0} for 10 dimensional multivariate Gaussian when n = 5 𝑛 5 n=5 .

	${𝒟}_{data}\left({\mathbb{P}}_0\right)$	${\delta }_{ϵ}$	$N$ for Extended SO	$N$ for Extended FAST
${\mu }_1\left({\delta }_{\hat{\theta }}\right)$	37.9161	$6.5766\times {10}^{-5}$	285601	70221
${\mu }_2$	11.0368	$2.2454\times {10}^{-4}$	83649	20707
${\mu }_3$	14.7391	$1.6850\times {10}^{-4}$	111465	27528

Table 2. Table 2. Comparisons among choices of ℙ 0 subscript ℙ 0 \mathbb{P}_{0} for 10 dimensional multivariate Gaussian when n = 10 𝑛 10 n=10 .

	${𝒟}_{data}\left({\mathbb{P}}_0\right)$	${\delta }_{ϵ}$	$N$ for Extended SO	$N$ for Extended FAST
${\mu }_1\left({\delta }_{\hat{\theta }}\right)$	5.2383	$4.6857\times {10}^{-4}$	40081	10026
${\mu }_2$	3.3139	$7.3298\times {10}^{-4}$	25621	6481
${\mu }_3$	3.7926	$6.4275\times {10}^{-4}$	29219	7363

Table 3. Table 3. Single linear CCP under Gaussian with unknown mean for different ϵ italic-ϵ \epsilon and d 𝑑 d .

	$ϵ=0.1$	$ϵ=0.1$	$ϵ=0.1$	$ϵ=0.05$	$ϵ=0.05$	$ϵ=0.05$
	$d=5$	$d=10$	$d=20$	$d=5$	$d=10$	$d=20$
$n$	50	80	200	50	80	200
$N_{exact}$	113	183	312	229	371	631
$N$	449	743	1016	1443	2349	3118
$\widehat{ϵ}$	0.0050	0.0041	0.0041	0.0015	0.0015	0.0014
$Q_{95}$	0.0136	0.0103	0.0088	0.0045	0.0037	0.0031
$f_{val}$	-0.7577	-0.7447	-0.7360	-0.7353	-0.7243	-0.7128

Table 4. Table 4. Comparisons for single linear CCP under Gaussian: ϵ = 0.05 italic-ϵ 0.05 \epsilon=0.05 , d = 10 𝑑 10 d=10 and β = 0.001 𝛽 0.001 \beta=0.001 .

	RO	Extended SO	Extended FAST
$n$	80	80	80
$N_{exact}$	NA	447	447
$N$	NA	2887	1079
$\widehat{ϵ}$	0.0180	0.0011	0.00069
$Q_{95}$	0.0272	0.0029	0.0019
$f_{val}$	-0.8008	-0.7212	-0.7093

Table 5. Table 5. Comparisons for single linear CCP under Exponential: ϵ = 0.05 italic-ϵ 0.05 \epsilon=0.05 , d = 10 𝑑 10 d=10 and β = 0.001 𝛽 0.001 \beta=0.001 .

	RO	Extended SO	Extended FAST
$n$	80	80	80
$N_{exact}$	NA	447	447
$N$	NA	2887	1079
$\widehat{ϵ}$	0.0045	0.0047	0.0016
$Q_{95}$	0.0094	0.0100	0.0050
$f_{val}$	-0.6978	-0.6981	-0.6701

Table 6. Table 6. Comparisons for Joint linear CCP under Gaussian: ϵ = 0.05 italic-ϵ 0.05 \epsilon=0.05 , m = 3 𝑚 3 m=3 , d = 10 𝑑 10 d=10 and β = 0.001 𝛽 0.001 \beta=0.001 .

	RO	Extended SO	Extended FAST
$n$	80	80	80
$N_{exact}$	NA	291	291
$N$	NA	2388	1214
$\widehat{ϵ}$	0.0003	0.0012	0.0226
$Q_{95}$	0.0007	0.0033	0.0564
$f_{val}$	-0.6448	-0.6626	-0.6466

Table 7. Table 7. Comparisons for quadratic objective with joint linear chance constraint under Gaussian: ϵ = 0.05 italic-ϵ 0.05 \epsilon=0.05 , m = 5 𝑚 5 m=5 , d = 10 𝑑 10 d=10 and β = 0.001 𝛽 0.001 \beta=0.001 .

	RO	Extended SO	Extended FAST
$n$	200	200	200
$N_{exact}$	NA	447	447
$N$	NA	3888	1384
$\widehat{ϵ}$	0	0.0006	0.0096
$Q_{95}$	0	0.0017	0.0253
$f_{val}$	-0.4800	-0.5547	-0.5476

Table 8. Table 8. Linear objective with quadratic chance constraint for different ϵ italic-ϵ \epsilon , m 𝑚 m and d 𝑑 d .

	$ϵ=0.1,d=5,m=5$	$ϵ=0.05,d=10,m=10$	$ϵ=0.05,d=15,m=15$
$n$	80	200	300
$N_{exact}$	113	371	504
$N$	154	334	422
$\widehat{ϵ}$	0.0092	0.0050	0.0048
$Q_{95}$	0.0263	0.0133	0.0128
$f_{val}$	-0.8393	-0.8576	-0.8672

Equations

$$\underset{x\in \mathcal{X}\subseteq {\mathbb{R}}^d}{\min }$$

$$\mathbb{P}(x\in {\mathcal{X}}_{\xi })\ge 1-ϵ,$$

$$\underset{x\in \mathcal{X}\subseteq {\mathbb{R}}^d}{\min }$$

$$x\in {\mathcal{X}}_{{\xi }_i},i=1,\dots ,n.$$

$$V(x,\mathbb{P})=\mathbb{P}(x\notin {\mathcal{X}}_{\xi })$$

$$V(x,\mathbb{P})\le ϵ.$$

$${\mathbb{P}}_{data}(V(\hat{x},\mathbb{P})\le ϵ)\ge 1-\alpha ,$$

$${\mathbb{P}}_{data}(V(\hat{x},\mathbb{P})>ϵ)\le \sum _{i=0}^{d-1}\left(\genfrac{}{}{0ex}{}{n}{i}\right){ϵ}^i(1-ϵ{)}^{n-i},$$

$$B(ϵ,d,n)=\sum _{i=0}^{d-1}\left(\genfrac{}{}{0ex}{}{n}{i}\right){ϵ}^i(1-ϵ{)}^{n-i}\le \alpha ,$$

$${\mathbb{P}}_{data}(\mathbb{P}\in {\mathcal{U}}_{data})\ge 1-\alpha ,$$

$$\underset{x\in \mathcal{X}\subseteq {\mathbb{R}}^d}{\min }$$

$$\underset{\mathbb{Q}\in {\mathcal{U}}_{data}}{\min }\mathbb{Q}(x\in {\mathcal{X}}_{\xi })\ge 1-ϵ,$$

$$\mathbb{P}(\hat{x}\in {\mathcal{X}}_{\xi })\ge \underset{\mathbb{Q}\in {\mathcal{U}}_{data}}{\min }\mathbb{Q}(\hat{x}\in {\mathcal{X}}_{\xi })\ge 1-ϵ$$

$${\mathbb{P}}_{data}(\mathbb{P}(\hat{x}\in {\mathcal{X}}_{\xi })\ge 1-ϵ)\ge {\mathbb{P}}_{data}(\mathbb{P}\in {\mathcal{U}}_{data})\ge 1-\alpha ,$$

$$\underset{x\in \mathcal{X}\subseteq {\mathbb{R}}^d}{\min }$$

$$x\in {\mathcal{X}}_{{\xi }_i^{MC}},i=1,\dots ,N.$$

$$N_{exact}(ϵ,\beta ,d)=\min \left\{n:\sum _{i=0}^{d-1}\left(\genfrac{}{}{0ex}{}{n}{i}\right){ϵ}^i(1-ϵ{)}^{n-i}\le \beta \right\}.$$

$$\underset{\mathbb{Q}\in {\mathcal{U}}_{data}}{\max }$$

$$V({\hat{x}}^{MC},{\mathbb{P}}_0)\le \delta .$$

$$\underset{\mathbb{Q}\in {\mathcal{U}}_{data},A\subset \mathcal{Y}}{\max }$$

$${\mathbb{P}}_0(A)\le \delta .$$

$$\underset{\mathbb{Q}\in {\mathcal{U}}_{data}}{\max }V({\hat{x}}^{MC},\mathbb{Q})\le M({\mathbb{P}}_0,{\mathcal{U}}_{data},V({\hat{x}}^{MC},{\mathbb{P}}_0)),$$

$$M({\mathbb{P}}_0,{\mathcal{U}}_{data},{\delta }_{ϵ})\le ϵ.$$

$${\mathbb{P}}_{data}(\mathbb{P}\in {\mathcal{U}}_{data})\ge 1-\alpha ,$$

$${\mathbb{P}}_{MC,0}(V({\hat{x}}^{MC},{\mathbb{P}}_0)>{\delta }_{ϵ})\le \beta$$

$$V({\hat{x}}^{MC},{\mathbb{P}}_0)\le {\delta }_{ϵ}⟹\underset{\mathbb{Q}\in {\mathcal{U}}_{data}}{\max }V({\hat{x}}^{MC},\mathbb{Q})\le M({\mathbb{P}}_0,{\mathcal{U}}_{data},{\delta }_{ϵ}).$$

$${\mathbb{P}}_{data}\left(\underset{\mathbb{Q}\in {\mathcal{U}}_{data}}{\max }V({\hat{x}}^{MC},\mathbb{Q})>ϵ\right)\le {\mathbb{P}}_{data}(V({\hat{x}}^{MC},{\mathbb{P}}_0)>{\delta }_{ϵ})\le \beta$$

$$\underset{\mathbb{Q}\in {\mathcal{U}}_{data}}{\max }V({\hat{x}}^{MC},\mathbb{Q})\le ϵ⟹V({\hat{x}}^{MC},\mathbb{P})\le ϵ.$$

$$\{\Xi :V({\hat{x}}^{MC},\mathbb{P})>ϵ\}\subseteq \{\Xi :\mathbb{P}\notin {\mathcal{U}}_{data}\}\cup \{\Xi :V({\hat{x}}^{MC},{\mathbb{P}}_0)>{\delta }_{ϵ}\}.$$

$$\underset{n\to \infty }{\mathrm{liminf}}{\mathbb{P}}_{data}(\mathbb{P}\in {\mathcal{U}}_{data})\ge 1-\alpha ,$$

$$\underset{n\to \infty }{\lim }{\mathbb{P}}_{data}({\mathbb{P}}_{{\theta }_{true}}\in {\mathcal{U}}_{data})=1-\alpha$$

$${\mathcal{U}}_{data}=\left\{\mathbb{Q}\in {\mathcal{P}}_{para}:d_{\varphi }({\mathbb{P}}_{{\hat{\theta }}_n},\mathbb{Q})\le \frac{{\varphi }^{\prime \prime }(1){\chi }_{1-\alpha ,D}^2}{2n}\right\},$$

$$d_{\varphi }({\mathbb{P}}_1,{\mathbb{P}}_2)=$$

Full text

chance constraint, distributionally robust optimization, scenario optimization, parametric uncertainty

††copyright: acmcopyright††journal: TOMACS††ccs: Mathematics of computing Probability and statistics††ccs: Computing methodologies Modeling and simulation††ccs: Simulation types and techniques Uncertainty quantification††ccs: Mathematical analysis Mathematical optimization

1. Introduction

We consider optimization problems in the form

[TABLE]

where $\mathbb{P}$ is a probability measure governing the random variable $\xi$ on some space $\mathcal{Y}$ and $\mathcal{X}_{\xi}\subseteq\mathcal{X}\subseteq\mathbb{R}^{d}$ is a set depending on $\xi$. Problem (1.1) enforces a solution $x$ to satisfy $x\in\mathcal{X}_{\xi}$ with high probability, namely at least $1-\epsilon$. This problem is often known as a probabilistically constrained or chance-constrained program (CCP) (Prékopa, 2003). It provides a natural framework for decision-making under stochastic resource capacity or risk tolerance, and has been applied in various domains such as production planning (Murr and Prékopa, 2000), inventory management (Lejeune and Ruszczynski, 2007), reservoir design (Prékopa and Szántai, 1978; Prékopa et al., 1978), communications (Shi et al., 2015), and ranking and selection (Hong et al., 2015).

We focus on the situations where $\mathbb{P}$ is unknown, but some i.i.d. data, say $\xi_{1},\ldots,\xi_{n}$, are available. One common approach to handle (1.1) in these situations is to use the so-called scenario optimization (SO) or constraint sampling (Campi and Garatti, 2011; Nemirovski and Shapiro, 2006). This replaces the unknown constraint in (1.1) with $x\in\mathcal{X}_{\xi_{i}},i=1,\ldots,n$, namely, by considering

[TABLE]

Note that CCP (1.1) is generally difficult to solve even when the set $\mathcal{X}_{\xi}$ is convex for any given $\xi$ and the distribution $\mathbb{P}$ is known (Prékopa, 2003). Thus, the sampled problem (1.2) offers a tractable approximation for the difficult CCP even in non-data-driven situations, assuming the capability to generate these samples.

Our goal is to find a good feasible solution for (1.1) by solving (1.2) under the availability of i.i.d. data described above. Intuitively, as the sample size $n$ increases, the number of constraints in (1.2) increases and one expects them to sufficiently populate the safety set $\{\xi:x\in\mathcal{X}_{\xi}\}$, thus ultimately give rise to a feasible solution for (1.1). To make this more precise, we first mention that because of the statistical noise from the data, one must settle for finding a solution that is feasible with a high confidence. More specifically, define, for any given solution $x$,

[TABLE]

to be the violation probability of $x$ under probability measure $\mathbb{P}$ that generates $\xi$. Obviously, $x$ is feasible for (1.1) if and only if

[TABLE]

We would like to obtain a solution, say $\hat{x}$, from the data such that

[TABLE]

where $\mathbb{P}_{data}$ is the distribution that generates the i.i.d. data $\xi_{i},i=1,\ldots,n$ (each sampled from $\mathbb{P}$), and $1-\alpha$ is a given confidence level (e.g., $\alpha=5\%$). In other words, we want $\hat{x}$ to satisfy the chance constraint in (1.1) with the prescribed confidence.

Under the convexity of $\mathcal{X}_{\xi}$ and mild additional assumptions (namely, that every instance of (1.2) has a feasible region with nonempty interior and a unique optimal solution), the seminal work (Campi and Garatti, 2008) provides a tight estimate on the required data size $n$ to guarantee (1.4). They show that a solution $\hat{x}$ obtained by solving (1.2) satisfies

[TABLE]

with equality held for the class of “fully-supported” optimization problems (Campi and Garatti, 2008). Thus, suppose we have a sample size $n$ large enough such that

[TABLE]

then from (1.5) we have $\mathbb{P}_{data}(V(\hat{x},\mathbb{P})>\epsilon)\leq\alpha$ or (1.4).

However, in small-sample situations in which the data size $n$ is not large enough to support (1.6), the feasibility guarantee described above may not hold. It can be shown (Campi and Garatti, 2008) that the minimum $n$ that achieves (1.4) is linear in $d$ and reciprocal in $\epsilon$, thus may impose challenges especially in high-dimensional and low-tolerance problems. Similar dependence on the key problem parameters also appears in other related methods such as (De Farias and Van Roy, 2004), which uses the Vapnik-Chervonenkis dimension to infer required sample sizes, the sampling-and-discarding approach in (Campi and Garatti, 2011), and the closely related approach using sample average approximation in (Luedtke and Ahmed, 2008). Several recent lines of techniques have been suggested to overcome these challenges and reduce sample size requirements, including the use of support rank and solution-dependent support constraints (Schildbach et al., 2013; Campi and Garatti, 2018), regularization (Campi and Carè, 2013), and sequential approaches (Carè et al., 2014; Calafiore et al., 2011; Chamanbaz et al., 2016; Calafiore, 2017).

In this paper, we offer a different path to alleviate the data size requirement than the above methods, when $\mathbb{P}$ possesses known parametric structures. Namely, we assume $\mathbb{P}\in\{\mathbb{P}_{\theta}\}_{\theta\in\Theta}$ for some parametric family of distribution, where $\mathbb{P}_{\theta}$ satisfies two basic requirements: It is estimatable, i.e., the unknown quantity or parameter $\theta$ can be estimated from data, and simulatable, i.e., given $\theta$, samples from $\mathbb{P}_{\theta}$ can be drawn using Monte Carlo methods. Under these presumptions, our approach turns the CCP (1.1), with an unknown parameter, into a CCP that has a definite parameter and a suitably re-adjusted tolerance level, which then allows us to generate enough Monte Carlo samples and consequently utilize the guarantee provided from (1.5). On a high level, this approach replaces the data size requirement in using (1.2) (or, in fact, any of its variant methods) with a Monte Carlo size requirement, the latter potentially more available given cheap modern computational power. Our methodological contributions consist of the development of procedures, related statistical results on their sample size requirement translations, and also showing some key differences between parametric and nonparametric regimes.

Our approach starts with a distributionally robust optimization (DRO) to incorporate the data-driven parametric uncertainty. The latter is a framework for decision-making under modeling uncertainty on the underlying probability distributions in stochastic problems. It advocates the search for decisions over the worst case, among all distributions contained in a so-called uncertainty set or ambiguity set (e.g., (Wiesemann et al., 2014; Delage and Ye, 2010; Goh and Sim, 2010)). In CCP, this entails a worst-case chance constraint over this set (e.g., (Hanasusanto et al., 2015; Zymler et al., 2013; Hanasusanto et al., 2017; Li et al., 2019; Jiang and Guan, 2016; Zhang et al., 2016; Hu and Hong, 2013; Cheng et al., 2014; Xie and Ahmed, 2018; Chen et al., 2010, 2018; Ji and Lejeune, 2018)). When the uncertainty set covers the true distribution with a high confidence (i.e., the set is a confidence region), then feasibility for the distributionally robust CCP converts into a confidence guarantee on the feasibility for the original CCP. We follow this viewpoint and utilize uncertainty sets in the form of a neighborhood ball surrounding a baseline distribution, where the ball size is measured by a statistical distance (e.g., (Ben-Tal et al., 2013; Petersen et al., 2000; Hansen and Sargent, 2008; Love and Bayraksan, 2015; Dupuis et al., 2016; Lam and Zhou, 2017; Lam, 2019; Hu and Hong, 2013; Gotoh et al., 2018; Duchi et al., 2016; Esfahani and Kuhn, 2015; Blanchet and Murthy, 2016; Blanchet et al., 2016; Gao and Kleywegt, 2016)). In the parametric case, a suitable choice of this distance (such as the $\phi$-divergence that we focus on) allows easy and meaningful calibration of the ball size from the data, so that the resulting DRO provides a provable feasibility conversion to the CCP.

Our next step is to combine this DRO with Monte Carlo sampling and scenario approximation. The definition of DRO means that there are many possible candidate distributions that can govern the truth, whereas the statistical guarantee for SO assumes a specific distribution that generates the data or Monte Carlo samples. To resolve this discrepancy, we select a generating distribution that draws the Monte Carlo samples, and develop a translation of the guarantee from a fixed distribution into one on the DRO. We highlight the benefits in using SO to handle this DRO, as opposed to other potential methods. While there exist many good results on tractable reformulations of DRO for chance constraints (e.g., (Hanasusanto et al., 2015, 2017; Ji and Lejeune, 2018; Li et al., 2019; Xie and Ahmed, 2018)), the reformulation tightness typically relies on using moment-based uncertainty sets and particular forms of the safety condition. Compared to moments, divergence-based uncertainty sets can be calibrated with data to consistently shrink to the true distribution. Importantly, in the parametric case, the calibration of divergence-based sets is especially convenient, and achieves a tight convergence rate by using maximum likelihood theory that efficiently captures parametric information. Our condition for applying SO to this DRO is at the same level of generality as applying SO to an unambiguous CCP, which, as mentioned before, only requires the convexity of $\mathcal{X}_{\xi}$ and mild conditions.

To exploit the full capability of our approach, we investigate the optimal choice of the generating distribution in relation to the target DRO, in the sense of requiring the least Monte Carlo size. We show that, if there is no ambiguity on the distribution (i.e., a standard CCP), or when the uncertainty set of a DRO is constructed via a divergence ball in the nonparametric space, the best generating distribution is, in a certain sense, the true or the baseline distribution at the center of the ball. However, if there is parametric information, the optimal choice of the generating distribution can deviate from the baseline distribution in a divergence-based DRO. We derive these results by casting the problem of selecting a generating distribution into a hypothesis testing problem, which connects the sampling efficiency of the generating distribution with the power of the test and the Neyman-Pearson lemma (Lehmann and Romano, 2006). The results on DRO in particular combine this Neyman-Pearson machinery with the established DRO reformulation of chance constraints in (Jiang and Guan, 2016; Hu and Hong, 2013), with the discrepancy between the best generating distribution and the baseline distribution in the parametric case stemming from the removal of the extremal distributions in the corresponding nonparametric uncertainty set. These connections among hypothesis testing, SO and DRO are, to our best knowledge, the first of its kind in the literature.

Finally, given the non-optimality of the baseline distribution of a divergence-based DRO in generating Monte Carlo samples, we further develop procedures to search over generating distributions that improve upon this baseline. On a high level, this can be achieved by increasing the sampling variability to incorporate the uncertainty of the distributional parameters (one may intuit this from the perspective of a posterior distribution in a Bayesian framework), which is implemented by utilizing suitable mixture distributions. We provide several classes of mixture distributions to attain such a variability enlargement, and study descent-type algorithms to search for good distributions in these classes. In the experiments, we show our methods can be combined with SO or other SO-based methods including FAST (Carè et al., 2014) to solve a variety of optimization problems and data distributions, some of which are not amenable to RO, especially when the objective function is non-linear or the feasible sets are jointly chance-constrained. Furthermore, we also demonstrate how to search for more judicious choices of generating distributions that can significantly reduce the required number of Monte Carlo samples.

We conclude this introduction by briefly discussing a few other lines of related literature. The first is the so-called robust Monte Carlo or robust simulation that, like us, also considers using Monte Carlo sampling together with DRO (Hu et al., 2012; Glasserman and Xu, 2014; Lam, 2016; Hu and Hong, 2015; Lam, 2018; Ghosh and Lam, 2019). However, this literature focuses on approximating DRO with stochasticity in the objective function, and does not study the chance constraint feasibility and SO that constitute our main focus. We also contrast our work with (Lam, 2019) that also considers likelihood theory and utilizes simulation in tackling uncertain constraints. The study (Lam, 2019) focuses on the nonparametric regime and uses the empirical likelihood to construct uncertainty sets. Unlike our work, there is no parametric information there that can be leveraged to overcome sample size requirements in SO. Moreover, the simulation used in (Lam, 2019) is for calibrating the uncertainty set, instead of drawing sample constraints. Next, (Erdoğan and Iyengar, 2006) considers a scenario approach to distributionally robust CCP with an uncertainty set based on the Prohorov distance. Like (De Farias and Van Roy, 2004), (Erdoğan and Iyengar, 2006) utilizes the Vapnik-Chervonenkis dimension in studying feasibility, in contrast to the convexity-based argument in (Campi and Garatti, 2008) that we utilize. More importantly, we aim to optimize the efficiency of Monte Carlo sampling in handling limited-data CCP, thus motivating us to study the choice of distance, calibration schemes, and selection of generating distributions that are different from (Erdoğan and Iyengar, 2006). Finally, a preliminary conference version of this work has appeared in (Lam and Li, 2018), which contains a basic introduction of our framework, without detailed investigation of the optimality of generating distributions, improvement strategies, and extensive numerical demonstrations.

To summarize, our main contributions of this paper are:

1. (1)

We propose a framework to obtain good feasible solutions in data-driven CCPs in small-sample situations, where the data size is insufficient to support the use of SO with valid statistical guarantees. Focusing on the parametric regime, our framework operates by setting up a DRO, with an uncertainty set constructed from parameter estimates using the data, that can in turn be tackled by using SO with Monte Carlo samples. In doing so, our framework effectively leverages the parametric information to convert the SO requirement on the data size into a requirement on the Monte Carlo size, the latter can be much more abundant given cheap modern computational power. The overview of this framework and the DRO construction are in Sections 2.1 and 2.3. 2. (2)

We investigate and present the Monte Carlo size requirements needed to give statistically feasible solutions to the divergence-based DRO used in our framework. This relies on developing an implementable mechanism to connect the sample size requirement for SO, which attempts to solve a CCP with a fixed underlying distribution, to the sample size requirement needed to solve a DRO, by selecting a suitable generating distribution to draw the Monte Carlo samples. This contribution is presented in Section 2.2. 3. (3)

We study the optimality of generating distributions, in a sense of minimizing the Monte Carlo effort that we will describe precisely. In particular, we show that the optimal generating distributions for an unambiguous CCP, and for a distributionally robust CCP with nonparametric divergence-based uncertainty sets, are simply their respective natural choices, namely the original underlying distribution and the baseline distribution (i.e., center of the divergence ball). In contrast, the optimal generating distribution for a distributionally robust CCP in the parametric case is more delicate, and the baseline distribution there can be readily dominated by other generating distributions. These results are derived by bridging the Neyman-Pearson lemma in statistical hypothesis testing with SO and DRO, which appears to be the first of its kind in the literature as far as we know. This contribution is presented in Section 3. 4. (4)

Motivated by the non-optimality of the baseline distribution, we propose several approaches to construct generating distributions that dominate the baseline distributions for parametric DRO, by using mixture schemes that, on a high level, enlarge the variability of the generating distributions. We show how to use descent-type search procedures to construct these distributions. This contribution is presented in Section 4.

Lastly, we also present in full detail our implementation algorithms in Section 5, numerically demonstrate our approach and compare with other methods in Section 6, and conclude in Section 7.

2. From Data-Driven DRO to Scenario Optimization

This section introduces our overall framework. Recall our goal as to find a good feasible solution $\hat{x}$ for (1.1), and suppose that we have an i.i.d. data size $n$ possibly less than the requirement shown in (1.6). As discussed in the introduction, we first formulate a DRO that incorporates the parametric estimation noise and subsequently allows us to resort to Monte Carlo sampling to obtain a feasible solution for (1.1). In the following, Section 2.1 first describes the basic guarantees from DRO. Section 2.2 investigates Monte Carlo sampling that provides guarantees on DRO. Section 2.3 discusses the choice of the uncertainty set.

2.1. Overview of Data-Driven DRO

For concreteness, suppose the unknown true distribution $\mathbb{P}\in\mathcal{P}$, the class of possible probability distributions for $\xi$ (to be specified later). Given the observed data $\xi_{1},...,\xi_{n}$, the basic steps in our data-driven DRO are:

• •

Step 1: Find a data-driven uncertainty set $\mathcal{U}_{data}=\mathcal{U}_{data}(\xi_{1},\ldots,\xi_{n})\subseteq\mathcal{P}$ such that

[TABLE]

where $\mathbb{P}_{data}$ denotes the measure generating the data $\xi_{i},i=1,\ldots,n$.

• •

Step 2: Given $\mathcal{U}_{data}$, set up the distributionally robust CCP:

[TABLE]

where the probability measure $\mathbb{Q}$ is the decision variable in the minimization in the constraint.

• •

Step 3: Find a solution $\hat{x}$ feasible for (2.2).

It is straightforward to see that $\hat{x}$ obtained from the above procedure is feasible for (1.1) with confidence at least $1-\alpha$: If $\mathbb{P}\in\mathcal{U}_{data}$, then any $\hat{x}$ feasible for (2.2) satisfies

[TABLE]

Thus

[TABLE]

which gives our conclusion.

2.2. Monte Carlo Sampling for DRO

To use the above procedure, we need to provide a way to construct the depicted $\mathcal{U}_{data}$ and to find a (confidently) feasible solution for (2.2). We postpone the set construction to the next subsection and focus on finding a feasible solution here. We resort to SO, via Monte Carlo sampling, to handle (2.2). Note that, unlike in the standard SO discussed in the introduction, the distribution $\mathbb{Q}$ here can be any candidate within the set $\mathcal{U}_{data}$. Thus, let us select a generating distribution, called $\mathbb{P}_{0}$ (which can depend on the data), to generate Monte Carlo samples $\xi_{i}^{MC},i=1,\ldots,N$, and solve

[TABLE]

For convenience, denote, for any $\epsilon,\beta>0$,

[TABLE]

From the result of (Campi and Garatti, 2008) discussed in the introduction, using $N_{exact}(\epsilon,\beta,d)$ or more Monte Carlo samples from $\mathbb{P}_{0}$ in (2.4) would give a solution $\hat{x}^{MC}$ that satisfies $V(\hat{x}^{MC},\mathbb{P}_{0})\leq\epsilon$ with confidence level $1-\beta$. This is not exactly the distributionally robust feasibility statement for problem (2.2). To address this discrepancy, we consider, conditional on the data $\xi_{1},\ldots,\xi_{n}$,

[TABLE]

This optimization problem serves to translate a guarantee on the violation probability under $\mathbb{P}_{0}$ to any $\mathbb{Q}$ in $\mathcal{U}_{data}$. If we can bound the optimal value in (2.6), then we can trace back the level of $\delta$ that is required to ensure a chance constraint validity of tolerance level $\epsilon$. However, the event involved in defining $V(\hat{x}^{MC},\mathbb{P}_{0})$ and $V(\hat{x}^{MC},\mathbb{Q})$, namely $\{\xi:\hat{x}^{MC}\notin\mathcal{X}_{\xi}\}$, can be challenging to handle in general. Thus, we relax (2.6) to

[TABLE]

where the decision variables now include the set $A$ in addition to the probability measure $\mathbb{Q}$. Conditional on the data $\xi_{1},\ldots,\xi_{n}$, the optimal value of optimization problem (2.7), which we denote $M(\mathbb{P}_{0},\mathcal{U}_{data},\delta)$, is clearly an upper bound for that of (2.6). In fact, it is also clear from (2.7) that $M(\mathbb{P}_{0},\mathcal{U}_{data},\delta)$ is non-decreasing in $\delta>0$ and

[TABLE]

by simply taking $A=\{\xi:\hat{x}^{MC}\notin\mathcal{X}_{\xi}\}$ and $\delta=V(\hat{x}^{MC},\mathbb{P}_{0})$ in (2.7). We have the following guarantee:

Theorem 2.2.1.

Given $\mathbb{P}_{0}$, $\mathcal{U}_{data}$ and $\epsilon>0$, suppose there exists $\delta_{\epsilon}>0$ small enough such that

[TABLE]

If we solve (2.4) with $N_{exact}(\delta_{\epsilon},\beta,d)$ number of samples drawn from $\mathbb{P}_{0}$, then the obtained solution $\hat{x}^{MC}$ would be feasible for (2.2) with confidence at least $1-\beta$. Furthermore, if

[TABLE]

where $\mathbb{P}_{data}$ is the measure governing the real-data generation under the true distribution $\mathbb{P}$, then the obtained solution $\hat{x}^{MC}$ would be feasible for (1.1) with confidence at least $1-\alpha-\beta$.

Proof.

By results in (Campi and Garatti, 2008), we know that by solving (2.4) with $N_{exact}(\delta_{\epsilon},\beta,d)$ number of samples from $\mathbb{P}_{0}$, the obtained solution $\hat{x}^{MC}$ would satisfy

[TABLE]

where $\mathbb{P}_{MC,0}$ is the measure with respect to the Monte Carlo samples drawn from $\mathbb{P}_{0}$. Moreover, based on the monotonicity property of $M(\cdot)$ and (2.8), we have

[TABLE]

Thus (2.9) implies that

[TABLE]

and hence $\hat{x}^{MC}$ is feasible for (2.2) with confidence at least $1-\beta$. Furthermore, if $\mathbb{P}\in\mathcal{U}_{data}$, then a $\hat{x}^{MC}$ feasible for (2.2) is also feasible for (1.1) since $\max\limits_{\mathbb{Q}\in\mathcal{U}_{data}}V(\hat{x}^{MC},\mathbb{Q})\geq V(\hat{x}^{MC},\mathbb{P})$ and hence

[TABLE]

Thus, if we denote ${\Xi}=\{\xi_{1},...,\xi_{n},\xi^{MC}_{1},...,\xi^{MC}_{N}\}$ to be entire sequence consisting of real data and the generated Monte Carlo samples, it then follows that

[TABLE]

It now follows by (2.10) and (2.11) that $\hat{x}^{MC}$ is feasible for (1.1) with probability at least $1-\alpha-\beta$. ∎

Theorem 2.2.1 can be cast in terms of asymptotic instead of finite-sample guarantees by following the same line of arguments. We summarize it as the following corollary.

Corollary 2.2.2.

In Theorem 2.2.1, if the condition $\mathbb{P}_{data}(\mathbb{P}\in\mathcal{U}_{data})\geq 1-\alpha$ is substituted by the asymptotic condition

[TABLE]

then the feasibility of $\hat{x}^{MC}$ in the last conclusion of Theorem 2.2.1 holds with confidence asymptotically tending to at least $1-\alpha-\beta$.

To summarize, in the presence of data insufficiency, if we choose $\mathcal{U}_{data}$ to satisfy the confidence property (2.1), and are able to evaluate the bounding function $M(\mathbb{P}_{0},\mathcal{U}_{data},\delta)$ that translates the violation probability under $\mathbb{P}_{0}$ to a worst-case violation probability over $\mathcal{U}_{data}$, then we can run SO with $N_{exact}(\delta_{\epsilon},\beta,d)$ Monte Carlo samples from $\mathbb{P}_{0}$ to obtain a solution for (1.1) with confidence $1-\alpha-\beta$.

We also note that the above scheme still holds if the $N_{exact}(\epsilon,\beta,d)$ in (2.5) is replaced by the sample size requirements of other variants of SO (e.g., FAST (Carè et al., 2014)) that are potentially smaller. This works as long as we stay with the same SO-based procedure in using the Monte Carlo samples. For clarity, throughout most of our exposition we will focus on the sample size requirement depicted in (2.5), but we will discuss other variants in our implementation and numerical sections.

Finally, let us take a step back and justify why we use SO to tackle (2.2), as opposed to other potential means. Indeed, as pointed out in the introduction, there exist many good results on tractable reformulations of DRO. As will be discussed in detail in the next subsection, in the present context we will choose an uncertainty set that can leverage parametric information efficiently. Sets based on the neighborhoods of distributions measured by $\phi$-divergences are particularly attractive choices, as they can be calibrated easily (both the ball center and the size) in a way that efficiently uses parametric information. The dependence on the parameter dimension in particular is reflected in the degree of freedom in the $\chi^{2}$-distribution used in the calibrating the ball size, which shrinks to zero at a canonical rate as the data size increases. Other sets, such as moment-based ones, though possibly amenable to tight tractable reformulations, do not enjoy these statistical properties in the parametric context. Thus, in view of tackling $\phi$-divergence-based DRO, SO appears to be a natural choice, and we have set up a framework to utilize it under conditions at the same level of generality as required for the unambiguous counterpart. Sections 3 and 4 will study this framework in further depth and enhance its efficiency. We caution, however, that the conservativeness in our proposed uncertainty set (which affects the optimality of the obtained solution) relies on the dimensionality of the distributional parameters. Our approach is expected to work well when this dimension is moderate, but not in high-dimensional problems where other approaches could be better choices.

2.3. Constructing Uncertainty Sets

In this section we discuss the construction of the uncertainty set $\mathcal{U}_{data}$, using the $\phi$-divergence approach (Ben-Tal et al., 2013). We assume the true distribution $\mathbb{P}$ of $\xi$ lies in a parametric family. We denote the true parameter as $\theta_{true}$. To highlight the parametric dependence, we call the true distribution $\mathbb{P}_{\theta_{true}}\in\mathcal{P}_{para}=\{\mathbb{P}_{\theta}\}_{\theta\in\Theta\subset\mathbb{R}^{D}}$ indexed by $\theta$, where $D$ is the dimension of parameter space. Given data $\xi_{1},\xi_{2},...,\xi_{n}$, we want to construct an uncertainty set $\mathcal{U}_{data}$ satisfying

[TABLE]

so that Corollary 2.2.2 applies. To do so, we first estimate $\theta_{true}$ from the data. There are various approaches to do so; here we apply the common maximum likelihood estimator (MLE) $\hat{\theta}_{n}$, and set $\mathcal{U}_{data}$ to be

[TABLE]

where $\chi^{2}_{1-\alpha,D}$ is the $1-\alpha$ quantile of $\chi^{2}_{D}$, the $\chi^{2}$-distribution with degree of freedom $D$, and $d_{\phi}(\cdot,\cdot)$ is the $\phi$-divergence between two probability measures, i.e., given a convex function $\phi:\mathbb{R}_{+}\to\mathbb{R}_{+}$, with $\phi(1)=0$, a distance between two probability measures $\mathbb{P}_{1}$ and $\mathbb{P}_{2}$ defined as

[TABLE]

assuming $\mathbb{P}_{2}$ is absolutely continuous with respect to $\mathbb{P}_{1}$ with Radon-Nikodym derivative $\frac{d\mathbb{P}_{2}}{d\mathbb{P}_{1}}$ on $\mathcal{Y}$. Moreover, we assume that $\phi$ is twice continuously differentiable with $\phi^{\prime\prime}(1)\neq 0$, and if necessary set the continuation of $\phi$ to $\mathbb{R}_{-}$ as $\phi(x)=+\infty$ for $x<0$. In (2.17), we call the center of the divergence ball, $\mathbb{P}_{\hat{\theta}_{n}}$, the baseline distribution.

To guarantee desirable asymptotic properties of our uncertainty set, we make the following assumption:

A1.

Let $\theta_{true}\in\Theta$ be the true parameter and let $\hat{\theta}_{n}$ be the MLE of $\theta_{true}$ estimated from $n$ i.i.d. data points. Then, as $n\rightarrow\infty$, $\hat{\theta}_{n}$ satisfies consistency and asymptotic normality condition:

[TABLE]

where $\mathcal{I}(\theta)$ is the Fisher information for the parametric family $\mathcal{P}_{para}$ with well-defined inverse that is continuous in the domain $\theta\in\Theta$.

Assumption A1 of MLE estimator is known to hold under various regularity conditions (Van der Vaart, 2000; Lehmann, 2004). We list a set of such conditions in Appendix A.

Under Assumption A1, it can be shown (Pardo, 2005; Van der Vaart, 2000) that $\mathcal{U}_{data}$ in (2.17) satisfies the confidence guarantee (2.16). Furthermore, since we can identify each $\mathbb{P}_{\theta}$ in $\mathcal{P}_{data}$ with $\theta$, we can equivalently view $\mathcal{U}_{data}$ as a subset of $\theta\in\Theta$, and write it as

[TABLE]

For convenience, we shall use the two definitions of $\mathcal{U}_{data}$ interchangeably depending on the context. It is also known that the asymptotic confidence properties of (2.17) or (2.20) are the same among different choices within the $\phi$-divergence class. These can be seen via a second order expansion of the $\phi$-divergences. Moreover, they are asymptotically equivalent to

[TABLE]

where $\mathcal{I}(\hat{\theta}_{n})$ is the estimated Fisher information, under the regularity conditions above (Pardo, 2005; Nielsen and Nock, 2014; Van der Vaart, 2000). In other words, under Assumption A1, both (2.20) and (2.21) satisfy

[TABLE]

Note that the convergence rate of (2.16) or (2.22) depends on the higher-order properties of the parametric model, which in turn can depend on the parameter dimension. Different from the sample size requirements in SO, this convergence rate is a consequence of MLE properties. Some details on finite-sample behaviors of MLE can be found in (Korostelev and Korosteleva, 2011).

The $\mathcal{U}_{data}$ discussed above is a set over the parametric class of distributions (or parameter values). Considering tractability, DRO over nonparametric space could be easier to handle than parametric, which suggests a relaxation of the parametric constraint to estimate the bounding function $M$. This also raises the question of whether one can possibly contain $\mathcal{U}_{data}$ in a nonparametric ball with a shrunk radius and subsequently obtain a better $M$. These would be the main topics of Sections 3 and 4.

3. Bounding Functions and Generating Distributions

Given the uncertainty set $\mathcal{U}_{data}$ in (2.20), we turn to the choice of the generating measure $\mathbb{P}_{0}$ and the bounding function $M(\mathbb{P}_{0},\mathcal{U}_{data},\delta)$ which, as we recall, is the optimal value of optimization problem (2.7). In the discussed parametric setup, the latter becomes

[TABLE]

From Theorem 2.2.1 and the fact that $M(\mathbb{P}_{0},\mathcal{U}_{data},\delta)$ is non-decreasing in $\delta$, we want to choose $\mathbb{P}_{0}$ that minimizes $M(\mathbb{P}_{0},\mathcal{U}_{data},\delta)$ so that we can take the maximum $\delta_{\epsilon}$ and subsequently achieve overall confident feasibility with the least Monte Carlo sample size. Note that $M(\mathbb{P}_{0},\mathcal{U}_{data},\delta)$ is a multi-input function depending on both $\mathbb{P}_{0}$ and $\delta$, and so a priori it is not clear that a uniform minimizer $\mathbb{P}_{0}$ can exist across all values of $\delta$ so that the described task is well-defined. It turns out that this is possible in some cases, which we shall investigate in detail. In the following, we discuss results along this line at three levels: The unambiguous case, namely when $\mathcal{U}_{data}$ in (3.1) is a singleton (Section 3.1), the case where $\mathcal{U}_{data}$ is nonparametric (Section 3.2), and the case where $\mathcal{U}_{data}$ is parametric (Section 3.3). The first two cases pave the way to the last one, which is most important to our development and also motivates Section 4. With these results in hand, we also discuss the possibility of using other statistical distances in our framework in Section 3.4.

3.1. Neyman-Pearson Connections and A Least Powerful Null Hypothesis

We first consider, for a given $\theta_{1}\in\mathcal{U}_{data}$, the optimization problem

[TABLE]

This problem can be viewed as choosing a most powerful decision rule in a statistical hypothesis test. More precisely, one can think of $A$ as a rejection region for a simple test with null hypothesis $\mathbb{P}_{0}$ and alternate hypothesis $\mathbb{P}_{\theta_{1}}$. Subject to a tolerance of $\delta$ Type-I error, optimization problem (3.2) looks for a decision rule that maximizes the power of the test. By the Neyman-Pearson lemma (Lehmann and Romano, 2006), under mild regularity conditions on the parametric family, the optimal set $A^{\star}_{0,\theta_{1},\delta}$ of (3.2) takes the form

[TABLE]

with $K_{0,\theta_{1},\delta}^{\star}$ chosen so that $\mathbb{P}_{0}(A^{\star}_{0,\theta_{1},\delta})=\delta$. Also, then, the optimal value of (3.2) is $\mathbb{P}_{\theta_{1}}(A^{\star}_{0,\theta_{1},\delta})$. Generalizing the above analysis to all $\theta\in\mathcal{U}_{data}$, we conclude that

[TABLE]

is the optimal value of (3.1). These observations will be useful for deriving our subsequent results.

Our goal is to choose $\mathbb{P}_{0}$ to minimize (3.4). To start our analysis, let us first consider the extreme case where the uncertainty set $\mathcal{U}_{data}$ consists of only one point $\mathbb{Q}$. In this case, we look for $\mathbb{P}_{0}$ that minimizes $M(\mathbb{P}_{0},\{\mathbb{Q}\},\delta)$, the optimal value of

[TABLE]

That is, for a given measure $\mathbb{Q}$, we seek for the maximum discrepancy between $\mathbb{Q}$ and $\mathbb{P}_{0}$ over all $\mathbb{P}_{0}$-measure sets that have $\delta$ or less content. This is similar to minimizing the total variation distance between $\mathbb{Q}$ and $\mathbb{P}_{0}$, and hints that the optimal choice of $\mathbb{P}_{0}$ is $\mathbb{Q}$. The following theorem, utilizing the Neyman-Pearson lemma depicted above, confirms this intuition. We remark that the assumptions of the theorem can be relaxed by using more general versions of the lemma, but the presented version suffices for most purposes and also the subsequent examples we will give.

Theorem 3.1.1.

Given a measure $\mathbb{Q}$ with continuous density on $\mathcal{X}$, among all $\mathbb{P}_{0}$ such that $\frac{d\mathbb{Q}}{d\mathbb{P}_{0}}$ exists and is continuous and positive almost surely, the minimum $M(\mathbb{P}_{0},\{\mathbb{Q}\},\delta)$ is obtained by choosing $\mathbb{P}_{0}=\mathbb{Q}$, giving $M(\mathbb{P}_{0},\{\mathbb{Q}\},\delta)=\delta$.

Proof.

Under the assumptions, by the Neyman-Pearson lemma, for a fixed measure $\mathbb{P}_{0}$, the set achieving the optimal value of (3.5) takes the form $A^{\star}=\{\xi\in\mathcal{Y}:\frac{d\mathbb{Q}}{d\mathbb{P}_{0}}(\xi)>K^{\star}\}$ for some $K^{\star}\geq 0$ with $\mathbb{P}_{0}(A^{\star})=\delta$. It then follows that

[TABLE]

Under the absolute continuity assumption, we define

[TABLE]

which can be seen to be a non-increasing function for $K\geq 1$ and a non-decreasing function for $K\leq 1$. To see this, take $K_{1}\geq K_{2}$, and we have

[TABLE]

Thus, when $K_{1}\geq K_{2}\geq 1$, we have $g(K_{2})\geq g(K_{1})$ because

[TABLE]

while when $1\geq K_{1}\geq K_{2}$, we have $g(K_{2})\leq g(K_{1})$ because

[TABLE]

Then, to identify the minimum of $g(K)$, we either decrease $K$ from 1 to 0 which gives

[TABLE]

by using the dominated convergence theorem (e.g., by considering the set $\{1>d\mathbb{Q}/d\mathbb{P}_{0}(\xi)>K\}$) or we increase $K$ from 1 to $\infty$ which gives

[TABLE]

by Fatou’s lemma. Observations (3.6) and (3.7) suggest that $g(K)\geq 0$ for all $K\geq 0$ and imply that $g(K^{\star})\geq 0$. Thus, we must have $M(\mathbb{P}_{0},\{\mathbb{Q}\},\delta)\geq\delta$. Note that this holds for any $\mathbb{P}_{0}$. Now, since choosing $\mathbb{P}_{0}=\mathbb{Q}$ gives $M(\mathbb{Q},\{\mathbb{Q}\},\delta)=\delta$, an optimal choice of $\mathbb{P}_{0}$ is $\mathbb{Q}$. ∎

Theorem 3.1.1 shows that under mild regularity conditions, in terms of choosing the generating distribution $\mathbb{P}_{0}$ and minimizing $M(\mathbb{P}_{0},\{\mathbb{Q}\},\delta)$, we cannot do better than simply choosing $\mathbb{Q}$ itself. This means that if we had known the true distribution was $\mathbb{Q}$, and without additional knowledge of the event of interest, the safest choice (in the minimax sense) for sampling would be $\mathbb{Q}$, a quite intuitive result. In the language of hypothesis testing, given the simple alternate hypothesis $\mathbb{Q}$, the null hypothesis $\mathbb{P}_{0}$ that provides the least power for the test, i.e., makes it most difficult to distinguish between the two hypotheses, is $\mathbb{Q}$.

3.2. Nonparametric DRO

Building on the discussion in Section 3.1, we now consider the choice of generating distribution $\mathbb{P}_{0}$ to minimize the bounding function obtained from (3.1). Before so, we first discuss the nonparametric case, where the analog of (3.1) is in the form:

[TABLE]

for some ball radius $\lambda>0$, where the decision variables are $\mathbb{Q}$ in the space of all distributions absolutely continuous with respect to $\mathbb{P}_{\hat{\theta}}$, and $A$.

We show that the above setting can be effectively reduced to the unambiguous case, i.e., when $\mathbb{Q}$ lies in a singleton discussed in Section 3.1. This comes from an established equivalence between a distributionally robust chance constraint and an unambiguous chance constraint evaluated by the center of the divergence ball, when the event $A$ is fixed (Hu and Hong, 2013; Jiang and Guan, 2016). In particular, suppose the stochasticity space is $\mathcal{Y}=\mathbb{R}^{k}$, and $\mathbb{P}_{\hat{\theta}}$ admits a density $p_{\hat{\theta}}$. Theorem 1 in (Jiang and Guan, 2016) shows that for any $A$,

[TABLE]

where $\epsilon^{\prime}=\epsilon^{\prime}(\epsilon,\lambda,\phi)>0$ can be explicitly determined by $\epsilon$, $\lambda$ and $\phi$ as

[TABLE]

with $\phi^{*}(t)=\sup_{x}\{tx-g(x)\}$ being the conjugate function of $\phi$ and $\underline{m}(\phi^{*})=\sup\{m\in\mathbb{R}:\phi^{*}\text{\ is a finite constant on\ }(-\infty,m]\}$, $\overline{m}(\phi^{*})=\inf\{m\in\mathbb{R}:\phi^{*}(m)=+\infty\}$, $\ell_{\phi}=\lim_{x\to+\infty}\phi(x)/x$, and $\pi=-\infty$ if $Leb\{[p_{\hat{\theta}}=0]\}=0$, [math] if $Leb\{[p_{\hat{\theta}}=0]\}>0$ and $Leb\{[p_{\hat{\theta}}=0]\setminus A\}=0$, and $1$ otherwise, where $Leb\{\cdot\}$ is the Lebesgue measure on $\mathbb{R}^{k}$.

The above equivalence can be used to obtain the following result.

Theorem 3.2.1.

Suppose $\mathcal{Y}=\mathbb{R}^{k}$ and $\mathbb{P}_{\hat{\theta}}$ admits a density. Among all $\mathbb{P}_{0}$ such that $\frac{d\mathbb{P}_{\hat{\theta}}}{d\mathbb{P}_{0}}$ exists and is continuous, positive almost surely, an optimal choice of $\mathbb{P}_{0}$ that minimizes $M(\mathbb{P}_{0},\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda\},\delta)$, namely the optimal value of (3.8), is the center of the $\phi$-divergence ball $\mathbb{P}_{\hat{\theta}}$. Moreover, this gives $M(\mathbb{P}_{0},\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda\},\delta)={\epsilon^{\prime}}^{-1}(\delta,\lambda,\phi)$, where ${\epsilon^{\prime}}^{-1}(\cdot,\lambda,\phi)$ is the inverse of the function $\epsilon^{\prime}=\epsilon^{\prime}(\epsilon,\lambda,\phi)$ defined in (3.10) with respect to $\epsilon$, given by

[TABLE]

Proof.

From Theorem 1 in (Jiang and Guan, 2016), we know that, for any $A\subset\mathcal{Y}$ and $0\leq\epsilon\leq 1$, (3.9) holds. We can rewrite the optimal value of problem (3.8) in the form:

[TABLE]

which, according to (3.9), has the same optimal value as

[TABLE]

Since, fixing $\phi$ and $\lambda$, $\epsilon^{\prime}$ is a non-decreasing function of $\epsilon$, we see that minimizing $\epsilon$ is equivalent to minimizing $\epsilon^{\prime}$. Denoting $\nu^{\star}$ as the optimal value of the optimization problem

[TABLE]

then the optimal value of (3.13) is ${\epsilon^{\prime}}^{-1}(\nu^{*},\lambda,\phi)$. Moreover, this is achievable by setting $\mathbb{P}_{0}=\mathbb{P}_{\hat{\theta}}$ that gives the optimal value $\nu^{*}=\delta$ to (3.14) by Theorem 3.1.1.

∎

An implication of Theorem 3.2.1 is that, by noting that a parametric divergence ball lies inside a corresponding nonparametric ball, we can compute a bound for $M$ to obtain a required Monte Carlo size, drawn from the baseline $\mathbb{P}_{\hat{\theta}}$, to get a feasible solution for the distributionally robust CCP (2.2) and subsequently the CCP (1.1). More precisely, recall the bounding function $M(\mathbb{P}_{0},\mathcal{U}_{data},\delta)=M(\mathbb{P}_{0},\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda,\mathbb{Q}\in\mathcal{P}_{para}\},\delta)$ with $\lambda=\phi^{\prime\prime}(1)\chi^{2}_{1-\alpha,D}/(2n)$, given by (3.1), as the optimal value of

[TABLE]

We have:

Corollary 3.2.2.

Given a data size $n$, suppose $\mathcal{Y}=\mathbb{R}^{k}$ and $\mathbb{P}_{\hat{\theta}}$ admits a density, where $\hat{\theta}$ is the MLE under Assumption A1. If we choose $\delta_{\epsilon}=\epsilon^{\prime}(\epsilon,\phi^{\prime\prime}(1)\chi^{2}_{1-\alpha,D}/(2n),\phi)$ and draw $N_{exact}(\delta_{\epsilon},\beta,d)$ Monte Carlo samples from the generating distribution $\mathbb{P}_{\hat{\theta}}$ to construct the sampled problem (2.4), then the obtained solution will be feasible for (1.1) with asymptotic confidence level at least $1-\alpha-\beta$.

Proof.

Note that a parametric divergence ball lies inside a corresponding nonparametric ball in the sense that

[TABLE]

Thus, by the definition of $M$, we have

[TABLE]

In particular,

[TABLE]

where the equality follows from Theorem 3.2.1. Thus, if we choose $\delta_{\epsilon}$ such that ${\epsilon^{\prime}}^{-1}(\delta_{\epsilon},\lambda,\phi)\leq\epsilon$, or $\delta_{\epsilon}=\epsilon^{\prime}(\epsilon,\lambda,\phi)$, where $\lambda=\phi^{\prime\prime}(1)\chi^{2}_{1-\alpha,D}/(2n)$ as presented in (2.17), and the generating distribution as $\mathbb{P}_{\hat{\theta}}$, then Corollary 2.2.2 guarantees that running SO on $N_{exact}(\delta_{\epsilon},\beta,d)$ Monte Carlo samples gives a feasible solution for (1.1) with confidence asymptotically at least $1-\alpha-\beta$. ∎

Corollary 3.2.2 thus provides an implementable procedure to handle (1.1) through (2.2).

3.3. Parametric DRO

Next we discuss further the choice of generating distributions in parametric DRO beyond $\mathbb{P}_{\hat{\theta}}$. While the ball center $\mathbb{P}_{\hat{\theta}}$ is a valid choice, the equivalence relation (3.9) does not apply when the divergence ball is in a parametric class, and the optimal choice of the generating distribution may no longer be $\mathbb{P}_{\hat{\theta}}$, as shown in the next result.

Theorem 3.3.1.

In terms of selecting a generating distribution $\mathbb{P}_{0}$ to minimize $M(\mathbb{P}_{0},\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda,\mathbb{Q}\in\mathcal{P}_{para}\},\delta)$, the optimal value of (3.15), the choice $\mathbb{P}_{\hat{\theta}}$ can be strictly dominated by other distributions.

Intuitively, Theorem 3.3.1 arises because the extreme distribution that achieves the equivalence relation (3.9) may not be in the considered parametric family. It implies more flexibility in choosing the generating measure $\mathbb{P}_{0}$, in the sense of requiring less Monte Carlo samples than using $\mathbb{P}_{\hat{\theta}}$.

From the standpoint of hypothesis testing in Section 3.1, the imposed minimax problem (3.15) in searching for the best $\mathbb{P}_{0}$ can be viewed as finding a simple null hypothesis that is uniformly least powerful across the uncertainty set. This question is related and appears more general than finding the least favorable or powerful prior in testing against composite null hypothesis (Lehmann and Romano, 2006). In the latter context, given a set $\Theta_{1}$, one aims to find a distribution $\mu^{\star}(d\theta_{0})$ such that $\Gamma(\mu^{\star})\leq\Gamma(\mu)$ for all distributions $\mu(d\theta_{0})$ on $\Theta_{0}$, where $\Gamma(\mu)$ is the optimal value of

[TABLE]

The distribution $\mu(d\theta_{0})$ is interpreted as a prior on a composite null hypothesis parametrized by $\theta_{0}$, and $\mu^{*}(d\theta_{0})$ is the least favorable prior. The difference between (3.16) and our formulation (3.15) lies in the restriction to measures of the form $\mathbb{P}_{0}=\int_{\Theta_{0}}\mathbb{P}_{\theta_{0}}\mu(d\theta_{0})$ for the former, leading to a smaller search space than ours. This mixture-type $\mathbb{P}_{0}$ and the Bayesian connection will partly motivate our investigation in Section 4.

To prove Theorem 3.3.1, we present a counter example and also some related discussion.

Example 3.3.2.

Consider the uncertainty set $\mathcal{U}_{data}=\{\mathbb{P}_{\theta},:-1\leq\theta\leq 1\}$ within Gaussian location family on $\mathbb{R}$ with $\mathbb{P}_{\theta}(dy)=\frac{1}{\sqrt{2\pi}}e^{-\frac{(y-\theta)^{2}}{2}}$. This can be thought of, e.g., as an uncertainty set based on the $\chi^{2}$-distance, the latter defined between two probability measures $\mathbb{P}_{1}$ and $\mathbb{P}_{2}$ as

[TABLE]

Note that the $\chi^{2}$-distance is in the family of $\phi$-divergences, by choosing $\phi=(x-1)^{2}$. We aim to find a generating distribution $\mathbb{P}_{0}$ to minimize $M(\mathbb{P},\{\mathbb{P}_{\theta}:\theta\in\mathcal{U}_{data}\},\delta)$, the optimal value of

[TABLE]

We consider several symmetric distributions as $\mathbb{P}_{0}$ (symmetry is reasonably conjectured as a good property since an imbalanced shift might increase the power for the alternative hypothesis on one side and the worst case overall). We list these symmetric distributions in increasing variability:

[TABLE]

Given $0\leq\theta\leq 1$, it can be shown by the Neyman-Pearson lemma that the rejection region $A^{\star}$ (i.e. the set giving the optimal value of (3.18) for a given $\theta$) for $\mathbb{P}_{0}^{1}$ has the form $\{y:y>c_{1}\}$, for $\mathbb{P}_{0}^{2}$ the form $\{y:\lvert y-2\theta\rvert\leq c_{2}\}$ and for $\mathbb{P}_{0}^{3}$ the form $\{y:\frac{e^{\theta y}}{e^{y}+e^{-y}}>c_{3}\}$, for some $c_{1}$, $c_{2}$ and $c_{3}$. Let $\delta=0.05$ be the tolerance level, it can be shown through numerical verification that

[TABLE]

Thus, the natural choice $\mathbb{P}_{\hat{\theta}}=\mathbb{P}_{0}^{1}$ based on relaxing to nonparametric DRO yields a bounding function $M(\cdot)$ that is outperformed by $\mathbb{P}_{0}^{2}$ or $\mathbb{P}_{0}^{3}$. Later in Section 4 we will see numerically how $\mathbb{P}_{0}^{2}$ and $\mathbb{P}_{0}^{3}$ can lead to a smaller sample size requirements.

Although Theorem 3.3.1 reveals room to search for the best generating distribution, the involved optimization, or even just finding an improved distribution over $\mathbb{P}_{\hat{\theta}}$, appears to be nontrivial. In particular, the maximization problem in (3.15) depends on the computation of $A^{\star}$ for each alternative of $\theta\in\mathcal{U}_{data}$. Section 4 discusses some approaches to search for improvements. We conclude the current section with some discussion on the choice of statistical distances used in the uncertainty set.

3.4. Choice of Statistical Distance

We have chosen to use $\phi$-divergence to construct our uncertainty set $\mathcal{U}_{data}$, and we have seen how this allows us to effectively translate sample size requirements from the data to Monte Carlo. Note that another common type of distance is the Wasserstein distance (e.g., (Esfahani and Kuhn, 2015; Blanchet and Murthy, 2016; Gao and Kleywegt, 2016)). If one can translate the violation probability under a generating distribution into the worst-case violation probability over a Wasserstein ball, then the same line of arguments in Section 2 applies to using SO on this DRO. Presuming that the size of a parametric Wasserstein-based confidence region can be properly calibrated from data, it is conceivable that the above can give rise to an alternate solution route. It is known (Theorem 3 in (Blanchet and Murthy, 2016)), under suitable regularity conditions, that one can equate a Wasserstein-ambiguous probability$\sup\limits_{d_{W}(\mathbb{Q},\mathbb{P}_{\hat{\theta}})\leq\lambda}\mathbb{Q}(\xi\in{A})$, where $d_{W}$ denotes a Wasserstein distance of order 1 and cost function $c$, and $A$ is an event, to $\mathbb{P}_{\hat{\theta}}(c(\xi,A)\leq 1/\nu^{*})$ where $\nu^{*}\geq 0$ is a dual multiplier for the associated optimization problem, and $c(\xi,A)$ denotes the cost-induced distance between a point $\xi$ and a set $A$. Thus, $M(\mathbb{P}_{0},\{\mathbb{Q}:d_{W}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda\},\delta)$ can be written as

[TABLE]

Compared to the evaluation of $M(\mathbb{P}_{0},\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda\},\delta)$ in Theorem 3.2.1, the tightening of the tolerance level from $\epsilon$ to $\epsilon^{\prime}$ is now replaced by the set inflation from $A$ to the $(1/\nu^{*})$-neighborhood of $A$ given by $\{\xi:c(\xi,A)\leq 1/\nu^{*}\}$. Note that, regardless of the distance used, one could reduce the conservativeness of our analysis by focusing on $A$ in the form $\{x\notin\mathcal{X}_{\xi}\}$, but this would require looking at the specific form of the safety set $\mathcal{X}_{\xi}$.

4. Improving Generating Distributions

This section discusses some approaches to search for better generating distributions beyond the baseline distribution in a divergence ball of DRO. Section 4.1 first states a general result to create better generating distributions. Section 4.2 then specializes to using a mixture distribution on $\theta$ to exploit this result. Sections 4.3 and 4.4 then provide two specific ways to construct these mixtures. Finally, Section 4.5 demonstrates some numerical comparisons in using these new mixing generating distributions and also simply using the baseline.

4.1. A Framework to Reduce Divergence Ball Size by Incorporating Parametric Information

The reason why the best choice of generating distribution $\mathbb{P}_{0}$ is not the baseline of the divergence ball, $\mathbb{P}_{\hat{\theta}}$, in minimizing $M(\mathbb{P}_{0},\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda,\mathbb{Q}\in\mathcal{P}_{para}\},\delta)$ is that the equivalence relation (3.9) does not hold when $\mathbb{Q}$ is restricted to a parametric class. In some sense the reduction to the unambiguous chance constraint in the right hand side of (3.9) is over-conservative as it does not account for parametric information. Suppose we would still like to use the analytically tractable relation (3.9), but at the same time be less conservative. Then, one approach is to find a new baseline distribution, say $\tilde{\mathbb{P}}$, such that the parametrically restricted divergence ball $\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda,\mathbb{Q}\in\mathcal{P}_{para}\}$ lies inside a new nonparametric divergence ball at the center $\tilde{\mathbb{P}}$, namely $\{\mathbb{Q}:d_{\phi}(\tilde{\mathbb{P}},\mathbb{Q})\leq\tilde{\lambda}\}$. If we can obtain a nonparametric ball size $\tilde{\lambda}$ such that $\tilde{\lambda}<\lambda$ and the set inclusion holds, then this new ball is also a valid uncertainty set, and, when simply setting the generating distribution as $\mathbb{P}_{0}=\tilde{\mathbb{P}}$ and applying Theorem 3.2.1, we have a smaller upper bound for $M(\mathbb{P}_{0},\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda,\mathbb{Q}\in\mathcal{P}_{para}\},\delta)$ than ${\epsilon^{\prime}}^{-1}(\delta,\lambda,\phi)$ obtained from using Theorem 3.2.1 directly with the parametric constraint relaxed.

To above mechanism can be executed as follows. Let $\mathcal{U}_{data}=\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda,\mathbb{Q}\in\mathcal{P}_{para}\}$. For any $\mathbb{P}_{0}$, let

[TABLE]

Then we clearly have

[TABLE]

since the right-hand-side set includes distributions outside of the parametric family as well.

Our goal is to find $\mathbb{P}_{0}$ to minimize $\mathcal{D}_{data}(\mathbb{P}_{0},\phi)$ or any upper bound of $\mathcal{D}_{data}(\mathbb{P}_{0},\phi)$ so that it is smaller than the ball size $\lambda$ appearing in the original parametric divergence ball $\mathcal{U}_{data}$. We state the implication of this as follows:

Theorem 4.1.1.

Suppose $\mathcal{Y}=\mathbb{R}^{k}$ and $\mathbb{P}_{\hat{\theta}}$ admits a density. Consider the parametric divergence ball $\mathcal{U}_{data}=\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda,\mathbb{Q}\in\mathcal{P}_{para}\}$. Suppose we can find $\mathbb{P}_{0}$ such that $\mathcal{D}_{data}(\mathbb{P}_{0},\phi)$ defined in (4.1) satisfies $\mathcal{D}_{data}(\mathbb{P}_{0},\phi)<\lambda$. Then we have

[TABLE]

and

[TABLE]

where ${\epsilon^{\prime}}^{-1}(\epsilon,\lambda,\phi)$ is defined in (3.11).

Proof.

By the definition of $\mathcal{D}_{data}(\mathbb{P}_{0},\phi)$, (4.2) holds. Together with the condition $\mathcal{D}_{data}(\mathbb{P}_{0},\phi)<\lambda$, we have the set inclusions

[TABLE]

The inequalities (4.3) then follow from the definition of $M$. The inequalities (4.4) in turn follow immediately from Theorem 3.2.1. ∎

Theorem 4.1.1 stipulates that choosing $\mathbb{P}_{0}$ depicted in the theorem as the generating distribution, and setting ${\epsilon^{\prime}}^{-1}(\delta,\mathcal{D}_{data}(\mathbb{P}_{0},\phi),\phi)$ as an upper bound for $M(\mathbb{P}_{0},\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda,\mathbb{Q}\in\mathcal{P}_{para}\},\delta)$ to obtain the required Monte Carlo size $N_{exact}(\delta_{\epsilon},\beta,d)$ implied by Corollary 2.2.2, will give a lighter Monte Carlo requirement than using the bound ${\epsilon^{\prime}}^{-1}(\delta,\lambda,\phi)$ directly obtained by relaxing the parametric constraint and using $\mathbb{P}_{\hat{\theta}}$ as the generating distribution as in Corollary 3.2.2.

4.2. Mixture as Generating Distribution

Since optimization (4.1) can be difficult to solve generally, we focus on finding improved generating distribution $\mathbb{P}_{0}$ so that the implication of Theorem 4.1.1 holds, instead of fully optimizing (4.1). In this and the next subsections, we design a search space $\mathcal{P}_{0}$ for $\mathbb{P}_{0}$ that allows the construction of tractable procedures to achieve such improvements, while at the same time ensures the obtained $\mathbb{P}_{0}$ are amenable to Monte Carlo simulation.

From now on we will focus on $\chi^{2}$-distance as our choice of $\phi$ for convenience (as will be seen). Suppose that $\mathbb{P}_{\theta}$ has density $p(y;\theta)$. We then set $\mathcal{P}_{0}$ to be the collection of distributions with densities in the form

[TABLE]

for some probability measure $\mu$ on $\Theta$. This class of distributions is easy to sample assuming $p(y;\theta)$ and $\mu$ are, as one can first sample $\theta\sim\mu(d\theta)$ and then $\xi\sim\mathbb{P}_{\theta}$ given $\theta$.

Searching for the best $p_{0}(y)$ requires minimizing $\mathcal{D}_{data}(\mathbb{P}_{0})$ over $\mathbb{P}_{0}\in\mathcal{P}_{0}$ (where for convenience we denote $\mathcal{D}_{data}(\mathbb{P}_{0})$ as $\mathcal{D}_{data}(\mathbb{P}_{0},\phi)$ with $\phi$ representing the $\chi^{2}$-distance). We first use (3.17) to write

[TABLE]

Denoting $\mathcal{P}(\Theta)$ as the space of probability measures on $\Theta$, we define the function $L:\mathcal{P}(\Theta)\times\Theta\rightarrow\mathbb{R}$ to be

[TABLE]

assuming the integral is well-defined for $\mathcal{P}(\Theta)\times\Theta$ and further define

[TABLE]

Thus (4.2) can be written as $\mathcal{D}_{data}(\mathbb{P}_{0})=l(\mu)-1$, and minimizing $\mathcal{D}_{data}(\mathbb{P}_{0})$ is equivalent to solving

[TABLE]

Optimization (4.10) has the following convexity property:

Lemma 4.2.1.

The outer minimization in problem (4.10) is convex.

Lemma 4.2.1 can be proved by direct verification, which is shown in Appendix C. Note also that, if $\mu$ is the point mass $\delta_{\theta}$ for $\theta\in\Theta$, then the mixture distribution would recover the parametric distribution $\mathbb{P}_{\theta}$. Hence the proposed family $\mathcal{P}_{0}$ includes $\{\mathbb{P}_{\theta}\}_{\theta\in\Theta}$, and in particular the original baseline distribution $\mathbb{P}_{\hat{\theta}}$. Although the outer minimization of (4.10) is a convex problem, computing $l(\mu)$ involves a non-convex optimization and is difficult in general. Our approach is to search for a descent direction for the convex function $l(\cdot)$ from $\delta_{\hat{\theta}}$. In the following, we will study two types of search directions, each using its own version of Danskin’s Theorem (Bertsekas, 1971; Bertsekas et al., 2003). To proceed, we introduce the following definition:

Definition 4.2.2.

Define $\Theta^{*}(\mu)$ to be the set of optimal points for the maximization problem in ${l}(\mu)=\sup\limits_{\theta\in\mathcal{U}_{data}}L(\mu,\theta)$ given $\mu\in\mathcal{P}(\Theta)$ :

[TABLE]

It can be shown that $\Theta^{*}(\mu)$ is non-empty and $\Theta^{*}(\mu)\subseteq\mathcal{U}_{data}$ because $\mathcal{U}_{data}$ is compact and ${L}(\mu,\theta)$ is continuous in $\theta$.

4.3. Mixing with a Proposed Distribution

We consider mixing distributions in the form $(1-t)\delta_{\hat{\theta}}+t\mu_{prop}$ for some proposed distribution $\mu_{prop}$, and look for a descent direction by varying $t$ from 0 to 1. We have the following result that is a consequence of Danskin’s Theorem that involves a one-sided derivative. We provide proofs both for this theorem and our following result in Appendix C.

Theorem 4.3.1.

Fix any $\mu_{1},\mu_{2}\in\mathcal{P}(\Theta)$ and $\theta\in\Theta$. Under the assumptions that $\psi(t)=L((1-t)\mu_{1}+t\mu_{2},\theta)$ is well defined for $0\leq t\leq 1$, we know that the function $g(y,t)$

[TABLE]

is integrable for $t\in[0,1]$. If we further assume that there exists a integrable function $g_{0}(y)$ such that

[TABLE]

then we have the right derivative of $\psi(t)$ at $t=0$ given by

[TABLE]

The quantity $\psi^{+}(0)$ is the directional derivative of $L(\mu_{1})$ in the direction $\mu_{2}-\mu_{1}$. Thus, to improve on $\mathcal{D}_{data}(\mathbb{P}_{\hat{\theta}})$, we can propose a mixing distribution $\mu_{prop}(d\theta^{\prime})$, and substitute $\mu_{1}=\delta_{\hat{\theta}}$ and $\mu_{2}=\mu_{prop}$ in (4.12) to check if

[TABLE]

which indicates a strict descent for $l(\cdot)$ from $\delta_{\hat{\theta}}$ to $\mu_{prop}$. In this case, it follows from the convexity of $l(\cdot)$ that we can find some $0<t\leq 1$ such that $l((1-t)\delta_{\hat{\theta}}+t\mu_{prop})<l(\delta_{\hat{\theta}})$, so that

[TABLE]

gives rise to $\mathcal{D}_{data}(\mathbb{P}_{0})<\mathcal{D}_{data}(\mathbb{P}_{\hat{\theta}})$. Finding such a $t$ can be done by a bisection search or enumerating $\mathcal{D}_{data}(\mathbb{P}_{0})$ on $p_{t}$ over a grid of $t$. Note that the above can be implemented only if (4.13) can be verified and also if $\mathcal{D}_{data}(\mathbb{P}_{0})$ is computable. We will show that both properties are satisfied for the case of multivariate Gaussian when $\mu_{prop}$ is properly chosen. In particular, we will identify general sufficient conditions for $\mu_{prop}$ to guarantee (4.13), and also find $\mu_{prop}$ such that the maximization involved in computing $\mathcal{D}_{data}(\mathbb{P}_{0})$ in (4.2) can be reduced to a one-dimensional problem.

Consider a multivariate Gaussian distribution with unknown mean $\Theta\subset\mathbb{R}^{D}$ in an open convex set with density

[TABLE]

where $\Sigma$ is a fixed positive semi-definite covariance matrix. Direct verification (in Appendix C) shows that

[TABLE]

and thus

[TABLE]

We propose the following $\mu_{prop}$. First, we call a distribution on $\Theta$ symmetrical around $\theta\in\Theta$ if its probability density or mass function has the same value for any $\theta_{1},\theta_{2}\in\Theta$ such that $\theta=\frac{\theta_{1}+\theta_{2}}{2}$.

Proposition 4.3.2.

Let $\mu_{prop}(d\theta^{\prime})$ be any symmetrical distribution around $\hat{\theta}$. Given $\theta\in\Theta^{*}(\delta_{\hat{\theta}})$, we define $Y_{\theta}=(\theta-\hat{\theta})^{\intercal}\Sigma^{-1}(\theta^{\prime}-\hat{\theta})$ with $\theta^{\prime}\sim\mu_{prop}(d\theta^{\prime})$. Suppose there exists an integrable random variable $Y$ under the measure $\mu_{prop}$ such that $e^{2Y_{\theta}}\leq Y$ for all $\theta\in\Theta^{*}(\delta_{\hat{\theta}})$. If, for each $\theta\in\Theta^{*}(\delta_{\hat{\theta}})$, $Y_{\theta}$ does not equal to 0 with probability 1, then (4.13) holds and the mixture distribution produced by $\mu_{prop}(d\theta)$ would result in a descent direction on $\mathcal{D}_{data}(\mathbb{P}_{\hat{\theta}})$.

One can check that any Gaussian distribution with mean $\hat{\theta}$ satisfies the conditions of Proposition 4.3.2, and so does any $\mu_{prop}(d\theta^{\prime})$ that is discrete, symmetrical around $\hat{\theta}$, whose outcome directions $\theta^{\prime}-\theta$ constitute a basis of $\mathbb{R}^{D}$. Alternately, we also consider the following continuous $\mu_{prop}$. We set $\theta^{\prime}\sim\hat{\theta}+\sqrt{\frac{\chi^{2}_{1-\alpha,D}}{n}}\cdot\Sigma^{1/2}\eta$ where $\eta$ is a random vector uniformly distributed on the surface of the $D$-dimension unit ball. Note that this $\theta^{\prime}$ can be efficiently simulated by sampling $D$ independent standard Gaussian random variables and scaling their norm to unit length to obtain $\eta$. While this $\mu_{prop}$ can be readily checked to satisfy the conditions in Proposition 4.3.2, we also provide an alternate proof on the validity of this $\mu_{prop}$ in achieving a descent direction in Lemma C.0.5 in the Appendix, as results proven therein provide important reference to calculations in numerical experiments regarding $\mu_{prop}$.

Next, we discuss the computation of $\mathcal{D}_{data}(\mathbb{P}_{0})$ for a given $\mathbb{P}_{0}$. First, we call a random variable $Y$ on $\mathcal{Y}\subset\mathbb{R}^{k}$ rotationally invariant if $Y{\buildrel\mathcal{D}\over{=}}Q^{\intercal}Y$ for any rotational matrix $Q\in\mathbb{R}^{k\times k}$. Using this notion, the following shows how one can reduce the $D$-dimensional maximization problem in the definition of $\mathcal{D}_{data}(\mathbb{P}_{0})$ into a one-dimensional problem.

Proposition 4.3.3.

Given a nominal distribution $Y\sim\mathbb{P}_{0}$ and a multivariate Gaussian family with known covariance $\Sigma$ denoted $\mathbb{P}_{\theta}=\mathcal{N}(\theta,\Sigma)$. If the nominal distribution $Y\sim\mathbb{P}_{0}$ satisfies the condition that the random variable $Z=\Sigma^{-1/2}(Y-\hat{\theta})$ is rotationally invariant, then for any $\theta_{1},\theta_{2}$ satisfying $(\theta_{1}-\hat{\theta})^{\intercal}\Sigma^{-1}(\theta_{1}-\hat{\theta})=(\theta_{2}-\hat{\theta})^{\intercal}\Sigma^{-1}(\theta_{2}-\hat{\theta})$, we have

[TABLE]

Thus, for $\mathcal{D}_{data}(\mathbb{P}_{0})=\max_{\theta\in\mathcal{U}_{data}}\chi^{2}(\mathbb{P}_{0},\mathbb{P}_{\theta})$ with $\mathcal{U}_{data}=\{\theta\in\Theta:(\theta-\hat{\theta})^{\intercal}\Sigma^{-1}(\theta-\hat{\theta})\leq\lambda\}$ as in (4.3), we have

[TABLE]

given any $\theta^{\star}$ satisfying $(\theta^{\star}-\hat{\theta})^{\intercal}\Sigma^{-1}(\theta^{\star}-\hat{\theta})=\lambda$.

Proposition 4.3.4.

Given $0\leq t\leq 1$ and $\mu_{prop}(d\theta)=\hat{\theta}+\sqrt{\frac{\chi^{2}_{1-\alpha,D}}{n}}\cdot\Sigma^{1/2}\eta$, where $\eta$ is a random vector uniformly distributed on the surface of the $D$-dimension unit ball, the nominal measure $\mathbb{P}_{t}$ with density

[TABLE]

satisfies the conditions in Proposition 4.3.3.

Therefore, in computing $\mathcal{D}_{data}(\mathbb{P}_{0})$ derived from the proposed distribution $\mu_{prop}(d\theta)=\hat{\theta}+\sqrt{\frac{\chi^{2}_{1-\alpha,D}}{n}}\cdot\Sigma^{1/2}\eta$, using Propositions 4.3.3 and 4.3.4 we can change the domain of the involved maximization from $\Theta\subset\mathbb{R}^{D}$ into $\mathbb{R}$, leading to a substantial reduction in the search space and a tractable problem.

4.4. Enlarging Mixture Variability

Our next proposal is to consider a continuous mixing distribution $\mu_{r}(d\theta^{\prime})$ on $\Theta$ where $r\geq 0$ controls the variability of the distribution, so that $r=0$ corresponds to $\delta_{\hat{\theta}}$. Here, we can parametrize the density of the generating distribution as

[TABLE]

and our search direction is along $r$ starting from $r=0$. We propose two possible ways to define $\mu_{r}(d\theta^{\prime})$. First is to let $\mu^{1}_{r}(d\theta^{\prime})$ follow the distribution of $\theta^{\prime}\sim\hat{\theta}+\Sigma^{\frac{1}{2}}\cdot\eta_{\sqrt{r}}$ where $\eta_{\sqrt{r}}$ is the uniform distribution inside the $D$-dimensional unit ball with radius $\sqrt{r}$. Second is to let $\mu^{2}_{r}(d\theta^{\prime})$ follow $\mathcal{N}(\hat{\theta},r\Sigma)$. The second approach in particular can be intuited as the posterior distribution of the parameter from a Bayesian perspective. In both cases, we notice that letting $r=0$ would recover the original baseline distribution $p(y;\hat{\theta})$.

To analyze these schemes, we abuse notation slightly and now define $L:\mathbb{R}^{+}\times\Theta\rightarrow\mathbb{R}$ to be

[TABLE]

and

[TABLE]

We show that increasing $r$ to positive values would produce a descent direction for $l(r)$ at $r=0$, when the underlying distribution is Gaussian. Recall that in this case $\Theta^{*}(\delta_{\hat{\theta}})$ can be expressed by (4.17). As $l(r)$ is not necessarily convex in this situation, we use a generalized version of Danskin’s Theorem (Clarke, 1975) for non-convex problems to get the following result:

Theorem 4.4.1.

With $l(r)$ and $L(r,\theta)$ defined in (4.21) and (4.22), and $p(y;\theta)$ multivariate Gaussian with mean $\theta$ and known positive definite covariance $\Sigma$, we have

[TABLE]

The proof is in Appendix C. With Theorem 4.4.1, we can show that both $\mu_{r}^{1}$ and $\mu_{r}^{2}$ proposed above are valid choices to produce descent directions. Moreover, we can also show that they allow tractable computation of $\mathcal{D}_{data}(\mathbb{P}_{0})$. These are depicted as follows.

Corollary 4.4.2.

Under the assumptions in Theorem 4.4.1, $l^{+}(0)<0$ for both $\mu_{r}^{1}$ and $\mu_{r}^{2}$.

Corollary 4.4.3.

*Given $r\geq 0$ and $\mu_{prop}$ being $\mu_{r}^{1}(d\theta)$ or $\mu_{r}^{2}(d\theta)$, the nominal measure $\mathbb{P}_{r}$ with density given by (4.20) satisfies the conditions in Proposition 4.3.3. *

The proofs of Corollaries 4.4.2 and 4.3.4 are in Appendix C.

4.5. Numerical Demonstrations

To confirm our findings in Section 4.3 and 4.4, we perform several numerical experiments. Consider $\mathbb{P}_{\theta}$ to be multivariate Gaussian $\mathcal{N}(\theta,I_{D})$ with $k=D=10$. We set $\epsilon=\alpha=0.05$ while $\beta=0.01$ and data size $n=10$ or 5. Notice in this case, the dimension $D$ is high but the available sample $n$ is low and we would actually need $N_{exact}=371$ data points to perform standard SO. Based on our discussion, we compare three choices of $\mu_{prop}$:

• •

$\mu_{1}=\delta_{\hat{\theta}}$, the point mass at $\hat{\theta}$.

• •

$\mu_{2}\sim\hat{\theta}+\sqrt{\frac{\chi^{2}_{1-\alpha,D}}{n}}\cdot\eta$, where $\eta$ is the uniform random vector on the surface of a $D$-dimension unit ball, discussed in Section 4.3.

• •

$\mu_{3}\sim\mathcal{N}(\hat{\theta},I_{D}/n)$, the Gaussian distribution with mean $\hat{\theta}$ and covariance matrix $I_{D}/n$, discussed in Section 4.4.

For $\mu_{1}$, $\mu_{2}$ and $\mu_{3}$, the calculation of $\mathcal{D}_{data}(\mathbb{P}_{0})$ is tractable. We leave the details in the Appendix as remarks following Lemma C.0.5 and summarize the results in Table LABEL:10dim1 and LABEL:10dim2. We use $N$ to denote the number of Monte Carlo samples needed. Moreover, we use both algorithms Extended SO and Extended FAST discussed in Section 5 for demonstration. As we can see, the decrease in $N$ under a better sampling distribution can be considerable, down to less than a third compared to using the baseline in some cases. Mixing with a proposed uniform distribution ($\mu_{2}$) appears to reduce $N$ more than applying a Gaussian mixture ($\mu_{3}$). As a side note, we also observe Extended FAST requires significantly less sample size than Extended SO in this example.

5. Procedural Description

This section presents our procedures to find solutions for CCP (1.1) using SO-based methods, when the direct use of data $\xi_{1},...,\xi_{n}$ from $\mathbb{P}$ is possibly insufficient to achieve feasibility with a given confidence. Algorithm 1, which we call “Extended SO”, first presents the basic and most easily applicable procedure arising from Corollary 3.2.2. Notice that, given an overall target confidence level, say $c$, we have flexibility in choosing $\alpha$ and $\beta$ such that $\alpha+\beta=c$. In our experiments, we simply choose $\alpha=\beta=\frac{c}{2}$. However, if the required confidence level is high, it is more beneficial to choose a relatively small $\beta$, since the required Monte Carlo sample size depends only logarithmically on $\beta$ (i.e., the required sample size for SO is of order $\log\frac{1}{\beta}$) (Campi and Garatti, 2008). On the other hand, as the confidence level $1-\alpha$ grows higher, the size of uncertainty set $\mathcal{U}_{data}$ would grow and cause the tolerance level $\epsilon$ for the SO (under the baseline $\mathbb{P}_{0}$) to decrease. Here, the dependence of Monte Carlo sample size on $\epsilon$ is less favorable, typically of order $\frac{1}{\epsilon}$ (Campi and Garatti, 2008).

There are several variants of Algorithm 1. First, we have discussed the use of plain SO and that the required sample size is (2.5), while on the other hand, as mentioned at the end of Section 2.2, we can use other variants of SO such as FAST that requires a smaller sample size for either the data or the Monte Carlo samples we generate. In the case of FAST, we would have $N_{exact}(\epsilon,\beta,d)=20d+\frac{1}{\epsilon}\log\frac{1}{\beta}$, as suggested by (Carè et al., 2014). Thus, a variant of Algorithm 1 is to replace $N_{exact}$ with this latter quantity, and replace (2.4) with the FAST procedure in (Carè et al., 2014) for the last step of Algorithm 1 (we call this algorithm “Extended FAST” which will also be used in the next section).

The explicit expression for $\epsilon^{\prime}(\epsilon,\lambda,\phi)$ for different $\phi,\epsilon$ and $\lambda$ can be found in (Jiang and Guan, 2016). For example, if we choose $\phi=(x-1)^{2}$ which corresponds to the $\chi^{2}$-distance, then for $\epsilon<1/2$, we have $\epsilon^{\prime}=\max\{0,\epsilon-\frac{\sqrt{\lambda^{2}+4\lambda(\epsilon-\epsilon^{2})}-(1-2\epsilon)\lambda}{2\lambda+2}\}$. We can also replace $\epsilon^{\prime}(\epsilon,\lambda,\phi)$ by any $\delta_{\epsilon}$ that achieves $M(\mathbb{P}_{\hat{\theta}},\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda\},\delta_{\epsilon})\leq\epsilon$. In Appendix B, we derive a self-contained easy upper bound for $M(\mathbb{P}_{\hat{\theta}},\{\mathbb{Q}:d_{\phi}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda\},\delta)$ in the case of $\chi^{2}$-distance and use it to find such a $\delta_{\epsilon}$. This easy computation of $\delta_{\epsilon}$ will also be used in our numerics in the next section.

Section 4.1 has investigated some proposals to improve the generating distributions. Algorithm 2 depicts these proposals in a general form. The main difference of Algorithm 2 compared to Algorithm 1 is the introduction of $\mathcal{D}_{data}(\mathbb{P}_{0},\phi)$ that one can attempt to minimize over a class of generating distribution $\mathbb{P}_{0}$ or evaluate for trial-and-error choices of $\mathbb{P}_{0}$, so that at the end we have $\mathcal{D}_{data}(\mathbb{P}_{0},\phi)<\phi^{\prime\prime}(1)\chi^{2}_{1-\alpha,D}/(2n)$. As discussed in Section 4.1, using this $\mathbb{P}_{0}$ allows us to obtain a smaller Monte Carlo size requirement than simple relaxation of the parametric constraint. Sections 4.3 and 4.4 describe the possibilities of achieving such a reduction, in the case of Gaussian underlying distributions and using $\chi^{2}$-distance. Note that, just like in Algorithm 1, we can consider other variants such as incorporating FAST and using alternate bounds for $M$ instead of $\epsilon^{\prime}$, by undertaking the same modifications as in Algorithm 1.

6. Numerical Experiments

This section presents some numerical examples to support our theoretical findings and illustrate the performance of our proposed procedures for data-driven CCPs. We focus on Algorithm 1 (Extended SO) and its FAST variant discussed in Section 5 (Extended FAST). We consider both single and joint CCPs (i.e., one and multiple inequalities respectively in the safety condition of the probability) as well as quadratic optimization problems. Moreover, we compare numerically with methods of robust optimization (RO) in (Marandi et al., 2017; Ben-Tal et al., 2009). The experimental outputs that we report include:

• •

Under each setting, we repeat the experiment 1000 times with new data generated each time. For the solution $\hat{x}$ obtained in each trial from a given algorithm, we evaluate the violation probability $V(\hat{x},\mathbb{P})$ under the true probability measure $\mathbb{P}$ (under $\theta_{true}$) either through exact calculation or Monte Carlo simulation with sample size 10000. Moreover, using the empirical distribution for the violation probabilities, we report $\hat{\epsilon}$ as the average violation probability $V(\hat{x},\mathbb{P})$ as well as $Q_{95}$, the 95-percentile. Finally, we report and compare “$f_{val}$”, the average objective value for the optimization problem across all 1000 runs.

• •

We fix $\alpha=0.05$ and $\beta=0.01$ across different values of $\epsilon$ and $d$. However, when we compare our methods with robust optimization approaches, we set $\alpha=0.05$ and $\beta=0.001$, since RO approaches essentially guarantee $\beta=0$. On the other hand, the sample size chosen for FAST is taken with default values $N_{1}=20d$ in stage 1 and $N_{2}=\frac{\log\beta-\log(B_{\epsilon}^{N_{1},d})}{\log(1-\epsilon)}$ in stage 2 as discussed in (Carè et al., 2014).

• •

For given $\epsilon$ and $d$, we denote $N_{exact}$ as the required sample size if we can directly sample from $\mathbb{P}$ and use standard SO. We denote $n$ as the available data size ($n<N_{exact}$) and $N$ as the Monte Carlo size needed for the our DRO-based methods. In DRO-based methods, we fix our generating distribution $\mathbb{P}_{0}$ as $\mathbb{P}_{\hat{\theta}}$ and use the $\chi^{2}$-distance across the experiments.

6.1. Single Linear Chance Constraint Problem

We first consider a single linear CCP

[TABLE]

where $x\in\mathbb{R}^{d}$ is the decision variable, $a,c\in\mathbb{R}^{d}$ and $b\in\mathbb{R}$ are fixed and $\xi\in\mathbb{R}^{d}$ is a random vector following some parametric distribution. We fix $a=[5,5,...,5]\in\mathbb{R}^{d}$, $b=5$ and $c=[-1,-1,...,-1]\in\mathbb{R}^{d}$ and the problem would have a non-empty feasible region with high probability for $\xi$ considered here. Moreover, a robustly feasible point for FAST (Carè et al., 2014) is chosen to be $\bar{x}=\mathbb{0}\in\mathbb{R}^{d}$ and an explicit $\mathcal{U}_{data}$ is constructed as (2.21) for our DRO.

6.1.1. Multivariate Gaussian

We conduct experiments when $\xi\sim\mathcal{N}(\theta,\Sigma)$ with fixed but a priori randomly generated positive definite covariance matrix $\Sigma\in\mathbb{R}^{d\times d}$ and unknown $\theta\in\mathbb{R}^{d}$. Due to the normality of $\xi$, for any given $\theta$, we can reformulate the chance constraint exactly as a second-order cone constraint, which can be robustified straightforwardly in the ambiguous chance constraint case. The underlying true parameter is taken to be $\theta_{true}=\mathbb{0}\in\mathbb{R}^{d}$ and the results are summarized in Table LABEL:single_results_gau and LABEL:single_compare.

6.1.2. Exponential Distribution

We conduct experiment when each coordinate $\xi_{i}$ of $\xi\in\mathbb{R}^{d}$ independently follows exponential distribution with rate $\lambda_{i}$. Since $\xi$ is no longer Gaussian and the domain of the moment generating moment function for exponential distribution depends on $\lambda=(\lambda_{1},\ldots,\lambda_{d})$, for convenience we use RO constructed from a convex approximation using Chebyshev’s inequality:

[TABLE]

which, combined with $\mathcal{U}_{data}$ as in (2.21), reduces the ambiguous chance constraint into a robust conic quadratic constraint

[TABLE]

The above can be tractably reformulated as in Section 5 of (Marandi et al., 2017) on problems in the form of 5(b), with $\Omega=(\min_{i}(\hat{\lambda}_{i})(1-\frac{\chi^{2}_{1-\alpha,d}}{n}))^{-1}$ where $\hat{\lambda}_{i}$ represents the MLE estimate of $\lambda_{i}$. Finally, the underlying true parameters are taken as $\lambda_{i}=1,\forall i$, and results are summarized in Table LABEL:single_compare_exp.

From the results of the experiments, we can see the vast majority of solutions produced by three methods satisfy statistical feasibility. In fact, all methods are conservative with respect to the violation probability $\epsilon$, although some are more conservative than the other. In particular, when $\xi$ is Gaussian, RO takes advantage of an exact formulation to produce less conservative solution with lower objective value (closer to the optimal value). This can be seen in Table 4, where $\hat{\epsilon}=0.018$ $f_{val}=-0.80$ for RO and $\hat{\epsilon}=0.0011$ $f_{val}=-0.72$ only for Extended SO. When $\xi$ is no longer Gaussian, RO appears to produce similar-quality solutions as Extended SO in terms of feasibility or optimality. For example in Table 5, we have $\hat{\epsilon}=0.0045$ $f_{val}=-0.6978$ for RO and $\hat{\epsilon}=0.0047$ $f_{val}=-0.6981$ for Extended SO. Note that while the validity of RO depends crucially on the applicability and accuracy of convex approximation, the validity of Extended SO or Extended FAST is not restricted by the distributions of $\xi$, and they also do not require intensive, case-specific analysis as RO. In general, we observe consistent performances of our methods in both experiments.

6.2. Joint Linear Chance Constraint Problem

Next, we consider a joint chance-constrained linear problem:

[TABLE]

where $x\in\mathbb{R}^{d}$ is the decision variable, $A\in\mathbb{R}^{m\times d},c\in\mathbb{R}^{d}$ and $b\in\mathbb{R}^{m}$ are fixed and $\Xi\in\mathbb{R}^{m\times d}$ is a random matrix following some parametric distribution. We set $c$, each row of $A$ and $b$ to be the same as in the single linear CCP. We treat $\Xi\in\mathbb{R}^{m\times d}$ as a matrix concatenated from a random vector $\xi\in\mathbb{R}^{md}\sim\mathcal{N}(\theta,\Sigma)$ with fixed but a priori randomly generated positive definite covariance matrix $\Sigma\in\mathbb{R}^{md\times md}$ and unknown $\theta\in\mathbb{R}^{md}$. To solve RO, we use Bonferroni’s inequality as in (Nemirovski and Shapiro, 2007) to first divide the violation probability $\epsilon$ uniformly across $m$ individual chance constraints and then follow the procedure as in single linear CCP. The results are summarized in Table LABEL:joint_compare_exp.

In this joint linear example, Extended FAST provides the least conservative solution in terms of the achieved tolerance level ($\hat{\epsilon}=0.0226$, which is closer to $0.05$, compared to $0.0003$ in RO and $0.0012$ in Extended SO), and Extended SO is the least conservative in terms of the objective value ($f_{val}=-0.6626$ compared to $-0.6448$ in RO and $-0.6466$ in Extended FAST). RO appears to be the most conservative in terms of both the achieved tolerance level and objective value. Note that this occurs even though the underlying randomness is Gaussian, which allows exact reformulation in the single chance constraint case. The conservative performance here is likely (and unsurprisingly) due to the crude Bonferroni’s correction. Note that other alternatives to using Bonferroni, if one considers tractable reformulation, would be to use moment-based DRO where tractability can be achieved (e.g., (Xie and Ahmed, 2018)). However, it is unclear if using moment-based DRO would be more or less conservative than using Bonferroni correction along with exact reformulation for the individualized constraints, which could comprise an interesting comparison for a future study. Nonetheless, our Extended SO/FAST, being purely sampled-based, avoids the additional conservativeness coming from the Bonferroni correction. However, we note that a large number of Monte Carlo samples are required due to the large size of $\mathcal{U}_{data}$ in this high-dimensional problem.

6.3. Non-Linear Chance Constrained Problems

In this section, we conduct numerical experiments for non-linear CCP. We consider two examples. First is a quadratic objective with joint linear chance constraints, and second is a linear objective with a quadratic chance constraint, similar as the QM problem considered in (Lejeune and Margot, 2016).

6.3.1. Quadratic Objective with Joint Linear Chance Constraint

We adopt the same setup (thus the robust feasibility condition remains the same) as in (6.2) except we modify the objective with a quadratic term

[TABLE]

for a fixed but a priori randomly generated positive definite matrix $H$. We use $\epsilon=0.05$. Results are summarized in Table LABEL:ccp_qp. As we can see, feasibility in terms of violation probability is satisfied by all methods, though RO suffers from higher conservativeness compared to Extended SO/FAST in terms of the objective value ($f_{val}=-0.48$ compared to $-0.5547$ and $-0.5476$ for Extended SO and FAST respectively). Like the previous example, this could be attributed to the Bonferroni correction used in the RO. Extended FAST gives the least conservative solution in terms of the tolerance level ($\hat{\epsilon}=0.0096$), using only one third of the samples compared to Extended SO (3888 vs 1384). On the other hand, Extended SO gives the least conservative solution in terms of the objective value ($f_{val}=-0.5547$).

6.3.2. Linear Objective with Quadratic Chance Constraint

We consider the following setup:

[TABLE]

We set $\Xi=\frac{1}{m}\sum_{i=1}^{m}\xi_{i}\xi_{i}^{T}$ and $\xi_{i}\in\mathbb{R}^{d}\sim\mathcal{N}(\theta,\Sigma)$ are i.i.d. with unknown $\theta$. We set $\theta_{true}=0\in\mathbb{R}^{d}$ and consequently $m\Xi$ follows a Wishart distribution $\mathcal{W}(\Sigma,m)$ with $m$ degrees of freedom and covariance matrix $\Sigma$ under $\mathbb{P}$. We use $\epsilon=0.05$. The RO formulation for this problem is not readily available while our sampling-based methods are still directly applicable. We thus focus on evaluating the performance of Extended FAST under different hyper-parameters. The results are summarized in Table LABEL:LPQCC. As we can see, the high dimensions of the problem do not affect the sample size requirement of Extended FAST dramatically, as it increases moderately form $N=154$ when $d=5$ to $N=334$ when $d=10$ and to $N=422$ when $d=15$. Moreover, the average optimal value $f_{val}$ is around $-0.85$ and feasibility is satisfied ($\hat{\epsilon}$ all within $0.05$), showing the consistent effectiveness of our method.

7. Conclusion

We consider data-driven chance constrained problems with limited data. In such situation, standard approaches in SO may not be able to generate statistically feasible solutions. We investigate an approach that uses divergence-based DRO to efficiently incorporate parametric information through a data-driven uncertainty set, and subsequently uses Monte Carlo sampling to generate enough samples to handle the distributionally robust chance constraint. In this way our framework translates the data size requirement in SO into a Monte Carlo requirement, the latter could be much more abundant thanks to cheap modern computational power.

To exploit the full capability of our framework, we have investigated the optimality of the generating distribution in drawing the Monte Carlo samples in the sense of minimizing its required sample size. We have shown that, while the optimal choice is the baseline distribution in the unambiguous and nonparametric DRO cases, this natural choice can be dominated by other distributions in the parametric DRO case. We proved this by connecting the Neyman-Pearson lemma in statistical hypothesis testing to DRO and SO, which comprises the first such results of its kind as far as we know. We then studied several ways to find better generating distributions by searching for mixtures that enhance distributional variability. Lastly, we showed some numerical results to demonstrate how our approach can give rise to feasible solutions in a wide range of settings where other methods such as RO cannot be utilized directly or give more conservative solutions.

Acknowledgements.

We gratefully acknowledge support from the National Science Foundation under grants CMMI-1542020, CAREER CMMI-1653339/1834710 and IIS-1849280.

Appendix A Regularity Conditions to Verify Assumption A1

We consider the following conditions:

(C1) $p(x,\theta_{1})=p(x,\theta_{2})$ for all $x$ implies $\theta_{1}=\theta_{2}$.

(C2) $\theta_{true}$ is an inner point of $\Theta\subseteq\mathbb{R}^{D}$.

(C3) The support of distribution $\{x:p(x,\theta)>0\}$ does not depend on $\theta$.

(C4) There exists a measurable function $L_{1}(x)$ such that $\mathbb{E}_{\theta_{true}}L_{1}^{2}<\infty$ and

[TABLE]

for all $\theta_{1},\theta_{2}$ in a neighborhood of $\theta_{true}$.

(C5) $I(\theta_{true})$ is non-singular.

(C6) The density family $\{{\mathbb{P}}_{\theta}\}_{\theta\in\Theta}$ is differentiable in quadratic mean at $\theta_{true}$, i.e., there exists a measurable function $L_{2}(x):\mathcal{X}\rightarrow\mathbb{R}^{D}$ such that for any $h\in\mathbb{R}^{D}$ that converges to 0,

[TABLE]

The consistency and asymptotic normality of MLE in Assumption A1 is guaranteed under conditions (C1)-(C6). See (Van der Vaart, 2000; Lehmann, 2004).

Appendix B Alternate Bounds Using $\chi^{2}$ Distance

Consider the $\chi^{2}$-based uncertainty set $\mathcal{U}_{data}=\{\mathbb{Q}\in\mathcal{P}_{para}:\chi^{2}(\mathbb{P}_{\hat{\theta}},\mathbb{Q})\leq\lambda\}$. Here we provide an alternate upper bound for the function $M(\mathbb{P}_{0},\mathcal{U}_{data},\delta)$, which we call $\tilde{M}(\mathbb{P}_{0},\mathcal{U}_{data},\delta)$. That is, we find $\tilde{M}(\mathbb{P}_{0},\mathcal{U}_{data},\delta)$ that satisfies

[TABLE]

For any $\mathbb{Q}$ absolutely continuous with respect to $\mathbb{P}_{0}$, we have

[TABLE]

where the fourth line follows from the Cauchy-Schwarz inequality. Thus, we can set

[TABLE]

which is non-decreasing in $\delta$. By (2.9), we can choose $\delta_{\epsilon}$ such that $\delta_{\epsilon}+\delta_{\epsilon}^{1/2}\left(\mathcal{D}_{data}(\mathbb{P}_{0})\right)^{1/2}\leq\epsilon$, or equivalently,

[TABLE]

by solving the quadratic equation. In the case where we relax the parametric constraint completely, we have $\mathcal{D}_{data}(\mathbb{P}_{0})=\lambda$. Compared to the bound obtained from Theorem 3.2.1 and Corollary 3.2.2, (B.2) is less tight, but the gap can be shown to asymptotically vanish when $\epsilon,\frac{\chi^{2}_{1-\alpha,D}}{n}\rightarrow 0$.

Appendix C Proofs and Other Technical Results

Proof of Lemma 4.2.1.

First, by definition $\mathcal{P}(\Theta)$ is a convex set and, for any $\mu_{1},\mu_{2}\in\mathcal{P}(\Theta)$ and $0<t<1$, we have

[TABLE]

Next, fixing $\theta\in\mathcal{U}_{data}$, the function $L(\cdot,\theta)$ is convex since:

[TABLE]

for any $0<t<1$ where the inequality follows from the convexity of the function $1/x$ for $x>0$. Thus, as the supremum of convex functions, $l(\mu)\triangleq\sup\limits_{\theta\in\mathcal{U}_{data}}L(\mu,\theta)$ is also convex. ∎

We provide a version of Danskin’ Theorem needed to prove Theorem 4.3.1. Alternately, one can also resort to a generalized version in (Clarke, 1975) by verifying the conditions there. Here we opt for the former and provide a self-contained proof, which mostly relies on the techniques from Proposition 4.5.1 of (Bertsekas et al., 2003) but with some slight modification to handle issues regarding the domain of the involved function. We have:

Lemma C.0.1.

Fix probability measures $\mu_{1},\mu_{2}\in\mathcal{P}(\Theta)$. Suppose $t_{k}\downarrow 0$ is a positive sequence such that $(1-t_{k})\mu_{1}+t_{k}\mu_{2}\in\mathcal{P}(\Theta)$ for all $k$ and $\theta_{k}\in\Theta^{*}((1-t_{k})\mu_{1}+t_{k}\mu_{2})$ is a sequence such that $\theta_{k}\rightarrow\theta_{0}$ for some $\theta_{0}\in\mathcal{U}_{data}$, then we have

[TABLE]

if we assume $L((1-t)\mu_{1}+t\mu_{2},\theta)$ is jointly continuous in $0\leq t\leq 1$ and $\theta\in\Theta$.

Proof.

It is known that if $f:\mathbb{I}\rightarrow\mathbb{R}$ is a convex function with $\mathbb{I}$ being an open interval containing some point $x$, we then have the following results (Bertsekas et al., 2003):

[TABLE]

[TABLE]

and

[TABLE]

In other words, these limits exist and satisfy the above relations for convex functions. Thus, if we define $f_{k}(t)={L}((1-t_{k})\mu_{1}+t_{k}\mu_{2}+t(\mu_{2}-\mu_{1}),\theta_{k})$, it follows from the convexity of $\mathcal{P}(\Theta)$ and ${L}(\cdot,\theta_{k})$ that $f_{k}(t)$ is convex and well-defined for $-t_{k}\leq t\leq 1-t_{k}$. Using the above results in (C.1), (C.2) and (C.3), we then have

[TABLE]

On the other hand, if we define $f_{0}(t)=L((1-t)\mu_{1}+t\mu_{2},\theta_{0})$, it also follows that $f_{0}(t)$ is convex and well-defined for $0\leq t\leq 1$. It follows from the convexity of $f_{0}(\cdot)$ as well as (C.1) that

[TABLE]

Then, it again follows from the convexity of $f_{0}(\cdot)$ that, given any $\tau>0$, we can find some $\eta>0$ such that

[TABLE]

for all $0<s<\eta$. It then follows from definitions and (C.6) that

[TABLE]

for all $0<s<\eta$. Fixing one such $s$, since the function $L((1-t)\mu_{1}+t\mu_{2},\theta)$ is jointly continuous in $0\leq t\leq 1$ and $\theta\in\Theta$, and the sequence satisfies $\theta_{k}\rightarrow\theta_{0}$, we have

[TABLE]

as long as we make $\eta>s>0$ small enough so that $\eta\leq 1-t_{k}$ for all $k$. Then, for $k$ large enough, we have

[TABLE]

Combining (C), (C) and (C.8), we have that, for $k$ large enough,

[TABLE]

Finally, since $\tau$ is arbitrary, we conclude that

[TABLE]

∎

We now prove the following version of Danskins’ Theorem:

Theorem C.0.2.

Fix $\mu_{1},\mu_{2}\in\mathcal{P}(\mathcal{U}_{data})$. Suppose $t_{k}\downarrow 0$ is a positive sequence such that $(1-t_{k})\mu_{1}+t_{k}\mu_{2}\in\mathcal{P}(\Theta)$ for all $k$ and $L((1-t)\mu_{1}+t\mu_{2},\theta)$ is jointly continuous in $0\leq t\leq 1$ and $\theta\in\Theta$. Then, if we let ${\psi}(t)={l}((1-t)\mu_{1}+t\mu_{2})$ for $0\leq t\leq 1$ with ${l}(\cdot)=\sup\limits_{\theta\in\mathcal{U}_{data}}L(\cdot,\theta)$ defined as (4.9), we have

[TABLE]

Proof.

For any $\theta_{0}\in\Theta^{*}(\mu_{1})$ and $\theta_{t}\in\Theta^{*}((1-t)\mu_{1}+t\mu_{2})$, we have

[TABLE]

Thus, by taking $t\downarrow 0$ and taking the supremum over all $\theta_{0}\in\Theta^{*}(\mu_{1})$, we have

[TABLE]

Notice that the existence of the several limits above follows from the convexity of related functions. To prove the reverse inequality, we consider a sequence $\{t_{k}\}$ with $0<t_{k}<1$ and $t_{k}\downarrow 0$. Then, we pick another sequence $\{\theta_{k}\}\subseteq\mathcal{U}_{data}$ with $\theta_{k}\in\Theta^{*}((1-t_{k})\mu_{1}+t_{k}\mu_{2})$ for all $k$. Since $\mathcal{U}_{data}$ is compact, there exist a subsequence of $\{\theta_{k}\}$ converge to some $\theta_{0}\in\mathcal{U}_{data}$. Without loss of generality, we drop the subsequence and simply assume $\theta_{k}\rightarrow\theta_{0}$. We first show $\theta_{0}\in\Theta^{*}(\mu_{1})$. To do this, pick any $\tilde{\theta}_{0}\in\Theta^{*}(\mu_{1})$. Since ${L}((1-t)\mu_{1}+t\mu_{2},\theta)$ is jointly continous in $t$ and $\theta$, we have

[TABLE]

where the inequality follows from the definition of $\theta_{k}$. Now, since $\tilde{\theta}_{0}\in\Theta^{*}(\mu_{1})$ and ${L}(\mu_{1},\theta_{0})\geq{L}(\mu_{1},\tilde{\theta}_{0})$, we must have

[TABLE]

Now, using the definition of $\Theta^{*}(\mu_{1})$, we can write

[TABLE]

Now we use Lemma C.0.1 to conclude that

[TABLE]

Finally, we combine (C.10) and (C) to conclude the proof

[TABLE]

∎

Proof of Theorem 4.3.1.

The result can be obtained from Leibniz’s integral rule (i.e. differentiation under the integral sign). See, for example, Theorem 2.27 in (Folland, 2013). ∎

Next we prove Proposition 4.3.2. For convenience, we note that (4.15) can be written in a compact form for exponential family (Nielsen and Nock, 2014):

[TABLE]

where $\langle a,b\rangle=a^{\intercal}b$ represents the usual inner product in the Euclidean space, and $t(\cdot),F(\cdot)$ and $k(\cdot)$ are known functions. In particular, we have

[TABLE]

To facilitate the calculation, we first introduce two lemmas involving the exponential parametric family based on (Nielsen and Nock, 2014).

Lemma C.0.3.

Pick $\theta_{1},\theta_{2}\in\Theta$. If $2\theta_{2}-\theta_{1}\in\Theta$, then we have

[TABLE]

In particular, if $F(\theta)=\frac{\theta^{\intercal}\Sigma^{-1}\theta}{2}$, then $\int_{\mathcal{Y}}\frac{(p(y;\theta_{2}))^{2}}{p(y;\theta_{1})}dy=e^{(\theta_{2}-\theta_{1})^{\intercal}\Sigma^{-1}(\theta_{2}-\theta_{1})}$.

Proof.

It follows from (C.13) that

[TABLE]

∎

Lemma C.0.4.

Pick $\theta_{1}$,$\theta_{2}$ and $\theta_{3}\in\Theta$. If $2\theta_{2}-2{\theta}_{1}+\theta_{3}\in\Theta$, then we have

[TABLE]

In particular, if $F(\theta)=\frac{\theta^{\intercal}\Sigma^{-1}\theta}{2}$, then $\int_{\mathcal{Y}}\frac{(p(y;\theta_{2}))^{2}p(y;\theta_{3})}{(p(y;\theta_{1}))^{2}}dy=e^{(\theta_{2}-\theta_{1})^{\intercal}\Sigma^{-1}(\theta_{2}-\theta_{1})+2(\theta_{2}-\theta_{1})\Sigma^{-1}(\theta_{3}-\theta_{1})}$.

Proof.

The proof follows from the same techniques as in Lemma C.0.3. ∎

Then (4.3) follows from (2.20), (C.14) and Lemma C.0.3 so that

[TABLE]

and (4.17) follows. We now prove Proposition 4.3.2:

Proof of Proposition 4.3.2.

Following Theorem 4.3.1, Lemma C.0.3 and Lemma C.0.4, we have

[TABLE]

Notice the second equality follows from Fubini’s theorem. The third equality follows from Lemma C.0.3 and Lemma C.0.4. The fourth equality follows from (4.17). Now, following the last line (C), for the search of descent direction, it is sufficient to prove

[TABLE]

However, since $\mu_{prop}(d\theta^{\prime})$ is a symmetrical distribution around $\hat{\theta}$, we know that

[TABLE]

for any $\theta\in\Theta^{*}(\delta_{\hat{\theta}})$. Then, it follows from Jensen’s inequality that

[TABLE]

Now suppose for the sake of contradiction that

[TABLE]

Then, let $\{\theta_{k}\}_{k}\subseteq\Theta^{*}(\delta_{\hat{\theta}})$ be a subsequence such that $\mathbb{E}_{\theta^{\prime}\sim\mu_{prop}}[e^{2(\theta_{k}-\hat{\theta})^{\intercal}\Sigma^{-1}(\theta^{\prime}-\hat{\theta})}]\rightarrow 1$. Due to the compactness of $\Theta^{*}(\delta_{\hat{\theta}})$, we can find a subsequence of $\{\theta_{k}\}_{k}$ converging to some $\theta_{0}\in\Theta^{*}(\delta_{\hat{\theta}})$. For convenience we drop the subsequence and suppose $\theta_{k}\rightarrow\theta_{0}.$ Then the existence of $Y$ allows us to use dominated convergence theorem:

[TABLE]

However, Jensen’s inequality would indicate that $\mathbb{E}[e^{2Y_{\theta_{0}}}]=1$ if and only $\mathbb{P}(Y_{\theta_{0}}=0)=1$, which contradicts our assumption. Thus, we know that

[TABLE]

as claimed. ∎

Proof of Proposition 4.3.3.

First we prove (4.18). Letting $c=\frac{1}{(2\pi)^{D}|\Sigma|}$, we know that

[TABLE]

where we denote $p_{Z}(\cdot)$ to be the density function of random variable $Z=\Sigma^{-1/2}(Y-\hat{\theta})$ with $Y\sim\mathbb{P}_{0}$ and the last two lines follow from a change of variable $z=\Sigma^{-1/2}(y-\hat{\theta})$. Now, since $\|\Sigma^{-1/2}(\theta_{1}-\hat{\theta})\|_{2}^{2}=\|\Sigma^{-1/2}(\theta_{2}-\hat{\theta})\|_{2}^{2}=r$ for some $r$ by assumption, it follows from (C) that $\chi^{2}(\mathbb{P}_{0},\mathbb{P}_{\theta_{1}})=\chi^{2}(\mathbb{P}_{0},\mathbb{P}_{\theta_{2}})$ if we can show

[TABLE]

However, since $p_{Z}(z)$ and $e^{-z^{T}z}$ are both rotationally invariant functions (i.e. $f(z)=f(Q^{\intercal}z)$ for all $z$ and rotational matrix $Q$, with $|Q|=1$), it can be shown that $\int\frac{e^{-z^{T}z}\cdot e^{-2z^{T}\nu}}{p_{Z}(z)}dz$ holds the same value for any $\nu$ such that $\|\nu\|_{2}^{2}=r$. Notice the rotational invariance of $p_{Z}(z)$ follows from the rotational invariance of $Z$. This proves (4.18). To prove (4.19), notice that for any $\theta\in\mathcal{U}_{data}$, we can find some $0\leq t\leq 1$ such that

[TABLE]

and hence $\chi^{2}(\mathbb{P}_{0},\mathbb{P}_{((1-t)\hat{\theta}+t\theta^{*})})=\chi^{2}(\mathbb{P}_{0},\mathbb{P}_{\theta})$ by (4.18). ∎

Proof of Proposition 4.3.4.

To check that $Y\sim\mathbb{P}_{t}$ with density

[TABLE]

leads to rotationally invariant $Z=\Sigma^{-1/2}(Y-\hat{\theta})$, simply notice that

[TABLE]

where $U_{t}$ is an independent Bernoulli variable with success rate $t$, $\eta$ is a random vector uniformly distributed on the surface of the $D$-dimensional unit ball and $X_{1},X_{2}$ are independent $\mathcal{N}(0,\Sigma)$. Then, it follows that

[TABLE]

where $Z_{1},Z_{2}$ are now independent $\mathcal{N}(0,I_{D})$. Consequently, the rotational invariance of $Z$ now follows from the rotational invariance of $Z_{1},Z_{2},\eta$ and their independence. ∎

Following the comments after Proposition 4.3.2, we show that $\theta\sim\mu_{prop}$ with $\theta\overset{\mathcal{D}}{=}\hat{\theta}+\sqrt{\frac{\chi^{2}_{1-\alpha,D}}{n}}\cdot\Sigma^{1/2}\cdot\eta$ provides a descent direction, with an alternate proof using the following lemma and the last line of (C).

Lemma C.0.5.

Fixing $\theta_{1}\in\Theta^{*}(\delta_{\hat{\theta}})$, we have

[TABLE]

for $\theta_{2}\sim\mu_{prop}(d\theta)$ where $\theta_{2}\overset{\mathcal{D}}{=}\hat{\theta}+\sqrt{\frac{\chi^{2}_{1-\alpha,D}}{n}}\cdot\Sigma^{1/2}\cdot\eta$ with $\eta$ following the uniform distribution on the surface of the $D$-dimensional unit ball.

Proof of Lemma C.0.5.

Let $u_{1}\in\mathbb{R}^{D}$ denote an arbitary point on the surface of $D$ dimensional unit ball ( $\|u_{1}\|_{2}^{2}=1$) and let $\eta=[\eta_{1},\eta_{2},...,\eta_{D}]$ be the random vector in $\mathbb{R}^{D}$ uniformly distributed on the surface of $D$ dimensional unit ball. Then we claim that $\frac{\mu_{1}^{T}\eta+1}{2}\sim Beta(\frac{D-1}{2},\frac{D-1}{2})$.

To show this, assume without loss of generality that $u_{1}=[1,0,...,0]\in\mathbb{R}^{D}$. Then for any $t\in[-1,1]$, it follows that $\mathbb{P}(u_{1}^{T}\eta\in dt)$ is proportional to the infinitesimal surface area on the ball corresponding to $\eta_{1}\in dt$, which is in turn proportional to the product of the sub-dimension $D-2$ surface area on the belt $x_{2}^{2}+x_{3}^{2}+...+x_{D}^{2}=1-t^{2}$ with the infinitesimal width of this belt. Specifically, the sub-dimension $D-2$ surface area around the belt is proportional to $(\sqrt{1-t^{2}})^{D-2}$. This follows from the fact that points of the form $[0,\sqrt{1-t^{2}},0,...,0],[0,0,\sqrt{1-t^{2}},0,...,0],...,[0,0,...,0,\sqrt{1-t^{2}}]$ are on this belt. Also, the width of this belt, according to the Pythagorean theorem, is $dt\cdot\sqrt{(\frac{d\sqrt{1-t^{2}}}{dt})^{2}+1}=\frac{dt}{\sqrt{1-t^{2}}}$. Thus,

[TABLE]

Now, we can substitute $\frac{t+1}{2}=s$ with $s\in[0,1]$ to get

[TABLE]

which can only be the density function for $Beta(\frac{D-1}{2},\frac{D-1}{2})$. It now follows from (Winkelbauer, 2012) that $\frac{u_{1}^{T}\eta+1}{2}$ has moment generating function

[TABLE]

for some $c>0$ where ${}_{1}F_{1}(\cdot,\cdot,\cdot)$ and ${}_{0}F_{1}(;,\cdot,\cdot)$ are the confluent hypergeometric function with identity ${}_{1}F_{1}(a,2a,x)=e^{x/2}_{0}F_{1}(;a+1/2,x^{2}/16)$ (see (Daalhuis, 2010)),

[TABLE]

with $(\gamma)_{k}=\frac{\Gamma(\gamma+k)}{\Gamma(\gamma)}$ being the Pochhammer symbol (Winkelbauer, 2012). To conclude the proof, denote $\rho_{n}=\sqrt{\log(1+\frac{\chi^{2}_{1-\alpha,D}}{n})}\cdot\sqrt{\frac{\chi^{2}_{1-\alpha,D}}{n}}$ and use (4.3), (4.17) and (C) to write

[TABLE]

∎

Remark*.*

Following Lemma C.0.5, we discuss the numerical calculations of $\mathcal{D}(\mathbb{P}_{0})$ following Proposition 4.3.4. We use $\mathcal{U}_{data}=\{\mathbb{P}_{\theta}:\|\theta-\hat{\theta}\|_{2}^{2}\leq\frac{\chi^{2}_{1-\alpha,D}}{n}\}$ where $p(y;\theta)=(2\pi)^{-\frac{D}{2}}e^{-\frac{1}{2}\|(y-\theta)\|_{2}^{2}}$. Then, for $\mu_{1}$, the nominal $p_{0}(y)$ is simply $p(y;\hat{\theta})$ and $\mathcal{D}_{data}(\mathbb{P}_{0})=\mathcal{D}_{data}(\mathbb{P}_{\hat{\theta}})=\max_{\theta\in\mathcal{U}}e^{\|\theta-\hat{\theta}\|_{2}^{2}}-1=e^{\frac{\chi^{2}_{1-\alpha,D}}{n}}-1$ according to (4.3) and Lemma C.0.3. For $\mu_{2}$, it can be shown that the nominal $\mathbb{P}_{0}$ follows $\mathcal{N}(\hat{\theta},(1+\frac{1}{n})\cdot I_{D})$, and a direct computation would show that $\mathcal{D}_{data}(\mathbb{P}_{0})=\max_{\theta\in\mathcal{U}}(\frac{(n+1)^{2}}{n(n+2)})^{\frac{d}{2}}e^{\frac{n}{n+2}\|\theta-\hat{\theta}\|_{2}^{2}}-1=(\frac{(n+1)^{2}}{n(n+2)})^{\frac{d}{2}}e^{\frac{n}{n+2}\frac{\chi^{2}_{1-\alpha,D}}{n}}-1$. Finally, for $\mu_{3}$, assume w.l.o.g that $\hat{\theta}=0$. Then we use the derivation in Lemma C.0.5 that $\frac{\mu_{1}^{T}\eta+1}{2}\sim Beta(\frac{D-1}{2},\frac{D-1}{2})$ for any $u_{1}$ on the $D$-dimensional unit ball surface to show that, for any $v\in\mathbb{R}^{D}$,

[TABLE]

and consequently

[TABLE]

Then, to calculate $\mathcal{D}_{data}(\mathbb{P}_{0})$, we note that

[TABLE]

Furthermore, through either direct verification or analysis similar to those in Lemma C.0.5, we note that $Y\sim\mathcal{N}(0,I_{D})$ shares the same distribution of $L\eta$ where $L\in\mathbb{R}^{+}$ and $\eta\in\mathbb{R}^{D}$ are two independent random variables with $L$ being the norm of $\mathcal{N}(0,I_{D})$ bearing density $f_{L}(l)=1_{\{l\geq 0\}}\frac{2^{1-\frac{D}{2}}}{\Gamma(\frac{D}{2})}l^{d-1}e^{-\frac{l^{2}}{2}}$ and $\eta$ being the random vector on the $D$-dimensional unit ball surface. Thus, it follows from (Remark) that (Remark) equals

[TABLE]

which is numerically tractable.

Proof of Theorem 4.4.1.

It follows from routine calculation that we can find a compact neighborhood of $r$ around 0 such that $\nabla_{r}L(r,\theta)$ exists and is continuous. Thus we can use the main theorem in (Clarke, 1975) to show that

[TABLE]

∎

To prove Corollary 4.4.2, we present two technical Lemmas C.0.6 and C.0.7.

Lemma C.0.6.

For any $\theta\in\Theta^{*}(\delta_{\hat{\theta}})$, \lim_{r\downarrow 0}\frac{1}{r}\Big{(}1-\mathbb{E}_{\theta^{\prime}\sim\mu^{1}_{r}}[e^{2(\theta-\hat{\theta})^{\intercal}\Sigma^{-1}(\theta^{\prime}-\hat{\theta})}]\Big{)} is a fixed negative value.

Proof of Lemma C.0.6.

For any $\theta\in\Theta^{*}(\delta_{\hat{\theta}})$, we have $\|\Sigma^{-1/2}(\theta-\hat{\theta})\|_{2}=\sqrt{\log(1+\frac{\chi^{2}_{1-\alpha,D}}{n})}$. Denote $\rho_{n}=\sqrt{\log(1+\frac{\chi^{2}_{1-\alpha,D}}{n})}$. Furthermore, under $\theta^{\prime}\sim\mu^{1}_{r}(d\theta^{\prime})$, we have $\Sigma^{-1/2}(\theta^{\prime}-\hat{\theta})\sim\eta_{\sqrt{r}}$, the uniform distribution inside the $D$-dimensional unit ball with radius $\sqrt{r}$, which can be viewed as the product of two independent random variables

[TABLE]

where $U$ is the uniform distribution on the surface of the $D$-dimensional unit ball and $R$ is the norm of the random vector ranged from [math] to $\sqrt{r}$. For any $0\leq s\leq\sqrt{r}$, since $\eta_{\sqrt{r}}$ follows a uniform distribution inside a $D$-dimensional unit ball, and the volume of a $D$-dimensional ball with radius $s$ is proportional to $s^{D}$, then $f_{R}(s)$, the density of $R$, must satisfy

[TABLE]

which is equivalent to saying

[TABLE]

Thus, we have that $\mathbb{E}[R^{2}]=c_{1}r$ for some $c_{1}>0$. Now we let $u_{1}=[1,0,...,0]\in\mathbb{R}^{D}$. We utilize the proof in Lemma C.0.5 as well as the independence of $R,U$ to show that

[TABLE]

Now it follows that

[TABLE]

∎

Lemma C.0.7.

For any $\theta\in\Theta^{*}(\delta_{\hat{\theta}})$, \lim_{r\downarrow 0}\frac{1}{r}\Big{(}1-\mathbb{E}_{\theta^{\prime}\sim\mu^{2}_{r}}[e^{2(\theta-\hat{\theta})^{\intercal}\Sigma^{-1}(\theta^{\prime}-\hat{\theta})}]\Big{)} is a fixed negative value.

Proof of Lemma C.0.7.

For any $\theta\in\Theta^{*}(\delta_{\hat{\theta}})$, we have $\|\Sigma^{-1/2}(\theta-\hat{\theta})\|_{2}=\sqrt{\log(1+\frac{\chi^{2}_{1-\alpha,D}}{n})}$. Denote $\rho_{n}=\sqrt{\log(1+\frac{\chi^{2}_{1-\alpha,D}}{n})}$. Furthermore, under $\theta^{\prime}\sim\mu^{2}_{r}(d\theta^{\prime})$, we have $\Sigma^{-1/2}(\theta^{\prime}-\hat{\theta})\sim\mathcal{N}(0,rI_{D})$. Using the moment generating function for Gaussian random variables, we have

[TABLE]

Now it follows that

[TABLE]

∎

Proof of Corollary 4.4.2.

Lemmas C.0.6 and C.0.7 combined with (4.23) indicate that increasing $r$ to positive value would produce a descent direction for $l(r)$ at $r=0$. ∎

Proof of Corollary 4.4.3.

We proceed the proof as in Proposition 4.3.4. The proof for the case of $\mu_{r}^{1}(d\theta)$ is entirely similar. For the proof of the case $\mu_{r}^{2}(d\theta)$, we simply notice that if $Y\sim\mathbb{P}_{t}$, then

[TABLE]

where $U_{t}$ is an independent Bernoulli variable with success rate $t$ and $X_{1},X_{2},X_{3}$ are independent $\mathcal{N}(0,\Sigma)$. Then, it follows that

[TABLE]

where $Z_{1},Z_{2},Z_{3}$ are now independent $\mathcal{N}(0,I_{D})$. Consequently, the rotational invariance of $Z$ now follows from the rotational invariance of $Z_{1},Z_{2},Z_{3}$ and their independence.

∎