{'fraction_non_alphanumeric': 0.046841315558008384, 'fraction_numerical': 0.015383819395340014, 'mean_word_length': 4.083194529674761, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Validating Sample Average Approximation Solutions with Negatively Dependent Batches

Sample-average approximations (SAA) are a practical means of finding approximate solutions of stochastic programming problems involving an extremely large (or infinite) number of scenarios. SAA can also be used to find estimates of a lower bound on the optimal objective value of the true problem which, when coupled with an upper bound, provides confidence intervals for the true optimal objective value and valuable information about the quality of the approximate solutions. Specifically, the lower bound can be estimated by solving multiple SAA problems (each obtained using a particular sampling method) and averaging the obtained objective values. State-of-the-art methods for lower-bound estimation generate batches of scenarios for the SAA problems independently. In this paper, we describe sampling methods that produce negatively dependent batches, thus reducing the variance of the sample-averaged lower bound estimator and increasing its usefulness in defining a confidence interval for the optimal objective value. We provide conditions under which the new sampling methods can reduce the variance of the lower bound estimator, and present computational results to verify that our scheme can reduce the variance significantly, by comparison with the traditional Latin hypercube approach.


Introduction
Stochastic programming provides a means for formulating and solving optimization problems that contain uncertainty in the model. When the number of possible scenarios for the uncertainty is extremely large or infinite, sample-average approximation (SAA) provides a means for finding approximate solutions for a reasonable expenditure of computational effort. In SAA, a finite number of scenarios is sampled randomly from the full set of possible scenarios, and an approximation to the full problem of reasonable size is formulated from the sampled scenarios and solved using standard algorithms for deterministic optimization (such as linear programming solvers). Solutions obtained from SAA procedures are typically suboptimal. Substantial research has been done in assessing the quality of an obtained solution (or a set of candidate solutions), and in understanding how different sampling methods affect the quality of the approximate solution.
Important information about a stochastic optimization problem is provided by the optimality gap (Mak et al. 1999, Bayraksan andMorton 2006), which provides a bound on the difference between the objective value for the computed SAA solution and the true optimal objective value. An estimate of the optimality gap can be computed using upper and lower bounds on the true optimal objective value. Mak et al. (1999) proves that the expected objective value of an SAA problem is a lower bound of the objective of the true solution, and that this expected value approaches the true optimal objective value as the number of scenarios increases. We can estimate this lower bound (together with confidence intervals) by solving multiple SAA problems, a task that can be implemented in parallel in an obvious way. An upper bound can be computed by taking a candidate solution x and evaluating the objective by sampling from the scenario set, typically taking a larger number of samples than were used to set up the SAA optimization problem for computing x.
Much work has been done to understand the quality of SAA solutions obtained from Monte Carlo (MC) sampling, Latin hypercube (LH) sampling (McKay et al. 1979), and other methods. MC generates independent identically distributed scenarios where the value of each variable is picked independently from its corresponding distribution. The simplicity of this method has made it an important practical tool; it has been the subject of much theoretical and empirical research. Many results about the asymptotic behavior of optimal solutions and values of MC have been obtained; see (Shapiro et al. 2009, Chapter 5) for a review. By contrast with MC, LH stratifies each dimension of the sample space in such a way that each strata has the same probability, then samples the scenarios so that each strata is represented in the scenario sample set. This procedure introduces a dependence between the different scenarios of an individual SAA problem. The sample space is in some sense "uniformly" covered on a per-variable basis, thus ensuring that there are no large unsampled regions and leading to improved performance. Linderoth et al. (2006) provides empirical results showing that the bias and variance of a lower bound obtained by solving multiple SAA problems constructed with LH sampling is considerably smaller than the statistics obtained from an MC-based procedure. Theoretical results about the asymptotic behavior of these estimates were provided later by Homem-de Mello (2008). Other results on the performance of LH have been obtained, including results on large deviations (Drew and Homem-de Mello 2005), rate of convergence to optimal values (Homem-de Mello 2008), and bias reduction of approximate optimal values (Freimer et al. 2012), all of which demonstrate the superiority of LH over MC. This favorable empirical and theoretical evidence makes LH the current state-of-the-art method for obtaining high-quality lower bounds and optimality gaps via SAA. In this paper, we build on the ideas behind the LH method to obtain LH variants with even better variance properties.
In the past, when solving a set of SAA problems to obtain a lower-bound estimate of the true optimal objective value, each batch of scenarios determining each SAA was chosen independently of the other batches. In this paper, we introduce two approaches to sampling -sliced Latin hypercube (SLH) sampling and sliced orthogonal-array Latin hypercube (SOLH) sampling -that yield better estimates of the lower bound by imposing negative dependencies between the batches in the different SAA approximations. These approaches not only stratify within each batch (as in LH) but also between all batches. The SLH approach is easy to implement, while the SOLH approach provides better variance reduction. With these methods, we can significantly reduce the variance of the lower bound estimator even if the size of each SAA problem or the number of SAA problems were kept the same, which can be especially useful when solving each SAA problem is time consuming or when computing resources are limited. We will provide theoretical results analyzing the variance reduction properties of both approaches, and present empirical results demonstrating their efficacy across a variety of stochastic programs studied in the literature. Sliced Latin hypercube sampling was introduced first in Qian (2012) and has proven to be useful for the collective evaluation of computer models and ensembles of computer models.
Here we briefly outline the rest of our paper. The next section begins with a brief description of how the optimality gap can be estimated ( § 2.1), a review of Latin hypercube sampling ( § 2.2) and an introduction of functional analysis of variance ( § 2.3). Section 3 focuses on sliced Latin hypercube sampling. It outlines the construction of dependent batches ( § 3.1), describes theoretical results of variance reduction based on a certain monotonicity condition ( § 3.2), applies the results to two-stage stochastic linear program ( § 3.3), and finally studies the relation between the lower bound estimation and numerical integration ( § 3.4). Section 4 reviews orthogonal arrays and introduces a method to incorporate these arrays into the sliced Latin hypercube sampling, which leads to stronger between-batch negative dependence. The next two sections deal with our computational experiments. Section 5 describes the setup and some of the implementation details, while Section 6 describes and analyzes the performance of the new sampling methods in the lower bound estimation problem. We end the paper in Section 7 with a summary of our results and a discussion of possible future research.

Stochastic Programs and Lower Bound Estimators
We consider a stochastic program of the form where X ⊂ R p is a compact feasible set, x ∈ X is a vector of decision variables, ξ = (ξ 1 , ξ 2 , . . . , ξ m ) is a vector of random variables, and G : R p × R m → R is a real-valued measurable function. Unless stated otherwise, we assume that ξ is a random vector with uniform distribution on (0, 1] m and that E is the expectation with respect to the distribution of ξ. If ξ has a different distribution on R m , we can transform it into a uniform random vector on (0, 1] m as long as ξ 1 , ξ 2 , . . . , ξ m are either (i) independent discrete or continuous random variables or (ii) dependent random variables which are absolutely continuous (Rosenblatt 1953). Problem (1) may be extremely challenging to solve directly, since the evaluation of g(x) involves a high-dimensional integration. We can replace (1) with the following approximation: where ξ 1 , ξ 2 , . . . , ξ n are scenarios sampled from the uniform distribution on (0, 1] m . The functionĝ n is called a sample-average approximation (SAA) to the objective g in (1). In this paper we will frequently use the term SAA problem to refer to equation (2). We use x * and v * to denote a solution and the optimal value of the true problem (1), respectively, while x * n and v * n denote the solution and optimal value of the SAA problem (2), respectively.
We introduce here some items of terminology that are used throughout the remainder of the paper. Let D denote an n × m matrix with ξ T i in (2) as its ith row. Hence, D represents a batch of scenarios that define an SAA problem. We will refer ξ i to as the ith scenario in D. We use D(:, k) to denote the kth column of D, and ξ ik to denote the (i, k) entry of D, that is, the kth entry in the ith scenario in D.
We can express the dependence of v * n in (2) on D explicitly by writing this quantity as v n (D), where v n : (0, 1] n×m → R. Given a distribution over the D matrices where ξ 1 , ξ 2 , . . . , ξ n are each drawn from the uniform (0, 1] m distribution but not necessarily independently, it is well known and easy to show that the expectation with respect to the D matrices E[v n (D)] ≤ v * giving us a lower bound of the true optimal value (Norkin et al. 1998, Mak et al. 1999). E[v n (D)] can be estimated as follows. Generate t independent batches D 1 , D 2 , . . . , D t and compute the optimal value min v n (D r ) by solving subproblem (2) for each D r , r = 1, 2, . . . , t. From Mak et al. (1999), a lower bound estimate of v * is (3)

Latin Hypercube Sampling
Latin hypercube sampling, which stratifies sample points along each dimension (McKay et al. 1979), has been used in numerical integration for many years. It has been shown that the variance of the mean output of a computer experiment under Latin hypercube sampling can be much lower than for experiments based on Monte Carlo methods (McKay et al. 1979, Stein 1987, Loh 1996. Let v M C n (D) and v LH n (D) denote the approximate optimal value when the ξ i in D are generated using Monte Carlo and Latin hypercube sampling, respectively. Homem-de Mello (2008) showed that the asymptotic variance of v LH n (D) is smaller than the variance of v M C n (D) under some fairly general conditions. Extensive numerical simulations have provided empirical demonstrations of the superiority of Latin hypercube sampling (Homem-de Mello 2008, Linderoth et al. 2006.
To define Latin hypercube sampling, we start with some useful notation. Given an integer p ≥ 1, we define Z p := {1, 2, . . . , p}. Given an integer a, the notation Z p ⊕ a denotes the set {a + 1, a + 2, . . . , a + p}. For a real number y, y denotes the smallest integer no less than y. A "uniform permutation on a set of p integers" is obtained by randomly taking a permutation on the set, with all p! permutations being equally probable.
We have the following definition.
Definition An n × m array A is a Latin hypercube if each column of A is a uniform permutation on Z n . Moreover, A is an ordinary Latin hypercube if all its columns are generated independently.
Using an ordinary Latin hypercube A, an n × m matrix D with scenarios ξ 1 , ξ 2 , . . . , ξ n that defines an SAA problem is obtained as follows (McKay et al. 1979): where all γ ik are U [0, 1) random variables, and the quantities a ik and the γ ik are mutually independent. We refer the matrix D thus constructed as an ordinary Latin hypercube design.
We now introduce a different way of looking at design matrices D ∈ (0, 1] n×m that will be useful when we discuss extensions to sliced Latin hypercube designs in later sections. We can write a design matrix D as where When D is an ordinary Latin hypercube design, B is an ordinary Latin hypercube and θ ik corresponds to γ ik in (4) for all i = 1, 2, . . . , n and k = 1, 2, . . . , m. By the properties of an ordinary Latin hypercube design, the entries in Θ are mutually independent, and Θ and B are independent. We refer B to as the underlying array of D.
The lower bound on v * can be estimated from (3) by taking t independently generated ordinary Latin hypercube designs D 1 , D 2 , . . . , D t (Linderoth et al. 2006). We denote this sampling scheme by ILH and denote the estimator obtained from (3) by L ILH n,t . To illustrate the limitations of the ILH scheme, Figure 1 displays three independent 3 × 3 ordinary Latin hypercube designs generated under ILH with n = t = m = 3. Scenarios from each three-dimensional design are denoted by the same symbol, and are projected onto each of the three bivariate planes. The dashed lines stratify each dimension into three partitions. For any design, each of the these three intervals will contain exactly one scenario in each dimension. This scheme covers the space more "uniformly" than three scenarios that are picked identically and independently from the uniform distribution, as happens in Monte Carlo schemes. However, the combination of points from all three designs does not cover the space particularly well, which gives some cause for concern, since all designs are being used in the lower bound estimation. Specifically, when we combine the three designs together (to give nine scenarios in total), it is usually not the case that each of the nine equally spaced intervals of (0,1] contains exactly one scenario in any dimension. This shortcoming provides the intuition behind the sliced Latin hypercube (SLH) design, which we will describe in the subsequent sections.

Functional ANOVA Decomposition
In order to understand the asympototic properties of estimators arising from Latin-Hypercube based samples, it is necessary to review the functional analysis of variance decomposition (Owen 1994), also known as functional ANOVA. Let D := {1, 2, . . . , m} represent the axes of (0, 1] m associated with an input vector ξ = (ξ 1 , ξ 2 , . . . , ξ m ) defined in Section 2.1. Let F k denote the uniform distribution for ξ k , with F := m k=1 F k . For u ⊆ D, let ξ u denote a vector consisting of ξ k for k ∈ u. Define where dF D−u = k / ∈u dF k integrates out all components except for those included in u, and v ⊂ u is a proper subset of u. Hence, we have that is the main effect function for factor k, and between factors k and l.

Sliced Latin Hypercube Sampling
Instead of generating D 1 , D 2 , . . . , D t independently for each SAA subproblem, we propose a new scheme called sliced Latin hypercube (SLH) sampling that generates a family of correlated designs. An SLH design (Qian 2012) is a nt × m matrix that can be partitioned into t separate LH designs, represented by the matrices D r , r = 1, 2, . . . , t, each having the same properties as ILH, with respect to SAA. SLH was originally introduced to aid in the collective evaluation of computer models, but here we demonstrate its utility in creating multiple SAA problems to solve.
Let L SLH n,t denote the lower bound estimator of v * under SLH. Because the individual designs D r , r = 1, 2, . . . , t are LH designs, we have that E(L SLH n,t ) = E(L ILH n,t ). Consider two distinct batches of scenarios D r and D s for any r, s = 1, 2, . . . , t and r = s. We will show that when v n (D) fulfills a certain monotonicity condition, v n (D r ) and v n (D s ) are negatively dependent under SLH. Compared with ILH, SLH introduces between-batch negative dependence while keeping the within-batch structure intact. As a result, we expect a lower-variance estimator: var(L SLH n,t ) ≤ var(L ILH n,t ).
Step 3. For r = 1, 2, . . . , t, obtain the (i, k)th entry of D r as follows: where the γ r,ik are U [0, 1) random variables that are mutually independent of the a r,ik .
By vertically stacking the matrices D 1 , D 2 , . . . , D t , we obtain the nt × m matrix representing the SLH design, as defined in Qian (2012).
As in (5), we can express each D r as where B r = (b r,ik ) n×m with b r,ik = nξ r,ik and Θ r = (θ r,ik ) n×m . We have the following proposition regarding properties of the SLH design, including dependence of the batches. (This result is closely related to (Qian 2012, Lemma 2).) Proposition 3.1 Consider the SLH design with D r , r = 1, 2, . . . , t constructed according to Algorithm 1, with B r and Θ r , r = 1, 2, . . . , t defined as in (9). The following are true.
(iv) For r, s = 1, 2, . . . , t with r = s, θ r,ik is dependent on θ s,ik if and only if B r,ik = B s,ik ; (v) The D r , r = 1, 2, . . . , t are ordinary Latin hypercube designs, but they are not independent.
Proof (i) According to (9) and Step 2 in the above construction, B r is the same as A r in Step 1 prior to the replacement step, and the result follows.
where the a r,ik are values in A r after the replacement in Step 2 of the construction above. By (8) and (9), we have According to Step 2 of the construction, for each r = 1, 2, . . . , t, i = 1, 2, . . . , n and k = 1, 2, . . . , m, the quantity t a r,ik t − a r,ik is randomly selected among Z t and is independent of B r . Since the γ r,ik are independent of the a r,ik , the claim is proved.
(iii) For each r = 1, . . . , t and k = 1, . . . , m, the quantities {t a r,1k t −a r,1k , t a r,2k t −a r,2k , . . . , t a r,nk t − a r,nk } are independently and randomly selected among n different Z t 's, respectively. Thus, the θ r,ik are mutually independent within Θ r . In other words, the t a r,ik t − a r,ik , i = 1, 2, . . . , n are mutually independent discrete uniform random variables on Z t , such that the θ r,ik are U [0, 1) random variables, by (10).
(iv) It suffices to show that t a r,1k t − a r,1k and t a r,nk t − a s,nk are dependent if and only if B r,ik = B s,ik . That is true because t a r,1k t − a r,1k and t a r,nk t − a s,nk are selected from the same Z t when B r,ik = B s,ik .
(v) The result follows directly from (i), (ii), (iv), and the definition of the ordinary Latin hypercube design. Figure 2 displays the bivariate projection of the three 3 × 3 ordinary Latin hypercube designs, each denoted by a different symbol, which are generated under an SLH scheme. For each design, each of the three equally spaced intervals of (0, 1] contains exactly one scenario in each dimension. In contrast to Figure 1, when we combine the three designs together, each of the nine equally spaced intervals of (0, 1] contains exactly one scenario in any one dimension.

Monotonicity Condition
We derive theoretical results to show var(L SLH n,t ) ≤ var(L ILH n,t ) under a monotonicity condition that will be defined shortly.
The monotonicity condition can be checked by directly studying the function v n (D), but it can also be shown to be satisfied if the stochastic program has certain nice properties. Later we will prove some sufficient conditions for two-stage stochastic linear programs and give a simple argument to show how some problems from the literature have the monotonicity property. Qian (2012) proves that the function values of any two scenarios in a sliced Latin hypercube design are negatively quadrant dependent. The next lemma extends this result, showing the function values of any two batches in a sliced Latin hypercube design are also negatively quadrant dependent, under the monotonicity assumption on v n (D).
Proof Given two ordinary Latin hypercubes B r and B s in (9), let D r = (B r −Θ r )/n and D s = (B s −Θ s )/n denote two slices generated by (8). Since the underlying Latin hypercubes are fixed for D r and D s , the only random parts in the definition are Θ r and Θ s .

Taking expectations of both sides gives
The last equality holds because B r and B s are independent, by Proposition 3.1 (i).
The next result is an immediate consequence of Lemma 3.2. It indicates that the variance of the lower-bound estimator could be reduced by using SLH, when v n (D) is monotonic. in (3) obtaind under ILH and SLH, respectively. Suppose v n (D) in monotonic, then for any n and t, we have that var(L SLH n,t ) ≤ var(L ILH n,t ).
Even if the monotonicity condition does not hold, we can prove similar statements about the asymptotic behavior of the variance of ILH and SLH.
Theorem 3.4 gives an asymptotic result that shows the same conclusion can be drawn as in Theorem 3.3 even if the monotonicity condition does not hold.
Theorem 3.4 Let D denote an n × m Latin hypercube design based on a Latin hypercube B such that where the last equality in (11) holds because batches are exchangeable. Let H(∆ 1 ; B 1 ) = v n (D 1 ) and where the second-last equality is from (Qian 2009, Lemma 1). The result follows by substituting (12) into (11). The functional ANOVA decomposition is properly defined because the quantities δ k are independent, according to Proposition 3.1 (iii).
The second equality holds because E[v n (D)|B] is the same regardless of the underlying ordinary Latin hypercube B when v n (D) is additive. The third equation holds due to the functional ANOVA decomposition on H(∆; B). We complete the proof by combining (13) with the result in (ii).

Two-Stage Linear Program
We now discuss the theoretical properties of SLH for two-stage stochastic linear programs. Consider problems of the form min where X is a polyhedron and and ξ := (h, q, T ). The problem has fixed recourse since the recourse matrix W does not depend on the random variable ξ. Defining G(x, ξ) to be the function c T x + Q(x, ξ), we see that (14) is a special case of (1). Let x = (x 1 , x 2 , . . . , x p ) and T = (T kj ). By (15), Q(x, ξ) is a decreasing function of any entry in h, for any x ∈ X. Furthermore, Q(x, ξ) is a decreasing function of any entry T kj in T if x j is nonpositive, and an increasing function of any entry T kj in T if x j is nonnegative. By LP duality, we have and hence Q(x, ξ) is an increasing function of any entry in q for any x ∈ X. We conclude that G(·, ξ) is monotonic in each component of ξ if the recourse matrix W is fixed.
Lemma 3.5 Let v n (D) in (2) denote the approximated optimal value of the two-stage stochastic program in (14) with fixed recourse. Then v n (D) is monotonic if (i) T is fixed or (ii) for every j = 1, 2, . . . , p, x j is always nonnegative or nonpositive, given any D.
Proof For any ordinary Latin hypercube design D, let x n be the arg min of the approximation problem, that is, where ξ i is the ith scenario of D. Without loss of generality, increase only ξ 1k to obtain a new design D * with scenarios ξ * 1 , ξ 2 , . . . , ξ n . Let x * n be the arg min of the approximation problem with design D * , that is, If either (i) or (ii) is satisfied, the value of x does not affect whether G(·, ξ 1 ) is increasing or decreasing in Example Consider the newsvendor problem from Freimer et al. (2012), which can be expressed as a twostage stochastic program. In the first stage, choose an order quantity x. After demand ξ has been realized, we decide how much of the available stock y to sell. Assume that demand is uniformly distributed on the interval (0, 1], and there is a shortage cost α ∈ (0, 1) and an overage cost 1 − α. The second stage problem is P : Since the first-stage cost is zero, the two-stage stochastic program is The optimal value is v * = α(1 − α)/2 with solution x * = α.
Based on a sample of n demands ξ 1 , ξ 2 , . . . , ξ n , the approximated optimal value is given by v n (D) = min where D = (ξ 1 , ξ 2 , . . . , ξ n ) T . The optimal solution x * n is the αn th order statistic of {ξ 1 , ξ 2 , . . . , ξ n }. Let α = .4, for which v * = 0.12. Table 1 gives means and standard errors for the estimators of L n,t = E[v n (D)] for several values of n and t. This table shows that SLH reduces the variance of L n,t significantly when compared with ILH. Analytically, we have v n (D) = H(∆; B) = n −2 where r * = αn and B is an arbitrary underlying Latin hypercube for D. We notice that for any B, H(∆; B) is decreasing in δ 1 , . . . , δ r * −1 , increasing in δ r * +1 , . . . , δ n and monotonic in δ r * (the direction depends on α). Thus, v n (D) is monotonic, which can alternatively be checked by applying Lemma 3.5. By Theorem 3.3, we should have var(L SLH n,t ) ≤ var(L ILH n,t ). We notice that var(L ILH n,t ) = O(t −1 ) while var(L SLH n,t ) = o(t −1 ), a fact that can be explained by Theorem 3.4 (iii), since the newsvendor problem only has one random variable and v n (D) is additive.

Discrete Random Variables
Theorems 3.3 and 3.4 confirm the effectiveness of sliced Latin hypercube designs in reducing the variance of lower-bound estimates in SAA. However, the assumptions in those theorems limit their applicability to fairly specialized problems. Theorem 3.3 does not apply to two-stage problem in (14) with random recourse. Theorem 3.4 does not apply when n → ∞, which is a more practical assumption than t → ∞. In this section, we consider the case in which ξ 1 , ξ 2 , . . . , ξ m are discrete random variables, as discussed by Homem-de Mello (2008). In fact, we assume ξ 1 , ξ 2 , . . . , ξ m to be independent discrete random variables. We plan to show that estimating the lower bound of v * is almost equivalent to a numerical integration problem. Several existing results regarding numerical integration provide us with tools for studying effects of different sampling schemes for lower bound estimation more generally.
Assumption 3.2 The feasible set X is compact and polyhedral; the function G(·, ξ) is convex piecewise linear for every value of ξ; and the distribution of ξ has finite support.
Assumption 3.1 holds under Latin hypercube sampling if E [G(x, ξ)] 2 < ∞ for every x ∈ X, by (Loh 1996, Theorem 3). Assumption 3.2 holds in practice for stochastic linear programs in which the random vector is discretized. Let S * denote the set of optimal solutions of (1). Based on both assumptions, Homem-de Mello (2008) shows the consistency of the approximated optimal solution in (2).
The following proposition connects two-stage stochastic linear program in (14) and numerical integration.
Proposition 3.7 Consider problem (14), and suppose that E [G(x,ξ)] 2 < ∞ for every x ∈ X, and that S * = {x * } is a singleton. Then for n sufficiently large, we have Proof Let x * n (D r ) denote an optimal solution to the SAA problem with scenarios given by D r for r = 1, 2, . . . , t, we have L n,t = t −1 t r=1 v n (D r ) = (nt) −1 t r=1 n i=1G (x * n (D r ),ξ r,i ). Because x * n (D r ) = x * w.p.1 for n sufficiently large by Proposition 3.6 andG(·,ξ) is continuous, the result follows according to the Continuous Convergence Theorem (Mann and Wald 1943).
Proposition 3.7 indicates that L n,t becomes an integral estimator ofG(x * , ξ r,i ) and that L n,t is directly estimating the true optimal value v * of (1), for n large enough. Results about the variance formula for L ILH n,t and L SLH n,t are given in the next result.
Theorem 3.8 Consider problem (14). With the same conditions in Proposition 3.7, we have that the following results hold as n, t → ∞: (ii) The result holds as a consequence of (Qian 2012, Theorem 2).
(iii) Since any batch sampled under SLH is statistically equivalent to an ordinary Latin hypercube design, by Proposition 3.1 (v), the variances var [v n (D i )] are the same as those under ILH. Due to exchangeability of batches and scenarios within the same batch, it suffices to show that under SLH, we have cov G (x * ,ξ 1,1 ),G(x * ,ξ 2,1 ) ≤ 0, whereξ r,i denote the ith scenario in D r . For 0 < z 1 , z 2 ≤ 1 and an integer p, define α p (z 1 , z 2 ) = 1, pz 1 = pz 1 0, o.w.
To prove the second claim, we define Again, by following the proof of (Stein 1987, Theorem 1), we have Because each Ξ k is finite, the number of non-zero values of ζ n i (ξ k ) is also finite. Hence, var(L ILH n,t ) = t −1 var[v n (D 1 )] = O(n −2 t −1 ). Under SLH, using the same arguments as (17) Let ν k denote the number of all possible distinct values for ξ k in Ξ k , and denote by p (i)k the probability of ξ taking its ith smallest possible value. To sample ξ k , we will first takeξ 1 ,ξ 2 , . . . ,ξ m in n equally spaced subintervals in (0, 1], respectively. If anyξ k taken from the same subinterval always leads to the same value for ξ k = F −1 k (ξ k ), then v n (D) will entirely depend on the underlying ordinary Latin hypercube B. The following proposition gives sufficient conditions for equal performance between ILH and SLH.
Proof Divide (0, 1] into the n equal subintervals as (0, 1/n], (2/n, 3/n], . . . , ((n − 1)/n, 1], and let I(i, n) = ((i − 1)/n, i/n], as before. Assume p (0)k = 0. For any k = 1, 2, . . . , m, when np (1)k is an integer, we have ξ k = F −1 k (ξ k ) equal to the smallest possible value in Ξ k , provided thatξ k is drawn from I(1, n), I(2, n), . . . , I(np (1)k , n). Similarly, when np (2)k is an integer, we have that ξ k = F −1 k (ξ k ) equal to the second smallest value in Ξ k , provided thatξ k is taken from I(np (1)k +1, n), I(np (1)k +2, n), . . . , I(np (2)k , n), and so on. As a result, ξ 1 , ξ 2 , . . . , ξ n is determined only by B. That is, given b ik ,ξ ik is taken from I(b ik , n) and ξ ik = F −1 k (ξ ik ) is the th smallest value in Ξ k , where is an integer and np ( −1)k + 1 ≤ b ik ≤ n ( )k . Under SLH, v n (D r ) and v n (D s ) would be independent for any r = s, because they depend on B r and B s , which are independent according to Proposition 3.1 (i). Theorem 3.8 (i) and (ii) indicate that ILH and SLH are equally effective, in general, for estimating v * when n, t → ∞. Proposition 3.9 gives a specific case in which SLH is exactly the same as ILH. Another type of sliced Latin hypercube design is introduced in § 4, which possesses similar (if not the same) within-batch structure as ILH and SLH, but much stronger between-batch negative dependence than SLH to further reduce the variance var(L n,t ) of the lower-bound estimator.

Sliced Orthogonal Array Based Latin Hypercube Sampling
Section 3 indicates that SLH may yield significant improvement in estimating L n,t when the monotonicity condition in is satisfied or whenG(x * ,ξ k ) is an additive function inξ k , under Assumptions 3.1 and 3.2. Both the monotonicity and additivity conditions emphasize individual random variables. As a consequence, the results in Theorems 3.3 and 3.8 are intuitive, because combining all batches under SLH (rather than ILH) gives a design with better stratification in each dimension; see Figures 1 and 2. On the other hand, the combined design under SLH does not possess better stratification when we consider groups of variables. Theorem 3.8 (i) and (ii) suggest that we would need a better sampling scheme than SLH if neither the monotonicity nor additivity conditions are satisfied.
We now introduce another sampling scheme with the following properties. First, the design for each SAA subproblem in this scheme is an ordinary Latin hypercube, as for ILH and SLH. Fixing this property between our different sampling schemes allows us to better study and understand the benefits of improving the between-batch stratification. Second, when we choose the number of batches t sufficiently large, the increased size of the combined design matrix achieves better stratification not just for each variable but also for every pair of variables. With the additional stratification, the var(L n,t ) can be further reduced, under the assumptions in Proposition 3.7.

Orthogonal Arrays
Orthogonal arrays originate from the pioneering work of Rao (1946Rao ( , 1947Rao ( , 1949. Patterson (1954) introduced lattice sampling based on randomized orthogonal arrays, which is found to be suitable for computer experiments and other related fields (Owen 1992). Tang (1993) and Owen (1994) independently studied the Monte Carlo variance of means over orthogonal-array-based Latin hypercube designs and randomized orthogonal arrays, respectively.
We introduce some properties and notation for orthogonal arrays that pertain to Latin hypercube designs. Let P denote a set of s symbols or levels. To be consistent with the notation for a Latin hypercube A in § 2.2, we denote levels by 1, 2, . . . , s. The following formal definition of orthogonal arrays is due to Hedayat et al. (1999).
Definition An N × m array C = (c ik ) with entries from P is said to be an orthogonal array with s levels, strength τ and index λ (for some τ satisfying 0 ≤ τ ≤ m) if every N × τ subarray of C contains each τ −tuple based on P exactly λ times as a row.
We denote an orthogonal array as OA (N, m, s, τ ), where N , m, s, τ , and λ are all integers, with N = λs τ . An ordinary Latin hypercube is an orthogonal array with τ = 1. Furthermore, combined hypercubes under ILH and SLH are special cases of orthogonal arrays with τ = 1. We summarize these observations in the following proposition.
Proposition 4.1 Let D 1 , D 2 , . . . , D t denote t Latin hypercube designs of dimension n × m, which will be used to solve t SAA problems. Let D = (ξ ik ) denote the N × m design by stacking D 1 , D 2 , . . . , D t row by row, such that N = nt. Let A = (a ik ) denote an N × m array such that a ik = nξ ik . Let B = (b ik ) denote an N × m array such that b ik = N ξ ik . Then A is an OA(N, m, n, 1) under ILH and SLH. Furthermore, B is an OA (N, m, N, 1) under SLH.
Proposition 4.1 reveals the fact that L n,t is constructed based on an OA(N, m, n, 1) and an OA(N, m, N, 1) under ILH and SLH, respectively. Theorems 3.3, 3.4, and 3.8 (iii) all indicate that OA(N, m, N, 1) is superior to OA(N, m, n, 1) in estimating L n,t . In other words, if two orthogonal arrays have the same values of N and m, and have τ = 1, the one with a larger value of s and a smaller value of λ would be more desirable. In the remainder of this section, we use τ = 2, but the same result might still apply when generalized to larger values of τ .
There exist many methods for constructing orthogonal arrays with τ = 2, most of which use Galois fields and finite geometry. We summarize some popular constructions and their restrictions for N , m, s, and λ in Table 2.  Bush (1952) OA(s 2 , m, s, 2) m ≤ s + 1 s is a prime power Bose and Bush (1952) OA(λs 2 , m, s, 2) m ≤ λs + 1 (s and λ are powers of the same prime) Addelman and Kempthorne (1961) OA(2s 2 , m, s, 2) m ≤ 2s + 1 (s is an odd prime power) Addelman and Kempthorne (1961) OA(2s 3 , m, s, 2) m ≤ 2s 2 + 2s + 1 (s is an odd prime power, 2 or 4) Table 3 presents two strength-two orthogonal arrays. Both have sixteen scenarios and five columns, the left one having four levels while the right one has just two levels. Intuitively, the left one seems preferable, as it includes all 16 possible level combinations in any two columns. Owen (1994) defines the concept of coincidence defect, which can be used to more formally justify the superiority of the left array which has more levels. An orthogonal array A with strength τ has coincidence defect if there exist two rows of A that agree in τ + 1 columns. The left array in Table 3 does not have coincidence defect because no two rows of A agree in more than a single column. The right one contains coincidence defects; for example, the second and the third rows agree in columns 2, 3, and 4.
Comparing orthogonal arrays can be difficult. As in Owen (1994), we define ω ij (u) for each u ⊆ D := {1, 2, . . . , m} as ω ij (u) := {k ∈ u|c ik = c jk }. We further define for u ⊆ D and r = 0, 1, . . . , |u| to be the number of pairs of ith and jth rows in an orthogonal array C (i can be the same as j) such that c i and c j agree on exactly r of the axes in u. For an OA(N, m, s, 2) without coincidence defect, M (u, 3) = N for any u ⊆ D. For the array on the right in Table 3, M (u, 3) may be much larger than 16 -for example, M ({1, 2, 3}, 3) = 3N . In general, we would like to select the orthogonal array C with no coincidence defect such that there are no duplicates in any three columns of C.
If we are forced to use orthogonal arrays with coincidence defect, we should pick the one with the smallest value of M (u, r). Discussion on the existence of coincidence defects for orthogonal arrays constructed using the methods of Table 2 can be found in Owen (1994).

Construction
Let SOLH denote the scheme that generates batches D 1 , D 2 , . . . , D t as slices of a Latin hypercube design based on sliced orthogonal arrays. The original purpose of SOLH was to share strength across all batches for numerical integration with higher accuracy. Given an OA(N, m + 1, t, 2) with N = nt, n = λt, and s = t

Theoretical Results
We provide a general variance formula of L n,t under SOLH introduced in the last section. The result is a direct consequence of (Hwang et al. 2013, Theorem 1).
Theorem 4.2 Consider problem (14), and suppose that the conditions in Proposition 3.7 hold. Based on an OA(N, m + 1, t, 2) with N = nt, n = λt, and s = t symbols, we have as s → ∞ that where u is a subset of D, andξ u containsξ k for k ∈ u.
Theorem 4.2 shows that the variance of L n,t can be reduced under SOLH after filtering out the bivariate interactions inG(x * ,ξ). This fact remains true for SOLH based on any strength-two orthogonal array even if coincidence defects are present, provided that s → ∞. According to the functional ANOVA decomposition and Theorem 3.8 (i), we have var(L SLH n,t ) = N −1 |u|≥2 var[G u (x * ,ξ u )] + o(N −1 ), where coefficients of variances due to interactions with more than two variables are all of order N −1 . As a result, SOLH based on an orthogonal array with M (u, |u|) > N , for some u ⊆ D and |u| ≥ 3, can inflate the variance due to higher-order interactions. This side effect is less significant if the bivariate interactions dominate higher order interactions, which is true in most problems.
In addition, we need to emphasize that each batch under SOLH is based on an OA(n, m, t, 1), which is a Latin hypercube design instead of an ordinary Latin hypercube. We conjecture that Assumption 3.1 still holds in this case, and the bivariate variance reduction is due to the between-batch negative dependence rather than the within-batch non-ordinary Latin hypercube structure. This conjecture is consistent with our computational observations, as shown below.

Practical Considerations
Using the SOLH method is not just a matter of picking the desired n, m values of the dimensions of the design matrices and the number of batches t, like it was in the SLH case. Picking the appropriate orthogonal array requires some understanding of the tradeoffs between the different choices.
Consider an OA(n 2 , m + 1, n, 2). The benefits of using these arrays include (i) each batch is based on an ordinary Latin hypercube, and (ii) there is no coincidence defect. However, we have to solve n batches of n scenarios each. We want n to be large enough to ensure that the quantity x * n converges and that the bias v * − E[v n (D)] is small, but solving n batches could be computationally infeasible.
One way around this computational hurdle is to relax the restriction that we use all the batches that are generated by the SOLH method. We can use just the first t batches, and we will call this the sliced partial orthogonal-aray-based Latin hypercube method (SPOLH). Let υ denote the fraction of the number of batches generated by SOLH. SPOLH can perform better than SLH since υ of the variance due to bivariate interactions can be removed. However, when υ is small (10%, say), it does not guarantee substantial variance reduction from SLH to SPOLH.
We can alternatively use orthogonal arrays with a higher index value (i.e. λ > 1) allowing us to have a t that is much smaller than n. In this case, there is coincidence defect and picking and generating the right orthogonal arrays with desirable n, m, t values with low M (u, r) can be difficult.

Experimental Setup
To illustrate the effectiveness of negatively dependent designs for estimation of the lower bound on the objective value in stochastic programs, we perform computational experiments using data sets from the literature. This section provides some details on our implementations and test problems.

Implementation
We estimate a single lower bound by choosing a design family, obtaining a set of sampled approximations using this design family, and finally solving these approximations. The performance of a sampling scheme is assessed by repeating this process to obtain a number of lower bound estimates. We then calculate the mean and variance of these estimates.
We augmented the SUTIL library (Czyzyk et al. 2000), a C and C++ based library for manipulating stochastic programming instances with a design-based sampling framework. This augmented SUTIL library can take an orthogonal array as input and generate an orthogonal-array-based design family as output. The library can generate other design families from scratch. We used a C library due to Owen (available at http://lib.stat.cmu.edu/) to generate the orthogonal arrays. The SUTIL software produces the extensive form (or deterministic equivalent) linear program of each sampled instance and outputs the linear program as an MPS file. These files are fed into the commercial linear programming solver Cplex v12.5 and solved using the barrier method. All the timed experiments were run on an Intel Xeon X5650 (24 cores @ 2.66Ghz) server with 128GB of RAM, and 16 of the 24 cores available on the machine were used. Other experiments were run on the HTCondor grid (Thain et al. 2005) of the Computer Sciences department at University of Wisconsin-Madison. They required more than a week of wall-clock time to complete.

Test Problems
Our five test problems were drawn from the stochastic programming literature. In fact, they are the same problems that were studied in Linderoth et al. (2006). These problems are specified in SMPS file format (a stochastic-programming extension of MPS), and have finite discrete distributions for their random variables.
20term, first described in Mak et al. (1999), is a model of a motor freight carrier's operations. The first stage variables represent positions of a fleet of vehicles at the beginning of a day. The secondstage variables determine the movements of the fleet through a network to satisfy point-to-point demands for shipments (with penalties for unsatisfied demands), and to end the day with the fleet back in their initial positions.
gbd is derived from an aircraft allocation problem originally described in the textbook of Dantzig (1963, Chapter 28). Four different types of aircraft are to be allocated to routes to maximize the profit under uncertain demand for each route. There are costs associated with using each aircraft type, and when the capacity does not meet demand.
LandS is a modification by Linderoth et al. (2006) of a simple problem in electrical investment planning described by Louveaux and Smeers (1988). The first-stage variables represent capacities of different new technologies, and the second-stage variables represent production of each of the different modes of electricity from each of the technologies.
ssn (Sen et al. 1994) is a problem from telecommunications network design. The owner of a network sells private-line services between pairs of nodes in the network. During the first stage of the problem, the owner has a budget for adding bandwidth (capacity) to the edges in the network. In the second stage, to satisfy uncertain demands for service between each pair of nodes, short routes between two nodes with sufficient total bandwidth must be identified. Unmet demands incur costs, and the goal of the problem is to minimize the expected costs.
storm is a cargo flight scheduling problem described by Mulvey and Ruszczyński (1993). The goal is to schedule cargo-carrying flights over a set of routes in a network, where the amount of cargo delivered to each node is uncertain. In the first stage, the number of planes of each type on each route is decided. In the second stage, the random variables are the demands on the amounts of cargo to be delivered between nodes. The goal is to minimize (1) the costs that comes from assigning the planes and balancing payloads, and (2) the penalties associated with unmet demands.
All of these problems fit our monotonicity assumption defined in § 3.2. The intuition for this property is straightforward -the objective in each problem is to minimize costs and penalties that come from unmet demands. Hence, as the demand increases, the objective value always increases.
We outline some facts about the random variables and optimal solution for each problem in Table 4. The random variables for all the problems are independent from each other. For the distribution column, we use 'uniform' to refer to a uniform distribution over all possible values, and 'irregular' for any other distribution.
Additionally, we note that the distribution of each of the random variables of LandS approximate a linear function as it has a uniform distribution over evenly spaced points. The upper-and lower-bound estimates (with 95% confidence intervals) are obtained from (Linderoth et al. 2006, Table 4) with n = 5000. We will refer to the quantities in this table in our discussion of the computational results.

Computational Results
In this section we will summarize the computational results and observations. For each combination of sampling scheme, number of scenarios, and number of batches, we perform the lower bound estimation process 100 times, and compute the mean and standard error of the 100 lower bounds obtained. Before we proceed with our observations, we should emphasize that while we use the same problems as in Linderoth et al. (2006), our results are not directly comparable. In particular, Linderoth et al. (2006) focus on results based on analysis of the objective values across 10 uncorrelated subproblems of a single replicate (or similarly, across 10 replicates of one subproblem each). We focus on the analysis of the mean SAA objective values between many (potentially correlated) subproblems across 100 replicates. Hence, the results in Linderoth et al. (2006) yield a measure of the quality of a single lower-bound estimate, while our results compare the quality of several different lower bound estimates.
All data used in this analysis -specifically, the objective values obtained by solving each LP corresponding to each batch -can be obtained at the web site for this paper at http://pages.cs.wisc.edu/ conghan/slhd/. Table 5 shows the wall-clock time required to solve a single sampled approximation for the five test problems using our timing setup as described in § 5.1. Even though the time to generate a batch of scenarios increases as the sampling method increases in complexity, this time is insignificant compared to the rest of the process and hence was not included in this table. The timings were obtained for each problem and number of scenarios by averaging across 16 subproblems constructed from all design methods we considered. The timing increases with the number of scenarios at a superlinear rate. Hence, when the underlying problem is difficult and the number of scenarios is sufficiently high, solving many batches could be extremely time consuming even with multiple machines.

Mean and Standard Error of the Different Sampling Methods
In our first set of experiments, we computed the means and standard errors of the lower bound estimates of the objective value over 100 replicates of each sampling method and each number of scenarios. We fixed the number of subproblems at t = 16, while varying the number of scenarios per batch according to n ∈ {128, 256, 512, 1024}. We tested four different sampling methods: Monte Carlo (MC), independent ordinary Latin hypercube (ILH), Sliced Latin hypercube (SLH) and Sliced Partial Orthogonal Array Based Latin hypercube based on Bush orthogonal arrays (BUSH), which are OA(n 2 , m + 1, n, 2). Since we are only using 16 out of n possible batches, we are discarding much of the orthogonal array, and therefore not achieving much two-dimensional stratification in our negatively dependent designs. Computational results are summarized in Table 6. We begin with some general observations about the mean of the lower bound estimates. As observed in other experiments (Linderoth et al. 2006, Freimer et al. 2012, the Monte Carlo method is significantly worse than Latin Hypercube-based methods in terms of the bias. As expected, the Latin Hypercube-based methods produce statistically-indistinguishable lower bounds. By comparing with the results of Table 4, for all problems except ssn, ILH attains a mean extremely close to the true mean when the number of scenarios is 1024, and is already very close with a smaller number of scenarios. ssn is known to be a challenging problem, requiring at least 512 scenarios to attain a reasonable estimate of the optimum even for the three Latin hypercube schemes. Increasing the number of scenarios beyond 1024 would continue to improve the quality of the estimates. Linderoth et al. (2006) provide a more detailed description of the behavior of SAA for the ssn problem.
We now turn to the standard error. By this measure, for every problem except ssn, we see that MC performs much worse than the other methods. It is only slightly worse for ssn. This table shows situations in which the negatively dependent designs of SLH and BUSH begin to distinguish themselves from ILH. For gbd and storm, we can see large improvements from ILH to SLH/BUSH. The improvement from ILH to SLH/BUSH is much less pronounced in ssn, being in general smaller than the improvement from MC to ILH, and in one case performing slightly worse than ILH. However, we should note that in all the problems, BUSH and SLH have roughly similar performance. This observation suggests that when we use too few subproblems, the effect of partial orthogonality on performance is not significant.
For 20term and LandS, the Latin hypercube schemes ILH, SLH, and BUSH perform similarly. In the case of 20term, this similarity is unsurprising. The benefits of SLH over ILH come from the increased stratification that comes from the sliced structure, but since each variable can only take on two values, each with probability 0.5, any stratification that divides the probability space into a multiple of two would perform equally well. In the case of LandS, we suspect that the similarity of performance is due to the smoothness of the distribution of the random variables. In fact, the cumulative distribution function is essentially linear. We have found that when we modify the distribution to be more irregular or to to be a uniform distribution over a much smaller set, it tends to drive up the standard error and to cause SLH to have a significantly smaller standard error than ILH.
We now consider a greater number of subproblems for each n and new alternatives for the underlying orthogonal arrays. In addition to the four different sampling methods considered earlier, we show two additional variants: Sliced Orthogonal Array Based Latin hypercube based on Bose-Bush orthogonal arrays (BB), and independent batches taken over several BB (INDBB). Bose-Bush orthogonal arrays have an OA(λs 2 , m + 1, s, 2) structure. We pick s to be equal to the number of subproblems, and define λ = n/s. With these choices, the sliced designs achieve full two-dimensional stratification, making BB an example of SOLH sampling.
We include INDBB in our experiments to help isolate the factors that lead to the stronger performance of BB. If each slice of the BB design has some special structure that leads to improved performance, then INDBB designs should perform better than ILH, and the performance of INDBB should be close to that of BB. However, if the performance gains of BB are primarily due to the better two-dimensional stratification, then we would expect INDBB to perform no better than ILH. The numerical results support the second claim.
Results are shown in Table 7. Many of the observations about MC/ILH/SLH/BUSH from Table 6 carry over. Also, comparing the standard error of MC/ILH/SLH/BUSH between the two tables, we notice a factor of √ 2 or 2 difference in standard error, depending on whether the number of subproblems was doubled or quadrupled. This factor is consistent with Theorem 3.4. The lower bound estimates are roughly the same in both tables, demonstrating that changing the number of subproblems does not affect the bias. Table 7 shows a considerable advantage for BB over ILH. In fact, BB is better than all other methods tested, except on gbd, where it performs similarly to SLH/BUSH. On LandS, BB performs about 5-10× better than the other sliced sampling methods. A possible explanation for this huge improvement is that the total number of random variables is just three, so having two-dimensional stratification would cover a large portion of the possible interactions between variables. In the case of ssn, the improvement from ILH to BB is comparable to the improvement from MC to ILH, a bigger factor than is observed for any other problem.
Finally, we turn our attention to running times. The standard error of the estimates between the results for 512 scenarios and 16 slices in Table 6 and 256 scenarios and 32 slices in Table 7 are similar. Since the timing scales superlinearly in the number of scenarios ( § 6.1), the amount of time it takes to solve ct sampled approximations of n scenarios sequentially can be substantially less than solving t sampled approximations of cn. Each batch could also be solved independently and in parallel. This suggests that if computing resources on each machine is limited and using SLH/SPOLH with cn scenarios and t batches is computationally infeasible, using SOLH with with n scenarios and ct batches can be an effective way of reducing the standard error. We conclude that for a fairly small number of subproblems, the sliced sampling methods perform at least as well as Latin hypercube sampling, and in fact show significant improvement in some cases. Once we increase the number of subproblems and exploit the full "orthogonality" property of the orthogonal arrays, we see a substantial improvement in all cases. Thus, if the computational budget will only allow a small number of batches for the given value of n for which a lower bound v n is being estimated, there is significant computational benefit to using the more sophisticated sampling methods introduced in this work.

Conclusions and Future Work
In this paper, we propose the use of two types of negatively dependent designs to improve the lower bound of the objective value. Sliced Latin hypercube sampling is easy to implement since SLH does not impose any restriction on the number of batches t and the number of scenarios in each batch. We introduce the concept of monotonicity for two-stage stochastic linear programs, and we provide a non-asymptotic result showing that SLH can be better than ILH for problems with this monotonicity property. On the other hand, we show that SLH is asymptotically equivalent to ILH if the distribution of the random vector has finite support and the approximate solutions converge. Our computational results supports the theory, showing that SLH performs no worse than ILH and in some cases performs significantly better than ILH.
To improve upon SLH, we consider sliced orthogonal array-based Latin hypercube sampling schemes, which achieve stronger negative dependence between batches. The choice of the underlying orthogonal array can make a huge difference in variance reduction. We provide empirical results showing that when we are able to exploit the full orthogonality of the underlying orthogonal array, using Bose-Bush orthogonal arrays (Bose and Bush 1952), the performance is significantly better than when we use only part of an orthogonal array.
Our work treats Latin hypercube sampling as the baseline method, in part because it was investigated in earlier work (Linderoth et al. 2006). Other sampling methods, such as U sampling (Tang 1993, Tang andQian 2010) and randomized quasi-Monte Carlo (Niederreiter 1992, Owen 1995, Homem-de Mello 2008 could have been used instead for constructing a single SAA problem. We can carry over the idea of negatively dependent designs to these advanced within-batch sampling techniques. For U sampling, a strength-two orthogonal array-based Latin hypercube design could be generated for each SAA problem. The t underlying strength-two orthogonal arrays can be obtained by slicing a larger strength-three orthogonal array. For randomized quasi-Monte Carlo, we can sample different batches based on the same low-discrepancy sequence such that batches are negatively dependent spontaneously. Comparing with Latin hypercube sampling, both U sampling and randomized quasi-Monte Carlo are extremely complicated to implement as they require more constraints on the selection of batch size n and the number of batches t. A potential research direction in the future for us is to study the theoretical and empirical performance of negatively dependent batches based on these advanced sampling methods.