# Two General Extension Algorithms of Latin Hypercube Sampling

Authors: Zhi-zhao Liu; Wei Li; Ming Yang
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

For reserving original sampling points to reduce the simulation runs, two general extension algorithms of Latin Hypercube Sampling (LHS) are proposed. The extension algorithms start with an original LHS of sizem and construct a new LHS of size m
+
n that contains the original points as many as possible. In order to get a strict LHS of larger size, some original points might be deleted. The relationship of original sampling points in the new LHS structure is shown by a simple undirected acyclic graph. The basic general extension algorithm is proposed to reserve the most original points, but it costs too much time. Therefore, a general extension algorithm based on greedy algorithm is proposed to reduce the extension time, which cannot guarantee to contain the most original points. These algorithms are illustrated by an example and applied to evaluating the sample means to demonstrate the effectiveness.


---

## Body

## 1. Introduction

Latin Hypercube Sampling (LHS) is one of the most popular sampling approaches, which is widely used in the fields of simulation experiment design [1], uncertainty analysis [2], adaptive metamodeling [3], reliability analysis [4], and probabilistic load flow analysis [5]. Compared with other random or stratified sampling algorithms, LHS has a better space filling effect, better robustness, and better convergence character. The extension of LHS is to obtain a LHS of a larger size that reserves the preexisting LHS (or the original LHS). There are at least two situations that need the extension of sampling, especially for time consuming simulation systems. One is sequential sampling for sequential analysis, adaptive metamodeling, and so on. The other is to consider the extension of LHS when the original LHS was subsequently determined to be too small and a new LHS of a larger size without original sampling points might be time consuming. But the LHS structure makes it difficult to increase the size based on an original LHS while simultaneously keeping the stratification properties of LHS.A special extension case is the integral-multiple extension where the new LHS is integral times the size of the original sampling. Tong [6] proposed integral-multiple extension algorithms for stratified sampling methods including LHS. Sallaberry et al. [7] gave a two-multiple extension algorithm of LHS with correlated variables. Later, two related techniques appeared in the papers named “nested Latin hypercube design” [8] and “nested orthogonal array-based Latin hypercube design” [9]. A nested Latin hypercube design with two layers is defined to be a Latin hypercube design that contains a smaller Latin hypercube design as a subset. A special integral-multiple extension method called k-extended LHS method was illustrated, where the new LHS contains k smaller LHSs [10]. Some related papers were produced by Vorechovsky [11–13], where the new sampling size is multiple times more than the original sampling size. The integral-multiple extension algorithms have a good feature that can obtain a strict LHS of larger size and simultaneously preserve all the original sampling points.In this study, we consider the general extension algorithm of LHS where the new sample size is more controllable and the algorithm can be applied more widely. Wang [14] and Blatman and Sudret [15] obtained an approximate LHS of a larger size, which might have two or more original points falling into the same variable interval. Wei [16] also proposed a general extension algorithm to get an approximate LHS, which might have no point falling into a variable interval. However, a sample is a LHS if (and only if) there is only one point in each variable interval. The approximate LHS does not satisfy the definition and is harmful to the extension with some criteria, such as correlated variables, maximizing minimum distance, orthogonal array, and so on.In this paper, we would like to obtain a strict extension of LHS (ELHS) rather than an approximate one and the new LHS contains original sampling points as many as possible. The extension algorithm includes two parts: the reservation of original sampling points and the generation of new ones. As the generation of new sampling points is almost the same as integral-multiple extension, the reservation of original sampling points is the main problem to discuss. The relationship of original sampling points in new LHS is expressed by a graph. Then, the reservation problem can be solved by graph theory.In Section2, the procedure of LHS and the mathematical description for extension problem of LHS are given, which is the basis of extension algorithm. In Section 3, a basic general extension algorithm of LHS (BGELHS) is proposed based on the graph theory, which reserves the most original sampling points. A general extension algorithm of LHS based on greedy algorithm (GGELHS) is intended to balance the simulation runs and the generation time of ELHS. In Section 4, the proposed extension algorithms are illustrated by an application, which is performed to evaluate sample means. Finally, the conclusion and some thoughts for the future research are given.


## 2. Latin Hypercube Sampling and Extension Problem

### 2.1. Procedure of LHS

Suppose that the input variables arex
1
,
x
2
,
…
,
x
r and the range of x
i is u
i
=
a
i
b
i, i
=
1,2
,
…
,
r. Then, LHS can be obtained as follows, which can ensure that each input variable has all portions among its range.(a)
Divide the rangeu
i into n equiprobable intervals u
i
1
,
u
i
2
,
…
,
u
i
n, i
=
1,2
,
…
,
r. So the intervals satisfy u
i
1
∪
u
i
2
,
…
,
∪
u
i
n
=
u
i, u
i
j
∩
u
i
k
=
⌀, and P
(
x
∈
u
i
j
)
=
1
/
n, where j
,
k
=
1,2
,
…
,
n.
(b)
For thejth interval of variable x
i, the cumulative probability can be obtained as(1)
P
r
o
b
j
i
=
j
-
1
n
+
r
j
i
n
,where r
j
i is a uniform random number ranging from 0 to 1. So all the probability values can be noted as P
r
o
b
=
(
P
r
o
b
j
i
)
n
×
r.
(c)
Transform the probability into the sample valuex
j
i by the inverse of the cumulative distribution function F
(
·
):(2)
x
j
i
=
F
i
-
1
P
r
o
b
j
i
.Then, the sample matrix is (3)
X
=
x
11
x
12
⋯
x
1
r
x
21
x
22
⋯
x
2
r
⋮
⋮
⋱
⋮
x
n
1
x
n
2
⋯
x
n
r
.
(d)
Then values of each variable are paired randomly or in some prescribed order with the n values of the other variables. Then the sample matrix of LHS can be written as(4)
X
′
=
x
11
′
x
12
′
⋯
x
1
r
′
x
21
′
x
22
′
⋯
x
2
r
′
⋮
⋮
⋱
⋮
x
n
1
′
x
n
2
′
⋯
x
n
r
′
,where each row is a sampling point.Figure1 shows an example of LHS in size of 5, where each interval of each input variable has one sampling point.Figure 1
An example of LHS.


### 2.2. Extension Problem of LHS

In this paper, the extension problem of LHS has two purposes. One is to construct a strict LHS of a larger size than original LHS. The other is to reserve the most original sampling points in the new LHS. Assume that there is an original LHS and the sample set isA
=
{
X
1
,
X
2
,
…
,
X
n
}, where X
i is a sampling point, i
=
1,2
,
…
,
n. The objective of extension is to get ELHS of size m
+
n noted as B
=
{
X
1
′
,
X
2
′
,
…
,
X
n
+
m
′
}, where m
≥
1 and m is an integer. B satisfies the constraint condition as shown in(5)
max
⁡
c
a
r
d
A
∩
B
,where c
a
r
d
(
·
) denotes the number of elements in set. It means that the objective is to obtain ELHS in size of n
+
m that contains the original LHS points as many as possible. Figure 2 shows the original LHS points in Figure 1 and the structure of ELHS in size of 7, which can be seen that the original sampling points do not accord with the new structure. So we need to delete original sampling points as little as possible and then create some new sampling points. The reserved points and the new ones are ELHS in size of 7.Figure 2
Distribution of generated sampling points in new structure.


## 2.1. Procedure of LHS

Suppose that the input variables arex
1
,
x
2
,
…
,
x
r and the range of x
i is u
i
=
a
i
b
i, i
=
1,2
,
…
,
r. Then, LHS can be obtained as follows, which can ensure that each input variable has all portions among its range.(a)
Divide the rangeu
i into n equiprobable intervals u
i
1
,
u
i
2
,
…
,
u
i
n, i
=
1,2
,
…
,
r. So the intervals satisfy u
i
1
∪
u
i
2
,
…
,
∪
u
i
n
=
u
i, u
i
j
∩
u
i
k
=
⌀, and P
(
x
∈
u
i
j
)
=
1
/
n, where j
,
k
=
1,2
,
…
,
n.
(b)
For thejth interval of variable x
i, the cumulative probability can be obtained as(1)
P
r
o
b
j
i
=
j
-
1
n
+
r
j
i
n
,where r
j
i is a uniform random number ranging from 0 to 1. So all the probability values can be noted as P
r
o
b
=
(
P
r
o
b
j
i
)
n
×
r.
(c)
Transform the probability into the sample valuex
j
i by the inverse of the cumulative distribution function F
(
·
):(2)
x
j
i
=
F
i
-
1
P
r
o
b
j
i
.Then, the sample matrix is (3)
X
=
x
11
x
12
⋯
x
1
r
x
21
x
22
⋯
x
2
r
⋮
⋮
⋱
⋮
x
n
1
x
n
2
⋯
x
n
r
.
(d)
Then values of each variable are paired randomly or in some prescribed order with the n values of the other variables. Then the sample matrix of LHS can be written as(4)
X
′
=
x
11
′
x
12
′
⋯
x
1
r
′
x
21
′
x
22
′
⋯
x
2
r
′
⋮
⋮
⋱
⋮
x
n
1
′
x
n
2
′
⋯
x
n
r
′
,where each row is a sampling point.Figure1 shows an example of LHS in size of 5, where each interval of each input variable has one sampling point.Figure 1
An example of LHS.


## 2.2. Extension Problem of LHS

In this paper, the extension problem of LHS has two purposes. One is to construct a strict LHS of a larger size than original LHS. The other is to reserve the most original sampling points in the new LHS. Assume that there is an original LHS and the sample set isA
=
{
X
1
,
X
2
,
…
,
X
n
}, where X
i is a sampling point, i
=
1,2
,
…
,
n. The objective of extension is to get ELHS of size m
+
n noted as B
=
{
X
1
′
,
X
2
′
,
…
,
X
n
+
m
′
}, where m
≥
1 and m is an integer. B satisfies the constraint condition as shown in(5)
max
⁡
c
a
r
d
A
∩
B
,where c
a
r
d
(
·
) denotes the number of elements in set. It means that the objective is to obtain ELHS in size of n
+
m that contains the original LHS points as many as possible. Figure 2 shows the original LHS points in Figure 1 and the structure of ELHS in size of 7, which can be seen that the original sampling points do not accord with the new structure. So we need to delete original sampling points as little as possible and then create some new sampling points. The reserved points and the new ones are ELHS in size of 7.Figure 2
Distribution of generated sampling points in new structure.


## 3. General Extension Algorithms of LHS

The general extension problem can be described as (5). In this section, two general extension algorithms of LHS are demonstrated to get ELHS. The problem of selecting reserved points is transformed into a problem to get the maximum independent set, which is a graph theory problem. The generation procedure is based on a subtraction rule [14]. These two stages construct BGELHS in Section 3.1. As BGELHS costs too much time to obtain ELHS, GGELHS is proposed in Section 3.2.


### 3.1. Basic General Extension Algorithm of LHS

Some original sampling points might be deleted before generating new sample. At the beginning of BGELHS, the original sample is transformed into a graph. Here the definition of simple graph is introduced.Definition A (see [17]).
A simple graphG
=
(
V
,
E
) is a graph where there is no loop or multiple edges, where V
(
G
) is a vertex set and E
(
G
) is an edge set. No loop means that there is no edge from a vertex to itself. No multiple edges means that there is 0 or 1 edge between two vertices.
Let the original sample matrix beX
n
×
r
(
o
l
d
) and note the sampling points as 1,2
,
…
,
n by the row orders. Transform the sampling points into vertices of a graph as follows: Every point is a vertex and links two vertices by a line if the corresponding two sampling points are in the same interval of the ELHS structure. Then we get a simple undirected graph. Order the sampling points in Figure 2 from left to right and the corresponding simple undirected graph is obtained as shown in Figure 3. The original problem in (5) is to delete the least vertices and their corresponding lines so that the remaining vertices have no line to link, which can be proved by reduction. Independent set, maximal independent set, and maximum independent set are introduced to solve this problem.Figure 3
The simple undirected acyclic graph of Figure2.Definition B (see [17]).
It is supposed that there is a simple graphG
=
(
V
,
E
) where V is the set of vertices and E is the set of edges. If a set of vertices I
⊆
V
(
G
) and there is no edge between any two vertices of I, the set I is called an independent set.Definition C (see [17]).
It is supposed that there is a simple graphG
=
(
V
,
E
) where V is the set of vertices and E is the set of edges. I and I
′ are independent sets. If ∀
v
i
∈
V and v
i
∉
I, I
′
=
I
∪
{
v
i
} is not an independent set. Then, the set I is the maximal independent set. If ∀
I
′, c
a
r
d
(
I
)
≥
c
a
r
d
(
I
′
). Then the set I is the maximum independent set. The number of elements in I is the independent number of graph G noted as α
(
G
) or α.
From the definitions, reserving the most original sampling points can be transformed into computing the maximum independent set; that is, delete the least vertices so that there is no edge between any two vertices. In graph theory, the simple graph can be expressed as an adjacent matrixA
n
×
n
=
(
a
i
j
), which shows the lines between vertices. If there is a line between vertex v
i and vertex v
j, a
i
j
=
1. Otherwise a
i
j
=
0. Obviously, the properties are given as follows:(a)
symmetrical characteristicA
n
×
n
=
A
n
×
n
T;
(b)
a
i
i
=
0, i
=
1,2
,
…
,
n;
(c)
ifA
n
×
n
=
0
→, the vertex set V is an independent set.So the reservation problem is transformed into a problem to get the maximum independent set; that is, delete the least vertices so that the corresponding adjacent matrix of remaining vertices is a zero matrix. Suppose that there is a LHS of size n, the sample matrix is X
n
×
r
(
o
l
d
)
=
{
x
i
j
(
o
l
d
)
}, the order matrix is L
n
×
r
(
o
l
d
)
=
{
l
i
j
(
o
l
d
)
} (l
i
j
(
o
l
d
) is the interval order number of x
i
j
(
o
l
d
)in the jth variable), the cumulative probability matrix is U
n
×
r
(
o
l
d
)
=
{
u
i
j
(
o
l
d
)
} (u
i
j
(
o
l
d
) is the cumulative probability of x
i
j
(
o
l
d
), which is the P
r
o
b
j
i in (1)), and the random matrix is R
n
×
r
(
o
l
d
)
=
{
r
i
j
(
o
l
d
)
} (r
i
j
(
o
l
d
) is the random value of x
i
j
(
o
l
d
), which is the r
j
i in (1)). The objective is to get ELHS of size n
+
m noted as X
(
n
+
m
)
×
r
(
n
e
w
). The detailed steps of BGELHS algorithm are as follows.
Step 1. Compute the positions of original sampling points in ELHS structure.(a)
TransformL
n
×
r
(
o
l
d
) into the new sample structure and get the order matrix L
n
×
r
(
1
) as shown in the following:(6)
l
i
,
j
1
=
u
i
,
j
o
l
d
·
n
+
m
+
1
,where i
=
1,2
,
…
,
n, j
=
1,2
,
…
,
r, and · is the floor function (f
l
o
o
r
(
x
)
=
x is the largest integer not greater than x).
(b)
TransformR
n
×
r
(
1
) into the new sample structure and get the random matrix R
n
×
r
(
1
) as shown in the following:(7)
r
i
,
j
1
=
u
i
,
j
o
l
d
·
n
+
m
-
l
i
,
j
1
-
1
,where i
=
1,2
,
…
,
n, j
=
1,2
,
…
,
r.
Step 2. Reserve the most original sampling points, noted as the sample matrix X
n
′
×
r
(
2
).(a)
Compute the adjacent matrixA
n
×
n of the original LHS based on L
n
×
r
(
1
) as shown in the following:(8)
a
i
j
=
1
,
If
∃
p
∈
1
n
,
p
∈
Z
,
l
i
p
1
=
l
j
p
1
0
,
If
∀
p
∈
1
n
,
p
∈
Z
,
l
i
p
1
≠
l
j
p
1
.
(b)
Choosek sampling points of original LHS to delete, k
=
0,1
,
2
,
…
n. Then we can get a combination matrix B
q
×
n
=
B
1
T
B
2
T
⋯
B
q
T
T, where q
=
C
n
0
+
C
n
1
+
C
n
2
+
⋯
+
C
n
n
=
2
n, B
i is a combination, and the number of elements in B
i is not more than that in B
j, 1
≤
i
<
j
≤
q (e.g., let r
=
2 and n
=
3, B
4
×
2
=
n
u
l
l
1
2
1
n
u
l
l
n
u
l
l
n
u
l
l
2
T, where n
u
l
l means there is no element). Let p
=
1 and delete the point(s) corresponding with the element(s) of B
p; then judge whether the remaining adjacent matrix A
′ is equal to 0
→. If A
′
=
0
→, record the order number of deleted point(s) as B
p
′ and the dimension number of A
′ is noted as n
′, that is, the number of remaining points. Otherwise, set p
=
p
+
1, delete the points corresponding with the elements of B
p, and repeat the Step 2-(b).
(c)
Delete the rows ofL
n
×
r
(
1
), R
n
×
r
(
1
), U
n
×
r
(
o
l
d
), and X
n
×
r
(
o
l
d
) according to the element of B
p
′. The remaining matrices are noted as L
n
′
×
r
(
2
), R
n
′
×
r
(
2
), U
n
′
×
r
(
2
), and X
n
′
×
r
(
2
).
Step 3. Generate the new sample matrix X
(
n
+
m
-
n
′
)
×
r
(
3
).(a)
Generate a new order matrixL
n
+
m
×
r
′ of size (
n
+
m
)
×
r and every column of L
n
+
m
×
r
′ is a random permutation of the integer number from 1 to n
+
m.
(b)
For every column, delete the element that is the same as the corresponding column inL
n
′
×
r
(
2
) and get the matrix L
(
n
+
m
-
n
′
)
×
r
(
3
).
(c)
Generate a random matrixR
(
n
+
m
-
n
′
)
×
r
(
3
) where r
i
j is a uniform random number ranging from 0 to 1.
(d)
Compute the cumulative probability matrixU
(
n
+
m
-
n
′
)
×
r
(
3
) as shown in the following:(9)
u
i
,
j
3
=
l
i
,
j
3
-
1
+
r
i
,
j
3
n
+
m
.
(e)
Compute the sample matrixX
(
n
+
m
-
n
′
)
×
r
(
3
) as shown in (2).
Step 4. Get the order matrix L
(
n
+
m
)
×
r
(
n
e
w
), the random matrix R
(
n
+
m
)
×
r
(
n
e
w
), the cumulative probability matrix U
(
n
+
m
)
×
r
(
n
e
w
), and the sample matrix X
(
n
+
m
)
×
r
(
n
e
w
) of ELHS:(10)
L
n
+
m
×
r
n
e
w
=
L
n
′
×
r
2
L
n
+
m
-
n
′
×
r
3
,
R
n
+
m
×
r
n
e
w
=
R
n
′
×
r
2
R
n
+
m
-
n
′
×
r
3
,
U
n
+
m
×
r
n
e
w
=
U
n
′
×
r
2
U
n
+
m
-
n
′
×
r
3
,
X
n
+
m
×
r
n
e
w
=
X
n
′
×
r
2
X
n
+
m
-
n
′
×
r
3
.
As taking (7) into (1) as (11), we can prove that the matrix X
(
n
+
m
)
×
r
(
n
e
w
) is a sample matrix of LHS. Equation (11) shows that the remaining original sampling points meet the basic principle of LHS. As the new sampling points have the same feature, X
(
n
+
m
)
×
r
(
n
e
w
) is a strict LHS:(11)
u
i
,
j
1
=
l
i
,
j
1
-
1
+
r
i
,
j
1
n
+
m
=
l
i
,
j
1
-
1
+
u
i
,
j
o
l
d
·
n
+
m
-
l
i
,
j
1
-
1
n
+
m
=
u
i
,
j
o
l
d
.Example 1.
BGELHS is illustrated by the original LHS in Figure1. Let x
1 have a uniform distribution on [
0,10
] and let x
2 have a triangular distribution on [
5,10
] with mode at 7.5. The distribution function F
(
·
) is(12)
F
1
x
1
=
0.1
×
x
1
,
F
2
x
2
=
x
2
-
5
2
12.5
,
5
≤
x
2
≤
7.5
1
-
10
-
x
2
2
12.5
,
7.5
<
x
2
≤
10
.
The inverse of the distribution functionF
(
·
) is(13)
x
1
=
10
×
F
1
x
1
,
x
2
=
5
+
12.5
×
F
2
x
2
,
0
≤
F
2
x
2
≤
0.5
10
-
12.5
×
1
-
F
2
x
2
,
0.5
<
F
2
x
2
≤
1
.
The sample matrix and related matrices are(14)
X
5
×
2
o
l
d
=
1.29
6.12
3.57
8.46
5.86
7.60
6.57
8.30
9.00
6.90
,
L
5
×
2
o
l
d
=
1
1
2
5
3
3
4
4
5
2
,
R
5
×
2
o
l
d
=
0.645
0.500
0.785
0.070
0.930
0.710
0.285
0.860
0.500
0.430
,
U
5
×
2
o
l
d
=
0.129
0.100
0.357
0.814
0.586
0.542
0.657
0.772
0.900
0.286
.
Then, the objective is to get ELHS of size 7, which contains the original sampling points as many as possible. Son
=
5, m
=
2.
From Step 1, the distribution of original sampling points in the new sample is(15)
L
5
×
2
1
=
U
5
×
2
o
l
d
×
5
+
2
+
I
5
×
2
=
0.903
0.700
2.499
5.698
4.102
3.794
4.599
5.404
6.300
2.002
+
I
5
×
2
=
0
0
2
5
4
3
4
5
6
2
+
I
5
×
2
=
1
1
3
6
5
4
5
6
7
3
,
R
5
×
2
1
=
U
5
×
2
o
l
d
·
5
+
2
-
L
5
×
2
1
-
I
5
×
2
=
0.903
0.700
0.499
0.698
0.102
0.794
0.599
0.404
0.300
0.002
.
From Step 2-(a), the adjacent matrix is computed:(16)
A
5
×
5
=
0
0
0
0
0
0
0
0
1
0
0
0
0
1
0
0
1
1
0
0
0
0
0
0
0
.
From Step 2-(b), all the combinations are noted asB:(17)
B
=
n
u
l
l
n
u
l
l
n
u
l
l
n
u
l
l
n
u
l
l
1
n
u
l
l
n
u
l
l
n
u
l
l
n
u
l
l
2
n
u
l
l
n
u
l
l
n
u
l
l
n
u
l
l
3
n
u
l
l
n
u
l
l
n
u
l
l
n
u
l
l
⋮
⋮
⋮
⋮
⋮
2
3
4
5
n
u
l
l
1
2
3
4
5
32
×
5
.
Delete the point(s) corresponding with the element(s) ofB
p, p
=
1,2
,
…
,
32. For the fifth iteration, p
=
5 and B
5
=
4
n
u
l
l
n
u
l
l
n
u
l
l
n
u
l
l. Delete the fourth point (i.e., delete the fourth row and fourth column of A
5
×
5) and the new adjacent matrix is A
4
×
4
′
=
0
→. Then, B
4
′
=
1
2
3
5 and n
′
=
4. From Step 2-(c), the addition sample is computed:(18)
L
4
×
2
2
=
1
1
3
6
5
4
7
3
,
R
4
×
2
2
=
0.903
0.700
0.499
0.698
0.102
0.794
0.300
0.002
,
U
4
×
2
2
=
0.129
0.100
0.357
0.814
0.586
0.542
0.900
0.286
,
X
4
×
2
2
=
1.29
6.12
3.57
8.46
5.86
7.60
9.00
6.90
.
From Step 3, an order matrix is original:(19)
L
7
×
2
′
=
6
3
7
5
1
2
4
4
7
1
6
2
3
5
T
.
Then(20)
L
3
×
2
3
=
6
7
2
2
4
5
,
R
3
×
2
3
=
0.800
0.916
0.679
0.734
0.743
0.392
,
U
3
×
2
3
=
L
3
×
2
3
-
1
→
+
R
3
×
2
3
7
=
0.829
0.988
0.240
0.248
0.535
0.628
,
X
3
×
2
3
=
8.29
9.61
2.40
6.76
5.35
7.84
.
Finally, the ELHS of size 7 is obtained, which is shown in Figure4:(21)
X
7
×
2
n
e
w
=
X
4
×
2
2
X
3
×
2
3
=
1.29
3.57
5.86
9.00
8.29
2.40
5.35
6.12
8.46
7.60
6.90
9.61
6.76
7.84
T
.Figure 4
ELHS of size 7.


### 3.2. General Extension Algorithm of LHS Based on Greedy Algorithm

In graph theory, obtaining the maximum independent set is a NP-hard problem. Step 2-(b) of BGELHS needs the most time and the time complexity of BGELHS isO
(
2
n
). It means that, with the increase of the number of sampling points, the time of BGELHS will be increased exponential. In this section, GGELHS is proposed to reduce the time of extension algorithm. First of all, the definition of degree is introduced.Definition D (see [17]).
In a simple graphG
=
(
V
,
E
), the number of edges associated with the vertex v is the degree of v noted as d
(
v
).
d
(
v
) indicates the number of vertices linked with the vertex v. The basic idea of greedy algorithm is to delete a point one time whose d
(
v
) is the maximum until d
(
v
) of any remaining vertex is 0. The steps of GGELHS are the same as those of BGELHS except for Step 2-(b). The details of Step 2-(b) are given as follows.
Step 2-(b): compute the degree of every vertexd
(
v
i
)
=
∑
j
=
1
n
a
i
j, j
=
1,2
,
…
,
n. Delete the vertex v
i whose degree is the maximum. Delete the corresponding row and column in the adjacent matrix A (A
n
×
n at the beginning) and note the remaining matrix as A
′. Record the order of v
i in B
p. If A
′
≠
0
→, A
=
A
′, and repeat Step 2-(b) to reconstruct A
′. Otherwise, the element of B
p is the order of points to delete. The dimension of A
′ is noted as n
′ and the order of points to remain is recorded in B
p
′.
In Step 2-(a) of GGELHS, when computing the adjacent matrixA
n
×
n, the time required is n
(
n
-
1
)
/
2. It is easy to know that the time complexities of other steps are not greater than O
(
n
). So the time complexity of GGELHS is O
(
n
2
), which is much less than BGELHS. But this algorithm gets the maximal independent set rather than the maximum independent set; that is, this algorithm cannot guarantee to reserve the most original sampling points.Example 2.
To compare BGELHS and GGELHS, we consider the same question of Example1. Since GGELHS is almost the same as BGELHS except for Step 2-(b), we mainly give Step 2-(b) of GGELHS.
From Step 2-(b), the degrees of five points areD
5
=
0
1
1
2
0
T. The fourth point has the maximum degree. So delete the fourth point, and the adjacent matrix is A
4
×
4
′
=
0
→. So, B
p
′
=
1
2
3
5, n
′
=
4.
BGELHS run five times in Step 2-(b) and GGELHS run only one time. As computing the degree of original sampling points costs a little, GGELHS is less time consuming and more practical. From the result above, we can know that Step 2-(b) gets the maximum independent set in this example, which means that this algorithm deletes the least original sampling points. But the greedy algorithm cannot guarantee to have the same effect for every circumstance.


## 3.1. Basic General Extension Algorithm of LHS

Some original sampling points might be deleted before generating new sample. At the beginning of BGELHS, the original sample is transformed into a graph. Here the definition of simple graph is introduced.Definition A (see [17]).
A simple graphG
=
(
V
,
E
) is a graph where there is no loop or multiple edges, where V
(
G
) is a vertex set and E
(
G
) is an edge set. No loop means that there is no edge from a vertex to itself. No multiple edges means that there is 0 or 1 edge between two vertices.
Let the original sample matrix beX
n
×
r
(
o
l
d
) and note the sampling points as 1,2
,
…
,
n by the row orders. Transform the sampling points into vertices of a graph as follows: Every point is a vertex and links two vertices by a line if the corresponding two sampling points are in the same interval of the ELHS structure. Then we get a simple undirected graph. Order the sampling points in Figure 2 from left to right and the corresponding simple undirected graph is obtained as shown in Figure 3. The original problem in (5) is to delete the least vertices and their corresponding lines so that the remaining vertices have no line to link, which can be proved by reduction. Independent set, maximal independent set, and maximum independent set are introduced to solve this problem.Figure 3
The simple undirected acyclic graph of Figure2.Definition B (see [17]).
It is supposed that there is a simple graphG
=
(
V
,
E
) where V is the set of vertices and E is the set of edges. If a set of vertices I
⊆
V
(
G
) and there is no edge between any two vertices of I, the set I is called an independent set.Definition C (see [17]).
It is supposed that there is a simple graphG
=
(
V
,
E
) where V is the set of vertices and E is the set of edges. I and I
′ are independent sets. If ∀
v
i
∈
V and v
i
∉
I, I
′
=
I
∪
{
v
i
} is not an independent set. Then, the set I is the maximal independent set. If ∀
I
′, c
a
r
d
(
I
)
≥
c
a
r
d
(
I
′
). Then the set I is the maximum independent set. The number of elements in I is the independent number of graph G noted as α
(
G
) or α.
From the definitions, reserving the most original sampling points can be transformed into computing the maximum independent set; that is, delete the least vertices so that there is no edge between any two vertices. In graph theory, the simple graph can be expressed as an adjacent matrixA
n
×
n
=
(
a
i
j
), which shows the lines between vertices. If there is a line between vertex v
i and vertex v
j, a
i
j
=
1. Otherwise a
i
j
=
0. Obviously, the properties are given as follows:(a)
symmetrical characteristicA
n
×
n
=
A
n
×
n
T;
(b)
a
i
i
=
0, i
=
1,2
,
…
,
n;
(c)
ifA
n
×
n
=
0
→, the vertex set V is an independent set.So the reservation problem is transformed into a problem to get the maximum independent set; that is, delete the least vertices so that the corresponding adjacent matrix of remaining vertices is a zero matrix. Suppose that there is a LHS of size n, the sample matrix is X
n
×
r
(
o
l
d
)
=
{
x
i
j
(
o
l
d
)
}, the order matrix is L
n
×
r
(
o
l
d
)
=
{
l
i
j
(
o
l
d
)
} (l
i
j
(
o
l
d
) is the interval order number of x
i
j
(
o
l
d
)in the jth variable), the cumulative probability matrix is U
n
×
r
(
o
l
d
)
=
{
u
i
j
(
o
l
d
)
} (u
i
j
(
o
l
d
) is the cumulative probability of x
i
j
(
o
l
d
), which is the P
r
o
b
j
i in (1)), and the random matrix is R
n
×
r
(
o
l
d
)
=
{
r
i
j
(
o
l
d
)
} (r
i
j
(
o
l
d
) is the random value of x
i
j
(
o
l
d
), which is the r
j
i in (1)). The objective is to get ELHS of size n
+
m noted as X
(
n
+
m
)
×
r
(
n
e
w
). The detailed steps of BGELHS algorithm are as follows.
Step 1. Compute the positions of original sampling points in ELHS structure.(a)
TransformL
n
×
r
(
o
l
d
) into the new sample structure and get the order matrix L
n
×
r
(
1
) as shown in the following:(6)
l
i
,
j
1
=
u
i
,
j
o
l
d
·
n
+
m
+
1
,where i
=
1,2
,
…
,
n, j
=
1,2
,
…
,
r, and · is the floor function (f
l
o
o
r
(
x
)
=
x is the largest integer not greater than x).
(b)
TransformR
n
×
r
(
1
) into the new sample structure and get the random matrix R
n
×
r
(
1
) as shown in the following:(7)
r
i
,
j
1
=
u
i
,
j
o
l
d
·
n
+
m
-
l
i
,
j
1
-
1
,where i
=
1,2
,
…
,
n, j
=
1,2
,
…
,
r.
Step 2. Reserve the most original sampling points, noted as the sample matrix X
n
′
×
r
(
2
).(a)
Compute the adjacent matrixA
n
×
n of the original LHS based on L
n
×
r
(
1
) as shown in the following:(8)
a
i
j
=
1
,
If
∃
p
∈
1
n
,
p
∈
Z
,
l
i
p
1
=
l
j
p
1
0
,
If
∀
p
∈
1
n
,
p
∈
Z
,
l
i
p
1
≠
l
j
p
1
.
(b)
Choosek sampling points of original LHS to delete, k
=
0,1
,
2
,
…
n. Then we can get a combination matrix B
q
×
n
=
B
1
T
B
2
T
⋯
B
q
T
T, where q
=
C
n
0
+
C
n
1
+
C
n
2
+
⋯
+
C
n
n
=
2
n, B
i is a combination, and the number of elements in B
i is not more than that in B
j, 1
≤
i
<
j
≤
q (e.g., let r
=
2 and n
=
3, B
4
×
2
=
n
u
l
l
1
2
1
n
u
l
l
n
u
l
l
n
u
l
l
2
T, where n
u
l
l means there is no element). Let p
=
1 and delete the point(s) corresponding with the element(s) of B
p; then judge whether the remaining adjacent matrix A
′ is equal to 0
→. If A
′
=
0
→, record the order number of deleted point(s) as B
p
′ and the dimension number of A
′ is noted as n
′, that is, the number of remaining points. Otherwise, set p
=
p
+
1, delete the points corresponding with the elements of B
p, and repeat the Step 2-(b).
(c)
Delete the rows ofL
n
×
r
(
1
), R
n
×
r
(
1
), U
n
×
r
(
o
l
d
), and X
n
×
r
(
o
l
d
) according to the element of B
p
′. The remaining matrices are noted as L
n
′
×
r
(
2
), R
n
′
×
r
(
2
), U
n
′
×
r
(
2
), and X
n
′
×
r
(
2
).
Step 3. Generate the new sample matrix X
(
n
+
m
-
n
′
)
×
r
(
3
).(a)
Generate a new order matrixL
n
+
m
×
r
′ of size (
n
+
m
)
×
r and every column of L
n
+
m
×
r
′ is a random permutation of the integer number from 1 to n
+
m.
(b)
For every column, delete the element that is the same as the corresponding column inL
n
′
×
r
(
2
) and get the matrix L
(
n
+
m
-
n
′
)
×
r
(
3
).
(c)
Generate a random matrixR
(
n
+
m
-
n
′
)
×
r
(
3
) where r
i
j is a uniform random number ranging from 0 to 1.
(d)
Compute the cumulative probability matrixU
(
n
+
m
-
n
′
)
×
r
(
3
) as shown in the following:(9)
u
i
,
j
3
=
l
i
,
j
3
-
1
+
r
i
,
j
3
n
+
m
.
(e)
Compute the sample matrixX
(
n
+
m
-
n
′
)
×
r
(
3
) as shown in (2).
Step 4. Get the order matrix L
(
n
+
m
)
×
r
(
n
e
w
), the random matrix R
(
n
+
m
)
×
r
(
n
e
w
), the cumulative probability matrix U
(
n
+
m
)
×
r
(
n
e
w
), and the sample matrix X
(
n
+
m
)
×
r
(
n
e
w
) of ELHS:(10)
L
n
+
m
×
r
n
e
w
=
L
n
′
×
r
2
L
n
+
m
-
n
′
×
r
3
,
R
n
+
m
×
r
n
e
w
=
R
n
′
×
r
2
R
n
+
m
-
n
′
×
r
3
,
U
n
+
m
×
r
n
e
w
=
U
n
′
×
r
2
U
n
+
m
-
n
′
×
r
3
,
X
n
+
m
×
r
n
e
w
=
X
n
′
×
r
2
X
n
+
m
-
n
′
×
r
3
.
As taking (7) into (1) as (11), we can prove that the matrix X
(
n
+
m
)
×
r
(
n
e
w
) is a sample matrix of LHS. Equation (11) shows that the remaining original sampling points meet the basic principle of LHS. As the new sampling points have the same feature, X
(
n
+
m
)
×
r
(
n
e
w
) is a strict LHS:(11)
u
i
,
j
1
=
l
i
,
j
1
-
1
+
r
i
,
j
1
n
+
m
=
l
i
,
j
1
-
1
+
u
i
,
j
o
l
d
·
n
+
m
-
l
i
,
j
1
-
1
n
+
m
=
u
i
,
j
o
l
d
.Example 1.
BGELHS is illustrated by the original LHS in Figure1. Let x
1 have a uniform distribution on [
0,10
] and let x
2 have a triangular distribution on [
5,10
] with mode at 7.5. The distribution function F
(
·
) is(12)
F
1
x
1
=
0.1
×
x
1
,
F
2
x
2
=
x
2
-
5
2
12.5
,
5
≤
x
2
≤
7.5
1
-
10
-
x
2
2
12.5
,
7.5
<
x
2
≤
10
.
The inverse of the distribution functionF
(
·
) is(13)
x
1
=
10
×
F
1
x
1
,
x
2
=
5
+
12.5
×
F
2
x
2
,
0
≤
F
2
x
2
≤
0.5
10
-
12.5
×
1
-
F
2
x
2
,
0.5
<
F
2
x
2
≤
1
.
The sample matrix and related matrices are(14)
X
5
×
2
o
l
d
=
1.29
6.12
3.57
8.46
5.86
7.60
6.57
8.30
9.00
6.90
,
L
5
×
2
o
l
d
=
1
1
2
5
3
3
4
4
5
2
,
R
5
×
2
o
l
d
=
0.645
0.500
0.785
0.070
0.930
0.710
0.285
0.860
0.500
0.430
,
U
5
×
2
o
l
d
=
0.129
0.100
0.357
0.814
0.586
0.542
0.657
0.772
0.900
0.286
.
Then, the objective is to get ELHS of size 7, which contains the original sampling points as many as possible. Son
=
5, m
=
2.
From Step 1, the distribution of original sampling points in the new sample is(15)
L
5
×
2
1
=
U
5
×
2
o
l
d
×
5
+
2
+
I
5
×
2
=
0.903
0.700
2.499
5.698
4.102
3.794
4.599
5.404
6.300
2.002
+
I
5
×
2
=
0
0
2
5
4
3
4
5
6
2
+
I
5
×
2
=
1
1
3
6
5
4
5
6
7
3
,
R
5
×
2
1
=
U
5
×
2
o
l
d
·
5
+
2
-
L
5
×
2
1
-
I
5
×
2
=
0.903
0.700
0.499
0.698
0.102
0.794
0.599
0.404
0.300
0.002
.
From Step 2-(a), the adjacent matrix is computed:(16)
A
5
×
5
=
0
0
0
0
0
0
0
0
1
0
0
0
0
1
0
0
1
1
0
0
0
0
0
0
0
.
From Step 2-(b), all the combinations are noted asB:(17)
B
=
n
u
l
l
n
u
l
l
n
u
l
l
n
u
l
l
n
u
l
l
1
n
u
l
l
n
u
l
l
n
u
l
l
n
u
l
l
2
n
u
l
l
n
u
l
l
n
u
l
l
n
u
l
l
3
n
u
l
l
n
u
l
l
n
u
l
l
n
u
l
l
⋮
⋮
⋮
⋮
⋮
2
3
4
5
n
u
l
l
1
2
3
4
5
32
×
5
.
Delete the point(s) corresponding with the element(s) ofB
p, p
=
1,2
,
…
,
32. For the fifth iteration, p
=
5 and B
5
=
4
n
u
l
l
n
u
l
l
n
u
l
l
n
u
l
l. Delete the fourth point (i.e., delete the fourth row and fourth column of A
5
×
5) and the new adjacent matrix is A
4
×
4
′
=
0
→. Then, B
4
′
=
1
2
3
5 and n
′
=
4. From Step 2-(c), the addition sample is computed:(18)
L
4
×
2
2
=
1
1
3
6
5
4
7
3
,
R
4
×
2
2
=
0.903
0.700
0.499
0.698
0.102
0.794
0.300
0.002
,
U
4
×
2
2
=
0.129
0.100
0.357
0.814
0.586
0.542
0.900
0.286
,
X
4
×
2
2
=
1.29
6.12
3.57
8.46
5.86
7.60
9.00
6.90
.
From Step 3, an order matrix is original:(19)
L
7
×
2
′
=
6
3
7
5
1
2
4
4
7
1
6
2
3
5
T
.
Then(20)
L
3
×
2
3
=
6
7
2
2
4
5
,
R
3
×
2
3
=
0.800
0.916
0.679
0.734
0.743
0.392
,
U
3
×
2
3
=
L
3
×
2
3
-
1
→
+
R
3
×
2
3
7
=
0.829
0.988
0.240
0.248
0.535
0.628
,
X
3
×
2
3
=
8.29
9.61
2.40
6.76
5.35
7.84
.
Finally, the ELHS of size 7 is obtained, which is shown in Figure4:(21)
X
7
×
2
n
e
w
=
X
4
×
2
2
X
3
×
2
3
=
1.29
3.57
5.86
9.00
8.29
2.40
5.35
6.12
8.46
7.60
6.90
9.61
6.76
7.84
T
.Figure 4
ELHS of size 7.


## 3.2. General Extension Algorithm of LHS Based on Greedy Algorithm

In graph theory, obtaining the maximum independent set is a NP-hard problem. Step 2-(b) of BGELHS needs the most time and the time complexity of BGELHS isO
(
2
n
). It means that, with the increase of the number of sampling points, the time of BGELHS will be increased exponential. In this section, GGELHS is proposed to reduce the time of extension algorithm. First of all, the definition of degree is introduced.Definition D (see [17]).
In a simple graphG
=
(
V
,
E
), the number of edges associated with the vertex v is the degree of v noted as d
(
v
).
d
(
v
) indicates the number of vertices linked with the vertex v. The basic idea of greedy algorithm is to delete a point one time whose d
(
v
) is the maximum until d
(
v
) of any remaining vertex is 0. The steps of GGELHS are the same as those of BGELHS except for Step 2-(b). The details of Step 2-(b) are given as follows.
Step 2-(b): compute the degree of every vertexd
(
v
i
)
=
∑
j
=
1
n
a
i
j, j
=
1,2
,
…
,
n. Delete the vertex v
i whose degree is the maximum. Delete the corresponding row and column in the adjacent matrix A (A
n
×
n at the beginning) and note the remaining matrix as A
′. Record the order of v
i in B
p. If A
′
≠
0
→, A
=
A
′, and repeat Step 2-(b) to reconstruct A
′. Otherwise, the element of B
p is the order of points to delete. The dimension of A
′ is noted as n
′ and the order of points to remain is recorded in B
p
′.
In Step 2-(a) of GGELHS, when computing the adjacent matrixA
n
×
n, the time required is n
(
n
-
1
)
/
2. It is easy to know that the time complexities of other steps are not greater than O
(
n
). So the time complexity of GGELHS is O
(
n
2
), which is much less than BGELHS. But this algorithm gets the maximal independent set rather than the maximum independent set; that is, this algorithm cannot guarantee to reserve the most original sampling points.Example 2.
To compare BGELHS and GGELHS, we consider the same question of Example1. Since GGELHS is almost the same as BGELHS except for Step 2-(b), we mainly give Step 2-(b) of GGELHS.
From Step 2-(b), the degrees of five points areD
5
=
0
1
1
2
0
T. The fourth point has the maximum degree. So delete the fourth point, and the adjacent matrix is A
4
×
4
′
=
0
→. So, B
p
′
=
1
2
3
5, n
′
=
4.
BGELHS run five times in Step 2-(b) and GGELHS run only one time. As computing the degree of original sampling points costs a little, GGELHS is less time consuming and more practical. From the result above, we can know that Step 2-(b) gets the maximum independent set in this example, which means that this algorithm deletes the least original sampling points. But the greedy algorithm cannot guarantee to have the same effect for every circumstance.


## 4. Application

BGELHS and GGELHS can be applied for many other applications, such as adaptive metamodel building and sequential experiments design and analysis. In Section4, pick up two test examples to evaluate the sample means [6]. The functions are shown in (22) where x
i
∈
[
0,1
], i
=
1,2
,
3,4, and x
i obeys uniform distribution:(22)
y
1
=
∑
i
=
1
4
x
i
+
∑
i
=
1
4
∑
j
=
i
+
1
4
x
i
x
j
,
y
2
=
e
∑
i
=
1
4
x
i
2
+
∑
i
=
1
4
sin
⁡
16
π
x
i
.Four algorithms are applied including Monte Carlo Sample (MCS), classical LHS (CLHS), BGELHS, and GGELHS. Set the maximum sample size noted asN
m
a
x. The steps are as follows.(a)
PickM samples (MCS or LHS) of size N.
(b)
Compute the outputY
i
j of samples, i
=
1,2
,
…
,
M, j
=
1,2
,
…
,
N.
(c)
Compute the average of every sample outputsY
-
i as shown in the following:(23)
Y
-
i
=
∑
j
=
1
N
Y
i
j
N
,
i
=
1,2
,
…
,
M
.
(d)
Compute the standard deviationσ of Y
-
i as shown in the following:(24)
σ
=
∑
i
=
1
M
Y
-
i
-
Y
-
2
M
-
1
,where Y
- is the mean of Y
-
i, i
=
1,2
,
…
,
M.
(e)
LetN
=
N
+
P. If N
>
N
m
a
x, stop. Otherwise, generate sample by the five algorithms and go to step b.In this application,M
=
10, N
=
10, N
m
a
x
=
100, and P
=
10 for MCS, CLHS, BGELHS, and GGELHS. Figure 5 is the result of these four algorithms. It can be concluded that the convergence rates of two extension algorithms are the same as CLHS and better than MCS. Therefore, the extension algorithms reserve the good convergence of CLHS.Figure 5
Change of standard deviationσ ((a) function 1, (b) function 2).
(a)
(b)Compared with CLHS, an outstanding advantage of extension algorithms is the less number of runs. Figure6 shows the number of function runs in this application. The run number of CLHS is (
10
+
20
+
30
+
⋯
+
100
)
×
10
=
5500. BGELHS saves 3302 runs for function 1 and saves 3325 runs for function 2. GGELHS saves 3263 runs for function 1 and saves 3259 runs for function 2. Therefore, BGELHS and GGELHS save more than half of the runs compared with CLHS.Figure 6
Number of runs ((a) function 1, (b) function 2).
(a)
(b)The total times of these extension algorithms are shown in Figure7. The times of MCS and CLHS are almost the same. The time of GGELHS is one and two orders of magnitude slower than that of MCS and CLHS because of Step 2. BGELHS cost the most time, which is a serious problem when applied. GGELHS is easy to apply in practice.Figure 7
Time cost ((a) function 1, (b) function 2).
(a)
(b)


## 5. Conclusion

A drawback to LHS is the extension problem and this paper mainly discusses how to obtain a new strict LHS that contains the original sampling points as many as possible. There are three problems to solve including how to express the relationship of original sampling points in the ELHS structure, how to reserve the most original sampling points and how to generate the new sampling points. In this paper, BGELHS and GGELHS are proposed to solve these problems. BGELHS can reserve the most original points in theory. However, it is difficult to apply in practice because of the problem of time consuming. GGELHS cannot guarantee to reserve the most original points, but the time complexity is much less than that of BGELHS. Therefore, BGELHS is suitable for the analysis of time consuming system when the running time of system is much more valuable than the generation time of ELHS. GGELHS can be applied for more situations. A simple example and application for the evaluation of sample means show the effectiveness of these algorithms.This paper illustrates the relationship of ELHS and graph theory. Some conclusions of graph theory can be applied for the general extension of LHS directly, which is the main contribution of this paper. The relationship of original sampling points in ELHS is expressed in the form of adjacent matrix, but the reverse does not hold. The reason is that some information of original sampling points is lost including the dimensions of variables, the value of sampling points, and so on. Therefore, we can improve BGELHS based on the lost information and some conclusions of graph theory, which is under investigation. On the other hand, LHS is often applied with some criteria for a better distribute or space-filling, such as correlated variables, maximizing minimum distance, and minimizing discrepancy. The general extension of LHS with these criteria is also the future work.



---

*Source: 450492-2015-07-29.xml*