Periodic version of the minimax distance criterion for Monte Carlo integration

Periodic version of the minimax distance criterion for Monte Carlo integration ELIÁŠ, J.; VOŘECHOVSKÝ, M.; SADÍLEK, V. Advances in Engineering Software 2020, vol. 149, November 2020, pp. 1−13 ISSN: 0965-9978 DOI: https://doi.org/10.1016/j.advengsoft.2020.102900 Accepted manuscript © 2020. This manuscript version is made available under the CC-BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-nd/4.0/ dspace.vutbr.cz Periodic version of the minimax distance criterion for Monte Carlo integration Jan Eli´aˇsa,∗, Miroslav Voˇrechovsk´ya, V´aclav Sad´ıleka aInstitute of Structural Mechanics, Faculty of Civil Engineering, Brno University of Technology, Czech Republic Abstract The selection of points for numerical integration of the Monte Carlo type, largely used in analysis of engineering problems, is developed. It is achieved by modification of the metric in the minimax optimality criterion. The standard minimax criterion ensures the design exhibits good space-filling property and therefore reduces the variance of the estimator of the integral. We, however, show that the points are not selected with the same probability over the space of sampling probabilities: some regions are overor under-sampled when designs are generated repetitively. This violation of statistical uniformity may lead to systematically biased integral estimators. We propose that periodic metric be considered for calculation of the minimax distance. Such periodic minimax criterion provides statistically uniform designs leading to unbiased integration results and also low estimator variance due to retained space-filling property. These conclusions are demonstrated by examples integrating analytical functions. The designs are constructed by two different algorithms: (i) a new time-stepping algorithm resembling a damped system of attracted particles developed here, and (ii) the heuristic swapping of coordinates. The designs constructed by the time-stepping algorithm are attached to the paper as a supplementary material. The computer code for construction of the designs is attached too. Keywords: numerical integration, design of experiments, space-filling designs, Latin Hypercube Sampling, Vorono¨ı tessellation, periodic space 1. Introduction Many engineering applications are inevitably accompanied by random nature of input variables. The randomness, arising from environment, material structure, construction process, maintenance or many other sources, needs to be properly taken into account. The resulting output of the engineering analysis is a random function and we are interested in questions which can be subsumed under the analysis of reliability, uncertainty and sensitivity. Answers to these questions are obtained by numerical integration in high dimensional parameter spaces. Such task is extremely time consuming and often becomes a bottleneck of the analysis. We focus on this numerical integration of functions of random vectors. Let us assume a real-valued function g(X) of a random vector X= (X1, X2, . . . , Xs) of selements (we assume the function is integrable and repeatable). This function can represent simple analytical formula that provides stress in elastic bent beam, a complex numerical finite element solution of a large bridge under severe condition, performance of a system or any other engineering problem. The statistical characteristics of the output are defined via s-dimensional integration E[S[g(X)]] = ∞ Z −∞ ··· ∞ Z −∞ S[g(x)]fX(x) dx1. . . dxs(1) where fXis the joint probability density function of the input random vector and the statistical characteristic of interest is determined by function S[·]. The mean value is obtained by taking S[g(X)] = g(X), polynomials of g(X) yield higher statistical moments provided the integrals exist. The integration is typically evaluated numerically as the weighted summation of integrand values at ndiscrete points. To decrease extreme computational demands, one can employ some type of surrogate model [3] or adaptive technique [23, 41]. We limit ourselves to the basic case where integration points are generated a-priori and their weights are equal to 1/n (the points are supposedly selected uniformly with respect to fX), the numerical integration can be written using the following point estimator (average) E[S[g(X)]] ≈1 n n X i=1 SgF−1 X(ui) (2) ∗Corresponding author Email address: [email protected] (Jan Eli´aˇs) Preprint submitted to Advances in Engineering Software September 14, 2020 with FXbeing the joint cumulative distribution function of vector X, and uirepresenting realizations of s independent variables with a uniform distribution spanning interval the from 0 to 1. The realizations of the original input vector xiare obtained in Eq. (2) with the help of the inverse transformation xi=F−1 X(ui). Therefore, the realizations uiare drawn from a unit hypercube U= [0,1]srepresenting the probability space. Other types of design domains are also sometimes required [11], but these are not exploited here. The selection of integration points ui(called design points here) can be performed deterministically (via various forms of grids or quasi Monte Carlo sequences) or with some degree of randomness. In the latter case, which is explored in this paper, the result of Eq. (2) is a random variable. Its mean value should coincide with the exact solution, i.e., Eq. (2) should provide an unbiased estimator of Eq. (1). Simultaneously, the variance of the estimator should be as low as possible. These characteristics are dictated to a large extent by process in which the realizations of uiare generated. The famous Koksma-Hlawka inequality [22, 29] postulates an upper limit for the integration error as a multiplication of the discrepancy between the set of design points ui(called a design Dhereinafter) and a term dependent solely on the integrated function. If we have a-priori knowledge about the function, we may optimize the location of points based on that knowledge [41]. However, we focus here on a situation in which the model-free (“universal”) design is created. In such a case, the Koksma-Hlawka [22, 17] inequality motivates us to use points that have low discrepancy – this property is referred to as sample uniformity. The research area known as design of experiments offers several criteria that should provide such designs based on a particular measure of discrepancy [29, 10, 40]. There is a group of point sets in the literature known as quasi Monte Carlo (QMC) sequences [14, 31, 29, 15, 45]. These sequences are designed to achieve fast decrease of discrepancy with an increasing number of points in the domain. QMC sequences are generally assumed to be the most efficient tool for numerical integration of smooth functions [44, 21, 8]. Similarly, the so-called low-dispersion sequences [30, 32] play a similar role in QMC methods for global optimization. Another group of design criteria delivering designs with high sample uniformity measure space-filling property [1, 47, 46]. This group comprises several distance-based criteria among which the most famous are maximin and minimax [19], Audze-Eglajs [1] and φp(known as phi criterion) [28]. For an overview of the design of experiments please see, e.g., the recent review paper [12]. Designs optimized with respect to discrepancy or space-filling criteria indeed exhibit reduced estimator variability. However, this is not a sufficient condition for unbiased Monte Carlo integration. Eq. (2) is derived from the assumption that every possible domain point has the same probability of being included in the design. This statistical uniformity is usually violated in standard space-filling criteria, as has been shown by authors for Audze-Eglajs in [9] and for maximin and φpcriteria in [49]. Voˇrechovsk´y and Eli´aˇs [49] also show the statistical non-uniformity of the widely used centered L2discrepancy [16] and recently proposed support points [24]. A simple remedy which is and easy to implement was also proposed in [9, 49] for distance-based criteria to ensure their statistical uniformity via the usage of periodic space. The same idea is expanded here for the minimax criterion and conceptually for its relaxed version. Besides designs from U, we will also consider Latin hypercube (LH) designs [6, 27] in the unit hypercube, i.e., the location of points will be limited to the finite set L∈U,L=(i−1/2)/n sfor i∈ {1,2, . . . , n}, and each coordinate value must be used exactly once for each dimension. The fundamental concepts of this paper can be applied to LH designs from the finite set Las well as to designs from the infinite set U, but LH designs are particularly attractive because they are non-collapsing and also represent each input variable in an optimal manner. 2. The standard minimax criterion φmM Let us consider a design Dwith npoints ui. The minimax optimal design, DmM, minimizes the maximum distance from an arbitrary point u∈Uto the nearest design point ui[19]. To find such a design, one needs to search for D= (u1,u2, . . . , un) so that the minimax measure φmM is minimized φmM(D) = max u∈Uhmin i=1...n d(u,ui)i(3) A more appropriate name for the criterion would therefore be the minimaximin criterion, as noted by Pronzato and M¨uller [35]. The function dcan be an arbitrary metric; typically Minkowski distance is used dr(ui,uj) = s X v=1 [∆v(ui,uj)]r!1/r (4a) where ∆v(ui,uj) = |ui v−uj v|and r≥1 (4b) with r= 1 (Manhattan distance), r= 2 (Euclidean distance) or r=∞(Chebyshev distance) [47]. Only Euclidean distance is considered further in this paper, as the application of other powers in Minkowski distance 2 Figure 1: Optimal designs with respect to the minimax criterion (top row, designs taken from https://spacefillingdesigns.nl) and the proposed periodic minimax criterion (bottom row) for s= 2 and n=7, 10 and 13. The red circles with centers at the design points have a radius equal to the φmM(D) or φmM(D) value of that particular design and demonstrate the design’s quality – the whole domain must be covered by the circles but overlapping should be minimized. is straightforward. Euclidean distance is isotropic, i.e., the distance between two points does not depend on the orientation of the coordinate axes. However, its usage in a higher dimensional hypercube is problematic due to the concentration of distances [39]. The search for DmM is also known as the facility location problem or the set covering location problem [43]. It has a number of practical applications such as the creation of an ozone monitoring network [37]. The seminal paper by Johnson et al. [19] showed that the minimax design is G-optimal for certain correlation structures, i.e., it minimizes the maximum Kriging variance. It is therefore often used to find designs optimal for the construction of emulators. Eq. (3) can be relaxed analogously to the way the maximin criterion was relaxed to formulate the φpcriterion. The relaxation replaces strict requirement of Eq. (3) to minimize only the maximum distance to minimization of an integral of distances raised to some power over the whole domain. A family of new criteria is produced by the relaxation [37] φmM[p,q](D) =   Z U"n X i=1 d(u,ui)−q#−p/q dU   1/p (5) with p, q > 0. The minimax criterion is recovered for p, q → ∞. Examples of minimax optimal LH designs in two dimensions with Euclidean distance for n= 7, 10 and 13 are shown in the top row of Fig. 1. These designs were taken from the database https://spacefillingdesigns.nl [47]. The sample uniformity (or space-filling property) guaranteed by the design criterion is achieved. However, it will be shown below that the sample uniformity is not maintained and therefore these designs result in biased integration. 3. The periodic minimax criterion φmM Inspired by our work on the Audze-Eglajs, maximin and φpcriteria [9, 49], we propose a periodic version of the minimax criterion. We start with a periodic Euclidean metric d d(ui,uj) = v u u t s X v=1 ∆v(ui,uj)2(6a) where ∆v(ui,uj) = min {∆v(ui,uj),1−∆v(ui,uj)}(6b) The periodic metric measures distance in the periodic space defined by the unit hypercube U. It considers all the periodic images of points at coordinates u+k, with kbeing a vector of arbitrary integers (Z) of size s, 3 Figure 2: Visualization of periodic Vorono¨ı tessellation on a unit square and its folding into a torus, n=50. but only the shortest distance among all the images is taken into account. This type of periodic space is often referred to as the minimum image convention. Plugging the periodic metric into the definition of the minimax distance criterion from Eq. (3) provides its periodic version φmM(D) = max u∈Uhmin i=1...n d(u,ui)i(7) which is called the periodic minimax criterion hereinafter. Examples of periodic minimax optimal LH designs are shown in the bottom row of Fig. 1, including the nearest periodic images of all the design points. One can see that the space-filling property (sample uniformity) of the design is preserved. To prove that the statistical uniformity is ensured, imagine an arbitrary shift of the whole design Dby a vector lof arbitrary real numbers (without a loss of generality limited to interval h0,1i) of size s. Whenever a point leaves the hypercube, it enters it again from the opposite side due to the periodicity. Thus the individual coordinates of the new points will be either uv+lvor uv+lv−1. For all possible variants, the quantity ∆ from Eq. (6b) is always preserved for every dimension, implying that the periodic distance between any two points is preserved as well. Therefore, φmM(D) is invariant with respect to an arbitrary shift. When a design is φmM optimal, so are all its shifted versions, and all of them have the same probability of being found in the optimization process (assuming the optimization process itself is unbiased). Therefore, all the locations share the same probability of being selected and the design is statistically uniform. If the metric imposed by Eq. (6) were used in the relaxed version of the minimax criterion from Eq. (5) then the periodic relaxed version would be obtained. It would only be approximative, since only the shortest distances are considered – not all periodic images. The error would be greater for low powers q, for which the ignored distances to other periodic images of uicontribute significantly. In any case, the periodic minimax distance from Eq. (7) is exact, since only the shortest distances contribute. 4. Evaluation of φmM and φmM The evaluation of the minimax criterion (Eq. 3) is relatively computationally expensive. If one has a design Dwith design points ui, one needs to solve the maximization problem by finding the location u∗∈Uwith the largest distance from the design points. In the literature one can find three approaches which can be used to determine the location of u∗:  Typically, the problem is approximated by limiting the candidates for u∗to a finite set of points chosen artificially (usually a regular grid). This is the most common method that allows φmM(D) to be estimated in reasonable time. However, to keep the error of such an estimate low, the grid needs to be dense and the number of candidate points in the grid grows fast with dimension s. The modification of this method for periodic space is easy: one can use the same grid and the same evaluation method; it suffices to replace the metric dwith the periodic version ¯ d.  Another option for the evaluation of φmM(D) is via what is known as Markov chain Monte Carlo optimization. It is a sequential construction of points placed randomly but at a decreasing distance from design points ui. As the distance decreases, the φmM value is estimated with diminishing error. The algorithm is thoroughly described in [13]. Pronzato [34] shows that this method is the most effective for higher dimensions (tested up to s= 5) as the computational time increases linearly with s. The periodic version is again straightforward; as before, it suffices to change the metric to ¯ d.  The last method uses Delaunay triangulation and/or Vorono¨ı tessellation. These dual methods essentially directly identify candidates for u∗. In the case of Delaunay triangulation, the φmM(D) is the radius of the largest hypersphere that is circumscribed to Delaunay simplices, while in the case of Vorono¨ı tessellation, the φmM(D) is the maximum distance of the Vorono¨ı vertex to the design points ui. This method yields 4 Figure 3: Clipped and periodic Vorono¨ı tessellations for s= 2 and n=16. exact results as the set of candidates for u∗determined by the tessellation or triangulation contains the true solution. Unfortunately, the computational complexity of Delaunay triangulation and Vorono¨ı tessellation becomes prohibitive for higher dimensions. In addition, the Vorono¨ı tessellation can be constructed only approximately, as is effectively done in [26]. However, the method then becomes essentially the same as the first method on this list. One also needs to properly analyze the boundary region of the tessellation/triangulation. Both the Delaunay triangulation and Vorono¨ı tessellation approaches do not provide candidates for u∗at the boundary of the hypercube, so they must be found by additional geometrical calculations [7]. As the boundary region occupies more volume for high s, such additional calculations may take a significant time. One can easily get rid of the problems with boundaries by extending npoints by their symmetric images with respect to all hypercube boundaries [35], i.e., by calculating triangulation on n(2s+ 1) points. Even if not all n symmetric images are considered, the computational complexity still grows fast with s. The modification for the periodic minimax distance is in this case more elaborative since the boundaries must be treated in a completely different, periodic way. One can imagine it as computing Vorono¨ı tessellation or Delaunay triangulation in a folded space (anuloid), see Fig. 2 or Fig. 3. The next section elucidates the developed algorithm for such a construction. Since the third method is the only exact one, we decided to employ it for the evaluation of both the minimax (φmM) and the periodic minimax (φmM) criteria. The Vorono¨ı tessellation with auxiliary points placed symmetrically (mirrored) with respect to all the boundaries (needed for the calculation of φmM) is often called clipped in the computational geometry literature [54]. There are also papers computing the periodic Vorono¨ı tessellation needed for the calculation of φmM. However, these algorithms are usually focused on centroidal Vorono¨ı tessellation in s= 2 or s= 3 dimensions and complex, often nonconvex regions as the application is the generation of meshes [53, 38, 42]. In our case, we would like to compute clipped and periodic tessellation in higher dimensions, but our region is a simple convex hypercube. For that reason, we have developed our own algorithm employing Q-hull code [2]. 4.1. Incrementally constructed clipped Vorono¨ı diagram The clipped version starts by placing auxiliary outside points symmetrically with respect to all the boundaries. One can mirror all the ninner points to make sure the correct solution will be obtained. However the calculation of the tessellation is then performed with n+2sn points, which can be unnecessarily time consuming. We therefore mirror only the N0points that are closest to each mirroring plane. The solution is then computed on n+ 2sN0points only, which can be much faster, though it may happen that the result is not fully correct. The incorrect solution is enriched by the subsequent mirroring of additional points in a loop until the correct solution is found. The correctness of the solution and additional auxiliary points are determined via Delaunay simplices. All the simplices that connect inner (original) and outer (auxiliary) points are selected. It is required that each inner point of such a simplex is symmetrically placed with respect to all the boundaries crossed by that simplex. If such a point is missing, it is added in the next step. If no auxiliary point is missing, we have the correct clipped Vorono¨ı tessellation. An example of this procedure is shown in Fig. 4 left with n= 16 points and s= 2 dimensions. The gray color marks all the simplices that cross boundaries but do not have all of their inner points mirrored. Four new auxiliary points are identified and the triangulation/tessellation is updated. The tessellation in the second step is already the correct solution. Based on trial–error analyses, the reasonable N0was estimated as N0= 3n`char, where the characteristic length `char =n−1/s represents the shortest distance in the regular orthogonal grid [39]. The procedure usually terminates in the second step, but sometimes more steps are needed. Though Q-hull software allows the incremental extension of sets of points and the updating of the solution, it has been found to be much quicker to delete the previous solution and compute the new one from scratch. 5 Figure 4: Incremental construction of clipped and periodic Vorono¨ı tessellations and Delaunay triangulation for s= 2 and n=16. 100101102103104 n 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 time fraction ti/tf s= 2 100101102103104 n s= 3 100101102103104 n s= 4 100101102103 n s= 5 clipped, N0= 1 clipped, N0= 3n`char periodic, N0= 1 periodic, N0= 3n`char 100101102103104 n 10−3 10−1 101 103 time tf[s] s= 5, clipped s= 5, periodic s= 2, cliped & periodic s= 5, clipped s= 5, periodic s= 2, cliped & periodic Figure 5: Performance of the incremental construction of clipped and periodic Vorono¨ı tessellations using N0= 1 and 3n`char initially reflected or periodically repeated points. The four graphs on the left report the computational time of the incremental solution, ti, divided by the computational time of the full solution, tf, measured when all the inner points are reflected or periodically repeated in all necessary directions. The right-hand graph shows real time tffor s= 2 and 5. The Vorono¨ı constructions for each case were repeated 50 times with different locations of design points. Solid lines with circular markers show averages based on those 50 repetitions, while shaded bands show averages ±standard deviations. 4.2. Incrementally constructed periodic Vorono¨ı diagram For periodic tessellation, the mirroring is replaced by periodic repetition of the N0points that are closest to the boundaries. Unfortunately, one also needs to add auxiliary periodic images of points in the “corner regions”, i.e., in all regions where two or more coordinates are above 1 or below 0. This is ignored in the first step, however, and the first Vorono¨ı tessellation is found for n(inner) plus 2sN0(outer) points. The solution must again be checked for correctness. Whenever any simplex crosses a boundary hyperplane, it is required that all its inner points be periodically repeated with respect to the crossed boundaries. Otherwise, the missing points are added in the next step. The situation is shown in Fig. 4 right, where some (not all) of these simplices are colored in gray. Simplices band csimply cross one boundary, and thus must be periodically repeated with respect to this boundary. Simplices a,dand erequire periodic repetition with respect to two boundaries. In the third step we find the correct solution. Once again, the algorithm seems to work well for N0= 3n`char, and the fastest Q-hull implementation compute a completely new solution in every step. Usually more than two steps are needed. A comparison of approximate computation times for the incremental algorithm (ti) and full solution (tf) in several dimensions is shown in Fig. 5. The full solution assumes the mirroring or periodic repetition of all the points directly (a total of n(2s+ 1) points for the clipped and n3spoints for the periodic tessellation). The results show that the full solution is only faster for extremely low n. 6 Figure 6: Illustration of the analogies between distance-based criteria and a system of particles. a) Mm step. b) mM step. c) Decision between the two types of steps in the initial rearrangement. 5. Construction of minimax optimal designs The search for an optimal minimax design is an extremely difficult problem. Even if we restrict ourselves to the LH sampling, the finite set of options has size (n!)s−1/(s−1)! (if we ignore permutation of coordinates). The exhaustive search method quickly becomes prohibitively time-consuming, even for s= 2. Instead, heuristic optimization methods that, unfortunately, often deliver only the local minimum are being used. Some possible optimization methods are listed in [34]. In this paper, one existing and one new method are employed. 5.1. Shuffling of coordinates based on simulated annealing In our previous papers [9, 49] focused on different design criteria we used a heuristic optimization method based on simulated annealing [20, 5]. A thorough description of the algorithm can be found in [50]. Starting with a random LH design, new trial designs are generated by swapping a randomly chosen coordinate vof a randomly chosen pair of points aand b:ua vub v. The swap is accepted based on the difference in the minimax distance value of the original and new designs ∆φmM =φmM(Di+1)−φmM(Di) with a probability based on the Boltzmann distribution P=   1 ∆φmM ≤0 1 1 + exp (∆φmM/t)∆φmM >0(8) For the periodic version, the criterion φmM is used. Temperature tgoverns the process and decreases throughout it from the initial to the final temperature in steps. At each temperature step, a user defined number of trial swaps is performed. Simulated annealing is a heuristic optimization technique in which φmM or φmM are typically calculated hundreds of thousands of times. In the case of the minimax or periodic minimax criteria, the computational burden associated with so many evaluations is huge, and the temperature evolution and number of trial swaps at each temperature step (known as the cooling schedule) needs to be drastically modified to reduce the number of evaluations of the criteria. The price paid for such a reduction is that only near-optimal solutions are found. 5.2. Time-stepping algorithm for the construction of distance-based designs The insufficient performance of the switching algorithm for the minimax criterion resulted in the development of a completely new optimization method. Voˇrechovsk´y et al. [51] argued that the distance-based criteria for a design consisting of npoints can be viewed in analogy with a system of nparticles in s-dimensional space: the particles, each having scoordinates, represent the design points. Various formulas for the optimality criteria involve terms for all pairs of points, and these terms can represent the pairwise potential energy. Given this analogy, designs can be constructed by minimizing the total potential energy of a system of particles via the dynamical simulation of an N-body system [25]. As argued in [51], the (periodic) minimax criterion provides potential that changes abruptly and is therefore not suitable for a fully dynamic solution. Instead, this paper proposes a simplified version of the dynamical algorithm that neglects particle inertia and is well suited for the maximin and minimax criteria. It can employ the intersite and periodic metric equally well. At the beginning, a random arrangement of particles is generated. The simulation proceeds in time steps in which only selected particles iare displaced from their current positions: ui(t+ ∆t) = ui(t) + si(t). The magnitude of the displacement (the Euclidean length of the steps s(t)) decreases with increasing step number (time t). When the miniMax (mM) criterion is to be minimized, up to (s+1) particles lying on the surface of the largest empty hypersphere (forming the critical simplex) are displaced towards the center of this hypersphere. In 7 such a new configuration, the criterion is recalculated by identifying the current critical simplex and displacing its points. The displacements depend only on the current positions; no information from the previous step is utilized (one can think of neglected inertia in a dynamical system of particles). Such a displacement is called the mM-step and is illustrated in Fig. 6b. Similarly, when the maximin (Mm) criterion is targeted, the individual steps proceed by selecting the nearest pair of particles and moving the two particles in opposite directions (see Fig. 6a for an illustration). Experience shows that when an algorithm solely employs a series of mM-steps, particle clusters typically appear as there are no repelling forces between the particles. Analogously, when the algorithm involves a sequence of Mm-steps only, clusters are eliminated but large empty regions may remain unfilled. Therefore, before imposing the target criterion, we recommend always performing initial rearrangement by alternating both types of steps. The type of step to be performed is selected based on the comparison of two current length parameters: dmin and dmax. The definitions of these lengths are: dmin =|`char −φMm|and dmax =|c φmM −`char|, where φMm is the Maximin criterion, i.e., the smallest pairwise distance, φmM is the radius of the largest empty hypersphere and c=p2(s+ 1)/s is the ratio between the circumradius and the edge length of a regular simplex (for the periodic version replace φmM and φMm with φmM and φMm). When dmin > dmax, the spreading of the cluster is applied (Mm-step), otherwise the shrinking mM-step is taken. Once the initial rearrangement is finished, the algorithm can proceed with a series of steps of a single type depending on the target criterion. In order to increase the robustness of the algorithm, we recommend always applying a small random perturbation to each point in a random direction after performing a certain number of steps. The magnitudes of these perturbations also decrease with the number of steps (analogously to the decreasing excitations in annealing due to decreasing temperature). Even though these perturbations usually increase the value of the criterion, they can help to prevent locking in a local minimum. Once the final configuration is reached, the coordinates of the points are taken as the design. It is worth noting that the periodic version of the time-stepping algorithm typically finds the solution faster than the intersite version. The intersite version is constrained by boundaries and more point displacements are needed to satisfy these constraints. In contrast, the periodic version is unconstrained and progressively develops towards the “closest” of all the shifts of the optimal pattern with the lowest minimax criterion. If LH design is desired, a postprocessing step must be taken in order to distribute the coordinates uniformly along each dimension. The predefined coordinates are selected as midpoints of nsegments of the unit length: (πi−1/2)/n, where the integers πiare the rank numbers of points along the current coordinate. We remark that the latinized sample obtained with “free points” may not necessarily be the best possible LH sample given the φp criterion. However, the application of latinization via post-processing in this manner is much more efficient than the combinatorial optimization that shuffles predefined coordinates based on the simulated annealing algorithm described in Sec. 5.1. 5.3. Comparison of SA and TS algorithm performance The designs generated by the simulated annealing (SA) and time-stepping (TS) algorithms are compared from the viewpoint of the following criteria: the (i) standard minimax, (ii) periodic minimax, (iii) periodic phi criterion with power p=s+ 1 and (iv) wrap around discrepancy. We intend to show via such comparison the versatility of the proposed periodic minimax criterion: it provides designs with relatively good quality also with respect to criteria for which they were not directly optimized. We also intend to compare performance of the φmM designs with random designs or designs optimized by other criteria from the literature. The designs under comparison are:  random LH designs as a reference solution.  a quasi Monte Carlo (QMC) sequence used widely for numerical integration. From the pool of available QMCs we have selected the Sobol’ sequence. The sequence is generated from the R package “fOptions” [36, 52] using the “runif.sobol” function.  optimal minimax and periodic minimax designs obtained by simulated annealing (denoted φmM LH SA and ¯ φmM LH SA, respectively). These designs are optimized in 2 and 3 dimensions only due to their heavy computational demands, and they are limited to the Lspace.  optimal minimax and periodic minimax designs obtained by a time-stepping algorithm without latinization (denoted φmM TS and ¯ φmM TS, respectively). The effectiveness of the time-stepping algorithm allows optimization even in dimension 5, but we only studied it up to n= 512 points.  latinized designs from the previous item (denoted φmM LH TS and ¯ φmM LH TS, respectively).  optimal periodic phi LH designs (labeled ¯ φs+1 LH) obtained by simulated annealing optimization with the distance exponent p=s+ 1. These designs were presented in [49] and are included here for comparative purposes. 8 12 3 4 5 6 7 8 9 10 11 12 mode number 0.00 0.01 0.02 95th percentile εk 0.000 0.005 0.010 0.015 median error εk s= 3, n= 1024 random LH Sobol’ scram. Sobol’ ¯ φs+1 LH ¯ φmM LH TS 1234567 8 910111213 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 mode number 0.00 0.01 0.02 0.03 95th percentile εk 0.000 0.005 0.010 0.015 0.020 median error εk s= 2, n= 1024 Figure 11: Integration error evaluated using Chebyshev polynomials for s= 2,3 integrated by a design of 1024 points optimized with respect to various criteria. We consider K=ds √neto be the maximum polynomial order for which integration error is measured. For higher K, polynomials have more changes of signs inside the interval h0,1ithan the total number of points used for integration. Therefore, one can hardly expect good performance for K > ds √ne. The performance of a design should not depend on the order of modes; we would like to see the same performance for arbitrary permutations of k. For that reason, we always integrate all permutations of kand report the statistical data. Only sorted vectors kare therefore used, i.e., k={1,1,2}actually represents all 3 different permutations {1,1,2},{2,1,1}and {1,2,1}. We can also order all the sorted kthat have no item greater than ds √nein the following way: ka≤kb⇔ ∀v:ka v≤kb v. Fig. 11 shows integration performance in dimensions two and three for different mode vectors kusing n= 1024 points. Separate graphs are shown for the median of the error (top graph) and the 95th percentile (bottom graph). We eliminated all the non-periodic designs from the comparison since we had already shown that they result in strong bias. We also omitted the φmM design generated by simulated annealing (¯ φmM LH SA), since for this large number of points shuffling optimization already suffers from an insufficient number of trial swaps. For the sake of clarity, we did not plot latinized φmM designs generated via the time-stepping process (¯ φmM LH TS) because they exhibit behavior almost distinguishable from raw time-stepping (¯ φmM TS) designs. The only difference is at modes where all the kitems are zero except the last one, for example (0,3), (0,8) or (0,0,7). For these single-degree mode vectors, the Chebyshev basis Tkis dependent on a single variable only, and therefore the LH property becomes extremely advantageous. These modes are located in Fig. 11 at the ticks of the horizontal axis labeled by the single nonzero integer in the vector k. At these particular locations, the latinized time-stepping generated φmM designs exhibit significantly lower integration error compared to the raw ones plotted. Let us now comment on the individual curves shown in Fig. 11:  The reference solution is again provided by the random LH designs. They exhibit large error and also large variance. However, the error suddenly drops for the single-degree mode vector thanks to the LH restriction. 15  The QMC Sobol’ sequences integrate with low errors for initial mode vectors, but the integration errors quickly increase and they often even exceed the random LH design error.  Numerical integration is often performed with re-randomized QMC sequences in the literature. For comparison purposed, we included scrambled Sobol’ sequence using two types of scrambling: Owen scrambling [33] and also the Faure-Tezuka [45] type of scrambling; see also [18]. The median performance is more stable, however, 95th percentile of error is relatively large which is a price paid for the re-randomization.  The periodic phi criterion (described in [49]) shows excellent performance. These LH designs have the lowest error from all of those presented, and also their variance seems to be the lowest.  Finally, the φmM designs developed here via time-stepping are just a little bit worse than the phi designs. If the latinized version of these designs were used, we would also obtain diminishing integration errors for the single-degree mode vectors. Along with the phi criterion, the φmM criterion (with a proper optimization method) is capable of delivering robust designs that integrate all the basis functions with low error and low error variability. They are outperformed only for low Chebyshev modes by QMC sequences because these sequences are purposely prepared in order to integrate very smooth functions well [8]. 9. Conclusions The designs optimized with respect to the standard (intersite) minimax distance are statistically nonuniform and thus yield systematically biased estimators when used for numerical Monte Carlo type integration. When the intersite distance is replaced by the periodic metric, the statistical uniformity is guaranteed and unbiased estimators are obtained. The reduction of estimator variance is achieved by excellent space-filling property and discrepancy (and thus also sample uniformity) of the minimax and periodic minimax optimal designs. The numerical integration is tested using three simple analytical functions and thousands of functions generated in the space of multidimensional Chebyshev polynomials. In all the tests, the periodic minimax designs exhibited low integration error. They are sometimes outperformed by quasi Monte Carlo sequences (represented by the Sobol’ sequence), but these sequences often exhibit extremely large errors (greater than those obtained from random LH designs). The excellent and stable integration ability (robustness) of the minimax designs makes them the best choice (along with the phi designs included for comparison). A disadvantage of the (periodic) minimax criterion is its extremely large computational complexity. In an attempt to speed up the calculations, we have developed a method employing the incremental construction of (periodic) Vorono¨ı tessellation. The computational time is however still excessive and prohibits the use of the otherwise reliable shuffling optimization method based on simulated annealing. Instead, a novel time-stepping algorithm has been developed to deliver φmM optimized designs in reasonable time. Even with these algorithms, one can hardly imagine the use of this design criterion in dimensions greater than five. Acknowledgement The authors acknowledge the financial support provided by the Czech Science Foundation under projects Nos. GA19-12197S and GC19-06684J. References [1] P. Audze and Vilnis Egl¯ajs. New approach for planning out of experiments. Problems of Dynamics and Strengths, 35:104–107, 1977. (in Russian). [2] C. Bradford Barber, David P. Dobkin, and Hannu Huhdanpaa. The quickhull algorithm for convex hulls. ACM Trans. Math. Softw., 22(4):469–483, 1996. ISSN 0098-3500. doi: 10.1145/235815.235821. [3] Mohamed Amine Bouhlel, John T. Hwang, Nathalie Bartoli, R´emi Lafage, Joseph Morlier, and Joaquim R.R.A. Martins. A Python surrogate modeling framework with derivatives. Advances in Engineering Software, 135:102662, 2019. ISSN 0965-9978. doi: 10.1016/j.advengsoft.2019.03.005. [4] Evgeny Burnaev, Ivan Panin, and Bruno Sudret. Efficient design of experiments for sensitivity analysis based on polynomial chaos expansions. Annals of Mathematics and Artificial Intelligence, 81(1):187–207, 2017. ISSN 1573-7470. doi: 10.1007/ s10472-017-9542-1. [5] V. ˇ Cern´y. Thermodynamical approach to the traveling salesman problem: An efficient simulation algorithm. Journal of Optimization Theory and Applications, 45(1):41–51, 1985. ISSN 1573-2878. doi: 10.1007/BF00940812. [6] W.J. Conover. On a better method for selecting input variables. unpublished Los Alamos National Laboratories manuscript, reproduced as Appendix A of “Latin Hypercube Sampling and the Propagation of Uncertainty in Analyses of Complex Systems” by J.C. Helton and F.J. Davis, Sandia National Laboratories report SAND2001-0417, 1975. [7] Jorge Cort´es and Francesco Bullo. Nonsmooth coordination and geometric optimization via distributed dynamical systems. SIAM Review, 51(1):163–189, 2009. ISSN 00361445, 10957200. [8] Josef Dick. Higher order scrambled digital nets achieve the optimal rate of the root mean square error for smooth integrands. The Annals of Statistics, 39(3):1372–1398, 2011. ISSN 00905364. doi: 10.1214/11-AOS880. 16 [9] Jan Eli´aˇs and M. Voˇrechovsk´y. Modification of the Audze–Egl¯ajs criterion to achieve a uniform distribution of sampling points. Advances in Engineering Software, 100:82–96, 2016. doi: 10.1016/j.advengsoft.2016.07.004. [10] Kai-Tai Fang and Y. Wang. Number-Theoretic Methods in Statistics. Chapman and Hall/CRC, London ; New York, 1st edition edition, 1993. ISBN 9780412465208. [11] Fabian Fuerle and Johann Sienz. Formulation of the Audze-Eglais uniform latin hypercube design of experiments for constrained design spaces. Advances in Engineering Software, 42(9):680–689, 2011. ISSN 0965-9978. doi: 10.1016/j.advengsoft.2011.05.004. [12] Sushant S. Garud, Iftekhar A. Karimi, and Markus Kraft. Design of computer experiments: A review. Computers & Chemical Engineering, 106:71–95, 2017. ISSN 0098-1354. doi: 10.1016/j.compchemeng.2017.05.010. [13] Arnaud Guyader, Nicolas Hengartner, and Eric Matzner-Løber. Simulation and estimation of extreme quantiles and extreme probabilities. Applied Mathematics & Optimization, 64(2):171–196, 2011. ISSN 1432-0606. doi: 10.1007/s00245-011-9135-z. [14] J.H. Halton. On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals. Numerische Mathematik, 2(1):84–90, 1960. ISSN 0029-599X. doi: 10.1007/BF01386213. [15] Fred J. Hickernell. Lattice rules: How well do they measure up? In Peter Hellekalek and Gerhard Larcher, editors, Random and Quasi-Random Point Sets, chapter Lattice Rules: How Well Do They Measure Up?, pages 109–166. Springer New York, New York, NY, 1998. ISBN 978-1-4612-1702-2. doi: 10.1007/978-1-4612-1702-2 3. [16] Fred J. Hickernell. A generalized discrepancy and quadrature error bound. Mathematics of Computation, 67(221):299–322, 1999. doi: 10.1090/S0025-5718-98-00894-1. [17] Edmund Hlawka. Funktionen von beschr¨ankter Variation in der Theorie der Gleichverteilung. Annali di Matematica Pura ed Applicata, 54(1):325–333, Dec 1961. ISSN 1618-1891. doi: 10.1007/BF02415361. [18] Hee Sun Hong and Fred J. Hickernell. Algorithm 823: Implementing scrambled digital sequences. ACM Transactions on Mathematical Software, 29(2):95–109, June 2003. ISSN 0098-3500. doi: 10.1145/779359.779360. [19] M.E. Johnson, L.M. Moore, and D. Ylvisaker. Minimax and maximin distance designs. Journal of Statistical Planning and Inference, 2(26):131–148, 1990. ISSN 0378-3758. doi: 10.1016/0378-3758(90)90122-B. [20] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. Optimization by simulated annealing. Science, 220(4598):671–680, 1983. ISSN 0036-8075. doi: 10.1126/science.220.4598.671. [21] Ronald Kleiss and Achilleas Lazopoulos. Error in Monte Carlo, quasi-error in Quasi-Monte Carlo. Computer Physics Communications, 175(2):93–115, 2006. ISSN 0010-4655. doi: 10.1016/j.cpc.2006.02.001. [22] J. F. Koksma. Een algemeene stelling uit de theorie der gelijkmatige verdeeling modulo 1. Mathematica B, 11:7–11, 1942/1943. [23] Wang-Sheng Liu, Sai Hung Cheung, and Wen-Jun Cao. An efficient surrogate-aided importance sampling framework for reliability analysis. Advances in Engineering Software, 135:102687, 2019. ISSN 0965-9978. doi: 10.1016/j.advengsoft.2019. 102687. [24] Simon Mak and V. Roshan Joseph. Support points. The Annals of Statistics, 46(6A):2562–2592, 12 2018. doi: 10.1214/ 17-AOS1629. [25] Jan Maˇsek and Miroslav Voˇrechovsk´y. Parallel implementation of hyper-dimensional dynamical particle system on CUDA. Advances in Engineering Software, 125:178–187, 2018. ISSN 0965-9978. doi: 10.1016/j.advengsoft.2018.03.009. [26] Jan Maˇsek and Miroslav Voˇrechovsk´y. Approximation of Vorono¨ı cell attributes using parallel solution. Advances in Engineering Software, 132:7–17, 2019. ISSN 0965-9978. doi: 10.1016/j.advengsoft.2019.03.012. [27] M. D. McKay, W. J. Conover, and R. J. Beckman. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 21:239–245, 1979. doi: 10.1080/00401706.1979.10489755. [28] Max D. Morris and Toby J. Mitchell. Exploratory designs for computational experiments. Journal of Statistical Planning and Inference, 43(3):381–402, 1995. ISSN 0378-3758. doi: 10.1016/0378-3758(94)00035-T. [29] H. Niederreiter. Random Number Generation and Quasi-Monte Carlo Methods. CBMS-NSF Regional Conference Series in Applied Mathematics. Society for Industrial and Applied Mathematics, Philadelphia, Pennsylvania, 1992. ISBN 978-0-898712957. [30] Harald Niederreiter. Point sets and sequences with small discrepancy. Monatshefte f¨ur Mathematik, 104:273–337, 1987. ISSN 0022-314X. doi: 10.1007/BF01294651. [31] Harald Niederreiter. A quasi-monte carlo method for the approximate computation. Journal of Number Theory, 30(1):51–70, 1988. ISSN 0022-314X. doi: 10.1016/0022-314X(88)90025-X. [32] Harald Niederreiter. Low-discrepancy and low-dispersion sequences. Journal of Number Theory, 30(1):51–70, 1988. ISSN 0022-314X. doi: 10.1016/0022-314X(88)90025-X. [33] Art B. Owen. Scrambling sobol’ and niederreiter-xing points. Journal of Complexity, 14(4):466–489, 1998. ISSN 0885-064X. doi: 10.1006/jcom.1998.0487. [34] Luc Pronzato. Minimax and maximin space-filling designs: some properties and methods for construction. Journal de la Societe Fran¸caise de Statistique, 158(1):7–36, 2017. [35] Luc Pronzato and Werner G. M¨uller. Design of computer experiments: space filling and beyond. Statistics and Computing, 22(3):681–701, 2012. ISSN 1573-1375. doi: 10.1007/s11222-011-9242-3. [36] R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2013. URL http://www.R-project.org/. [37] J. Andrew Royle and Doug Nychka. An algorithm for the construction of spatial coverage designs with implementation in SPLUS. Comput. Geosci., 24(5):479–488, 1998. ISSN 0098-3004. doi: 10.1016/S0098-3004(98)00020-X. [38] Chris H. Rycroft. VORO++: A three-dimensional Voronoi cell library in C++. Chaos: An Interdisciplinary Journal of Nonlinear Science, 19(4):041111, 2009. doi: 10.1063/1.3215722. [39] V´aclav Sad´ılek and Miroslav Voˇrechovsk´y. Evaluation of pairwise distances among points forming a regular orthogonal grid in a hypercube. Journal of Civil Engineering and Management, pages 410–423, 2019. ISSN 1392-3730. doi: 10.3846/jcem.2018.5189. [40] Thomas J. Santner, Brian J. Williams, and William I. Notz. The Design and Analysis of Computer Experiments. Springer Series in Statistics. Springer-Verlag New York, 2003. ISBN 978-0-387-95420-2. doi: 10.1007/978-1-4757-3799-8. [41] Michael D. Shields. Adaptive Monte Carlo analysis for strongly nonlinear stochastic systems. Reliability Engineering & System Safety, 175:207–224, 2018. ISSN 0951-8320. doi: 10.1016/j.ress.2018.03.018. [42] Liang Shuai, Xiaohu Guo, and Miao Jin. GPU-based computation of discrete periodic centroidal Voronoi tessellation in hyperbolic space. Computer-Aided Design, 45(2):463–472, 2013. ISSN 0010-4485. doi: 10.1016/j.cad.2012.10.029. [43] Matthias H. Y. Tan. Minimax designs for finite design regions. Technometrics, 55(3):346–358, 2013. doi: 10.1080/00401706. 2013.804439. [44] Stefano Tarantola, William Becker, and Dirk Zeitz. A comparison of two sampling methods for global sensitivity analysis. Computer Physics Communications, 183(5):1061–1072, 2012. ISSN 0010-4655. doi: 10.1016/j.cpc.2011.12.015. [45] Shu Tezuka and Henri Faure. i-binomial scrambling of digital nets and sequences. Journal of Complexity, 19(6):744–757, December 2003. ISSN 0885-064X. doi: 10.1016/S0885-064X(03)00035-9. [46] William Notz Thomas J. Santner, Brian J. Williams. The Design and Analysis of Computer Experiments. Springer Series in 17 Statistics. Springer-Verlag New York, 2003. ISBN 978-0-387-95420-2. doi: 10.1007/978-1-4757-3799-8. [47] Edwin R. van Dam. Two-dimensional minimax Latin hypercube designs. Discrete Applied Mathematics, 156(18):3483–3493, 2008. ISSN 0166-218X. doi: 10.1016/j.dam.2008.02.009. [48] M. Voˇrechovsk´y. Hierarchical refinement of Latin Hypercube Samples. Computer-Aided Civil and Infrastructure Engineering, 30(5):394–411, 2015. ISSN 1467-8667. doi: 10.1111/mice.12088. [49] M. Voˇrechovsk´y and Jan Eli´aˇs. Modification of the maximin and φp(phi) criteria to achieve statistically uniform distribution of sampling points. Technometrics, -(-):–, 2020. doi: 10.1080/00401706.2019.1639550. in press. [50] M. Voˇrechovsk´y and D. Nov´ak. Correlation control in small sample Monte Carlo type simulations I: A Simulated Annealing approach. Probabilistic Engineering Mechanics, 24(3):452–462, 2009. ISSN 0266-8920. doi: 10.1016/j.probengmech.2009.01. 004. [51] Miroslav Voˇrechovsk´y, Jan Maˇsek, and Jan Eli´aˇs. Distance-based optimal sampling in a hypercube: Analogies to N-body systems. Advances in Engineering Software, 137:102709, 2019. ISSN 0951-8320. doi: 10.1016/j.advengsoft.2019.102709. [52] Diethelm Wuertz, Tobias Setz, and Yohan Chalabi. fOptions: Rmetrics – Pricing and Evaluating Basic Options, 2017. R package, version 3.0.1. [53] Dong-Ming Yan, Kai Wang, Bruno L´evy, and Laurent Alonso. Computing 2D Periodic Centroidal Voronoi Tessellation. In 8th International Symposium on Voronoi Diagrams in Science and Engineering - ISVD2011, Qingdao, China, 2011. doi: 10.1109/ISVD.2011.31. [54] Dong-Ming Yan, Wenping Wang, Bruno L´evy, and Yang Liu. Efficient computation of clipped Voronoi diagram for mesh generation. Computer-Aided Design, 45(4):843–852, 2013. ISSN 0010-4485. doi: 10.1016/j.cad.2011.09.004. Geometric Modeling and Processing 2010. 18