Space-Filling designs

Description

This pacakge offers a comprehensive suite of functions to construct various types of space-filling designs, including Latin hypercube designs, clustering-based designs, maximin designs, maximum projection designs, and uniform designs (Joseph 2016). It also offers the option to optimize designs based on user-defined criteria.

Author(s)

Shangkun Wang, V. Roshan Joseph

Maintainer: Shangkun Wang <shangkunwang01@gmail.com>

References

Wang, Shangkun, Xie, Weijun and V. Roshan Joseph. SFDesign: An R package for Space-Filling Designs.

Joseph, V. R. (2016). Space-filling designs for computer experiments: A review. Quality Engineering, 28(1), 28-35.

Clustering error

Description

This function computes the clustering error.

Usage

cluster.error(design, X = NULL, alpha = 1)

Arguments

Details

cluster.error computes the clustering error. The clustering error for a design D=[\bm x_1, \dots, \bm x_n]^T is defined as \frac{1}{N}\sum_{i=1}^n\sum_{\bm x\in{V_i}}\|\bm x - \bm x_i\|^\alpha, where V_i is the Voronoi cell of each design point \bm x_i for i=1,\dots,n, N is the size of X. When \alpha=2, we obtain K-means and when \alpha=1, we obtain K-medians.

Value

clustering error of the design.

Examples

n = 20
p = 3
D = randomLHD(n, p)
cluster.error(D)

Designs generated by clustering algorithms

Description

This function is for producing designs by minimizing the clustering error.

Usage

clustering.design(
  n,
  p,
  X = NULL,
  D.ini = NULL,
  multi.start = 1,
  alpha = 1,
  Lloyd.iter.max = 100,
  cen.iter.max = 10,
  Lloyd.tol = 1e-04,
  cen.tol = 1e-04
)

Arguments

Details

clustering.design produces a design by clustering algorithms. It minimize the clustering error (see cluster.error) by Lloyd's algorithm. When \alpha > 2, accelerated gradient descent is used to find the center for each cluster (Mak, S. and Joseph, V. R. 2018). When \alpha \leq 2: Weizfeld algorithm is used to find the center for each cluster. Let \bm{x}^{(0)}_i=\bm{x}_i denote the initial position of the ith center and and let \bm{S}_i represent the points within its Voronoi cell. The center is then updated as: \bm{x}^{(k+1)}_i = \left.\left(\sum_{\bm{s}\in\bm{S}_i}\frac{\bm{s}}{\|\bm{s}-\bm{x}^{(k)}_i\|_2^{2-\alpha}}\right)\middle/ \left(\sum_{\bm{s}\in\bm{S}_i}\frac{1}{\|\bm{s}-\bm{x}^{(k)}_i\|_2^{2-\alpha}}\right)\right. \quad \text{for } k=0, 1, \dots.

Value

References

Mak, S. and Joseph, V. R. (2018), “Minimax and minimax projection designs using clustering,” Journal of Computational and Graphical Statistics, 27, 166–178.

Examples

n = 20
p = 3
X = spacefillr::generate_sobol_set(1e5*p, p)
D = clustering.design(n, p, X)

Continuous optimization of a design

Description

This function does continuous optimization of an existing design based on a specified criterion. It has an option to run simulated annealing after the continuous optimization.

Usage

continuous.optim(
  D.ini,
  objective,
  gradient = NULL,
  iteration = 10,
  sa = FALSE,
  sa.objective = NULL
)

Arguments

Details

continuous.optim optimizes an existing design based on a specified criterion. It is a wrapper for the L-BFGS-B function from the nloptr packakge (Johnson 2008) and/or GenSA function in GenSA package (Xiang, Gubian, Suomela and Hoeng 2013).

Value

the optimized design.

References

Johnson, S. G. (2008), The NLopt nonlinear-optimization package, available at https://github.com/stevengj/nlopt. Xiang Y, Gubian S, Suomela B, Hoeng (2013). "Generalized Simulated Annealing for Efficient Global Optimization: the GenSA Package for R". The R Journal Volume 5/1, June 2013.

Examples

# Below is an example showing how to create functions needed to generate MaxPro design manually by
# continuous.optim without using the maxpro.optim function in the package.
compute.distance.matrix <- function(A){
   log_prod_metric = function(x, y) 2 * sum(log(abs(x-y)))
   return (c(proxy::dist(A, log_prod_metric)))
}
optim.obj = function(x){
  D = matrix(x, nrow=n, ncol=p)
  d = exp(compute.distance.matrix(D))
  d_matrix = matrix(0, n, n)
  d_matrix[lower.tri(d_matrix)] = d
  d_matrix = d_matrix + t(d_matrix)
  fn = sum(1/d)
  lfn = log(fn)
  I = diag(n)
  diag(d_matrix) = rep(1,n)
  A = B = D
  for(j in 1:p)
  {
    A = t(outer(D[,j], D[,j], "-"))
    diag(A) = rep(1, n)
    B[, j] = diag((1/A - I) %*% (1/d_matrix - I))
  }
  grad = 2 * B / fn
  return(list("objective"=lfn, "gradient"=grad))
}
n = 20
p = 3
D.ini = maxproLHD(n, p)$design
D = continuous.optim(D.ini, optim.obj)

Generate a Latin-hypercube design (LHD) based on a custom criterion

Description

This function generates a LHD by minimizing a user-specified design criterion.

Usage

customLHD(
  compute.distance.matrix,
  compute.criterion,
  update.distance.matrix,
  n,
  p,
  design = NULL,
  max.sa.iter = 1e+06,
  temp = 0,
  decay = 0.95,
  no.update.iter.max = 400,
  num.passes = 10,
  max.det.iter = 1e+06,
  method = "full",
  scaled = TRUE
)

Arguments

Details

customLHD generates a LHD by minimizing a user-specified design criterion.

• If method='sa', then a simulated annealing algorithm is used to optimize the LHD. To custom the optimization process, you can change the default values for max.sa.iter, temp, decay, no.update.iter.max. In this optimization step, two design points are randomly chosen and their coordinate along one dimension are swaped. If the new design improves the criterion, then it is accepted; otherwise, it is accepted with some probability.

• If method='deterministic', then a deterministic swap algorithm is used to optimize the LHD. To custom the optimization process, you can change the default values for num.passes, max.det.iter. In this optimization step, we swap the coordinates of all pairs of design points (start with design point 1 with design point 2, then 1 with 3, ... 1 with n, then 2 with 3 until n-1 with n). Only accept the change if the swap leads to an improvement.

• If method='full', then optimization goes through the above two stages.

Value

Examples

# Below is an example showing how to create functions needed to generate
# MaxPro LHD manually by customLHD without using the maxproLHD function in
# the package.
compute.distance.matrix <- function(A){
   s = 2
   log_prod_metric = function(x, y) s * sum(log(abs(x-y)))
   return (c(proxy::dist(A, log_prod_metric)))
}
compute.criterion <- function(n, p, d) {
  s = 2
  dim <- as.integer(n * (n - 1) / 2)
  # Find the minimum distance
  Dmin <- min(d)
  # Compute the exponential summation
  avgdist <- sum(exp(Dmin - d))
  # Apply the logarithmic transformation and scaling
  avgdist <- log(avgdist) - Dmin
  avgdist <- exp((avgdist - log(dim)) * (p * s) ^ (-1))
  return(avgdist)
}

update.distance.matrix <- function(A, col, selrow1, selrow2, d) {
  s = 2
  n = nrow(A)
  # transform from c++ idx to r idx
  selrow1 = selrow1 + 1
  selrow2 = selrow2 + 1
  col = col + 1
  # A is the updated matrix
  row1 <- min(selrow1, selrow2)
  row2 <- max(selrow1, selrow2)

  compute_position <- function(row, h, n) {
    n*(h-1) - h*(h-1)/2 + row-h
  }

  # Update for rows less than row1
  if (row1 > 1) {
    for (h in 1:(row1-1)) {
      position1 <- compute_position(row1, h, n)
      position2 <- compute_position(row2, h, n)
      d[position1] <- d[position1] + s * log(abs(A[row1, col] - A[h, col])) -
          s * log(abs(A[row2, col] - A[h, col]))
      d[position2] <- d[position2] + s * log(abs(A[row2, col] - A[h, col])) -
          s * log(abs(A[row1, col] - A[h, col]))
    }
  }

  # Update for rows between row1 and row2
  if ((row2-row1) > 1){
    for (h in (row1+1):(row2-1)) {
      position1 <- compute_position(h, row1, n)
      position2 <- compute_position(row2, h, n)
      d[position1] <- d[position1] + s * log(abs(A[row1, col] - A[h, col])) -
          s * log(abs(A[row2, col] - A[h, col]))
      d[position2] <- d[position2] + s * log(abs(A[row2, col] - A[h, col])) -
          s * log(abs(A[row1, col] - A[h, col]))
    }
  }

  # Update for rows greater than row2
  if (row2 < n) {
    for (h in (row2+1):n) {
      position1 <- compute_position(h, row1, n)
      position2 <- compute_position(h, row2, n)
      d[position1] <- d[position1] + s * log(abs(A[row1, col] - A[h, col])) -
          s * log(abs(A[row2, col] - A[h, col]))
      d[position2] <- d[position2] + s * log(abs(A[row2, col] - A[h, col])) -
          s * log(abs(A[row1, col] - A[h, col]))
    }
  }
  return (d)
}

n = 6
p = 2
# Find an appropriate initial temperature
crit1 = 1 / (n-1)
crit2 = (1 / ((n-1)^(p-1) * (n-2))) ^ (1/p)
delta = crit2 - crit1
temp = - delta / log(0.99)
result_custom = customLHD(compute.distance.matrix,
function(d) compute.criterion(n, p, d),
update.distance.matrix, n, p, temp = temp)

Full factorial design

Description

This function generates a full factorial design.

Usage

full.factorial(p, level)

Arguments

Details

full.factorial generates a p dimensional full factorial design.

Value

a full factorial design matrix (scaled to [0, 1])

Examples

p = 3
level = 3
D = full.factorial(p, level)

Augment a design by adding new design points that minimize the reciprocal distance criterion greedily

Description

This function augments a design by adding new design points one-at-a-time that minimize the reciprocal distance criterion.

Usage

maximin.augment(n, p, D.ini, candidate = NULL, r = 2 * p)

Arguments

Details

maximin.augment augments a design by adding new design points that minimize the reciprocal distance criterion (see maximin.crit) greedily. In each iteration, the new design points is selected as the one from the candidate points that has the smallest sum of reciprocal distance to the existing design, that is, \bm x_{k+1} = \arg\min_{\bm x} \sum_{i=1}^k \frac{1}{\|\bm x - \bm x_i\|^r}.

Value

the augmented design.

Examples

n.ini = 10
n = 20
p = 3
D.ini = maximinLHD(n.ini, p)$design
D = maximin.augment(n, p, D.ini)

Maximin criterion

Description

This function calculates the maximin distance or the average reciprocal distance of a design.

Usage

maximin.crit(design, r = 2 * ncol(design), surrogate = FALSE)

Arguments

Details

maximin.crit calculates the maximin distance or the average reciprocal distance of a design. The maximin distance for a design D=[\bm x_1, \dots, \bm x_n]^T is defined as \phi_{Mm} = \min_{i\neq j} \|\bm x_i- \bm x_j\|_2. In optimization, the average reciprocal distance is usually used (Morris and Mitchell, 1995):

\phi_{\text{rec}} = \left(\frac{2}{n(n-1)} \sum_{i\neq j}\frac{1}{\|\bm{x}_i-\bm{x}_j\|_2^r}\right)^{1/r}.

The r is a power parameter and when it is large enough, the reciprocal distance is similar to the exact maximin distance.

Value

the maximin distance or reciprocal distance of the design.

References

Morris, M. D. and Mitchell, T. J. (1995), “Exploratory designs for computational experiments,” Journal of statistical planning and inference, 43, 381–402.

Examples

n = 20
p = 3
D = randomLHD(n, p)
maximin.crit(D)

Optimize a design based on maximin or reciprocal distance criterion

Description

This function optimizes a design by continuous optimization based on reciprocal distance criterion. A simulated annealing step can be enabled in the end to directly optimize the maximin distance criterion.

Usage

maximin.optim(D.ini, iteration = 10, sa = FALSE, find.best.ini = FALSE)

Arguments

Details

maximin.optim optimizes a design by L-BFGS-B algorithm (Liu and Nocedal 1989) based on the reciprocal distance criterion. A simulated annealing step can be enabled in the end to directly optimize the maximin distance criterion. Optimization detail can be found in continuous.optim. We also provide the option to try other initial designs generated internally besides the D.ini provided by the user (see argument find.best.ini).

Value

References

Liu, D. C., & Nocedal, J. (1989). On the limited memory BFGS method for large scale optimization. Mathematical programming, 45(1), 503-528.

Examples

n = 20
p = 3
D = maximinLHD(n, p)$design
D = maximin.optim(D, sa=FALSE)$design
# D = maximin.optim(D, sa=TRUE)$design # Let sa=TRUE only when the n and p is not large.

Sequentially remove design points from a design while maintaining low reciprocal distance criterion as possible

Description

This function sequentially removes design points one-at-a-time from a design while maintaining low reciprocal distance criterion as possible.

Usage

maximin.remove(D, n.remove, r = 2 * p)

Arguments

Details

maximin.remove sequentially removes design points from a design in a greedy way while maintaining low reciprocal distance criterion (see maximin.crit) as possible. In each iteration, the design point with the largest sum of reciprocal distances with the other design points is removed, that is, k^* = \arg\max_k \sum_{i\neq k}\frac{1}{\|\bm x_k - \bm x_i\|^r}.

Value

the updated design.

Examples

n = 20
p = 3
n.remove =  5
D = maximinLHD(n, p)$design
D = maximin.remove(D, n.remove)

Generate a maximin Latin-hypercube design (LHD)

Description

This function generates a LHD with large maximin distance.

Usage

maximinLHD(
  n,
  p,
  design = NULL,
  power = 2 * p,
  max.sa.iter = 1e+06,
  temp = 0,
  decay = 0.95,
  no.update.iter.max = 100,
  num.passes = 10,
  max.det.iter = 1e+06,
  method = "full",
  scaled = TRUE
)

Arguments

Details

maximinLHD generates a LHD or optimize an existing LHD to achieve large maximin distance by optimizing the reciprocal distance (see maximin.crit). The optimization details can be found in customLHD.

Value

Examples

# We show three different ways to use this function.
n = 20
p = 3
D.random = randomLHD(n, p)
# optimize over a random LHD by SA
D = maximinLHD(n, p, D.random, method='sa')
# optimize over a random LHD by deterministic swapping
D = maximinLHD(n, p, D.random, method='deterministic')
# Directly generate an optimized LHD for maximin criterion which goes
# through the above two optimization stages.
D = maximinLHD(n, p)

Maximum projection (MaxPro) criterion

Description

This function calculates the MaxPro criterion of a design.

Usage

maxpro.crit(design, delta = 0)

Arguments

Details

maxpro.crit calculates the MaxPro criterion of a design. The MaxPro criterion for a design D=[\bm x_1, \dots, \bm x_n]^T is defined as

\left\{\frac{1}{{n\choose 2}}\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}\frac{1}{\prod_{l=1}^p(x_{il}-x_{jl})^2+ \delta}\right\}^{1/p},

where p is the dimension of the design (Joseph, V. R., Gul, E., & Ba, S. 2015).

Value

the MaxPro criterion of the design.

References

Joseph, V. R., Gul, E., & Ba, S. (2015). Maximum projection designs for computer experiments. Biometrika, 102(2), 371-380.

Examples

n = 20
p = 3
D = randomLHD(n, p)
maxpro.crit(D)

Optimize a design based on the maximum projection criterion

Description

This function optimizes a design by continuous optimization based on maximum projection criterion (Joseph, V. R., Gul, E., & Ba, S. 2015).

Usage

maxpro.optim(D.ini, iteration = 10)

Arguments

Details

maxpro.optim optimizes a design by L-BFGS-B algorithm (Liu and Nocedal 1989) based on the maximum projection criterion maxpro.crit.

Value

References

Liu, D. C., & Nocedal, J. (1989). On the limited memory BFGS method for large scale optimization. Mathematical programming, 45(1), 503-528.

Joseph, V. R., Gul, E., & Ba, S. (2015). Maximum projection designs for computer experiments. Biometrika, 102(2), 371-380.

Examples

n = 20
p = 3
D = maxproLHD(n, p)$design
D = maxpro.optim(D)$design

Sequentially remove design points from a design while maintaining low maximum projection criterion as possible

Description

This function sequentially removes design points one-at-a-time from a design while maintaining low maximum projection criterion as possible.

Usage

maxpro.remove(D, n.remove, delta = 0)

Arguments

Details

maxpro.remove sequentially removes design points from a design while maintaining low maximum projection criterion (see maxpro.crit) as possible. The maximum projection criterion is modified to include a small delta term:

\phi_{\text{maxpro}}(\bm{D}_n) = \left\{\frac{1}{{n\choose 2}}\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}\frac{1}{\prod_{l=1}^p(x_{il}-x_{jl})^2 +\delta}\right\}^{1/p}.

The index of the point to remove is k^* =\arg\min_k \sum_{i \neq k}\frac{1}{\prod_{l=1}^p(x_{il}-x_{kl})^2 +\delta}.

Value

the updated design.

Examples

n = 20
p = 3
n.remove =  5
D = maxproLHD(n, p)$design
D = maxpro.remove(D, n.remove)

Generate a MaxPro Latin-hypercube design

Description

This function generates a MaxPro Latin-hypercube design.

Usage

maxproLHD(
  n,
  p,
  design = NULL,
  max.sa.iter = 1e+06,
  temp = 0,
  decay = 0.95,
  no.update.iter.max = 400,
  num.passes = 10,
  max.det.iter = 1e+06,
  method = "full",
  scaled = TRUE
)

Arguments

Details

maxproLHD generates a MaxPro Latin-hypercube design (Joseph, V. R., Gul, E., & Ba, S. 2015). The major difference with the MaxPro packages is that we have a deterministic swap algorithm, which can be enabled by setting method="deterministic" or method="full". For optimization details, see the detail section in customLHD.

Value

References

Joseph, V. R., Gul, E., & Ba, S. (2015). Maximum projection designs for computer experiments. Biometrika, 102(2), 371-380.

Examples

n = 20
p = 3
D = maxproLHD(n, p)

Random Latin hypercube design

Description

This function generates a random Latin hypercube design.

Usage

randomLHD(n, p)

Arguments

Details

randomLHD generates a random Latin hypercube design.

Value

a random Latin hypercube design.

Examples

n = 20
p = 3
D = randomLHD(n, p)

Uniform criterion

Description

This function calculates the wrap-around discrepancy of a design.

Usage

uniform.crit(design)

Arguments

Details

uniform.crit calculates the wrap-around discrepancy of a design. The wrap-around discrepancy for a design D=[\bm x_1, \dots, \bm x_n]^T is defined as (Hickernell, 1998):

\phi_{wa} = -\left(\frac{4}{3}\right)^p + \frac{1}{n^2}\sum_{i,j=1}^n\prod_{k=1}^p\left[\frac{3}{2} - |x_{ik}-x_{jk}|(1-|x_{ik}-x_{jk}|)\right].

Value

wrap-around discrepancy of the design

References

Hickernell, F. (1998), “A generalized discrepancy and quadrature error bound,” Mathematics of computation, 67, 299–322.

Examples

n = 20
p = 3
D = randomLHD(n, p)
uniform.crit(D)

Generate a uniform design for discrete factors with different number of levels

Description

This function generates a uniform design for discrete factors with different number of levels.

Usage

uniform.discrete(
  t,
  p,
  levels,
  design = NULL,
  max.sa.iter = 1e+06,
  temp = 0,
  decay = 0.95,
  no.update.iter.max = 400,
  num.passes = 10,
  max.det.iter = 1e+06,
  method = "full",
  scaled = TRUE
)

Arguments

Details

uniform.discrete generates a uniform design of discrete factors with different number of levels by minimizing the wrap-around discrepancy criterion (see uniform.crit).

Value

Examples

p = 5
levels = c(3, 4, 6, 2, 3)
t = 1
D = uniform.discrete(t, p, levels)

Optimize a design based on the wrap-around discrepancy

Description

This function optimizes a design through continuous optimization of the wrap-around discrepancy.

Usage

uniform.optim(D.ini, iteration = 10)

Arguments

Details

uniform.optim optimizes a design through continuous optimization of the wrap-around discrepancy (see uniform.crit) by L-BFGS-B algorithm (Liu and Nocedal 1989).

Value

References

Liu, D. C., & Nocedal, J. (1989). On the limited memory BFGS method for large scale optimization. Mathematical programming, 45(1), 503-528.

Examples

n = 20
p = 3
D = uniformLHD(n, p)$design
D = uniform.optim(D)$design

Generate a uniform Latin-hypercube design (LHD)

Description

This function generates a uniform LHD by minimizing the wrap-around discrepancy.

Usage

uniformLHD(
  n,
  p,
  design = NULL,
  max.sa.iter = 1e+06,
  temp = 0,
  decay = 0.95,
  no.update.iter.max = 400,
  num.passes = 10,
  max.det.iter = 1e+06,
  method = "full",
  scaled = TRUE
)

Arguments

Details

uniformLHD generates a uniform LHD minimizing wrap-around discrepancy (see uniform.crit). The optimization details can be found in customLHD.

Value

Examples

n = 20
p = 3
D = uniformLHD(n, p)