Efficient Scenario Generation for Heavy-Tailed Chance Constrained Optimization

1. Introduction

In this paper, we consider the following family of chance constrained optimization problems:

where $x\in {\mathbb{R}}^{d_x}$ is a d_x-dimensional decision vector and L is a d_l-dimensional random vector in ${\mathbb{R}}^{d_l}$. The elements of L are often referred to as risk factors; the function $\varphi :{\mathbb{R}}^{d_x}\times {\mathbb{R}}^{d_l}\to \mathbb{R}$ is often assumed to be convex in x and often models a cost constraint; the parameter $\delta >0$ is the risk level of the tolerance. Our framework encompasses the joint chance constraint of the form $\text{P}({\varphi }_j(x,L)>0, \exists j\in \{1,\dots ,n\})\le \delta$, by setting $\varphi (x,L)={\text{max}}_{j=1,\dots ,n}{\varphi }_j(x,L)$.

Chance constrained optimization problems have a rich history in operations research. Introduced by Charnes et al. (1958), chance constrained optimization formulations have proved to be versatile in modeling and decision making in a wide range of settings. For example, Prekopa (1970) used these types of formulations in the context of production planning. The work of Bonami and Lejeune (2009) illustrates how to take advantage of chance constrained optimization formulations in the context of portfolio selection. In the context of power and energy control the use of chance constrained optimization is illustrated in Andrieu et al. (2010). These are just examples of the wide range of applications that have benefited (and continue to benefit) from chance constrained optimization formulations and tools.

Consequently, there has been a significant amount of research effort devoted to the solution of chance constrained optimization problems. Unfortunately, however, these types of problems are provably NP-hard in the worst case (Luedtke et al. 2010). As a consequence, much of the methodological effort has been placed into developing (a) solutions in the case of specific models; (b) convex and, more generally, tractable relaxations; (c) combinatorial optimization tools; and (d) Monte Carlo sampling schemes. Of course, hybrid approaches are also developed. For example, as a combination of type b and type d approaches, Hong et al. (2011) show that the solution to a chance constraint optimization problem can be approximated by optimization problems with constraints represented as the difference of two convex functions. In turn, this is further approximated by solving a sequence of convex optimization problems, each of which can be solved by a gradient based Monte Carlo method. Another example is Peña-Ordieres et al. (2020), which combines relaxations of type b with sample-average approximation associated with type d methods. In addition to the aforementioned types, Hong et al. (2021) provides an upper bound for the chance constraint optimization problem using a robust optimization with a data-driven uncertainty set, achieving a dimension independent sample complexity.

Examples of type a approaches include the study of Gaussian or elliptical distributions when ϕ is affine both in L and x. In this case, the problem admits a conic programming formulation, which can be efficiently solved (Lagoa et al. 2005). Type b approaches include Hillier (1967); Seppälä (1971); Ben-Tal and Nemirovski (2000, 2002); Prékopa (2003); Bertsimas and Sim (2004); Nemirovski and Shapiro (2006a); Chen et al. (2010); and Tong et al. (2022). These approaches usually integrate probabilistic inequalities such as Chebyshev’s bound, Bonferroni’s bound, Bernstein’s approximations, or large deviation principles to construct tractable analytical approximations. Type c methods are based on branch and bounding algorithms, which connect squarely with the class of tools studied in areas such as integer programming (Ahmed and Shapiro 2008, Luedtke et al. 2010, Küçükyavuz 2012, Luedtke 2014, Zhang et al. 2014, Lejeune and Margot 2016). Type d methods include the sample gradient method, the sample average approximation, and the scenario approach. The sample gradient method is usually combined with a smooth approximation (see Hong et al. (2011) for example). The sample average approximations studied by Luedtke and Ahmed (2008) and Barrera et al. (2016), although simplifying the constraint’s probabilistic structure via replacing the population distribution by sampled empirical distribution, are nevertheless hard to solve due to nonconvex feasible regions. The method we consider in this paper is the scenario approach. The scenario approach is introduced and studied in Calafiore and Campi (2005) and is further developed in a series of papers, including Calafiore and Campi (2006); Nemirovski and Shapiro (2006b).

The scenario approach is the most popular generic method for (approximately) solving chance constrained optimization. The idea is to sample a number N of scenarios (each scenario consists of a sample of L) and enforce the constraint in all of these scenarios. The intuition is that if for any scenario, say $L^{(i)}$, the constraint $\varphi (L^{(i)},x)<0$ is convex in x, and $\delta >0$ is small, we expect that by suitably choosing N the constrained regions can be relaxed by enforcing $\varphi (L^{(i)},x)<0$ for all $i=1,\dots ,N$, leading to a good and, in some sense, tractable (if N is of moderate size) approximation of the chance constrained region. Of course, this intuition is correct only when $\delta >0$ is small and we expect the choice of N to be largely influenced by this asymptotic regime.

By choosing N sufficiently large, the scenario approach allows obtaining both upper and lower bounds which become asymptotically tighter as $\delta \to 0$. In a celebrated paper, Calafiore and Campi (2006) provide rigorous support for this claim. In particular, given a confidence level $\beta \in (0,1)$, if $N\ge (2/\delta )\times \log (1/\beta )+2d+(2d/\delta )\times \log (2/\delta )$, with probability at least $1-\beta$, the optimal solution of the scenario approach relaxation is feasible for the original chance constrained problem and, therefore, an upper bound to the problem is obtained. Unfortunately, the required sample size of N grows with $(1/\delta )\times \log (1/\delta )$ as δ becomes small, limiting the scope of the scenario methods in applications.

Many applications of chance constraint optimization require a very small δ. For example, in the 5G ultra-reliable communication system design, the failure probability δ is no larger than ${10}^{-5}$ (Alsenwi et al. 2019); for fixed income portfolio optimization, an investment grade portfolio has a historical default rate of ${10}^{-4}$, reported by Frank (2008).

Motivated by this, Nemirovski and Shapiro (2006b) developed a method that lowers the required sample size to the order of $\log (1/\delta )$, making additional assumptions on the function $\varphi$ (which is taken to be biaffine), and the risk factors L, which are to be assumed light-tailed. Specifically, the moment generating function $E[\exp (sL)]$ is assumed to be finite in a neighborhood of the origin. No guarantee is given in terms of how far the upper bound is from the optimal value function of the problem as $\delta \to 0$.

In the present paper, we focus on improving the scalability of N in terms of $1/\delta$ for the practically important case of heavy-tailed risk factors. Heavy-tailed distributions appear in a wide range of applications in science, engineering, and business (Wierman and Zwart 2012, Embrechts et al. 2013) but, in some aspects, are not as well understood as light-tails. One reason is that techniques from convex duality cannot be applied as the moment generating function of L does not exist in a neighborhood of zero. In addition, probabilistic inequalities, exploited in Nemirovski and Shapiro (2006b), do not hold in this setting. Only very recently, a versatile algorithm for heavy-tailed rare event simulation has been developed in Chen et al. (2019).

The main contribution of our paper is to develop an algorithm that has a sample complexity N uniformly bounded in the risk tolerance parameter, assuming a versatile class of heavy-tailed distributions for L. Specifically, we shall assume that L follows a semiparametric class of models known as multivariate regular variation, which is quite standard in multivariate heavy-tail modeling (Embrechts et al. 2013, Resnick 2013). Moreover, our estimator is shown to be within a constant factor to the solution to (${\text{CCP}}_{\delta }$) with high probability, uniformly as $\delta \to 0$. We are not aware of other approaches that provide a uniform performance guarantee of this type.

The main idea of our algorithm is to construct an analytically tractable event $C_{\delta }$ that uniformly contains the violation events $\{l\in {\mathbb{R}}^{d_l}|\varphi (x,l)>0\}$ for all x plausible to be feasible. In view of the reformulation of the probabilistic constraint in (${\text{CCP}}_{\delta }$) as $\text{P}(\varphi (x,L)>0|L\in C_{\delta })\le (\delta /\text{P}(L\in C_{\delta }))$, the problem (${\text{CCP}}_{\delta }$) can be solved by the scenario approach where L is sampled from the conditional distribution given $L\in C_{\delta }$. The risk tolerance parameter is adjusted to $\delta /\text{P}(L\in C_{\delta })$. The primary challenge is to construct $C_{\delta }$ as tight as possible so that the new risk tolerance parameter $\delta /\text{P}(L\in C_{\delta })$ is bounded. (This is at the heart of Property 1 defined later. This property is facilitated in the heavy-tailed setting if we assume that $\varphi (x,L)$ has appropriate scaling properties, similar to the distribution of L uniformly over a suitable compact set of decisions.)

We illustrate our assumptions and our framework with a risk problem of independent interest. This problem consists in computing a collective salvage fund in a network of financial entities whose liabilities and payments are settled in an optimal way using the Eisenberg-Noe model (Eisenberg and Noe 2001). The salvage fund is computed to minimize its size to guarantee a probability of collective default after settlements of less than a small prescribed margin. For the sake of demonstrating the broad applicability of our method, we also present a portfolio optimization problem with value-at-risk constraints as an additional running example.

The rest of the paper is organized as follows. In Section 2, we introduce the portfolio optimization problem and the minimal salvage fund problem as particular applications of chance constraint optimization. We use both problems as running examples to provide a concrete and intuitive explanation for the concepts we introduce throughout the paper. In Section 3, we provide a brief review of the scenario approach in Calafiore and Campi (2006). The ideas behind our main algorithmic contributions are given in Section 4, where we introduce its intuition, rooted in ideas originating from rare event simulation. Our algorithm requires the construction of several auxiliary functions and sets, and we summarized the explicit expressions of the sets for the running examples in Table 1. How to do this for a more general setting is detailed in Section 5, in which we also present several additional technical assumptions required by our constructions. In Section 5, we also explain that our procedure results in an estimate that is within a constant factor of the optimal solution of the underlying chance constrained problem with high probability as $\delta \to 0$. In Section 6, we show that the assumptions imposed are valid in our motivating example (as well as a second example with quadratic cost structure inside the probabilistic constraint). Numerical results for the examples are provided in Section 7. Throughout our discussion, in each section we present a series of results that summarize the main ideas of our constructions.

Table 1. Examples of $O_{\delta }$ and $C_{\delta }$ That Satisfy Property 1 When L Is Multivariate Regularly Varying

To keep the discussion fluid, we present the corresponding proofs in Appendix A unless otherwise indicated. In Appendix B, we introduce an importance sampling algorithm to sample from a parametric family of regularly varying distribution. We present additional numerical experiments in Appendix C.

1.1. Notations

In the rest of this paper, ${\mathbb{R}}_{+}=[0,+\infty )$ is the set of nonnegative real numbers, ${\mathbb{R}}_{++}=(0,+\infty )$ is the set of positive real numbers, and $\overline{\mathbb{R}}=[-\infty ,+\infty ]$ is the extended real line. A column vector with zeros is denoted by $\mathbf{0}$, and a column vector with ones is denoted by $\mathbf{1}$. For any matrix Q, the transpose of Q is denoted by $Q^{\top }$; the Frobenius norm of Q is denoted by $\Vert Q{\Vert }_F$. The identity matrix is denoted by I. For $\alpha \in \mathbb{R}$ and $x\in {\mathbb{R}}^d$, we use $\alpha ·x$ to denote the scalar multiplication of x with α. For two column vectors $x,y\in {\mathbb{R}}^d$, we say $x⪯y$ if and only if $y-x\in {\mathbb{R}}_{+}^d$. For a column vector $x\in {\mathbb{R}}^d$ and a scalar $\alpha \in \mathbb{R}$, we say that $x⪯\alpha$ if and only if $x⪯\alpha ·\mathbf{1}$. For $\alpha \in \mathbb{R}$ and $E\subseteq {\mathbb{R}}^d$, we define $\alpha ·E=\{\alpha ·x|x\in E\}$. The optimal value of an optimization problem $(\text{Prob})$ is denoted by $\text{Val}(\text{Prob})$. For any real-valued random variable X with probability measure $\text{P}$, define the inverse tail distribution function ${\overline{F}}_X^{-1}:[0,1]\to \overline{\mathbb{R}}$ as ${\overline{F}}_X^{-1}(\delta )=\inf \{x\in \mathbb{R}|\text{P}(x>X)\le \delta \}$. We also use Landau’s notation. In particular, if $f(·)$ and $g(·)$ are nonnegative real valued functions, we write $f(t)=O(g(t))$ if $f(t)\le c_0\times g(t))$ for some $c_0\in (0,\infty )$ and $f(t)=\mathrm{\Omega }(g(t))$ if $f(t)\ge g(t))/c_0$ for some $c_0\in (0,\infty )$.

2. Running Examples

2.1. Portfolio Optimization with Value-at-Risk Constraint

We first introduce a portfolio optimization problem. Suppose that there are d assets to invest. If we invest a dollar in the jth asset, the investment has mean return μ_j and a nonnegative random loss L_j. Let $x=(x_1,\dots ,x_d)$ represent the amount of dollars invested in different assets, and let $\mu =({\mu }_1,\dots ,{\mu }_d)$ and $L=(L_1,\dots ,L_d)$. We assume that L follows a multivariate heavy-tailed distribution.

A precise definition of this concept is rather involved and will be given in Section 5. Intuitively, $\text{P}(\Vert L{\Vert }_2>x)$ follows a power law, and the direction $L/\Vert L{\Vert }_2$ is assumed to converge in a suitable sense on the unit sphere, conditioned on the event that $\Vert L{\Vert }_2$ is large.

The portfolio manager’s goal is to maximize the mean return of the portfolio, which is equal to ${\mu }^{\top }x$, with a portfolio risk constraint prescribed by a risk measure called value-at-risk (VaR). The VaR at level $1-\delta \in (0,1)$ for a random variable X is defined as

For a given number $\eta >0$, we formulate the following portfolio optimization problem.

Using the definition of VaR and the fact that the cumulative distribution function is right continuous, we conclude that ${\text{VaR}}_{1-\delta }(x^{\top }L)\le \eta$ is equivalent to $\text{P}(x^{\top }L-\eta >0)\le \delta$. To facilitate the technical exposition, we apply the change of variable $x_j\mapsto 1/x_j$ to homogenize the constraint function, yielding the following equivalent chance constrained optimization problem in standard form:

where $\varphi (x,l)={\sum }_{j=1}^d(l_j/x_j)-\eta$. Despite the nonlinear objective, Calafiore and Campi (2005, section 4.3) show that it admits an epigraphic reformulation with a linear objective so that the standard scenario approach is applicable.

2.2. Minimal Salvage Fund

In this section, we use chance-constrained optimization to determine the minimal total salvage fund required for a reinsurance company to control its default probability, where the policy holders have complex liability structures.

A reinsurance policy is a contract sold to insurance companies for transferring the financial risk exposure and smoothing the cash flow. In a certain type of reinsurance contract, the reinsurance company is responsible to pay a fixed percentage of the net liability for its clients (in this paper we assume the percentage is 100% for simplicity). Therefore, the total amount of net liability is the minimal amount of salvage fund required for the reinsurance company to avoid default. However, calculating the distribution of the minimal salvage fund is nontrivial because the clients may also have insurance contracts with each other.

Suppose that the reinsurance company has d clients, each is an entity or an insurance firm. Let $L=(L_1,\dots ,L_d)\in {\mathbb{R}}_{+}^d$ denote the vector of incurred losses by each firm, where L_j denotes the total incurred loss that entity j is responsible to pay. We assume that L follows a multivariate heavy-tailed distribution as in the previous example. Let $Q=(Q_{i,j}:i,j\in \{1,\dots ,d\})$ be a deterministic matrix where $Q_{i,j}$ denotes the amount of money received by entity j when entity i pays one dollar. We assume that $Q_{i,j}\ge 0$ and ${\sum }_{j=1}^d{Q}_{i,j}<1$.

Let $x=(x_1,\dots ,x_d)$ denote the total amount that the salvage fund allocated to each entity, and $y^{\ast }=(y_1^{\ast },\dots ,y_d^{\ast })$ denote the amount of the final settlement. The amount of final settlement is determined by the following optimization problem:

In words, the system maximizes the payments subject to the constraint that nobody pays more than what they have (in the final settlement), and nobody pays more than what they owe. Notice that $y^{\ast }=y^{\ast }(x,L)$ is also a random variable (the randomness comes from L) satisfying $\mathbf{0}⪯y^{\ast }⪯L$.

Suppose that entity j bankrupts if the deficit $L_j-y_j^{\ast }\ge m_j$, where $m\in {\mathbb{R}}_{+}^d$ is a given vector. We are interested in finding the minimal amount of salvage fund that ensures no bankruptcy happens with probability at least $1-\delta .$ The problem can be formulated as a chance constraint programming problem as follows:

Now we write Problem (2) into standard form. Notice that $L-y^{\ast }(x,L)⪯m$ if and only if $\varphi (x,L)\le 0$, where $\varphi (x,L)$ is defined as follows:

Therefore, Problem (2) is equivalent to

3. Review of Scenario Approach

As mentioned in the Introduction, a popular approach to solve the chance constraint problem proceeds by using the scenario approach developed by Calafiore and Campi (2006). They suggest to approximate the probabilistic constraint $\text{P}(\varphi (x,L)>0)\le \delta$ by N sampled constraints $\varphi (x,L^{(i)})\le 0$ for $i=1,\dots ,N$, where $\{L^{(1)},\dots ,L^{(N)}\}$ are independent samples. Instead of solving the original chance constraint problem (${\text{CCP}}_{\delta }$), which is usually intractable, we turn to solve the following optimization problem:

The total sample size N should be large enough to ensure the feasible solution to the sampled problem (${\text{SP}}_N$) is also a feasible solution to the original problem (${\text{CCP}}_{\delta }$) with a high confidence level. According to Calafiore and Campi (2006), for any given confidence level parameter $\beta \in (0,1)$, if

then any feasible solution to the sampled optimization problem (${\text{SP}}_N$) is also a feasible solution to (${\text{CCP}}_{\delta }$) with probability at least $1-\beta$. However, when δ is small, the total number of sampled constraints is of order $\mathrm{\Omega }((1/\delta )\log (1/\delta ))$, which could be a problem for implementation. For example, as we shall see in Section 7, when $\beta ={10}^{-5}$, d = 15 and $\delta ={10}^{-3}$, the number of sampled constraints N is required to be larger than $2\times {10}^5$. In contrast, our method only requires sampling $2\times {10}^3$ constraints.

4. General Algorithmic Idea

To facilitate the development of our algorithm, we introduce some additional notation and a desired technical (property Property 1).

In our setting, a property is an intermediate assumption that facilitates the construction of an efficient scenario approach algorithm. We shall impose the technical property for now, and in Section 5, we will provide assumptions based on more direct model primitives, providing easy-to-verify sufficient conditions for the properties to hold.

We exploit key intuition borrowed from rare event simulation. A common technique exploited, for example, in Chen et al. (2019), is the construction of a so-called super set, which contains the rare event of interest. The super set should be analytically tractable and be constructed with a probability that is of the same order as that of the rare event of interest. If the conditional distribution given being in the super set is accessible, this can be used as an efficient sampling scheme. The first part of this section simply articulates the elements involved in setting the stage for constructing such a set in the outcome space of L. Later, in Section 5, we will impose assumptions in order to ensure that the probability of the superset, which eventually we will denote by $C_{\delta }$ is suitably controlled as $\delta \to 0$. Simply collecting the elements necessary to construct $C_{\delta }$ requires introducing some super sets involving the decision space, since the optimal decision is unknown.

Let $F_{\delta }\subseteq {\mathbb{R}}^{d_x}$ denote the feasible region of the chance constraint optimization problem (${\text{CCP}}_{\delta }$), that is,

Here, the subscript δ is involved to emphasize that the feasible region $F_{\delta }$ is parametrized by the risk level δ. For any fixed $x\in {\mathbb{R}}^{d_x}$, let $V_x≔\{L\in {\mathbb{R}}^{d_l}|\varphi (x,L)>0\}$ denote the violation event at x.

Property 1.

For any $\delta >0$, there exist a set $O_{\delta }\subseteq {\mathbb{R}}^{d_x}$, and an event $C_{\delta }\subseteq {\mathbb{R}}^{d_l}$ that satisfy the following statements.

• (a) The feasible set $F_{\delta }$ is a subset of $O_{\delta }$.

• (b) The event $C_{\delta }$ contains the violation event V_x for any $x\in O_{\delta }$.

• (c) There exist a constant M > 0 independent of δ such that $\text{P}(L\in C_{\delta })\le M·\delta$.

To visualize our intent with Property 1, keep in mind a feature that is often present in heavy-tailed rare-event simulation. In particular, if L is a one-dimensional random variable with, for example, power-law tail decay, then $\text{P}(L>b)\le M\times \text{P}(L>b/2)$ for some $M<\infty$ for all b. For example, if $\text{P}(L>b)=b^{-\alpha },b\ge 1$ we can take $M=2^{\alpha }$. In simple terms, “proportional enlargements” translate into “proportional likelihoods.” This sort of feature can be used to motivate the intent of Property 1 and the selection of event $C_{\delta }$, as it suggests the violation event V_x exhibits “proportional enlargements” when $\delta \to 0$. Specifically, suppose that for some specific $x\prime \in F_{\delta }$, we have that $V_{x\prime }=[b,\infty )$ and the safety constraint is active. That is, $P(L>b)=\delta$ and suppose that the enlarged region is of the form $C_{\delta }=[b/2,\infty )$. Then, if L is regularly varying we will have that (c) in Property 1 holds for all $\delta >0$ (which corresponds to all b large). Generally speaking, if $x\in F_{\delta }$, the set V_x is the set of “bad” outcomes for such a decision. One can imagine that in situations of interest, as we will illustrate, the set of all possible bad outcomes, which is ${\cup }_{x\in F_{\delta }}{V}_x$, can be conveniently enclosed by a region which is a “proportional enlargement” of the set of bad outcomes of a suitable feasible decision (as illustrated in the previous one-dimensional situation). Property 1 implies that the likelihood of the set of all bad outcomes is proportional to the constraint parameter δ. In our algorithms, knowing the constant M will not be relevant, we just need to know that M exists. The sets $O_{\delta }$ and $C_{\delta }$ are auxiliary sets introduced to enclose the set of all possible bad outcomes. We will explain how to construct these sets in examples later.

In the rest of this paper, we will refer to $O_{\delta }$ as the outer approximation set and $C_{\delta }$ as the uniform conditional event. A graphical illustration of $O_{\delta }$ and $C_{\delta }$ is shown in Figure 1.

As hinted in our earlier discussion that motivates Property 1, we shall focus on the case that L follows a multivariate regularly varying distribution (i.e. a multidimensional version of a power-law-type distribution). The definition of multivariate regular variation is provided in Section 5. In this case, to illuminate how the sets $O_{\delta }$ and $C_{\delta }$ can be constructed for different problems, we provide explicit expressions of them for two running examples in Table 1. In Section 5, we will illustrate how to construct $O_{\delta }$ and $C_{\delta }$ for more general settings.

Now, given $O_{\delta }$ and $C_{\delta }$ that satisfies Property 1, we define the conditionally sampled problem (${\text{CSP}}_{\delta ,N\prime }$):

Here, $L_{\delta }^{(i)}$ are independent and identically distributed (i.i.d.) samples generated from the conditional distribution $(L|L\in C_{\delta })$.

We now present our main result of this section in Lemma 1, which validates (${\text{CSP}}_{\delta ,N\prime }$) is an effective and sample efficient scenario approximation by incorporating (Calafiore and Campi 2006, theorem 2) and Property 1. The proof of Lemma 1 will be presented in Section 4.1.

Lemma 1.

Suppose that Property 1 is imposed and let $\beta >0$ be a given confidence level.

1. Let ${\delta }^{\prime }=\delta /\text{P}(L\in C_{\delta })\ge 1/M$ and $N\prime$ be any integer that satisfies

   $$N^{\prime }\ge \frac{2}{{\delta }^{\prime }}\log \frac{1}{\beta }+2d+\frac{2d}{{\delta }^{\prime }}\log \frac{2}{{\delta }^{\prime }}.$$(5)

   With probability at least $1-\beta$, if the conditionally sampled problem (${\text{CSP}}_{\delta ,N\prime }$) is feasible, then its optimal solution $x_N^{\ast }\in F_{\delta }$ and $\text{Val}({\text{CSP}}_{\delta ,N\prime })\ge \text{Val}({\text{CCP}}_{\delta }).$

2. Let $N\prime$ be any integer such that $N\prime \le \beta {\delta }^{-1}\text{P}(L\in C_{\delta })$. Assume that the chance constraint problem (${\text{CCP}}_{\delta }$) is feasible. Then, with probability at least $1-\beta$, (${\text{CSP}}_{\delta ,N\prime }$) is feasible and $\text{Val}({\text{CCP}}_{\delta })\ge \text{Val}({\text{CSP}}_{\delta ,N\prime })$.

Remark 1

(Size of Conditionally Sampled Problem). The lower bound of the sample size given in (5) is not greater than $2\mathrm{M\hspace{0.17em}log}\left(\frac{1}{\beta }\right)+2d+2d\mathrm{M\hspace{0.17em}log}(2M)$, which is independent of δ. Therefore, Lemma 1 shows that the chance constraint problem (${\text{CCP}}_{\delta }$) can be approximated by (${\text{CSP}}_{\delta ,N\prime }$) with sample complexity bounded uniformly as $\delta \to 0$, as long as Property 1 is satisfied.

Remark 2

(Feasibility of Conditionally Sampled Problem). In Lemma 1, part 1, the conditionally sampled problem (${\text{CSP}}_{\delta ,N\prime }$) is feasible with high probability if there exists small δ such that (${\text{CCP}}_{\delta }$) is feasible. In particular, we claim that

$$\text{P}(({\text{CSP}}_{\delta ,N\prime }) \text{is} \text{feasible})\ge {(1-{\delta }_{\text{min}}/\delta )}^{N\prime },$$(6)

where ${\delta }_{\text{min}}=\inf \{\delta \in {\mathbb{R}}_{++}:\text{Change} \text{Constrained} \text{problem} ({\text{CCP}}_{\delta }) \text{is} \text{feasible}\}.$ Recall from Remark 1 that the $N\prime$ can be chosen to be independent of δ; thus when ${\delta }_{\text{min}}$ is small, we have (${\text{CSP}}_{\delta ,N\prime }$) is feasible with high probability. For example, ${\delta }_{\text{min}}=0$ for both the minimal salvage fund problem and the portfolio optimization problem, which implies (${\text{CSP}}_{\delta ,N\prime }$) is almost surely feasible for these two examples.

We next prove (6). For arbitrarily small $ϵ>0$, the feasible region $F_{{\delta }_{\text{min}}+ϵ}$ for problem (${\text{CCP}}_{{\delta }_{\text{min}}+ϵ}$) is nonempty, and thus we can pick $x\in F_{{\delta }_{\text{min}}+ϵ}$ such that $\text{P}(L\in V_x)\le {\delta }_{\text{min}}+ϵ$. If $\delta >{\delta }_{\text{min}}+ϵ$, then $V_x\subseteq C_{\delta }$ and thus $\text{P}(L\in V_x|L\in C_{\delta })\le ({\delta }_{\text{min}}+ϵ)/\text{P}(L\in C_{\delta })\le ({\delta }_{\text{min}}+ϵ)/\delta$. Therefore, by the independence of samples,

$$\text{P}(x \text{is} \text{feasible} \text{for} ({\text{CSP}}_{\delta ,N\prime }))\ge \text{P}(L\notin V_x|L\in C_{\delta }){)}^{N\prime }\ge {(1-({\delta }_{\text{min}}+ϵ)/\delta )}^{N\prime }.$$

Letting $ϵ\to 0$, we conclude that (6) holds.

Remark 3

(Efficient Sampling Algorithm). Efficiently generating samples of $(L|L\in C_{\delta })$ when $\delta \to 0$ requires rare event simulation techniques. For example, when L is light-tailed, exponential tilting can be applied to achieve O(1) sample complexity uniformly in δ; when L is heavy-tailed, with the help of specific problem structure, one can apply importance sampling (Blanchet and Liu 2010) or Markov chain Monte Carlo (Gudmundsson and Hult 2014) to design an efficient sampling scheme. The specific structure of our salvage fund example results in $C_{\delta }$ being the complement of a box, which makes the sampling very tractable if the element of L are independent.

Even if the aforementioned rare event simulation techniques are hard to apply in practice, we can still apply a simple acceptance-rejection procedure to sample the conditional distribution $(L|L\in C_{\delta })$. It costs $O(1/\delta )$ samples of L on average to get one sample of $(L|L\in C_{\delta })$ because $\text{P}(L\in C_{\delta })=O(\delta )$. Consequently, the total complexity for generating $L_{\delta }^{(i)},i=1,\dots ,N\prime$ and solving $({\text{CSP}}_{\delta ,N\prime })$ is $O(1/\delta )$, which is still much more efficient than the scenario approach in Calafiore and Campi (2006), because it requires computational complexity $O({((1/\delta )\log (1/\delta ))}^3)$ for solving a linear programming problem with $O((1/\delta )\log (1/\delta ))$ sampled constraints by the interior point method.

Although Property 1 seems to be restrictive at first glance, we are still able to construct the sets $O_{\delta }$ and $C_{\delta }$ for a rich class of functions $\varphi (x,L)$, including the constraint function for the minimal salvage fund problem. As we shall see in the proof of Lemma 1, once $O_{\delta }$ and $C_{\delta }$ are constructed the sampled problem (${\text{CSP}}_{\delta ,N\prime }$) is a tractable approximation to the problem (${\text{CCP}}_{\delta }$). We explain how to construct the sets $O_{\delta }$ and $C_{\delta }$ in the next section under some additional assumptions. These assumptions relate in particular to the distribution of L. It turns out that, if L is heavy-tailed, the construction of $O_{\delta }$ and $C_{\delta }$ becomes tractable.

4.1. Proof of Lemma 1

If Property 1 is satisfied, (${\text{CCP}}_{\delta }$) is equivalent to

Let ${\delta }^{\prime }≔\delta /\text{P}(L\in C_{\delta })\ge 1/M$ denote the risk level in the equivalent problem (7). The sampled optimization problem related to Problem (7) is given by

where the $L_{\delta }^{(i)}$ are independently sampled from $\text{P}(·|L\in C_{\delta })$. Notice that

According to Calafiore and Campi (2006, corollary 1 and theorem 2), with probability at least $1-\beta$, if the sampled problem (${\text{CSP}}_{\delta ,N\prime }$) is feasible, then the optimal solution to problem (${\text{CSP}}_{\delta ,N\prime }$) is feasible to the chance constraint problem (7). Because (7) and (${\text{CCP}}_{\delta }$) are equivalent, the optimal solution to problem (${\text{CSP}}_{\delta ,N\prime }$) is also feasible to (${\text{CCP}}_{\delta }$). The proof of the first part of the lemma is complete.

Now we turn to prove the second part of the lemma. The equivalence between (${\text{CCP}}_{\delta }$) and (7) is still valid, so it is sufficient to compare the optimal values of (7) and (${\text{CSP}}_{\delta ,N\prime }$). By applying Calafiore and Campi (2006, theorem 2) again, we have with probability at least $1-\beta$ (${\text{CSP}}_{\delta ,N\prime }$) is feasible and the value of (${\text{CSP}}_{\delta ,N\prime }$) is no larger than the optimal value of

The proof is complete by using $1-{(1-\beta )}^{1/N\prime }\ge \beta /N\prime \ge \frac{\delta }{\text{P}(L\in C_{\delta })}$. Therefore, using $\text{Val}$ for “value of,” $\text{Val}(8)\le \text{Val}(7)=\text{Val}({\text{CCP}}_{\delta })$.

5. Constructing Outer Approximations and Summary of the Algorithm

In this section, we come full circle with the intuition borrowed from rare event simulation explained at the beginning of Section 4. The scale-free properties of heavy-tailed distributions (to be reviewed momentarily) coupled with natural (polynomial) growth conditions (like the linear loss) given by the structure of the optimization problem, provide the necessary ingredients to show that the set $C_{\delta }$ has a probability that is of order $O(\delta )$.

In the discussion immediately following Property 1, we imagined that the uniform conditional set $L\in C_{\delta }$ was of the form $L\in [b/2,\infty )$ for $b\to \infty$ as $\delta \to 0$. However, Property 1 can still be enforced if this statement applies to L² or any power of L. This is because power law-type decay (and more generally regular variation) is preserved under power transformations. We will provide assumptions that will enforce that regular variation properties can be applied when estimating the likelihood of the uniform conditional event.

We assume that the distribution of L is of multivariate regular variation. A definition that we now review. For background, we refer to Resnick (2013). Let ${\mathcal{M}}_{+}({\overline{\mathbb{R}}}^{d_l}\backslash \{\mathbf{0}\})$ denote all Radon measures on the space ${\overline{\mathbb{R}}}^{d_l}\backslash \{\mathbf{0}\}$ (recall that a measure is Radon if it assigns finite mass to all compact sets). If ${\mu }_n(·),\mu (·)\in {\mathcal{M}}_{+}({\overline{\mathbb{R}}}^{d_l}\backslash \{\mathbf{0}\})$, then μ_n converges to μ vaguely, denoted by ${\mu }_n\stackrel{v}{\to }\mu$, if for all compactly supported continuous functions $f:{\overline{\mathbb{R}}}^{d_l}\backslash \{\mathbf{0}\}\to {\mathbb{R}}_{+}$,

L is multivariate regularly varying with limit measure $\mu (·)\in {\mathcal{M}}_{+}({\overline{\mathbb{R}}}^{d_l}\backslash \{\mathbf{0}\})$ if

Assumption 1.

L is multivariate regularly varying with limit measure $\mu (·)\in {\mathcal{M}}_{+}({\overline{\mathbb{R}}}^{d_l}\backslash \{\mathbf{0}\})$.

Here are some intuitions behind the definition of multivariate regularly varying. Suppose that L is written in terms of polar coordinates, with R being the radius and Θ being a random variable taking values on the unit sphere. The radius $R=\Vert L{\Vert }_2$ has a one-dimensional regularly varying tail (i.e., we can write $\text{P}(R>x)=L(x)x^{-\alpha }$ for a slowly varying function L and $\alpha >0$). The angle Θ, conditioned on R being large, converges weakly (as $R\to \infty$) to a limiting random variable. The distribution of this limit can be expressed in terms of the measure μ. For a recent application of multivariate regular variation in operations research, see Kley et al. (2016).

In this section, we present two methods for the construction of $O_{\delta }$ and $C_{\delta }$ satisfying Property 1. We mostly focus on our “scaling method” which is presented in Section 5.1, which is facilitated precisely by the scale-free property that we will impose on L. After showing the construction of the outer sets under the scaling method, we summarize the algorithm at the end of Section 5.1. We supply a lower bound guaranteeing a constant approximation for the output of the algorithm in Section 5.2. Our second method for outer approximation constructions is summarized in Section 5.3. This method is simpler to apply because is based on linear approximations; however, it is less general because it assumes that $\varphi (x,L)$ is jointly convex.

5.1. Scaling Method

We start by analyzing the feasible region $F_{\delta }$ when $\delta \to 0$. Intuitively, if the violation probability $\text{P}(\varphi (x,L)>0)$ has a strictly positive lower bound in any compact set, then $F_{\delta }$ will ultimately be disjoint with the compact set when $\delta \to 0$. Thus, the set $F_{\delta }$ is expelled to infinity when $\delta \to 0$ in this case. $F_{\delta }$ is moving toward the direction that $\varphi (x,L)$ becomes small such that the violation probability becomes smaller. For instance, if x is one dimensional and $\varphi (x,L)$ is increasing in x, then $F_{\delta }$ is moving toward the negative direction. Consider the portfolio optimization problem as another example, in which ${\text{min}}_{j=1}^d{x}_j\to +\infty$ as $\delta \to 0$.

With such intuition, now we begin to construct the outer approximation set $O_{\delta }$. To this end, we need to introduce an auxiliary function which we shall call a level function. We assume the existence of a level function in Assumption 2, and the level function needs to be explicitly computable to construct the outer approximation set $O_{\delta }$.

Definition 1.

We say that $\pi :{\mathbb{R}}^{d_x}\to [0,+\infty ]$ is a level function if

1. For any $\alpha \ge 0$ and $x\in {\mathbb{R}}^{d_x}$, we have $\pi (\alpha ·x)=\alpha ·\pi (x)$,

2. The function π(x) is coersive, i.e., ${\lim }_{\delta \to 0}{\inf }_{x\in F_{\delta }}\pi (x)+\infty$.

For a given level function π, we define its unit level set as $\mathrm{\Pi }=\{x\in {\mathbb{R}}^{d_x}|\pi (x)=1\}$.

Assumption 2.

There exists an explicitly computable level function π and its unit level set Π.

The unit level set Π is used to characterize the moving direction of $F_{\delta }$ as $\delta \to 0$. The shape of Π is chosen in accordance with the moving direction of $F_{\delta }$ to reduce the size of $O_{\delta }$ to achieve better sample complexity. The outer approximation set $O_{\delta }$ is constructed as

where ${\alpha }_{\delta }$ characterizes the scaling rate of $F_{\delta }$. We will explain how to choose ${\alpha }_{\delta }$ in the proof of Lemma 2. Here are several examples of the level functions and unit level sets:

• Suppose that $\varphi (x,L)=-\Vert x\Vert {}^2-L$, then the level function π can be chosen as the Euclidean norm and Π can be chosen as the unit sphere in ${\mathbb{R}}^{d_x}$.

• For the portfolio optimization problem, the level function can be chosen as $\pi (x)={\text{min}}_{j=1}^d{x}_j+\infty ·I(x\notin {\mathbb{R}}_{++}^{d_x})$ in accordance with our intuition that ${\text{min}}_{j=1}^d{x}_j\to \infty$, and the unit level set can be chosen as $\mathrm{\Pi }=\{x\in {\mathbb{R}}^{d_x}|{\text{min}}_{j=1}^d{x}_j=1\}$.

To analyze the asymptotic shape of the uniform conditional event $C_{\delta }$, we connect the asymptotic distribution of L to the asymptotic distribution of $\varphi (x,L)$. Keep in mind that we wish to preserve the scaling property of the tail of L, when considering $\varphi (x,L)$, so that Property 1 can be ensured. We pick a continuous nondecreasing function $h:{\mathbb{R}}_{++}\to {\mathbb{R}}_{++}$ such that ${\lim }_{\alpha \to +\infty }h(\alpha )=+\infty$ to characterize the scaling rate of L. In addition, we pick another positive function $r:{\mathbb{R}}_{++}\to {\mathbb{R}}_{++}$ to characterize the scaling rate of $\varphi (\alpha ·x,h(\alpha )·L)$. Intuitively, the scaling function $r(·)$ and $h(·)$ should ensure the condition that the collection of probability measures of ${\{\frac{1}{r(\alpha )}\varphi (\alpha ·x,h(\alpha )·L)\}}_{\alpha \ge 1}$ is tight. For the minimal salvage fund problem with fixed δ, as the deficit $\varphi (x,L)$ is asymptotically linear with respect to the salvage fund x and the loss L, we can simply pick $r(\alpha )=h(\alpha )=\alpha$ in this problem. We next introduce two auxiliary functions ${\mathrm{\Psi }}_{+}$ and ${\mathrm{\Psi }}_{-}$.

Definition 2.

Let ${\mathrm{\Psi }}_{+}:{\mathbb{R}}^{d_l}\to \mathbb{R},\hspace{0.17em}{\mathrm{\Psi }}_{-}:{\mathbb{R}}^{d_l}\to \mathbb{R}$ be two Borel measurable functions. We say ${\mathrm{\Psi }}_{+}$ (respectively, c) is the asymptotic uniform upper (respectively, lower) bound of $\frac{1}{r(\alpha )}\varphi (\alpha ·x,h(\alpha )·l)$ over the unit level set $x\in \mathrm{\Pi }$ if for any compact set $K\subseteq {\mathbb{R}}^{d_l}$,

$$\underset{\alpha \to \infty }{\lim \inf }\underset{l\in K}{\inf }\left({\mathrm{\Psi }}_{+}(l)-\underset{x\in \mathrm{\Pi }}{\sup }\left[\frac{1}{r(\alpha )}\varphi (\alpha ·x,h(\alpha )·l)\right]\right)\ge 0,$$(9a)

$$\underset{\alpha \to \infty }{\lim \sup }\underset{l\in K}{\sup }\left({\mathrm{\Psi }}_{-}(l)-\underset{x\in \mathrm{\Pi }}{\inf }\left[\frac{1}{r(\alpha )}\varphi (\alpha ·x,h(\alpha )·l)\right]\right)\le 0.$$(9b)

We would like to have lower and upper bounds ${\mathrm{\Psi }}_{-}$ and ${\mathrm{\Psi }}_{+}$ so that $\varphi$ is of the order $r(\alpha )$ for every decision x, and in every direction l of the random vector L, whenever the norm of the latter is large (i.e., of the size $h(\alpha )$). A stronger assumption would have been to require an actual limit (rather than a liminf and a limsup), but this is not needed to provide big-O bounds for the complexity of our algorithm.

The functions ${\mathrm{\Psi }}_{+}$ and ${\mathrm{\Psi }}_{-}$ are used to define the event $C_{\epsilon ,-}$ and $C_{\epsilon ,+}$, which serve as the inner and outer approximation of the event ${\cup }_{x\in \mathrm{\Pi }}{V}_x$, where $V_x=\{l\in {\mathbb{R}}^{d_l}|\varphi (x,l)>0\}$ is the violation event at x.

Definition 3.

For $\epsilon >0$, let $C_{\epsilon ,+}$ (respectively, $C_{\epsilon ,-}$) be the ε-outer (respectively, inner) approximation event:

$$C_{\epsilon ,+}≔\{l\in {\mathbb{R}}^{d_l}|{\mathrm{\Psi }}_{+}(l)\ge -\epsilon \},$$(10a)

$$C_{\epsilon ,-}≔\{l\in {\mathbb{R}}^{d_l}|{\mathrm{\Psi }}_{-}(l)\ge +\epsilon \}.$$(10b)

We now define $O_{\delta }≔{\cup }_{\alpha \ge {\alpha }_{\delta }}\alpha ·\mathrm{\Pi }$. The following property ensures that the shape of Π is appropriate and ${\alpha }_{\delta }$ is large enough; hence, $O_{\delta }$ is an outer approximation of $F_{\delta }$.

Property 2.

There exist δ₀ such that for any $\delta <{\delta }_0$, we have an explicitly computable constant ${\alpha }_{\delta }$ that satisfies

$$\text{P}(\Vert L{\Vert }_2>h({\alpha }_{\delta }))=O(\delta )\hspace{1em}\hspace{1em}\mathit{\text{and}}\hspace{1em}\hspace{1em}{F}_{\delta }\subseteq {\displaystyle \underset{\alpha \ge {\alpha }_{\delta }}{\cup }\alpha }·\mathrm{\Pi }≕O_{\delta }.$$

If the violation probability is easy to analyze, we will directly derive the expression of ${\alpha }_{\delta }$ and verify Property 2. Otherwise, we resort to Lemma 2, which provides a sufficient condition of Property 2 by analyzing the asymptotic probability of the violation event as $\delta \to 0$.

Lemma 2.

Suppose that Assumptions 1 and 2 hold. If there exists an asymptotic uniform lower bound function ${\mathrm{\Psi }}_{-}(·)$ as given in (9b) and $\epsilon >0$ such that $\mu (C_{\epsilon ,-})>0$, then Property 2 is satisfied.

The high-level idea of Lemma 2 is to show that $F_{\delta }$ is disjoint from $\alpha ·\mathrm{\Pi }$ for small α. To this end, the asymptotic scaling (9b) is used to demonstrate that the violation probability is no less than $\text{P}(L\in h(\alpha )·C_{\epsilon ,-})$, which is approximately equal to $\text{P}(\Vert L{\Vert }_2>h({\alpha }_{\delta }))·\mu (C_{\epsilon ,-})$ according to the regularly varying property of L. The detailed proof of Lemma 2 is deferred to Appendix A.1.

We impose the following Assumption 3 on the asymptotic uniform upper bound ${\mathrm{\Psi }}_{+}(·)$ so that we can use the multivariate regular variation of L to estimate $\text{P}(L\in \alpha ·C_{\epsilon ,+})$ for large scaling factor α.

Assumption 3.

There exist an event $S\subseteq {\mathbb{R}}^{d_l}$ with $\mu (S^c)<\infty$ such that

$$S\subseteq \alpha ·S,\hspace{1em}\hspace{1em}{\mathrm{\Psi }}_{+}(l)\le {\mathrm{\Psi }}_{+}(\alpha ·l),\hspace{1em}\hspace{1em}\forall l\in S,\hspace{0.17em}\alpha \ge 1.$$

In addition, there exist some $\epsilon >0$ such that $C_{\epsilon ,+}$ is bounded away from the origin, that is, ${\inf }_{l\in C_{\epsilon ,+}}\Vert l{\Vert }_2>0.$

Moreover, both S and $C_{\epsilon ,+}$ have explicit expressions.

For the minimal salvage fund problem, because the deficit function $\varphi (x,L)$ is coordinate-wise nondecreasing with respect to the loss vector L, it is reasonable to assume that its asymptotic bound ${\mathrm{\Psi }}_{+}(·)$ is also coordinate-wise nondecreasing. For this example, the closed form expression of ${\mathrm{\Psi }}_{+}(·)$ and the detailed verification of all the assumptions are deferred to Proposition 2. Our next result summarizes the construction of the outer approximation sets.

Theorem 1.

Suppose that Property 2 and Assumption 3 are imposed. Then there exist ${\delta }_0>0$ such that the following sets

$$O_{\delta }={\displaystyle \underset{\alpha \ge {\alpha }_{\delta }}{\cup }\alpha }·\mathrm{\Pi },\hspace{1em}\hspace{1em}{C}_{\delta }=h({\alpha }_{\delta })·(C_{\epsilon ,+}\cup K^c\cup S^c)$$(11)

satisfy Property 1 for all $\delta <{\delta }_0$. Here, S is given in Assumption 3 and K is a ball in ${\mathbb{R}}^{d_l}$ with $\mu (K^c)<\infty$.

The main idea for proving Theorem 1 is as follows: If L lies in a “well-behaved” compact region, then by applying Assumption 3 and the asymptotic uniform lower bound (9b), the violation events ${\cup }_{x\in O_{\delta }}{V}_x$ is uniformly enclosed in $h({\alpha }_{\delta })·C_{\epsilon ,+}$. Otherwise L lies in the “ill-behaved” region $h({\alpha }_{\delta })·(K^c\cup S^c)$. Combining these two cases inspires that definition of $C_{\delta }≔h({\alpha }_{\delta })·(C_{\epsilon ,+}\cup K^c\cup S^c)$, and the probability of $L\in C_{\delta }$ is $O(\delta )$ because of $\text{P}(\Vert L{\Vert }_2>h({\alpha }_{\delta }))=O(\delta )$ and L is regularly varying. The detailed proof of Theorem 1 is deferred to Appendix A.1.

With the aid of Lemma 1 and Theorem 1, we provide an algorithm for approximating (${\text{CCP}}_{\delta }$) in which the sampled optimization problem is bounded in $1/\delta$.

Algorithm 1

(Scenario Approach with Optimal Scenario Generation)

input: Risk tolerance parameter δ, confidence level β, and all the elements and constants appearing in Property 2 and Assumption 3, including level function π or unit level set Π, constant ${\alpha }_{\delta }$, scaling function h, and explicit expression of $C_{\epsilon ,+},K$ and S.

• 1 Compute the expression of sets $O_{\delta }$ and $C_{\delta }$ by (11);

• 2 Compute required number of samples $N\prime$ by (5);

• 3 for $i=1,\dots ,N\prime$ do

• 4  Sample $L_{\delta }^{(i)}$ using acceptance-rejection or importance sampling.

• 5 end

• 6 Solve the conditionally sampled problem (${\text{CSP}}_{\delta ,N\prime }$).

5.2. Constant Approximation Guarantee

In Section 5.2, our objective is to show that the output of the previous algorithm is guaranteed to be within a constant factor of the optimal solution to (${\text{CCP}}_{\delta }$) with high probability, uniformly in δ.

We shall work under the setting of Theorem 1, so we enforce Property 2 and Assumptions 3. We want to show that there exist some constant $\mathrm{\Lambda }>1$ independent of δ, such that $\text{Val}({\text{CCP}}_{\delta })\le \text{Val}({\text{CSP}}_{\delta ,N\prime })\le \mathrm{\Lambda }\times \text{Val}({\text{CCP}}_{\delta })$ with high probability. This indicates that our result guarantees a constant approximation to (${\text{CCP}}_{\delta }$) for regularly varying distributions (under our assumptions) in O(1) sample complexity when $\delta \to 0$ with high probability.

Note that $({\text{CSP}}_{\delta ,N\prime })\le \mathrm{\Lambda }\times \text{Val}({\text{CCP}}_{\delta })$ is meaningful only if $\text{Val}({\text{CCP}}_{\delta })>0$. We assume that the outer approximation set is good enough such that the following natural assumption is valid.

Assumption 4.

There exist $\delta >0$ such that ${\text{min}}_{x\in O_{\delta }}{c}^{\top }x>0$.

The previous assumption will typically hold if c has strictly positive entries. Theorem 1 and the form of $O_{\delta }$ guarantee that the norm of the optimal solution of (${\text{CSP}}_{\delta ,N\prime }$) grows in proportion to ${\alpha }_{\delta }$, so we also assume the following scaling property for $\varphi (x,l)$.

Assumption 5.

There exist a function ${\varphi }_{\lim }:({\mathbb{R}}^{d_x}\backslash \{\mathbf{0}\})\times ({\mathbb{R}}^{d_l}\backslash \{\mathbf{0}\})\to \mathbb{R}$ such that for every compact set $E\subseteq {\mathbb{R}}^{d_l}\backslash \{\mathbf{0}\}$, we have

$$\underset{\alpha \to \infty }{\lim }\underset{l\in E}{\sup }\left|\frac{1}{r(\alpha )}\varphi (\alpha ·x,h(\alpha )·l)-{\varphi }_{\lim }(x,l)\right|=0.$$

In addition, ${\varphi }_{\lim }(x,l)$ is continuous in one.

Assumption 5 is satisfied by both running examples. For the portfolio optimization problem, we have $\varphi (x,l)={\sum }_{j=1}^d(L_j/x_j)-\eta$, thus ${\varphi }_{\lim }(x,l)=\varphi (x,l)$. For the minimal salvage fund problem, we have ${\varphi }_{\lim }(x,l)=\varphi (x,l)-m$ such that $|{\alpha }^{-1}\varphi (\alpha ·x,\alpha ·l)-{\varphi }_{\lim }(x,l)|\le {\alpha }^{-1}m$ and $|{\varphi }_{\lim }(x,l)-{\varphi }_{\lim }(x,l\prime )|\le \Vert l-l\prime {\Vert }_1$.

We define the following optimization problem, which will serve as an asymptotic upper bound of $({\text{CSP}}_{\delta ,N\prime })$ in stochastic order when $\delta \to 0$:

where $L_{\stackrel{(i)}{\lim }}$ are i.i.d. samples from a random variable $L_{\lim }$, whose distribution is characterized by $\text{P}(L_{\lim }\in (C_{\epsilon ,+}\cup K^c\cup S^c))=1$ and $\text{P}(L_{\lim }\in E)=\mu (E)/\mu (C_{\epsilon ,+}\cup K^c\cup S^c)$ for all measurable set $E\subseteq C_{\epsilon ,+}\cup K^c\cup S^c$.

Theorem 2.

Let $\beta >0$ be a given confidence level and $N\prime$ be a fixed integer that satisfies (5). If Assumptions 4 and 5 are enforced, and (${\text{CSP}}_{\lim ,N\prime }$) satisfies Slater’s condition with probability one, then there exist ${\delta }_0>0$ and $\mathrm{\Lambda }>0$ such that

$$\text{P}(\text{Val}({\text{CCP}}_{\delta })\le \text{Val}({\text{CSP}}_{\delta ,N\prime })\le \mathrm{\Lambda }\times \text{Val}({\text{CCP}}_{\delta }))\ge 1-2\beta ,\hspace{1em}\hspace{1em}\forall \delta <{\delta }_0.$$

In Theorem 2, the Slater’s condition (see section 5.2.3 in Boyd and Vandenberghe (2004) for reference) can be verified directly on the problem (${\text{CSP}}_{\lim ,N\prime }$). This condition is satisfied in the salvage fund problem by standard linear programming duality. We also remark that Assumptions 4 and 5 only require the existence rather than the explicit knowledge of $\{\delta |{\text{min}}_{x\in O_{\delta }}{c}^{\top }x>0\}$ and function ${\varphi }_{\lim }$.

5.3. Linear Approximation Method

Suppose that the constraint function $\varphi (x,l)$ is jointly convex in (x, l), and L is multivariate regularly varying. We will develop a simpler method in this section to construct the outer approximation set $O_{\delta }$ and the uniform conditional event $C_{\delta }$.

We first introduce a crucial assumption in the construction of $O_{\delta }$ and $C_{\delta }$.

Assumption 6.

There exist a convex piecewise linear function ${\varphi }_{-}(x,l):{\mathbb{R}}^{d_x}\times {\mathbb{R}}^{d_l}\to \mathbb{R}$ of the form

$${\varphi }_{-}(x,l)=\underset{j=1,\dots ,N}{\text{max}}{a}_j^{\top }l+b_j^{\top }x+c_j,\hspace{1em}\hspace{1em}{a}_j\in {\mathbb{R}}^{d_l},b_j\in {\mathbb{R}}^{d_x} \text{and} c_j\in \mathbb{R} \text{for} j=1,\dots ,N.$$

such that

1. The inequality ${\phi }_{-}(x,l)\le \phi (x,l), \text{holds\hspace{0.17em}for\hspace{0.17em}all\hspace{0.17em}}(x,l)\in {\mathbb{R}}^{d_x}\times {\mathbb{R}}^{d_l}$;

2. There exist some constant $C\in {\mathbb{R}}_{+}$ such that $\varphi (x,l)\le 0$ if ${\varphi }_{-}(x,l)\le -C$.

If $\varphi (x,l)$ itself is a piecewise affine function, then Assumption 6 is satisfied by simply taking ${\varphi }_{-}(x,l)=\varphi (x,l)$. For general jointly convex functions, the following lemma verifies Assumption 6 if $\varphi (x,l)$ has a compact zero sublevel set.

Lemma 3.

If the constraint function $\varphi (x,L):{\mathbb{R}}^{d_x}\times {\mathbb{R}}^{d_l}\to \mathbb{R}$ is convex and twice continuously differentiable, and it has a compact zero sublevel set $Z_{\varphi }≔\{(x,l)\in {\mathbb{R}}^{d_x}\times {\mathbb{R}}^{d_l}|\varphi (x,l)\le 0\}$, then Assumption 6 is satisfied.

With Assumption 6 enforced, we are now ready to provide our main result in this section to fully summarize the construction of $O_{\delta }$ and $C_{\delta }$.

Theorem 3.

If Assumptions 1 and 6 hold, we can construct $O_{\delta }$ and $C_{\delta }$ that satisfy Property 1 as

$$O_{\delta }≔{\displaystyle \underset{j=1}{\overset{N}{\cap }}\{x\in {\mathbb{R}}^{d_x}|b_j^{\top }x+c_j+{\overline{F}}_{a_j^{\top }L}^{-1}(}\delta )\le 0\},\hspace{1em}\hspace{1em}{C}_{\delta }≔{\displaystyle \underset{j=1}{\overset{N}{\cup }}\{L\in {\mathbb{R}}^{d_l}|a_j^{\top }L+C>{\overline{F}}_{a_j^{\top }L}^{-1}(}\delta )\},$$

where ${\overline{F}}_{a_j^{\top }L}^{-1}(\delta )=\inf \{x\in \mathbb{R}|\text{P}(x>a_j^{\top }L)\le \delta \}$.

6. Verifying the Assumptions in Examples

In this section, we verify the elements required to apply our algorithm. We provide explicit expressions for sets $O_{\delta }$ and $C_{\delta }$ in the statement of the propositions. The detailed verification process and the steps for constructing sets $O_{\delta }$ and $C_{\delta }$ are presented as the proofs in Appendix A.2.

6.1. Portfolio Optimization with VaR Constraint

In this section, we will verify that Theorem 1 is applicable to an equivalent form of the portfolio optimization problem (2).

Proposition 1.

The portfolio optimization problem (1) satisfies all assumptions required by Theorem 1, such that the sets $O_{\delta }$ and $C_{\delta }$ admit the following explicit expressions:

$$O_{\delta }=\{x\in {\mathbb{R}}_{++}^d|\eta ·x⪰{\overline{F}}_{{\mathbf{1}}^{\top }L}^{-1}(\delta )\},\hspace{1em}\hspace{1em}{C}_{\delta }=\{l\in {\mathbb{R}}_{++}^d|2·{\mathbf{1}}^{\top }l\ge {\overline{F}}_{{\mathbf{1}}^{\top }L}^{-1}(\delta )\}.$$

6.2. Minimal Salvage Fund

The key observation to solve the minimal salvage fund problem (3) is the following lemma, which provides a closed form piecewise linear expression for the constraint function $\varphi (x,L)$.

Lemma 4.

In the minimal salvage fund problem (3), we have

$$\varphi (x,L)=\underset{j=1,\dots ,d}{\text{max}}\hspace{0.17em}{L}_j-{\mathbf{e}}_j^{\top }{(I-Q^{\top })}^{-1}x-m_j,$$

where ${\mathbf{e}}_j$ denote the unit vector on the jth coordinate.

Now we prove that Theorem 3 is applicable to the minimal salvage fund problem (4).

Proposition 2.

The minimal salvage fund problem (3) satisfies all assumptions required by Theorem 3, such that the sets $O_{\delta }$ and $C_{\delta }$ admit the following explicit expressions:

$$O_{\delta }={\displaystyle \underset{j=1}{\overset{d}{\cap }}\{x\in {\mathbb{R}}^d|{\overline{F}}_{L_j}^{-1}(}\delta )\le {\mathbf{e}}_j^{\top }{(I-Q^{\top })}^{-1}x+m_j\},\hspace{1em}\hspace{1em}{C}_{\delta }={\displaystyle \underset{j=1}{\overset{d}{\cup }}\{l\in {\mathbb{R}}^d|l_j>{\overline{F}}_{L_j}^{-1}(}\delta )\}.$$

6.3. Quadratic Model

In this section, we consider a model with a quadratic control term in x as an additional example. Suppose that the constraint function $\varphi (x,l):{\mathbb{R}}^{d_x}\times {\mathbb{R}}^{d_l}\to \mathbb{R}$ is defined as

where $Q\in {\mathbb{R}}^{d_x\times d_x}$ is a symmetric matrix and $A\in {\mathbb{R}}^{d_x\times d_l}$ is a matrix with $\text{rank}(A)=d_x$; that is, there exists $\sigma >0$ such that $\Vert A^{\top }x{\Vert }_2\ge \sigma \Vert x{\Vert }_2$.

Proposition 3.

Consider the chance constraint optimization model with constraint function defined as (12).

1. If Q is a positive semidefinite matrix and L has a positive density, there exist some δ such that the problem is infeasible.

2. If Q has a negative eigenvalue and L is multivariate regularly varying, the model satisfies all the assumptions required by Theorem 1.

7. Numerical Experiments

To empirically study the computational complexity and compare the quality of the solutions, in this section, we conduct numerical experiments for two scenario generation algorithms:

1. The efficient scenario generation approach proposed in this paper (abbreviated as Eff-Sc)

2. The scenario approach in Calafiore and Campi (2006) (abbreviated as CC-Sc)

In Section 7.1, we present the results for the portfolio optimization problem. In Section 7.2, we present the results for the minimal salvage fund problem. The numerical experiment is conducted using a Laptop with a 2.2-GHz Intel Core i7 CPU, and the sampled linear programming problem is solved using CVXPY (Diamond and Boyd 2016) with the MOSEK solver (MOSEK ApS 2020).

7.1. Portfolio Optimization with VaR Constraint

First, we present the parameter selection and the implement details for the numerical experiment of portfolio optimization problem (1). Suppose that there are d = 10 assets to invest, and the parameters of the problem are chosen as follows:

• The mean return vector is $\mu =(1.0,1.5,2.0,2.5,3,1.6,1.2,1.1,1.8,2.2)$.

• The random variable L_j are i.i.d. with Pareto cumulative distribution function $\text{P}(L_j>l)=({\mathcal{\ell }}_j/l)$, for $l\ge {\mathcal{\ell }}_j$.

• The parameter $\mathcal{\ell }=({\mathcal{\ell }}_1,\dots ,{\mathcal{\ell }}_d)=(2.1,1.3,1.6,2.5,2.7,1.3,1.9,1.5,2.2,2.3)$.

• The loss threshold η = 1,000.

Now we explain the implementation detail of Eff-Sc. Recall the expression of $O_{\delta }$ and $C_{\delta }$ from Proposition 1, which involves the analytically unknown quantity ${\overline{F}}_{{\mathbf{1}}^{\top }L}^{-1}(\delta )$. Because quantile estimation is much more computationally efficient than solving the sampled optimization problem, we generate samples of L to estimate a confidence interval of ${\overline{F}}_{{\mathbf{1}}^{\top }L}^{-1}(\delta )$ with large enough confidence level $1-o(\beta )$, and we denote the resulting confidence interval by $(\widehat{\text{LB}},\widehat{\text{UB}})$. We replace the expressions of $O_{\delta }$ and $C_{\delta }$ by their sampled version conservative approximations, that is,

The value of $\text{P}(L\in C_{\delta })$ is also estimated using the generated samples. We compute the required number of samples $N\prime$ using Lemma 1, and the samples of $L_{\delta }$ is generated via acceptance-rejection.

In Figure 2, we compare the efficiency between Eff-Sc and CC-Sc. Figure 2(a) presents the required number of samples for both algorithms, in which one can quickly remark that Eff-Sc requires significantly fewer samples than CC-Sc, especially for the problems with small δ. In Figure 2(b), we compare the running time for both models. Whereas Eff-Sc costs slightly more time for δ around 0.1 due to the overhead cost of computing $O_{\delta }$ and $C_{\delta }$, the computational time stays nearly constant uniformly in δ, indicating that Eff-Sc is a substantially more efficient algorithm than CC-Sc.

Finally, we compare Eff-Sc and CC-Sc for the optimal values of the sampled problems and the violation probabilities of the optimal solutions. Because both methods require generating random samples, the generated solutions are also random. Thus, the optimal values and the violation probabilities are also random. To compare the distributions of the random quantities, we conduct 10³ independent experiments. In each experiment, we execute both algorithms and get two solutions, then we evaluate the solutions’ violation probabilities using 10⁶ samples of L. We use boxplots (McGill et al. 1978) to depict the samples’ distribution through their quantiles. A boxplot is constructed of two parts: a box and a set of whiskers. The box is drawn from the 25% quantile to the 75% quantile, with a horizontal line drawn in the middle to denote the median. Two whiskers indicate 5% and 95% quantiles, respectively, and the scatters represent all the rest sample points beyond the whiskers.

In Figure 3, we present (a) the optimal values and (b) the violation probabilities. One can quickly remark from Figure 3(a) that the optimal value of Eff-Sc is stochastically larger than the optimal value of CC-Sc, whereas Figure 3(b) indicates that the optimal solutions produced by both methods are feasible for all the 10³ experiments. Overall, with both methods successfully and conservatively approximating the probabilistic constraint, Eff-Sc is more computationally efficient and less conservative, producing solutions with better objective values than its counterpart.

7.2. Minimal Salvage Fund

In this section, we conduct a numerical experiment for the minimal salvage fund problem (3). In the experiment, we pick $d\in \{10,15,20\}$ to test the performance of the problem in different dimensions.

For each fixed d, the parameters of Problem (3) are chosen as follows:

• The matrix $Q=(Q_{i,j}:i,j\in \{1,\dots ,d\})$ where $Q_{i,j}=1/d$ if $i\ne j$ and otherwise $Q_{i,j}=0$.

• The vector $m=(m_j:j\in \{1,\dots ,d\})$ where m_j = 10 for each j.

• The random variables L_j are i.i.d. with Pareto cumulative distribution function $\text{P}(L_j>l)=(1/l)$, for $l\ge 1$.

Recall the explicit expressions for sets $O_{\delta }$ and $C_{\delta }$ from Proposition 2. To solve the conditionally sampled problem (${\text{CSP}}_{\delta ,N\prime }$), it remains to sample $L_{\delta }^{(i)}$ and compute $N\prime$, the required number of samples. When δ is small, when $\delta \le {10}^{-3}$, solving the optimization problem (${\text{CSP}}_{\delta ,N\prime }$) costs much more time than simulating $L_{\delta }^{(i)}$, despite that a simple acceptance rejection scheme is applied to sample $L_{\delta }^{(i)}$ in our experiments. We fix the confidence level parameter $\beta ={10}^{-5}$ and set ${\delta }^{\prime }=\delta /P(L\in C_{\delta })\ge d^{-1}$, and then we can compute $N\prime$ by the first part of Lemma 1.

Similar to Figure 2 of the portfolio optimization problem, we compare the efficiency between Eff-Sc and CC-Sc for different d and δ in Figure 4, in terms of (a) the required number of samples and (b) the CPU time for solving the sampled approximation problem. We observe that the Eff-Sc has uniformly smaller sample complexity and computational complexity than CC-Sc, where the superiority becomes significant for small δ. In particular, the required number of samples and the used CPU time are bounded for Eff-Sc, whereas they quickly deteriorate for CC-Sc when δ becomes smaller. It is also worth noting that Eff-Sc is consistently more efficient than CC-Sc for all the tested dimensions.

Finally, we compare optimal values of the sampled problems and violation probabilities of the optimal solutions in Figure 5. We present in Figure 5(a) the optimal values and in Figure 5(b) the violation probabilities, with fixed dimension d = 15 (we provide additional results for d = 5 and d = 10 in Appendix C.3). One can quickly remark from Figure 5(a) that the optimal value of Eff-Sc is stochastically smaller than the optimal value of CC-Sc, whereas Figure 5(b) indicates that the optimal solutions produced by both methods are feasible for all the 10³ experiments. Therefore, we are able to draw the same conclusion as we have from the portfolio optimization experiment: Eff-Sc efficiently produces less conservative solutions.