1. Introduction to pyLHD

pyLHD is a python implementation of the R package LHD by Hongzhi Wang, Qian Xiao, Abhyuday Mandal. As of now, only the algebraic construction of Latin hypercube designs (LHD) are implemented in this package. For search algorithms to construct LHDs such as: Simulated annealing, particle swarm optimization, and genetic algorithms refer to the R package.

In section 2 algebraic construction methods for LHDs are discussed

To evalute the generated LHDs we consider the following criteria

Maximin distance Criterion

Let \(X\) denote an LHD matrix. Define the \(L_q\)-distance between two run \(x_i\) and \(x_j\) of \(X\) as \(d_q(x_i,x_j) = \left( \sum_{k=1}^m |x_{ik}-x_{jk}|^q \right)^{1/q}\) where \(q\) is an integer. Define the \(L_q\)-distance of design \(X\) as \(d_q(X) = \min \{ d_q(x_i,x_j), 1 \leq i\leq j \leq n \}\). If \(q=1\), we are considering the Manhattan \((L_1)\) distance. If \(q=2\), the Euclidean \((L_2)\) distance is considered. A design \(X\) is called a maximim \(L_q\)-distance if it has the unique largest \(d_q(X)\) value.

Morris and Mitch (1995) and Jin et al. (2005) proposed the \(\phi_p\) criterion which is defined as \[ \phi_p = \left( \sum_{i=1}^{n-1} \sum_{j=i+1}^n d_q (x_i,x_j)^{-p} \right)^{1/p} \]

The \(\phi_p\) criterion is asymptotically equivalent to the Maximin distance criterion as \(p \rightarrow \infty\). In practice \(p=15\) often suffices.

Maximum Projection Criterion

Joseph et al (2015) proposed the maximum projection LHDs that consider designs’ space-filling properties in all possible dimensional spaces. Such designs minimize the maximum projection criterion, which is defined as

\[ \underset{X}{\min} \psi(X) = \left( \frac{1}{{n \choose 2}} \sum_{i=1}^{n-1} \sum_{j=i+1}^n \frac{1}{ \prod_{l=1}^k (x_{il}-x_{jl})^2} \right)^{1/k} \]

We can wee that any two design points should be apart from each other in any projection to minimize the value of \(\psi(x)\)

Orthogonality Criteria

Two major correlation-based criteria to measure designs’ orthogonality is the average absolute correlation criterion and the maximum absolute correlation

\[ ave(|q|) = \frac{2 \sum_{i=1}^{k-1} \sum_{j=i+1}^k |q_{ij}|}{k(k-1)} \quad \text{and} \quad \max |q| = \underset{i,j}{\max} |q_{ij}| \]

where \(q_{ij}\) is the correlation between the \(i\)th and \(j\)th columns of the design matrix \(X\). Orthogonal design have \(ave(|q|)=0\) and \(\max|q|=0\), which may not exist for all design sizes. Designs with smaller \(ave(|q|)\) or \(\max|q|\) are generally preferred in practice.

import pyLHD

Lets start by generating a random centered LHD with 5 rows and 3 columns

X = pyLHD.LatinHypercube(size = (5,3),scramble=False)
X

array([[0.1, 0.5, 0.5],
       [0.5, 0.1, 0.9],
       [0.9, 0.3, 0.1],
       [0.3, 0.9, 0.7],
       [0.7, 0.7, 0.3]])

We evaluate the above design with the different optimamlity criteria described earlier:

The maximin distance criterion (Manhattan)

pyLHD.phi_p(X,p=15,q=1) # using default parameters

1.312360649138747

The maximin distance criterion (Euclidean)

pyLHD.phi_p(X,p=10,q=2) # different p used than above

2.210980529249712

The average absolute correlation

pyLHD.AvgAbsCor(X)

0.3333333333333334

The maximum absolute correlation

pyLHD.MaxAbsCor(X)

0.6

The maximum projection criterion

pyLHD.MaxProCriterion(X)

10.757357196557185

We can apply Williams transformation on X defined as: \[ W(x) = \begin{cases} 2x & 0 \leq x \leq N/2 \\ 2(N-x)-1 & N/2 \leq x < N \end{cases} \]

W_x = pyLHD.WilliamsTransform(X)
W_x

array([[0. , 0.8, 0.8],
       [0.8, 0. , 1.6],
       [1.6, 0.4, 0. ],
       [0.4, 1.6, 1.2],
       [1.2, 1.2, 0.4]])

Lets evaluate the new transformed design

pyLHD.phi_p(W_x)

0.6561803245693735

The \(\phi_p\) value of transformed \(W_x\) is smaller than the original design \(X\)

2. Algebraic Construction Functions

The algebraic construction methods are demonstrated in the table below

	Ye98	Cioppa07	Sun10	Tang93	Lin09	Butler01
Run # \(n\)	\(2^m +1\)	\(2^m +1\)	\(r2^{m +1}\) or \(r2^{m +1} +1\)	\(n\)	\(n^2\)	\(n\)
Factor # \(k\)	\(2m-2\)	\(m + {m-1 \choose 2}\)	\(2^c\)	\(m\)	\(2fp\)	\(k \leq n-1\)
Note	\(m\) is a positive integer \(m\geq 2\)	\(m\) is a positive integer \(m\geq 2\)	\(r\) and \(c\) are positive integers	\(n\) and \(m\) are from \(OA(n,m,s,r)\)	\(n^2,2f\) and \(p\) are from \(OA(n^2,2f,n,2)\) and \(OLHD(n,p)\)	\(n\) is an odd prime number

For theoretical details on the construction methods, a good overview is Section 4.2: Algebraic Constuctions for Orthogonal LHDs from Musings about Constructions of Efficient Latin Hypercube Designs with Flexible Run-sizes

We start by implementing Ye 1998 construction, the resulting desig will have \(2^m+1\) runs and \(2m-2\) factors

Ye98 = pyLHD.OLHD_Ye98(m=4)
Ye98

array([[ 2., -6., -5., -1.,  8.,  3.],
       [ 6.,  2., -7., -8., -1., -4.],
       [ 7., -5.,  6., -4., -3.,  8.],
       [ 5.,  7.,  2., -3.,  4., -1.],
       [ 3., -4., -1.,  5.,  7., -2.],
       [ 4.,  3., -8.,  7., -5.,  6.],
       [ 8., -1.,  4.,  6., -2., -7.],
       [ 1.,  8.,  3.,  2.,  6.,  5.],
       [ 0.,  0.,  0.,  0.,  0.,  0.],
       [-2.,  6.,  5.,  1., -8., -3.],
       [-6., -2.,  7.,  8.,  1.,  4.],
       [-7.,  5., -6.,  4.,  3., -8.],
       [-5., -7., -2.,  3., -4.,  1.],
       [-3.,  4.,  1., -5., -7.,  2.],
       [-4., -3.,  8., -7.,  5., -6.],
       [-8.,  1., -4., -6.,  2.,  7.],
       [-1., -8., -3., -2., -6., -5.]])

pyLHD.MaxAbsCor(Ye98) # column-wise correlation are 0

0.0

Cioppa and Lucas 2007 construction, the resulting design will be a \(2^m+1\) by \(m+ {m-1 \choose 2}\) orthogonal LHD. Note \(m \geq 2\)

Cioppa07 = pyLHD.OLHD_Cioppa07(m=3)
Cioppa07

array([[ 1., -2., -4.,  3.],
       [ 2.,  1., -3., -4.],
       [ 3., -4.,  2., -1.],
       [ 4.,  3.,  1.,  2.],
       [ 0.,  0.,  0.,  0.],
       [-1.,  2.,  4., -3.],
       [-2., -1.,  3.,  4.],
       [-3.,  4., -2.,  1.],
       [-4., -3., -1., -2.]])

pyLHD.MaxAbsCor(Cioppa07) # column-wise correlation are 0

0.0

Sun et al. 2010 construction, the resulting design will be \(r2^{c+1}\) by \(2^c\) if type=‘even’. If type=‘odd’ the resulting design will be \(r2^{c+1} + 1\) by \(2^c\), where \(r\) and \(c\) are positive integers.

Sun10_odd = pyLHD.OLHD_Sun10(C=2,r=2,type='odd')
Sun10_odd

array([[ 1.,  2.,  3.,  4.],
       [ 2., -1., -4.,  3.],
       [ 3.,  4., -1., -2.],
       [ 4., -3.,  2., -1.],
       [ 5.,  6.,  7.,  8.],
       [ 6., -5., -8.,  7.],
       [ 7.,  8., -5., -6.],
       [ 8., -7.,  6., -5.],
       [ 0.,  0.,  0.,  0.],
       [-1., -2., -3., -4.],
       [-2.,  1.,  4., -3.],
       [-3., -4.,  1.,  2.],
       [-4.,  3., -2.,  1.],
       [-5., -6., -7., -8.],
       [-6.,  5.,  8., -7.],
       [-7., -8.,  5.,  6.],
       [-8.,  7., -6.,  5.]])

Sun10_even = pyLHD.OLHD_Sun10(C=2,r=2,type='even')
Sun10_even

array([[ 0.5,  1.5,  2.5,  3.5],
       [ 1.5, -0.5, -3.5,  2.5],
       [ 2.5,  3.5, -0.5, -1.5],
       [ 3.5, -2.5,  1.5, -0.5],
       [ 4.5,  5.5,  6.5,  7.5],
       [ 5.5, -4.5, -7.5,  6.5],
       [ 6.5,  7.5, -4.5, -5.5],
       [ 7.5, -6.5,  5.5, -4.5],
       [-0.5, -1.5, -2.5, -3.5],
       [-1.5,  0.5,  3.5, -2.5],
       [-2.5, -3.5,  0.5,  1.5],
       [-3.5,  2.5, -1.5,  0.5],
       [-4.5, -5.5, -6.5, -7.5],
       [-5.5,  4.5,  7.5, -6.5],
       [-6.5, -7.5,  4.5,  5.5],
       [-7.5,  6.5, -5.5,  4.5]])

Line et al. 2009 construction, the resulting design will be \(n^2\) by \(2fp\). This is obtained by using a \(n\) by \(p\) orthogonal LHD with a \(n^2\) by \(2f\) strength 2 and level \(n\) orthogonal array.

Start by generating an orthogonal LHD

OLHD_example = pyLHD.OLHD_Cioppa07(m=2)

Next, create an orthogonal array with 25 rows, 6 columns, 5 levels, and strength 2 OA(25,6,5,2)

import numpy as np

OA_example = np.array([[2,2,2,2,2,1],[2,1,5,4,3,5],
                      [3,2,1,5,4,5],[1,5,4,3,2,5],
                      [4,1,3,5,2,3],[1,2,3,4,5,2],
                      [1,3,5,2,4,3],[1,1,1,1,1,1],
                      [4,3,2,1,5,5],[5,5,5,5,5,1],
                      [4,4,4,4,4,1],[3,1,4,2,5,4],
                      [3,3,3,3,3,1],[3,5,2,4,1,3],
                      [3,4,5,1,2,2],[5,4,3,2,1,5],
                      [2,3,4,5,1,2],[2,5,3,1,4,4],
                      [1,4,2,5,3,4],[4,2,5,3,1,4],
                      [2,4,1,3,5,3],[5,3,1,4,2,4],
                      [5,2,4,1,3,3],[5,1,2,3,4,2],
                      [4,5,1,2,3,2] ])

Now using Lin at al. 2009 construction, we couple OLHD and OA to obtain

Lin09 = pyLHD.OLHD_Lin09(OLHD=OLHD_example,OA=OA_example)
Lin09

array([[ 12.,  -8.,  12.,  -8.,   7.,  -9.,   6.,  -4.,   6.,  -4.,  -9.,
         -7.],
       [  7.,  -9.,  -7.,   9., -10.,  -2.,  -9.,  -7.,   9.,   7.,  -5.,
         -1.],
       [ 10.,   2.,  -9.,  -7., -11.,   3.,   5.,   1.,  -7.,   9.,  -3.,
        -11.],
       [ -9.,  -7.,  -1.,   5.,  -8., -12.,  -7.,   9.,   2., -10.,  -4.,
         -6.],
       [  4.,   6., -10.,  -2.,   2., -10.,  -8., -12.,  -5.,  -1.,   1.,
         -5.],
       [ 11.,  -3.,  -5.,  -1.,   8.,  12.,   3.,  11.,  10.,   2.,   4.,
          6.],
       [  1.,  -5.,   8.,  12.,  -1.,   5.,  -2.,  10.,   4.,   6.,   2.,
        -10.],
       [  6.,  -4.,   6.,  -4.,   6.,  -4., -12.,   8., -12.,   8., -12.,
          8.],
       [ -1.,   5.,   7.,  -9., -12.,   8.,   2., -10.,  -9.,  -7.,  -6.,
          4.],
       [-12.,   8., -12.,   8.,   3.,  11.,  -6.,   4.,  -6.,   4., -11.,
          3.],
       [ -6.,   4.,  -6.,   4.,   4.,   6.,  12.,  -8.,  12.,  -8.,  -8.,
        -12.],
       [  5.,   1.,   9.,   7.,  -7.,   9., -10.,  -2.,   7.,  -9.,   9.,
          7.],
       [  0.,   0.,   0.,   0.,   5.,   1.,   0.,   0.,   0.,   0., -10.,
         -2.],
       [-10.,  -2.,  -3., -11.,   1.,  -5.,  -5.,  -1.,  11.,  -3.,  -2.,
         10.],
       [ -5.,  -1.,   3.,  11.,  12.,  -8.,  10.,   2., -11.,   3.,   6.,
         -4.],
       [ -7.,   9.,  10.,   2.,  -9.,  -7.,   9.,   7.,   5.,   1.,  -7.,
          9.],
       [  2., -10., -11.,   3.,  11.,  -3.,   1.,  -5.,  -3., -11.,   3.,
         11.],
       [ -8., -12.,   5.,   1.,  -6.,   4.,  -4.,  -6., -10.,  -2.,  12.,
         -8.],
       [ -4.,  -6.,  -8., -12.,  -5.,  -1.,   8.,  12.,  -4.,  -6.,  10.,
          2.],
       [  9.,   7.,  -2.,  10.,  -4.,  -6.,   7.,  -9.,  -1.,   5.,   8.,
         12.],
       [ -3., -11.,   1.,  -5.,  -2.,  10.,  11.,  -3.,  -2.,  10.,  -1.,
          5.],
       [ -2.,  10.,  -4.,  -6.,  -3., -11.,  -1.,   5.,   8.,  12.,  11.,
         -3.],
       [  8.,  12.,   4.,   6.,   0.,   0.,   4.,   6.,  -8., -12.,   0.,
          0.],
       [  3.,  11.,   2., -10.,   9.,   7., -11.,   3.,   1.,  -5.,   7.,
         -9.],
       [-11.,   3.,  11.,  -3.,  10.,   2.,  -3., -11.,   3.,  11.,   5.,
          1.]])

We can convert an orthogonal array into a LHD using the function OA2LHD. Consider the earlier OA_example with 25 rows and 6 columns.

pyLHD.OA2LHD(OA_example)

array([[10,  8,  7,  8,  6,  1],
       [ 9,  1, 22, 17, 14, 24],
       [11,  6,  4, 23, 18, 22],
       [ 4, 21, 17, 11,  9, 21],
       [19,  2, 13, 25,  7, 13],
       [ 5,  9, 11, 16, 23,  8],
       [ 3, 11, 24, 10, 17, 14],
       [ 2,  5,  1,  4,  2,  2],
       [20, 13,  9,  5, 22, 23],
       [25, 25, 25, 22, 21,  5],
       [17, 18, 19, 20, 20,  4],
       [13,  4, 16,  7, 24, 20],
       [15, 12, 15, 14, 13,  3],
       [14, 22,  8, 18,  5, 11],
       [12, 19, 23,  2,  8,  6],
       [23, 16, 14,  6,  4, 25],
       [ 7, 15, 20, 21,  3, 10],
       [ 6, 24, 12,  1, 19, 19],
       [ 1, 20,  6, 24, 15, 18],
       [16, 10, 21, 12,  1, 17],
       [ 8, 17,  2, 13, 25, 12],
       [22, 14,  5, 19, 10, 16],
       [21,  7, 18,  3, 11, 15],
       [24,  3, 10, 15, 16,  9],
       [18, 23,  3,  9, 12,  7]])

Lastly, we consider Butler 2001 construction by generating a \(n\) by \(k\) OLHD

Butler01 = pyLHD.OLHD_Butler01(size = (11,5))
Butler01 

array([[ 5.,  1.,  3.,  2.,  4.],
       [ 8.,  7., 11.,  3.,  2.],
       [ 3., 10.,  4.,  7.,  1.],
       [10.,  4.,  5., 11.,  3.],
       [ 1.,  3., 10.,  8.,  5.],
       [11.,  9.,  2.,  4.,  7.],
       [ 2.,  8.,  7.,  1.,  9.],
       [ 9.,  2.,  8.,  5., 11.],
       [ 4.,  5.,  1.,  9., 10.],
       [ 7., 11.,  9., 10.,  8.],
       [ 6.,  6.,  6.,  6.,  6.]])