- Email: [email protected]
Simplex based space filling designs
Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
Contents lists available at ScienceDirect
Journal of Statistical Planning and Inference journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / j s p i
Simplex based space filling designs ¨ Thomas Muhlenst a¨ dt∗ ,1 Fakulta¨ t Statistik, Technische Universita¨ t Dortmund, 44221 Dortmund, Germany
A R T I C L E
I N F O
A B S T R A C T
Article history: Received 19 June 2008 Received in revised form 2 June 2009 Accepted 4 August 2009 Available online 11 August 2009 Keywords: Triangulation Maximin design Uniform design Computer experiment
Space filling designs are important for deterministic computer experiments. Even a single experiment can be very time consuming and can have many input parameters. Furthermore the underlying function generating the output is often nonlinear. Thus, the computer experiment has to be designed carefully. There exist many design criteria, which can be numerically optimized. Here, a method is developed, which does not need algorithmic optimization. A mesh of nearly regular simplices is constructed and the vertices of the simplices are used as potential design points. The extraction of a design from these meshes is very fast and easy to implement once the underlying mesh has been constructed. The extracted designs are highly competitive regarding the maximin design criterion and it is easy to extract designs for nonstandard design spaces. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Frequently used approaches to obtain a space filling design D with n design points x1 , . . . , xn in k dimensions are maximin designs and uniform designs, see Santner et al. (2003) and Fang et al. (2006), respectively. The maximin approach consists of maximizing the minimum distance between two arbitrary points of the design: Mm(D) := min(xi − xj 2 ). ij
(1)
Johnson et al. (1990)argued that maximin designs have some asymptotic optimality property in a Bayesian context. As a maximin design is not necessarily uniquely defined, Morris and Mitchell (1995) suggest the following criterion: ⎛
p (D) := ⎝
m
⎞(1/p) −p Jj dj ⎠
(2)
j=1
with (d1 < · · · < dm ) being the sorted list of the point to point distances of the design and (J1 , . . . , Jm ) being the frequency list of the distances with m ⱕ ( 2n ). This criterion needs to be minimized numerically. p is only an approximation for a theoretical generalization of the maximin design and the parameter p influences how good the approximation is. For large p the approximation improves but p will tend to have many local minima. Hence, for increasing p it will be more difficult to find a global minimum numerically. Morris and Mitchell (1995) suggest to use p ∈ {20, . . . , 50}, therefore we use p = 40 in this article to have a good approximation. ∗ Tel.: +49 231 7557227; fax: +49 231 7555305. E-mail addresses: [email protected], [email protected] URL: http://www.statistik.tu-dortmund.de. 1 The financial support of the DFG (research training group Statistical Modelling) is gratefully acknowledged. 0378-3758/$ - see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2009.08.002
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
586
Another commonly used approach to create space filling designs are uniform designs which are trying to minimize a test statistic of a test on a U[0, 1]k -distribution. In this context special test statistics have been constructed, among those the wrap around criterion, which is defined by k n n k
4 1 3 (WD(D))2 = − + 2 − |xdi − xji |(1 − |xdi − xji |) . 3 2 n
(3)
d=1 j=1 i=1
See Fang et al. (2006) for a more general description of the wrap around discrepancy. Although there are some special cases for the WD criterion in which analytical optimization is possible (Fang et al., 2005), in general all criteria described above need to be optimized numerically. This can be done within the class of latin hypercube designs proposed by McKay et al. (1979) and discussed in Santner et al. (2003). A latin hypercube in k dimensions with sample size n is a Rn×k matrix where each column is a permutation of the numbers 1, . . . , n. This class is known to have good projection properties. Furthermore, a suitable optimization algorithm is simulated annealing, which has been proposed by Kirkpatrick et al. (1983) and since then has been modified for use in many different areas. Morris and Mitchell (1995) adapted it to search for optimal Latin hypercube designs according to the -criterion. Grosso et al. (2009) constructed maximin designs within the class of latin hypercubes by using iterated local search heuristics. In contrast to these approaches we will create space filling designs without algorithmic optimization. This is motivated by a simple model for deterministic computer experiments. Throughout the literature a broad variety of methods for modeling a deterministic computer experiment based on space filling designs exist, see Fang et al. (2006) for an overview. In general, such a model should fulfill the interpolation constraint: the model has to reproduce the observed data. A very simple way of analyzing computer experiments is to fit a polygon. In one dimension a polygon is uniquely defined, in higher dimensions, however, one has to determine a triangulation (see Section 2 for a definition of triangulations) which can then be used for fitting the polygon. There are several ways of defining a triangulation and the Delaunay triangulation is one with good properties, explained in Section 2. The prediction ability of a polygon model depends on the simplices within the triangulation. If there is no knowledge about the function generating the output of the computer experiment, an intuitive solution is to construct simplices which are very similar to each other according to size and shape. Thus, the problem is to determine a triangulation consisting of almost regular simplices. The paper is structured as follows: In Section 2 we give some definitions of concepts not frequently encountered in statistics and will present an algorithm which constructs a mesh of simplices. In Section 3 projection properties of designs extracted from these meshes are considered and in Section 4 we will compare the obtained designs with maximin and uniform designs. 2. The NRSM algorithm As the use of triangulations is not frequently encountered in statistics, a short definition is given. Definition 1. Consider a point set {x1 , . . . , xn } with xi ∈ Rk , i = 1, . . . , n. A triangulation T consisting of N simplices S1 , . . . , SN is j j partitioning the convex hull of the point set such that the vertices x0 , . . . , xk of simplex Sj are part of the point set {x1 , . . . , xn } and all simplices are disjunct besides their boundaries. It is easy to see that a triangulation of a point set with n > k always exists. In one dimension the triangulation is uniquely defined, in higher dimensions, however, if n > k + 1 there exist more than one triangulation, which rises the question, which triangulation should be chosen. A triangulation often used is the Delaunay triangulation, which has several interesting properties (Rajan, 1994). Here we will use a definition in terms of the circum sphere: In general, for each (nonempty) simplex S, a circum sphere C can be calculated, which has the vertices of the simplex on its boundary. Under weak constraints to the point set the Delaunay triangulation can uniquely be defined as follows. Definition 2. Let {x1 , . . . , xn } be a point set with xi ∈ Rk , i = 1, . . . , n. A triangulation T is the Delaunay triangulation if for each simplex Sj ∈ T, j = 1, . . . , N with its corresponding circum sphere Cj there are no points p ∈ {x1 , . . . , xn } inside the circum sphere Cj and the only points on the boundary of Cj are the vertices of Sj . The circum sphere of a simplex has the interesting property, that the distance from the center of the circum sphere to the vertices of the simplex is equal for all vertices of the simplex. This property will be used later. In the Appendix the calculation of the circum sphere for a given simplex is explained. In Fig. 1 the circum sphere property is explained graphically. Literally we use the Delaunay triangulation as it tends to form as few `sliver-like' simplices as possible. Now we draw our attention to regular tessellations. A regular tessellation is a covering of a subset of Rk (or of Rk itself) which consists of regular polytopes, which are all equal to each other. In R2 there exist point sets P such that the Delaunay triangulation of P forms a regular covering (tessellation) of the design space with simplices. But this result does not generalize to higher dimensions, e.g. in R3 the only covering consisting of regular polytopes is provided by regular cubes (Coxeter, 1963, p. 68) and our results suggest that there is no covering of general subsets of Rk for k > 2 by regular simplices.
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
587
Fig. 1. For four two-dimensional points there are two possible triangulations. The left hand one fulfils the circum sphere criterion and is hence the Delaunay triangulation, whereas for the right hand triangulation, the circum spheres are not empty.
As a possible way out, coverings with nearly regular simplices regarding size and shape of the simplices are calculated, by that Nearly Regular Simplex Mesh (NRSM) Algorithm, which is introduced in the following. Starting with a regular starting simplex a covering is constructed. The vertices of the simplices are then interpreted as potential design points in Rk .
2.1. Algorithm The algorithm used here needs two files: In P the vertices are saved. In S the Delaunay triangulation is stored, i.e. the indices of the vertices belonging to each simplex are stored line by line. Thus, with n vertices it is P ∈ Rn×k and with N simplices it is S ∈ NN×k+1 . At the beginning P contains only the vertices of the starting simplex and S contains one row with the indices 1, . . . , k + 1. The algorithm uses the Delaunay triangulation in order to determine which vertices are combined to a simplex. Beginning with the starting simplex the algorithm presented here suggests a new point xnew (i, j) based on simplex Si and vertex (i) xj of simplex Si . If there does not exist a point in a region around xnew (i, j) defined by {x ∈ Rk : x − xnew (i, j)2 ⱕ L(i, j)} with L(i, j) > 0, then the algorithm adds the new point to the existing points. Otherwise the suggested point is deleted. After adding a point, the Delaunay triangulation S is updated and sorted in two steps. First, each row is sorted in decreasing order. Then the rows are sorted according to the first column in ascending order, whereas ties are broken according to the next columns. The sorting is necessary in order to ensure that the mesh is constructed from inside out. The sorting is illustrated below: ⎛
3 ⎜3 ⎜ ⎜5 ⎜ ⎝1 2
1 4 1 3 6
⎞ ⎛ 4 4 ⎜6 6⎟ ⎟ ⎜ ⎜ 2⎟ ⎟ → ⎜5 ⎝3 2⎠ 5
6
3 4 2 2 5
⎞ ⎛ 1 3 ⎜4 3⎟ ⎟ ⎜ ⎜ 1⎟ ⎟ → ⎜5 ⎝6 1⎠ 2
6
2 3 2 4 5
⎞ 1 1⎟ ⎟ 1⎟ ⎟. 3⎠ 2
2.1.1. NRSM algorithm for the mesh construction (1) i := 1, j := 0, N ∈ N. (i) (i) (i) (i) (i) (2) From simplex Si with vertices x0 , . . . , xk and xj ∈ {x0 , . . . , xk } calculate xnew (i, j) and L(i, j). (3) If ∀x ∈ P : xnew (i, j) − x2 > L(i, j): Add xnew (i, j) to P and update the Delaunay triangulation. Sort each row of S in decreasing order. Sort the rows of S by the first column in increasing order and break ties according to the remaining columns. (4) If j = k then i := i + 1, j := 0 else j := j + 1. (5) If #P < N then go to 2. else stop. The NRSM algorithm is visualized in Fig. 2. The properties of this mesh are mainly determined by the choice of xnew (i, j) and L(i, j). (i) xnew (i, j) determines shape and size of the new simplex {xm |m = 0, . . . , k, m j} ∪ xnew (i, j). L(i, j) determines the region in which no other points are accepted.
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
588
Fig. 2. Schematic representation of the first four steps of the NRSM algorithm for two dimensions. The circles represent the area around a new vertex, in which no other vertices are allowed. The solid triangle is the simplex which is currently processed and the arrow indicates, which vertex is reflected. The number within each triangle is the row index of the simplex in the array S.
2.2. Starting simplex The starting simplex, which is needed for the initialization of the algorithm is constructed to be regular. We choose the edge length for the starting simplex to be 1 and such that all points of the starting simplex possess nonnegative coordinates. Then for dimensions k = 1, 2, 3 regular simplices are provided by (1)
(1)
S1 = {s0 = 0, s1 = 1},
(2) (2) (2) S2 = {s0 = 0, s1 = [1, 0] , s2 = [0.5, 3/4] }, (3) (3) (3) (3) S3 = {s0 = 0, s1 = [1, 0, 0] , s2 = [0.5, 3/4, 0] , s3 = [0.5, 1/12, 2/3] }.
(4)
For dimension k > 3 a starting simplex can be obtained by the following lemma which can be proved by complete induction. (k)
(k)
(k)
Lemma 1. The (k + 1) points s0 , . . . , sk generating a k-dimensional starting simplex with edge length 1 can be given by s0 = 0 and
for d > 0, sd = (sd,1 , . . . , sd,k ) with (k)
(k)
(k)
⎧ 1 ⎪ , i < d, ⎪ ⎪ ⎪ 2i(i + 1) ⎨ (k) sd,i = d + 1 ⎪ ⎪ , i = d, ⎪ ⎪ ⎩ 2d 0, i > d.
(5)
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
589
2.3. Choice of xnew (i, j) and L(i, j) The proposal function for xnew should have a self-regularizing property: from an existing simplex Si a new simplex is constructed by omitting one vertex and adding a new vertex. The new vertex together with the remaining vertices of the old simplex form a new simplex. The proposal function should find xnew in a way that the new simplex is more regular with respect to size and shape than the old simplex. Generally a proposal function should be `feasible' in the following sense. Definition 3. Given a simplex Sold with vertices {xd |d = 0, . . . , k}, a proposal function constructs a new simplex Snew by removing one vertex xj of Sold and adding one new vertex xnew (i, j) xj : Snew = {xnew (i, j), xd , d = 0, . . . , k, d j}. The proposal function is called feasible if it fulfills: for a regular simplex Sold with edge length 1 the simplex Snew is also regular with edge length 1. 2.3.1. Proposal function 1 The most simple feasible proposal function is defined by ⎛ 2 xnew (i, j) := ⎝ k
⎞
k
xd ⎠ − xj
(6)
d=0,d j
with {xd |d = 0, . . . , k} being the vertices of Si . Geometrically this means a point reflection of xj to the mean of the vertices {xd |d = 0, . . . , k, d j}. However, it is not a good choice as a simplex with poor shape and size properties will not be corrected. The new simplex has exactly the same content and the same point to point distances as the old one, which is not desired if the old simplex is poorly shaped or does not have approximately the same content as the starting simplex. 2.3.2. Proposal function 2 A better way which will be used in the following is one which works independently of the vertex xj , which is omitted. The new vertex is presented by a vector with end point v and origin at the circum center x¯ of the vertices {xd |d = 0, . . . , k, d j}. The computation of the circum center is explained in the Appendix. The vector is orthogonal to the affine hyperplane defined by {xd |d = 0, . . . , k, d j}. The length of the vector is constructed such that the distance between v and each of the vertices {xd |d = 0, . . . , k, d j} equals 1. The vector product is used to compute the orthogonal vector: Definition 4. Given k − 1 points x1 , . . . , xk−1 in Rk , the vector product VP(x1 , . . . , xk−1 ) given by VP(x1 , . . . , xk−1 ) :=
k
ed det(x1 , . . . , xk−1 , ed )
(7)
d=1
calculates a vector, which is orthogonal to the hyperplane defined by the origin and x1 , . . . , xk−1 . It is uniquely defined up to its sign. When applying the vector product, the simplex under consideration has to be moved such that a vertex xd xj equals the origin. Here we used the first vertex x0 for the translation (or, of j = 0, the second vertex x1 ): x˜ d := xd − x0 , d = 1, . . . , k. The vector product is applied to x˜ d , d = 1, . . . , k, d j: vˆ := VP(x˜ 1 , . . . , x˜ j−1 , x˜ j+1 , . . . , x˜ k ). In order to have that vˆ is orthogonal to the hyperplane with origin in the circum center x¯ of the vertices, the vertices are transformed once again: xˆ d := x˜ d − x¯ , d = 1, . . . , k. Now the length of vˆ can be scaled by a constant b such that Euclidean distance ˆ b ∈ R: between vˆ and xˆ d , d = 1, . . . , k, d j equals 1. Rewrite v := bv, ˆ 22 = 1 ⇔ b = ± xˆ d − bv
1 − xˆ d 22 ˆ 22 v
.
Considering to reverse the translations done, the point xnew is given by xnew = ±bvˆ + x¯ + x0 . In order to get the `correct' orientation of xnew , it has to be checked which one of the two possibilities is farther away from the old vertex xj . Fig. 3 shows how this proposal function works.
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
590
Fig. 3. Illustration of proposal function 2 for different simplices. The new vertex is constructed such that the point to point distances (denoted by | · |) from the new vertex to the already existing vertices always equals 1. By this, regular simplices produce regular simplices (left hand plot). For the middle and right hand plot, the old simplices are not regular as well as the new ones, but it can be achieved, that all new edges have the same length.
2.3.3. Choice of L(i, j) For the first proposal function a simple choice for L(i, j) is L(i, j) := c, c < 1 constant. However, as this proposal function has not delivered satisfying results in general, in the context of the second proposal function L(i, j) = 1 has been chosen. By this, the algorithm is constructing point sets for which the minimum distance among all pairs of points is 1. 2.4. Extracting a design from a mesh To extract a design from a mesh, first center the points in P with respect to the mesh mean: P¯ := {x− x¯ |x ∈ P, x¯ := (1/#P) x∈P x}. In order to choose a design search all points x ∈ P¯ with x∞ < a ∈ R+ . Then vary a to get the desired amount of design points and scale the resulting design points to the unit cube. In order to further improve the properties of the design, it is possible to scale every column of the design in a way that the maximum and the minimum of this column are 1 and 0, respectively. It is important that #P ?n. The NRSM algorithm constructs a sphere like structure. To prevent that the `corners' of a design are empty a sufficiently large number of points is required. The radius of the sphere can be approximated by r := maxx∈P¯ (x2 ). The design space is a cube√with side length 2a centered around the origin, for which the maximum L2 norm of any point within √ the cube is provided by a k. Then the design with largest sample size n extracted from P¯ should have a side length less than r/ k. 3. Modification of the NRSM algorithm The designs resulting from the NRSM algorithm are constructed such that they are space filling in Rk . However, this does not necessarily imply that it is also a good space filling design in lower dimensions. Fig. 4 shows a scatter plot matrix of two designs with sample size n = 50 and dimension k = 5. The plots in the upper right part show an NRSM algorithm based design, the plots in the lower left part show a random latin hypercube design. Obviously the random latin hypercube design fills out the two-dimensional design space better than the simplex based design although the corresponding 40 -value of the simplex based design is much lower. For this reason we suggest two methods to improve the projection properties of the mesh, which can be used in a combined manner or separately. The first one consists of a rotation of the mesh prior to extracting designs, which does not change the Euclidian point to point distances in Rk , but the properties of the mesh when projecting it to lower dimensions can be improved. The second method uses rank version designs extracted from simplex meshes which means to transfer the designs into the class of latin hypercubes. 3.1. Rotated meshes A real square matrix R satisfying RR = I (Orthogonality) and det(R) = 1
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
591
Fig. 4. Scatter plot matrix of two n = 50, k = 5 designs. The upper right scatter plots represent a design based on simplex meshes, the lower left ones a random latin hypercube design.
is a rotation matrix, see Meyer (2000). Rotation matrices can be constructed stepwise. Define a plane rotation matrix regarding dimension i and j as ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ Ri,j () = ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
1
⎞
... ..
⎟ ⎟ ⎟ ⎟ ⎟ ←i ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ←j ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
. cos
sin 1 ..
− sin
. cos 1 ..
↑ i
↑ j
. 1
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
592
Fig. 5. Scatter plot matrix of two n = 50, k = 5 designs. The upper right scatter plots represent a design based on simplex meshes, the lower left ones a design extracted from a rotated simplex mesh, to which the rank transformation has been applied.
with all nondiagonal entries equal to 0 except the two entries at position (i, j) and (j, i). A general rotation matrix can then be defined as the product of plane rotation matrices: R( = 1 , . . . , k(k−1)/2 ) := R1,2 (1 ) · · · R1,k (k−1 )R2,3 (k ) · · · R2,k (2k−3 ) · · · Rk−1,k (k(k−1)/2 ). The choice of the vector ∈ Rk(k−1)/2 is crucial for the improvement of the mesh. To compare rotated meshes a criterion is needed. We choose a maximin approach which can be applied to each column of P:
() = max {p ((PR())j )}, j=1, ...,k
where (PR())j is the j-th column of PR(). Minimizing this criterion means to search for the rotation which has the best one ∗ dimensional projection properties. Among all values {0 ⱕ ⱕ 2} := {0 ⱕ 1 ⱕ 2, . . . , 0 ⱕ k(k−1)/2 ⱕ 2} the value which minimizes () is searched for numerically. 3.2. Rank version of NRSM based designs As the class of latin hypercubes is very attractive due to its projection properties, the designs based on NRS meshes can be transferred into this class by applying the rank transformation to each column of the design and scaling the resulting design to
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
593
the unit cube. By numerical comparison, we show that these rank transformed designs still have good space filling properties regarding the maximin criterion but have better WD values than designs without rank transformation. If D ∈ [0, 1]n×k is a design extracted from a simplex mesh, for each column D.,d of D the n entries of the column can be ranked with numbers 1, . . . , n. Hence, if for entry (i, d) rank(Di,d ) = m ∈ {1, . . . , n}, this means, that there are m − 1 entries j in column d with Dj,d < Di,d and n − m entries j with Dj,d > Di,d . By applying the rank transformation to each column of the design D, the order of the one dimensional entries is preserved but the one-dimensional point to point distances are normalized. In Fig. 5 a scatter plot matrix is shown, which has the two dimensional projections of the simplex based design from Fig. 4 (from an unrotated mesh) at the upper right part and on the lower left part the two dimensional projections of a simplex based design from a rotated mesh with the same dimension and sample size, to which the rank transformation has been applied. Although the corresponding -criterion values are comparable, the rank design extracted from the rotated mesh has better projection properties. 4. Comparison and discussion The space filling designs constructed here will now be compared to designs constructed according to the maximin criteria
p (D ) and Mm (DMm ) and the WD criterion (DWD ). The criteria p and WD need to be minimized to get optimal designs. For
these two criteria the minimization has been done by simulated annealing within the class of latin hypercubes. The computation times for the numerically optimized designs are between 14.9 and 279.7 s (-criterion) and 17.8 and 766.8 s for the WD-criterion on a 2 GHz processor using R.2.8.1 under Windows XP. The extraction of the simplex based designs in all cases needed < 1 s, independent from dimension and sample size. The Mm criterion needs to be maximized. Here we used, where available, results explained by Grosso et al. (2009) and listed at www.spacefillingdesigns.nl. For these designs no computation times are reported. Let DR , DS be designs based on rotated meshes with and without rank transformation, respectively. As a minimum requirement would be that the constructed designs perform better than random latin hypercube designs, for each combination of dimension and sample size 1000 random latin hypercubes Drandom have been generated and the mean of the resulting criterion values is also included in the comparison. The results are shown in Table 1 (p ), Table 2 (Mm) and Table 3 (WD). For all three tables, the first and second column show the criterion values for the DS and DR designs. The third column reports the mean results of the respective criteria for 1000 random latin hypercube designs. The last column always shows criterion values for numerically optimized designs. The results from Tables 1 and 2 show that the simplex based design from rotated meshes without rank transformation (DS ) are highly competitive to numerically optimized maximin latin hypercube designs. For the 40 criterion in Table 1, all simplex based designs with and without rank transformation perform better than the latin hypercubes optimized by simulated annealing and as the mean of the random latin hypercube designs. The numerically optimized designs from Table 2 are constructed by Grosso et al. (2009) by using an iterated local search heuristic. The local search is applied to a slightly perturbed version of the latin hypercube from the last local search. The designs extracted from a rotated mesh without rank transformation perform slightly worse than the results reported by Grosso et al. (2009). The rank versions of the simplex based designs lose some efficiency w.r.t. the Mm criterion but still perform much better than the random latin hypercube designs. The results for the uniform designs with respect to the WD criterion in Table 3 show a different behavior. Here the simplex based designs without rank transformation perform worse than the mean of the criterion values for the 1000 random latin hypercube designs, indicating that the designs without rank transformation are not attractive as uniform designs (and hence are not recommendable for integral estimation). But when applying the rank transformation to the simplex based designs, the rank versions are performing only slightly worse than the numerically optimized but much better than the mean of the random latin hypercube designs. Table 1 Numerical results for the 40 -criterion.
40 (DS )
40 (DR )
40 (Drandom )
40 (D )@
3
30 60
3.214 3.976
3.394 4.582
10.570 18.140
3.787 5.269
4
40 80
2.429 3.072
2.723 3.293
7.368 11.630
2.768 3.563
5
50 100
2.058 2.467
2.237 2.557
5.556 7.916
2.327 2.993
6
60 120
1.829 2.144
1.966 2.193
4.444 5.919
2.088 2.430
7
70 140
1.642 1.883
1.723 1.919
3.634 4.668
1.794 2.185
8
80 160
1.484 1.736
1.598 1.727
3.113 3.797
1.619 1.904
k
n
The numerically optimized design has been searched for by simulated annealing within the class of latin hypercubes.
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
594 Table 2 Numerical results for the Mm-criterion. k
Mm(DS )
Mm(DR )
Mm(Drandom )
Mm(DMm )
3
n 30 60
0.3436 0.2830
0.2966 0.2249
0.1028 0.06046
0.3600 0.2800
4
40 80
0.4534 0.3715
0.3795 0.3157
0.1451 0.0937
0.4763 0.3899
5
50 100
0.5435 0.4659
0.4600 0.4194
0.1908 0.1339
0.5894 0.4949
6
60 120
0.6147 0.5402
0.5127 0.4855
0.2344 0.1776
0.6879 –
7
70 140
0.7004 0.6244
0.6021 0.5504
0.2855 0.2235
0.7798 –
8
80 160
0.7736 0.6781
0.6474 0.6201
0.3305 0.2712
0.8673 –
Where available, numerically optimized designs have been taken from www.spacefillingdesigns.nl.
Table 3 Numerical results for the uniform WD-criterion. k
WD(DS )
WD(DR )
WD(Drandom )
WD(DWD )@
3
n 30 60
0.01706 0.00709
0.00458 0.00156
0.00606 0.00244
0.00418 0.00144
4
40 80
0.02055 0.00949
0.00830 0.00293
0.01024 0.00457
0.00725 0.00262
5
50 100
0.03752 0.01483
0.01309 0.00553
0.01732 0.00808
0.01274 0.00485
6
60 120
0.05431 0.02527
0.02191 0.01010
0.02865 0.01371
0.02063 0.00851
7
70 140
0.12332 0.05574
0.04109 0.01889
0.04656 0.02254
0.03178 0.01447
8
80 160
0.17347 0.09463
0.06892 0.03039
0.07446 0.03628
0.05354 0.02618
The numerically optimized design has been searched for by simulated annealing within the class of latin hypercubes.
4.1. Numerical issues The computationally most intensive part of the algorithm is the calculation of the Delaunay triangulations. For a fixed number n of vertices, the number of simplices grows rapidly with the dimension, see Klee (1980). For dimensions larger than 7, this practically prohibits the calculation of a sufficient number of points on a personal computer. Also on servers with sufficient memory this limits the dimension up to which a sufficiently large number of points can be constructed. Here the quickhull algorithm (Barber et al., 1996) has been used to compute the Delaunay triangulation. The NRSM algorithm presented here has been implemented in `R' and Matlab. 5. Summary The designs suggested here are interesting alternatives to the fast construction of good space filling designs. The procedure to construct a design according to the methodology explained above looks as follows: First a simplex mesh with point set P ∈ Rn×k is constructed by the NRSM algorithm and the points in P are rotated such that good projection properties are achieved. These two steps only have to be done once for each dimension k. Hence it is acceptable that the construction and rotation of the point sets is time consuming. The rotated point sets for dimensions k = 3, . . . , 8 can be downloaded from www.statistik.tu-dortmund.de/ muehlenstaedt.html. The rotated point set PR can then be used repeatedly to extract designs from it. If a latin hypercube is desired, the rank transformation can be applied to the extracted design. Once the underlying simplex meshes are found and stored as a library file, the construction of a design is very simple and fast (usually < 1 s, independent from dimension and sample size). The designs we get by the methods presented here are generally competitive regarding designs optimized with respect to the maximin design criteria and Mm. The rank version of the simplex based designs are also nearly as good as numerically optimized designs regarding the WD criterion. In contrast to the simplex based designs, for all three criteria the optimization process can be very difficult and time expensive. Another advantage is that it is also very easy to construct space filling designs for nonstandard
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
595
design spaces, i.e. design spaces which are not rectangles. We only have to modify the way design points are extracted from the mesh. With other methods like uniform designs, this would be a much more difficult task. Appendix A. Calculation of the circum sphere Okabe et al. (2000, p. 80) state that in two dimensions, the circum circle with center c and radius r for a triangle can be calculated based on the following equation: 2 2 i=1 c i 2 x2 i=1 1,i 2 i=1 x22,i 2 x2 i=1 3,i
c1
c2
x1,1
x1,2
x2,1
x2,2
x3,1
x3,2
1 1 = 0. 1 1
(8)
Here (xj,1 , xj,2 )T , j = 1, 2, 3 are the coordinates of the vertices of the triangle and (c1 , c2 )T is the circum center of the circum circle of the triangle which is searched for. This equation can easily be extended to k > 2 dimensions: k 2 i=1 ci k 2 i=1 x1,i .. . k x2
c1
...
ck
x1,1
...
x1,k
xk+1,1
i=1 k+1,i
...
xk+1,k
1 1 = 0, .. . 1
(9)
where again (xj,1 , . . . , xj,k )T , j = 1, . . . , k + 1 are the coordinates of the vertices of the simplex and (c1 , . . . , ck )T are the coordinates of the circum center of the circum sphere of the simplex. In order to solve this equation cofactor expansion along the first row can be used. Define k 2 x1,1 . . . x1,j−1 x1,j+1 . . . x1,k 1 i=1 x1,i k 2 x x . . . x x . . . x 1 2,1 2,j−1 2,j+1 2,k i=1 2,i 2+j Dj := (−1) (10) .. .. . . k x2 x ... x x ... x 1 i=1 k+11,i
k+1,1
k+1,j−1
k+1,j+1
k+1,k
and x1,1 x2,1 := . .. xk+1,1
...
x1,k
...
x2,k
. . . xk+1,k
1 1 .. , . 1
k 2 i=1 x1,i k 2 i=1 x2,i := .. . k x2
i=1 k+1,i
x1,1
...
x2,1
...
xk+1,1
...
x1,k x2,k . .. . xk+1,k
(11)
With the cofactor expansion equation (8) becomes k
ci2 −
i=1
k
Dj cj + = 0.
(12)
j=1
Completing the squares yields cj = Dj /2. The radius r of the circum sphere can be calculated as r=
k 2 j=1 Dj
− 4
4 2
(13)
or alternatively just by calculating the distance of (c1 , . . . , ck )T to one of the vertices of the simplex. As the circum center is the point, which has the same distance to all vertices of the simplex, it is arbitrary, which vertex is chosen. For the second proposal function of Section 2.3, the circum center of the facet of a simplex is needed, i.e. for k vertices x0 , . . . , xj−1 , xj+1 , . . . , xk of simplex S with vertices x0 , . . . , xk . This can be deduced from the circum center of the complete simplex by (orthogonally) projecting the circum center c of S to the hyperplane defined by the k points x0 , . . . , xj−1 , xj+1 , . . . , xk , which can be done by using the Moore Penrose inverse.
596
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
References Barber, C., Dobkin, D., Huhdanpaa, H., 1996. The quickhull algorithm for convex hulls. ACM Transactions on Mathematical Software 22 (4), 469–483. Coxeter, H., 1963. Regular Polytopes. The Macmillan Company, New York. Fang, K., Tang, Y., Yin, J., 2005. Lower bounds for wrap-around l2 -discrepancy and construction of symmetrical uniform designs. Journal of Complexity 21, 757–771. Fang, K.-T., Li, R., Sudjianto, A., 2006. Design and Modeling for Computer Experiments. Computer Science and Data Analysis Series. Chapman & Hall, CRC, Boca Raton. Grosso, A., Jamali, A., Locatelli, M., 2009. Finding maximim latin hypercube designs by iterated local search heuristics. European Journal of Operational Research 197, 541–547. Johnson, M., Moore, L., Ylvisaker, D., 1990. Minimax and maximin distance designs. Journal of Statistical Planning and Inference 26, 131–148. Kirkpatrick, S., Gelatt, C., Vecchi, M., 1983. Optimization by simulated annealing. Science 220 (4598), 671–680. Klee, V., 1980. On the complexity of d-dimensional Voronoi diagrams. Archiv der Mathematik 34, 75–80. McKay, M., Beckman, R., Conover, W., 1979. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (2), 239–245. Meyer, C., 2000. Matrix Analysis and Applied Linear Algebra. SIAM, Philadelphia. Morris, M., Mitchell, T., 1995. Exploratory designs for computational experiments. Journal of Statistical Planning and Inference 43, 381–402. Okabe, A., Boots, B., Sugihara, K., Chiu, S., 2000. Spatial Tessellations, Concepts and Applications of Voronoi Diagrams. Wiley Series in Probability and Statistics. Wiley, Chichester. d Rajan, V., 1994. Optimality of the Delaunay triangulation in R . Discrete & Computational Geometry 12, 189–202. Santner, T., Williams, B., Notz, W., 2003. The Design and Analysis of Computer Experiments. Springer Series in Statistics. Springer, New York.
Contents lists available at ScienceDirect
Journal of Statistical Planning and Inference journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / j s p i
Simplex based space filling designs ¨ Thomas Muhlenst a¨ dt∗ ,1 Fakulta¨ t Statistik, Technische Universita¨ t Dortmund, 44221 Dortmund, Germany
A R T I C L E
I N F O
A B S T R A C T
Article history: Received 19 June 2008 Received in revised form 2 June 2009 Accepted 4 August 2009 Available online 11 August 2009 Keywords: Triangulation Maximin design Uniform design Computer experiment
Space filling designs are important for deterministic computer experiments. Even a single experiment can be very time consuming and can have many input parameters. Furthermore the underlying function generating the output is often nonlinear. Thus, the computer experiment has to be designed carefully. There exist many design criteria, which can be numerically optimized. Here, a method is developed, which does not need algorithmic optimization. A mesh of nearly regular simplices is constructed and the vertices of the simplices are used as potential design points. The extraction of a design from these meshes is very fast and easy to implement once the underlying mesh has been constructed. The extracted designs are highly competitive regarding the maximin design criterion and it is easy to extract designs for nonstandard design spaces. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Frequently used approaches to obtain a space filling design D with n design points x1 , . . . , xn in k dimensions are maximin designs and uniform designs, see Santner et al. (2003) and Fang et al. (2006), respectively. The maximin approach consists of maximizing the minimum distance between two arbitrary points of the design: Mm(D) := min(xi − xj 2 ). ij
(1)
Johnson et al. (1990)argued that maximin designs have some asymptotic optimality property in a Bayesian context. As a maximin design is not necessarily uniquely defined, Morris and Mitchell (1995) suggest the following criterion: ⎛
p (D) := ⎝
m
⎞(1/p) −p Jj dj ⎠
(2)
j=1
with (d1 < · · · < dm ) being the sorted list of the point to point distances of the design and (J1 , . . . , Jm ) being the frequency list of the distances with m ⱕ ( 2n ). This criterion needs to be minimized numerically. p is only an approximation for a theoretical generalization of the maximin design and the parameter p influences how good the approximation is. For large p the approximation improves but p will tend to have many local minima. Hence, for increasing p it will be more difficult to find a global minimum numerically. Morris and Mitchell (1995) suggest to use p ∈ {20, . . . , 50}, therefore we use p = 40 in this article to have a good approximation. ∗ Tel.: +49 231 7557227; fax: +49 231 7555305. E-mail addresses: [email protected], [email protected] URL: http://www.statistik.tu-dortmund.de. 1 The financial support of the DFG (research training group Statistical Modelling) is gratefully acknowledged. 0378-3758/$ - see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2009.08.002
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
586
Another commonly used approach to create space filling designs are uniform designs which are trying to minimize a test statistic of a test on a U[0, 1]k -distribution. In this context special test statistics have been constructed, among those the wrap around criterion, which is defined by k n n k
4 1 3 (WD(D))2 = − + 2 − |xdi − xji |(1 − |xdi − xji |) . 3 2 n
(3)
d=1 j=1 i=1
See Fang et al. (2006) for a more general description of the wrap around discrepancy. Although there are some special cases for the WD criterion in which analytical optimization is possible (Fang et al., 2005), in general all criteria described above need to be optimized numerically. This can be done within the class of latin hypercube designs proposed by McKay et al. (1979) and discussed in Santner et al. (2003). A latin hypercube in k dimensions with sample size n is a Rn×k matrix where each column is a permutation of the numbers 1, . . . , n. This class is known to have good projection properties. Furthermore, a suitable optimization algorithm is simulated annealing, which has been proposed by Kirkpatrick et al. (1983) and since then has been modified for use in many different areas. Morris and Mitchell (1995) adapted it to search for optimal Latin hypercube designs according to the -criterion. Grosso et al. (2009) constructed maximin designs within the class of latin hypercubes by using iterated local search heuristics. In contrast to these approaches we will create space filling designs without algorithmic optimization. This is motivated by a simple model for deterministic computer experiments. Throughout the literature a broad variety of methods for modeling a deterministic computer experiment based on space filling designs exist, see Fang et al. (2006) for an overview. In general, such a model should fulfill the interpolation constraint: the model has to reproduce the observed data. A very simple way of analyzing computer experiments is to fit a polygon. In one dimension a polygon is uniquely defined, in higher dimensions, however, one has to determine a triangulation (see Section 2 for a definition of triangulations) which can then be used for fitting the polygon. There are several ways of defining a triangulation and the Delaunay triangulation is one with good properties, explained in Section 2. The prediction ability of a polygon model depends on the simplices within the triangulation. If there is no knowledge about the function generating the output of the computer experiment, an intuitive solution is to construct simplices which are very similar to each other according to size and shape. Thus, the problem is to determine a triangulation consisting of almost regular simplices. The paper is structured as follows: In Section 2 we give some definitions of concepts not frequently encountered in statistics and will present an algorithm which constructs a mesh of simplices. In Section 3 projection properties of designs extracted from these meshes are considered and in Section 4 we will compare the obtained designs with maximin and uniform designs. 2. The NRSM algorithm As the use of triangulations is not frequently encountered in statistics, a short definition is given. Definition 1. Consider a point set {x1 , . . . , xn } with xi ∈ Rk , i = 1, . . . , n. A triangulation T consisting of N simplices S1 , . . . , SN is j j partitioning the convex hull of the point set such that the vertices x0 , . . . , xk of simplex Sj are part of the point set {x1 , . . . , xn } and all simplices are disjunct besides their boundaries. It is easy to see that a triangulation of a point set with n > k always exists. In one dimension the triangulation is uniquely defined, in higher dimensions, however, if n > k + 1 there exist more than one triangulation, which rises the question, which triangulation should be chosen. A triangulation often used is the Delaunay triangulation, which has several interesting properties (Rajan, 1994). Here we will use a definition in terms of the circum sphere: In general, for each (nonempty) simplex S, a circum sphere C can be calculated, which has the vertices of the simplex on its boundary. Under weak constraints to the point set the Delaunay triangulation can uniquely be defined as follows. Definition 2. Let {x1 , . . . , xn } be a point set with xi ∈ Rk , i = 1, . . . , n. A triangulation T is the Delaunay triangulation if for each simplex Sj ∈ T, j = 1, . . . , N with its corresponding circum sphere Cj there are no points p ∈ {x1 , . . . , xn } inside the circum sphere Cj and the only points on the boundary of Cj are the vertices of Sj . The circum sphere of a simplex has the interesting property, that the distance from the center of the circum sphere to the vertices of the simplex is equal for all vertices of the simplex. This property will be used later. In the Appendix the calculation of the circum sphere for a given simplex is explained. In Fig. 1 the circum sphere property is explained graphically. Literally we use the Delaunay triangulation as it tends to form as few `sliver-like' simplices as possible. Now we draw our attention to regular tessellations. A regular tessellation is a covering of a subset of Rk (or of Rk itself) which consists of regular polytopes, which are all equal to each other. In R2 there exist point sets P such that the Delaunay triangulation of P forms a regular covering (tessellation) of the design space with simplices. But this result does not generalize to higher dimensions, e.g. in R3 the only covering consisting of regular polytopes is provided by regular cubes (Coxeter, 1963, p. 68) and our results suggest that there is no covering of general subsets of Rk for k > 2 by regular simplices.
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
587
Fig. 1. For four two-dimensional points there are two possible triangulations. The left hand one fulfils the circum sphere criterion and is hence the Delaunay triangulation, whereas for the right hand triangulation, the circum spheres are not empty.
As a possible way out, coverings with nearly regular simplices regarding size and shape of the simplices are calculated, by that Nearly Regular Simplex Mesh (NRSM) Algorithm, which is introduced in the following. Starting with a regular starting simplex a covering is constructed. The vertices of the simplices are then interpreted as potential design points in Rk .
2.1. Algorithm The algorithm used here needs two files: In P the vertices are saved. In S the Delaunay triangulation is stored, i.e. the indices of the vertices belonging to each simplex are stored line by line. Thus, with n vertices it is P ∈ Rn×k and with N simplices it is S ∈ NN×k+1 . At the beginning P contains only the vertices of the starting simplex and S contains one row with the indices 1, . . . , k + 1. The algorithm uses the Delaunay triangulation in order to determine which vertices are combined to a simplex. Beginning with the starting simplex the algorithm presented here suggests a new point xnew (i, j) based on simplex Si and vertex (i) xj of simplex Si . If there does not exist a point in a region around xnew (i, j) defined by {x ∈ Rk : x − xnew (i, j)2 ⱕ L(i, j)} with L(i, j) > 0, then the algorithm adds the new point to the existing points. Otherwise the suggested point is deleted. After adding a point, the Delaunay triangulation S is updated and sorted in two steps. First, each row is sorted in decreasing order. Then the rows are sorted according to the first column in ascending order, whereas ties are broken according to the next columns. The sorting is necessary in order to ensure that the mesh is constructed from inside out. The sorting is illustrated below: ⎛
3 ⎜3 ⎜ ⎜5 ⎜ ⎝1 2
1 4 1 3 6
⎞ ⎛ 4 4 ⎜6 6⎟ ⎟ ⎜ ⎜ 2⎟ ⎟ → ⎜5 ⎝3 2⎠ 5
6
3 4 2 2 5
⎞ ⎛ 1 3 ⎜4 3⎟ ⎟ ⎜ ⎜ 1⎟ ⎟ → ⎜5 ⎝6 1⎠ 2
6
2 3 2 4 5
⎞ 1 1⎟ ⎟ 1⎟ ⎟. 3⎠ 2
2.1.1. NRSM algorithm for the mesh construction (1) i := 1, j := 0, N ∈ N. (i) (i) (i) (i) (i) (2) From simplex Si with vertices x0 , . . . , xk and xj ∈ {x0 , . . . , xk } calculate xnew (i, j) and L(i, j). (3) If ∀x ∈ P : xnew (i, j) − x2 > L(i, j): Add xnew (i, j) to P and update the Delaunay triangulation. Sort each row of S in decreasing order. Sort the rows of S by the first column in increasing order and break ties according to the remaining columns. (4) If j = k then i := i + 1, j := 0 else j := j + 1. (5) If #P < N then go to 2. else stop. The NRSM algorithm is visualized in Fig. 2. The properties of this mesh are mainly determined by the choice of xnew (i, j) and L(i, j). (i) xnew (i, j) determines shape and size of the new simplex {xm |m = 0, . . . , k, m j} ∪ xnew (i, j). L(i, j) determines the region in which no other points are accepted.
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
588
Fig. 2. Schematic representation of the first four steps of the NRSM algorithm for two dimensions. The circles represent the area around a new vertex, in which no other vertices are allowed. The solid triangle is the simplex which is currently processed and the arrow indicates, which vertex is reflected. The number within each triangle is the row index of the simplex in the array S.
2.2. Starting simplex The starting simplex, which is needed for the initialization of the algorithm is constructed to be regular. We choose the edge length for the starting simplex to be 1 and such that all points of the starting simplex possess nonnegative coordinates. Then for dimensions k = 1, 2, 3 regular simplices are provided by (1)
(1)
S1 = {s0 = 0, s1 = 1},
(2) (2) (2) S2 = {s0 = 0, s1 = [1, 0] , s2 = [0.5, 3/4] }, (3) (3) (3) (3) S3 = {s0 = 0, s1 = [1, 0, 0] , s2 = [0.5, 3/4, 0] , s3 = [0.5, 1/12, 2/3] }.
(4)
For dimension k > 3 a starting simplex can be obtained by the following lemma which can be proved by complete induction. (k)
(k)
(k)
Lemma 1. The (k + 1) points s0 , . . . , sk generating a k-dimensional starting simplex with edge length 1 can be given by s0 = 0 and
for d > 0, sd = (sd,1 , . . . , sd,k ) with (k)
(k)
(k)
⎧ 1 ⎪ , i < d, ⎪ ⎪ ⎪ 2i(i + 1) ⎨ (k) sd,i = d + 1 ⎪ ⎪ , i = d, ⎪ ⎪ ⎩ 2d 0, i > d.
(5)
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
589
2.3. Choice of xnew (i, j) and L(i, j) The proposal function for xnew should have a self-regularizing property: from an existing simplex Si a new simplex is constructed by omitting one vertex and adding a new vertex. The new vertex together with the remaining vertices of the old simplex form a new simplex. The proposal function should find xnew in a way that the new simplex is more regular with respect to size and shape than the old simplex. Generally a proposal function should be `feasible' in the following sense. Definition 3. Given a simplex Sold with vertices {xd |d = 0, . . . , k}, a proposal function constructs a new simplex Snew by removing one vertex xj of Sold and adding one new vertex xnew (i, j) xj : Snew = {xnew (i, j), xd , d = 0, . . . , k, d j}. The proposal function is called feasible if it fulfills: for a regular simplex Sold with edge length 1 the simplex Snew is also regular with edge length 1. 2.3.1. Proposal function 1 The most simple feasible proposal function is defined by ⎛ 2 xnew (i, j) := ⎝ k
⎞
k
xd ⎠ − xj
(6)
d=0,d j
with {xd |d = 0, . . . , k} being the vertices of Si . Geometrically this means a point reflection of xj to the mean of the vertices {xd |d = 0, . . . , k, d j}. However, it is not a good choice as a simplex with poor shape and size properties will not be corrected. The new simplex has exactly the same content and the same point to point distances as the old one, which is not desired if the old simplex is poorly shaped or does not have approximately the same content as the starting simplex. 2.3.2. Proposal function 2 A better way which will be used in the following is one which works independently of the vertex xj , which is omitted. The new vertex is presented by a vector with end point v and origin at the circum center x¯ of the vertices {xd |d = 0, . . . , k, d j}. The computation of the circum center is explained in the Appendix. The vector is orthogonal to the affine hyperplane defined by {xd |d = 0, . . . , k, d j}. The length of the vector is constructed such that the distance between v and each of the vertices {xd |d = 0, . . . , k, d j} equals 1. The vector product is used to compute the orthogonal vector: Definition 4. Given k − 1 points x1 , . . . , xk−1 in Rk , the vector product VP(x1 , . . . , xk−1 ) given by VP(x1 , . . . , xk−1 ) :=
k
ed det(x1 , . . . , xk−1 , ed )
(7)
d=1
calculates a vector, which is orthogonal to the hyperplane defined by the origin and x1 , . . . , xk−1 . It is uniquely defined up to its sign. When applying the vector product, the simplex under consideration has to be moved such that a vertex xd xj equals the origin. Here we used the first vertex x0 for the translation (or, of j = 0, the second vertex x1 ): x˜ d := xd − x0 , d = 1, . . . , k. The vector product is applied to x˜ d , d = 1, . . . , k, d j: vˆ := VP(x˜ 1 , . . . , x˜ j−1 , x˜ j+1 , . . . , x˜ k ). In order to have that vˆ is orthogonal to the hyperplane with origin in the circum center x¯ of the vertices, the vertices are transformed once again: xˆ d := x˜ d − x¯ , d = 1, . . . , k. Now the length of vˆ can be scaled by a constant b such that Euclidean distance ˆ b ∈ R: between vˆ and xˆ d , d = 1, . . . , k, d j equals 1. Rewrite v := bv, ˆ 22 = 1 ⇔ b = ± xˆ d − bv
1 − xˆ d 22 ˆ 22 v
.
Considering to reverse the translations done, the point xnew is given by xnew = ±bvˆ + x¯ + x0 . In order to get the `correct' orientation of xnew , it has to be checked which one of the two possibilities is farther away from the old vertex xj . Fig. 3 shows how this proposal function works.
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
590
Fig. 3. Illustration of proposal function 2 for different simplices. The new vertex is constructed such that the point to point distances (denoted by | · |) from the new vertex to the already existing vertices always equals 1. By this, regular simplices produce regular simplices (left hand plot). For the middle and right hand plot, the old simplices are not regular as well as the new ones, but it can be achieved, that all new edges have the same length.
2.3.3. Choice of L(i, j) For the first proposal function a simple choice for L(i, j) is L(i, j) := c, c < 1 constant. However, as this proposal function has not delivered satisfying results in general, in the context of the second proposal function L(i, j) = 1 has been chosen. By this, the algorithm is constructing point sets for which the minimum distance among all pairs of points is 1. 2.4. Extracting a design from a mesh To extract a design from a mesh, first center the points in P with respect to the mesh mean: P¯ := {x− x¯ |x ∈ P, x¯ := (1/#P) x∈P x}. In order to choose a design search all points x ∈ P¯ with x∞ < a ∈ R+ . Then vary a to get the desired amount of design points and scale the resulting design points to the unit cube. In order to further improve the properties of the design, it is possible to scale every column of the design in a way that the maximum and the minimum of this column are 1 and 0, respectively. It is important that #P ?n. The NRSM algorithm constructs a sphere like structure. To prevent that the `corners' of a design are empty a sufficiently large number of points is required. The radius of the sphere can be approximated by r := maxx∈P¯ (x2 ). The design space is a cube√with side length 2a centered around the origin, for which the maximum L2 norm of any point within √ the cube is provided by a k. Then the design with largest sample size n extracted from P¯ should have a side length less than r/ k. 3. Modification of the NRSM algorithm The designs resulting from the NRSM algorithm are constructed such that they are space filling in Rk . However, this does not necessarily imply that it is also a good space filling design in lower dimensions. Fig. 4 shows a scatter plot matrix of two designs with sample size n = 50 and dimension k = 5. The plots in the upper right part show an NRSM algorithm based design, the plots in the lower left part show a random latin hypercube design. Obviously the random latin hypercube design fills out the two-dimensional design space better than the simplex based design although the corresponding 40 -value of the simplex based design is much lower. For this reason we suggest two methods to improve the projection properties of the mesh, which can be used in a combined manner or separately. The first one consists of a rotation of the mesh prior to extracting designs, which does not change the Euclidian point to point distances in Rk , but the properties of the mesh when projecting it to lower dimensions can be improved. The second method uses rank version designs extracted from simplex meshes which means to transfer the designs into the class of latin hypercubes. 3.1. Rotated meshes A real square matrix R satisfying RR = I (Orthogonality) and det(R) = 1
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
591
Fig. 4. Scatter plot matrix of two n = 50, k = 5 designs. The upper right scatter plots represent a design based on simplex meshes, the lower left ones a random latin hypercube design.
is a rotation matrix, see Meyer (2000). Rotation matrices can be constructed stepwise. Define a plane rotation matrix regarding dimension i and j as ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ Ri,j () = ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
1
⎞
... ..
⎟ ⎟ ⎟ ⎟ ⎟ ←i ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ←j ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
. cos
sin 1 ..
− sin
. cos 1 ..
↑ i
↑ j
. 1
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
592
Fig. 5. Scatter plot matrix of two n = 50, k = 5 designs. The upper right scatter plots represent a design based on simplex meshes, the lower left ones a design extracted from a rotated simplex mesh, to which the rank transformation has been applied.
with all nondiagonal entries equal to 0 except the two entries at position (i, j) and (j, i). A general rotation matrix can then be defined as the product of plane rotation matrices: R( = 1 , . . . , k(k−1)/2 ) := R1,2 (1 ) · · · R1,k (k−1 )R2,3 (k ) · · · R2,k (2k−3 ) · · · Rk−1,k (k(k−1)/2 ). The choice of the vector ∈ Rk(k−1)/2 is crucial for the improvement of the mesh. To compare rotated meshes a criterion is needed. We choose a maximin approach which can be applied to each column of P:
() = max {p ((PR())j )}, j=1, ...,k
where (PR())j is the j-th column of PR(). Minimizing this criterion means to search for the rotation which has the best one ∗ dimensional projection properties. Among all values {0 ⱕ ⱕ 2} := {0 ⱕ 1 ⱕ 2, . . . , 0 ⱕ k(k−1)/2 ⱕ 2} the value which minimizes () is searched for numerically. 3.2. Rank version of NRSM based designs As the class of latin hypercubes is very attractive due to its projection properties, the designs based on NRS meshes can be transferred into this class by applying the rank transformation to each column of the design and scaling the resulting design to
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
593
the unit cube. By numerical comparison, we show that these rank transformed designs still have good space filling properties regarding the maximin criterion but have better WD values than designs without rank transformation. If D ∈ [0, 1]n×k is a design extracted from a simplex mesh, for each column D.,d of D the n entries of the column can be ranked with numbers 1, . . . , n. Hence, if for entry (i, d) rank(Di,d ) = m ∈ {1, . . . , n}, this means, that there are m − 1 entries j in column d with Dj,d < Di,d and n − m entries j with Dj,d > Di,d . By applying the rank transformation to each column of the design D, the order of the one dimensional entries is preserved but the one-dimensional point to point distances are normalized. In Fig. 5 a scatter plot matrix is shown, which has the two dimensional projections of the simplex based design from Fig. 4 (from an unrotated mesh) at the upper right part and on the lower left part the two dimensional projections of a simplex based design from a rotated mesh with the same dimension and sample size, to which the rank transformation has been applied. Although the corresponding -criterion values are comparable, the rank design extracted from the rotated mesh has better projection properties. 4. Comparison and discussion The space filling designs constructed here will now be compared to designs constructed according to the maximin criteria
p (D ) and Mm (DMm ) and the WD criterion (DWD ). The criteria p and WD need to be minimized to get optimal designs. For
these two criteria the minimization has been done by simulated annealing within the class of latin hypercubes. The computation times for the numerically optimized designs are between 14.9 and 279.7 s (-criterion) and 17.8 and 766.8 s for the WD-criterion on a 2 GHz processor using R.2.8.1 under Windows XP. The extraction of the simplex based designs in all cases needed < 1 s, independent from dimension and sample size. The Mm criterion needs to be maximized. Here we used, where available, results explained by Grosso et al. (2009) and listed at www.spacefillingdesigns.nl. For these designs no computation times are reported. Let DR , DS be designs based on rotated meshes with and without rank transformation, respectively. As a minimum requirement would be that the constructed designs perform better than random latin hypercube designs, for each combination of dimension and sample size 1000 random latin hypercubes Drandom have been generated and the mean of the resulting criterion values is also included in the comparison. The results are shown in Table 1 (p ), Table 2 (Mm) and Table 3 (WD). For all three tables, the first and second column show the criterion values for the DS and DR designs. The third column reports the mean results of the respective criteria for 1000 random latin hypercube designs. The last column always shows criterion values for numerically optimized designs. The results from Tables 1 and 2 show that the simplex based design from rotated meshes without rank transformation (DS ) are highly competitive to numerically optimized maximin latin hypercube designs. For the 40 criterion in Table 1, all simplex based designs with and without rank transformation perform better than the latin hypercubes optimized by simulated annealing and as the mean of the random latin hypercube designs. The numerically optimized designs from Table 2 are constructed by Grosso et al. (2009) by using an iterated local search heuristic. The local search is applied to a slightly perturbed version of the latin hypercube from the last local search. The designs extracted from a rotated mesh without rank transformation perform slightly worse than the results reported by Grosso et al. (2009). The rank versions of the simplex based designs lose some efficiency w.r.t. the Mm criterion but still perform much better than the random latin hypercube designs. The results for the uniform designs with respect to the WD criterion in Table 3 show a different behavior. Here the simplex based designs without rank transformation perform worse than the mean of the criterion values for the 1000 random latin hypercube designs, indicating that the designs without rank transformation are not attractive as uniform designs (and hence are not recommendable for integral estimation). But when applying the rank transformation to the simplex based designs, the rank versions are performing only slightly worse than the numerically optimized but much better than the mean of the random latin hypercube designs. Table 1 Numerical results for the 40 -criterion.
40 (DS )
40 (DR )
40 (Drandom )
40 (D )@
3
30 60
3.214 3.976
3.394 4.582
10.570 18.140
3.787 5.269
4
40 80
2.429 3.072
2.723 3.293
7.368 11.630
2.768 3.563
5
50 100
2.058 2.467
2.237 2.557
5.556 7.916
2.327 2.993
6
60 120
1.829 2.144
1.966 2.193
4.444 5.919
2.088 2.430
7
70 140
1.642 1.883
1.723 1.919
3.634 4.668
1.794 2.185
8
80 160
1.484 1.736
1.598 1.727
3.113 3.797
1.619 1.904
k
n
The numerically optimized design has been searched for by simulated annealing within the class of latin hypercubes.
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
594 Table 2 Numerical results for the Mm-criterion. k
Mm(DS )
Mm(DR )
Mm(Drandom )
Mm(DMm )
3
n 30 60
0.3436 0.2830
0.2966 0.2249
0.1028 0.06046
0.3600 0.2800
4
40 80
0.4534 0.3715
0.3795 0.3157
0.1451 0.0937
0.4763 0.3899
5
50 100
0.5435 0.4659
0.4600 0.4194
0.1908 0.1339
0.5894 0.4949
6
60 120
0.6147 0.5402
0.5127 0.4855
0.2344 0.1776
0.6879 –
7
70 140
0.7004 0.6244
0.6021 0.5504
0.2855 0.2235
0.7798 –
8
80 160
0.7736 0.6781
0.6474 0.6201
0.3305 0.2712
0.8673 –
Where available, numerically optimized designs have been taken from www.spacefillingdesigns.nl.
Table 3 Numerical results for the uniform WD-criterion. k
WD(DS )
WD(DR )
WD(Drandom )
WD(DWD )@
3
n 30 60
0.01706 0.00709
0.00458 0.00156
0.00606 0.00244
0.00418 0.00144
4
40 80
0.02055 0.00949
0.00830 0.00293
0.01024 0.00457
0.00725 0.00262
5
50 100
0.03752 0.01483
0.01309 0.00553
0.01732 0.00808
0.01274 0.00485
6
60 120
0.05431 0.02527
0.02191 0.01010
0.02865 0.01371
0.02063 0.00851
7
70 140
0.12332 0.05574
0.04109 0.01889
0.04656 0.02254
0.03178 0.01447
8
80 160
0.17347 0.09463
0.06892 0.03039
0.07446 0.03628
0.05354 0.02618
The numerically optimized design has been searched for by simulated annealing within the class of latin hypercubes.
4.1. Numerical issues The computationally most intensive part of the algorithm is the calculation of the Delaunay triangulations. For a fixed number n of vertices, the number of simplices grows rapidly with the dimension, see Klee (1980). For dimensions larger than 7, this practically prohibits the calculation of a sufficient number of points on a personal computer. Also on servers with sufficient memory this limits the dimension up to which a sufficiently large number of points can be constructed. Here the quickhull algorithm (Barber et al., 1996) has been used to compute the Delaunay triangulation. The NRSM algorithm presented here has been implemented in `R' and Matlab. 5. Summary The designs suggested here are interesting alternatives to the fast construction of good space filling designs. The procedure to construct a design according to the methodology explained above looks as follows: First a simplex mesh with point set P ∈ Rn×k is constructed by the NRSM algorithm and the points in P are rotated such that good projection properties are achieved. These two steps only have to be done once for each dimension k. Hence it is acceptable that the construction and rotation of the point sets is time consuming. The rotated point sets for dimensions k = 3, . . . , 8 can be downloaded from www.statistik.tu-dortmund.de/ muehlenstaedt.html. The rotated point set PR can then be used repeatedly to extract designs from it. If a latin hypercube is desired, the rank transformation can be applied to the extracted design. Once the underlying simplex meshes are found and stored as a library file, the construction of a design is very simple and fast (usually < 1 s, independent from dimension and sample size). The designs we get by the methods presented here are generally competitive regarding designs optimized with respect to the maximin design criteria and Mm. The rank version of the simplex based designs are also nearly as good as numerically optimized designs regarding the WD criterion. In contrast to the simplex based designs, for all three criteria the optimization process can be very difficult and time expensive. Another advantage is that it is also very easy to construct space filling designs for nonstandard
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
595
design spaces, i.e. design spaces which are not rectangles. We only have to modify the way design points are extracted from the mesh. With other methods like uniform designs, this would be a much more difficult task. Appendix A. Calculation of the circum sphere Okabe et al. (2000, p. 80) state that in two dimensions, the circum circle with center c and radius r for a triangle can be calculated based on the following equation: 2 2 i=1 c i 2 x2 i=1 1,i 2 i=1 x22,i 2 x2 i=1 3,i
c1
c2
x1,1
x1,2
x2,1
x2,2
x3,1
x3,2
1 1 = 0. 1 1
(8)
Here (xj,1 , xj,2 )T , j = 1, 2, 3 are the coordinates of the vertices of the triangle and (c1 , c2 )T is the circum center of the circum circle of the triangle which is searched for. This equation can easily be extended to k > 2 dimensions: k 2 i=1 ci k 2 i=1 x1,i .. . k x2
c1
...
ck
x1,1
...
x1,k
xk+1,1
i=1 k+1,i
...
xk+1,k
1 1 = 0, .. . 1
(9)
where again (xj,1 , . . . , xj,k )T , j = 1, . . . , k + 1 are the coordinates of the vertices of the simplex and (c1 , . . . , ck )T are the coordinates of the circum center of the circum sphere of the simplex. In order to solve this equation cofactor expansion along the first row can be used. Define k 2 x1,1 . . . x1,j−1 x1,j+1 . . . x1,k 1 i=1 x1,i k 2 x x . . . x x . . . x 1 2,1 2,j−1 2,j+1 2,k i=1 2,i 2+j Dj := (−1) (10) .. .. . . k x2 x ... x x ... x 1 i=1 k+11,i
k+1,1
k+1,j−1
k+1,j+1
k+1,k
and x1,1 x2,1 := . .. xk+1,1
...
x1,k
...
x2,k
. . . xk+1,k
1 1 .. , . 1
k 2 i=1 x1,i k 2 i=1 x2,i := .. . k x2
i=1 k+1,i
x1,1
...
x2,1
...
xk+1,1
...
x1,k x2,k . .. . xk+1,k
(11)
With the cofactor expansion equation (8) becomes k
ci2 −
i=1
k
Dj cj + = 0.
(12)
j=1
Completing the squares yields cj = Dj /2. The radius r of the circum sphere can be calculated as r=
k 2 j=1 Dj
− 4
4 2
(13)
or alternatively just by calculating the distance of (c1 , . . . , ck )T to one of the vertices of the simplex. As the circum center is the point, which has the same distance to all vertices of the simplex, it is arbitrary, which vertex is chosen. For the second proposal function of Section 2.3, the circum center of the facet of a simplex is needed, i.e. for k vertices x0 , . . . , xj−1 , xj+1 , . . . , xk of simplex S with vertices x0 , . . . , xk . This can be deduced from the circum center of the complete simplex by (orthogonally) projecting the circum center c of S to the hyperplane defined by the k points x0 , . . . , xj−1 , xj+1 , . . . , xk , which can be done by using the Moore Penrose inverse.
596
¨ T. Muhlenst a¨ dt / Journal of Statistical Planning and Inference 140 (2010) 585 -- 596
References Barber, C., Dobkin, D., Huhdanpaa, H., 1996. The quickhull algorithm for convex hulls. ACM Transactions on Mathematical Software 22 (4), 469–483. Coxeter, H., 1963. Regular Polytopes. The Macmillan Company, New York. Fang, K., Tang, Y., Yin, J., 2005. Lower bounds for wrap-around l2 -discrepancy and construction of symmetrical uniform designs. Journal of Complexity 21, 757–771. Fang, K.-T., Li, R., Sudjianto, A., 2006. Design and Modeling for Computer Experiments. Computer Science and Data Analysis Series. Chapman & Hall, CRC, Boca Raton. Grosso, A., Jamali, A., Locatelli, M., 2009. Finding maximim latin hypercube designs by iterated local search heuristics. European Journal of Operational Research 197, 541–547. Johnson, M., Moore, L., Ylvisaker, D., 1990. Minimax and maximin distance designs. Journal of Statistical Planning and Inference 26, 131–148. Kirkpatrick, S., Gelatt, C., Vecchi, M., 1983. Optimization by simulated annealing. Science 220 (4598), 671–680. Klee, V., 1980. On the complexity of d-dimensional Voronoi diagrams. Archiv der Mathematik 34, 75–80. McKay, M., Beckman, R., Conover, W., 1979. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (2), 239–245. Meyer, C., 2000. Matrix Analysis and Applied Linear Algebra. SIAM, Philadelphia. Morris, M., Mitchell, T., 1995. Exploratory designs for computational experiments. Journal of Statistical Planning and Inference 43, 381–402. Okabe, A., Boots, B., Sugihara, K., Chiu, S., 2000. Spatial Tessellations, Concepts and Applications of Voronoi Diagrams. Wiley Series in Probability and Statistics. Wiley, Chichester. d Rajan, V., 1994. Optimality of the Delaunay triangulation in R . Discrete & Computational Geometry 12, 189–202. Santner, T., Williams, B., Notz, W., 2003. The Design and Analysis of Computer Experiments. Springer Series in Statistics. Springer, New York.