Package {deepICMGP}

Active Learning Cohn for deepICMGP

Description

Select the next design point by minimizing the Active Learning Cohn (ALC) criterion implemented in deepICMGP over a set of candidate inputs. The criterion is averaged over reference inputs and retained MCMC samples.

Usage

ALC_deepICMGP(object, Xcand = NULL, Xref = NULL)

Arguments

Details

To select the next design point that provides information about all outputs jointly, ALC uses a determinant-based measure of the predictive uncertainty remaining after adding a candidate. This measure is averaged over reference inputs and retained MCMC samples. The candidate with the smallest value is selected.

The returned x_next is the row of Xcand corresponding to the lowest ALC value. Candidate inputs and the selected point are returned in their original units.

Value

A list containing:

• value: Numeric vector of ALC values at the candidate points, in the same row order as Xcand.

• Xcand: Matrix of candidate inputs supplied by the user or generated by the function, in original input units.

• x_next: One-row matrix containing the candidate with the smallest ALC value, in original input units.

• time: Computation time in seconds.

Code provenance

This implementation adapts code from the deepgp package and extends it to deep intrinsic coregionalization models. See Sauer, Gramacy, and Higdon (2023) for the underlying deep GP methodology, and https://CRAN.R-project.org/package=deepgp for the source package.

References

Chang, C.-Y. and Sung, C.-L. (2026). Deep Intrinsic Coregionalization Multi-Output Gaussian Process Surrogate with Active Learning. International Journal for Uncertainty Quantification, 16(4), 29–52. doi:10.1615/Int.J.UncertaintyQuantification.2026062570.

Sauer, A., Gramacy, R. B., and Higdon, D. (2023). Active Learning for Deep Gaussian Process Surrogates. Technometrics, 65(1), 4–18. doi:10.1080/00401706.2021.2008505.

See Also

pred_deepICMGP, trim_deepICMGP

Examples

library(lhs)

# Branin Function: 2d input with 3 outputs
f1 <- function(x) {
  x1 <- x[, 1] * 15 - 5
  x2 <- x[, 2] * 15
  10 * sqrt(
    (-1.275 * (1.2 * x1 + 0.4)^2 / pi^2 +
       5 * (1.2 * x1 + 0.4) / pi + (1.2 * x2 + 0.4) - 6)^2 +
      (10 - 5 / (4 * pi)) * cos(1.2 * x1 + 0.4) + 10
  ) + 2 * (1.2 * x1 + 1.9) - 3 * (3 * (1.2 * x2 + 2.4) - 1) -
    1 - 3 * x2 + 1
}
f2 <- function(x) {
  x1 <- x[, 1] * 15 - 5
  x2 <- x[, 2] * 15
  10 * sqrt(
    (-1.275 * (x1 + 2)^2 / pi^2 + 5 * (x1 + 2) / pi + (x2 + 2) - 6)^2 +
      (10 - 5 / (4 * pi)) * cos(x1 + 2) + 10
  ) + 2 * (x1 - 0.5) - 3 * (3 * x2 - 1) - 1
}
f3 <- function(x) {
  x1 <- x[, 1] * 15 - 5
  x2 <- x[, 2] * 15
  (-1.275 * x1^2 / pi^2 + 5 * x1 / pi + x2 - 6)^2 +
    (10 - 5 / (4 * pi)) * cos(x1) + 10
}
response <- function(x) cbind(f1(x), f2(x), f3(x))
set.seed(2025)
X <- maximinLHS(30, 2)
Y <- response(X)

Xcand <- maximinLHS(50, 2)

# Fit a deepICMGP model with 500 MCMC samples and set nugget to 1e-8 for deterministic case.
# For actual analysis, use a longer MCMC chain and assess convergence.
fit <- fit_deepICMGP(X, Y, nmcmc = 500, true_g = 1e-8)

# Compute ALC and select the next design point from Xcand
alc <- ALC_deepICMGP(fit, Xcand = Xcand)
alc$x_next


# Optionally, supply a separate reference grid.
Xref <- maximinLHS(200, 2)
alc <- ALC_deepICMGP(fit, Xcand = Xcand, Xref = Xref)
alc$x_next

# If candidates are not given, 100 candidate inputs will be generated automatically.
alc <- ALC_deepICMGP(fit)
nrow(alc$Xcand)
alc$x_next

Fit a two-layer multi-output Gaussian process

Description

Fit a deep intrinsic coregionalization Gaussian process (deepICMGP) to multiple outputs observed at shared input locations. The sampler updates kernel parameters and a latent Gaussian process layer. The coregionalization matrices describe dependence among outputs and among latent coordinates.

Usage

fit_deepICMGP(
  x,
  y,
  nmcmc = 5000,
  D = max(ncol(x), ncol(y)),
  trace = TRUE,
  true_g = NULL,
  settings = NULL,
  eps = sqrt(.Machine$double.eps)
)

Arguments

Details

The model uses squared-exponential kernels with a scalar lengthscale parameter for each layer. The outer layer maps the D-dimensional latent values to the q outputs. The prior settings apply to the model after the inputs are scaled and each output is standardized to mean zero and standard deviation one.

Metropolis-Hastings updates are used for the lengthscale parameters and, when not fixed, the nugget. The latent layer is updated by elliptical slice sampling. The proposal for a positive parameter a is uniform on [settings$l * a / settings$u, settings$u * a / settings$l]. Defaults are l = 1 and u = 2. Gamma prior densities are evaluated at the parameter minus eps, using shape and rate:

• alpha$g = 1.5, beta$g = 3.9.

• alpha$theta_y = 1.5, beta$theta_y = 3.9 / 6.

• alpha$theta_w = 1.5, beta$theta_w = 3.9 / 4.

where g is the nugget parameter, theta_y and theta_w are the lengthscale parameters for outer and inner-layer, respectively.

For example, settings = list(l = 1, u = 1.5, beta = list(g = 5)) overrides only those entries. Coregionalization matrices are estimated during MCMC sampling, rather than sampled directly from their posterior distributions.

Initial latent values use the first D scaled input columns when D <= ncol(x), or reuse the scaled input matrix otherwise. The inner-layer prior mean is zero, and the inner kernel scale is one. Input and output transformations are learned from training data only and retained for prediction and active learning.

One can inspect trace plots and convergence before using the samples. We suggest removing burn-in with trim_deepICMGP before prediction or active learning.

Value

A list of class deepICMGP. Let nt be the number of MCMC samples.

• x, y: the scaled input and standardized response matrices used by the sampler.

• scaling: training transformations x_min, x_range, y_mean, and y_sd.

• nmcmc, settings, eps: stored MCMC count and settings.

• g, theta_y, theta_w: length-nt vectors of nugget and lengthscale values.

• B_y: q \times q \times nt array of outer-layer coregionalization estimates.

• B_w: D \times D \times nt array of inner-layer coregionalization estimates.

• w: list of nt latent matrices, each n \times D.

• ll: length-nt vector of outer-layer log-likelihood values, up to omitted constants.

• time: elapsed computation time in seconds.

Code provenance

This implementation adapts code from the deepgp package and extends it to deep intrinsic coregionalization models. See Sauer, Gramacy, and Higdon (2023) for the underlying deep GP methodology, and https://CRAN.R-project.org/package=deepgp for the source package.

References

Chang, C.-Y. and Sung, C.-L. (2026). Deep Intrinsic Coregionalization Multi-Output Gaussian Process Surrogate with Active Learning. International Journal for Uncertainty Quantification, 16(4), 29–52. doi:10.1615/Int.J.UncertaintyQuantification.2026062570.

Sauer, A., Gramacy, R. B., and Higdon, D. (2023). Active Learning for Deep Gaussian Process Surrogates. Technometrics, 65(1), 4–18. doi:10.1080/00401706.2021.2008505.

See Also

trim_deepICMGP, pred_deepICMGP, ALC_deepICMGP

Examples

library(lhs)

# Forrester Function: 1d input with 2 outputs
f1 <- function(x) {
  (6 * x - 2)^2 * sin(12 * x - 4)
}
f2 <- function(x) {
  0.5 * f1(x) + 10 * (x - 0.5) + 5
}
response <- function(x) cbind(f1(x), f2(x))
set.seed(2025)
X <- maximinLHS(9, 1)
Y <- response(X)

# Fit a deepICMGP model with 2000 MCMC samples and set nugget to 1e-8 for deterministic case.
fit <- fit_deepICMGP(X, Y, nmcmc = 2000, true_g = 1e-8)

# Check the last coregionalization matrices for both layers.
fit$B_w[, , 2000]
fit$B_y[, , 2000]

# Trace plots for checking convergence of lengthscale parameters.
## Inner-layer lengthscale
ts.plot(fit$theta_w)

## Outer-layer lengthscale
ts.plot(fit$theta_y)


# Branin Function: 2d input with 3 outputs
f1 <- function(x) {
  x1 <- x[, 1] * 15 - 5
  x2 <- x[, 2] * 15
  10 * sqrt(
    (-1.275 * (1.2 * x1 + 0.4)^2 / pi^2 +
       5 * (1.2 * x1 + 0.4) / pi + (1.2 * x2 + 0.4) - 6)^2 +
      (10 - 5 / (4 * pi)) * cos(1.2 * x1 + 0.4) + 10
  ) + 2 * (1.2 * x1 + 1.9) - 3 * (3 * (1.2 * x2 + 2.4) - 1) -
    1 - 3 * x2 + 1
}
f2 <- function(x) {
  x1 <- x[, 1] * 15 - 5
  x2 <- x[, 2] * 15
  10 * sqrt(
    (-1.275 * (x1 + 2)^2 / pi^2 + 5 * (x1 + 2) / pi + (x2 + 2) - 6)^2 +
      (10 - 5 / (4 * pi)) * cos(x1 + 2) + 10
  ) + 2 * (x1 - 0.5) - 3 * (3 * x2 - 1) - 1
}
f3 <- function(x) {
  x1 <- x[, 1] * 15 - 5
  x2 <- x[, 2] * 15
  (-1.275 * x1^2 / pi^2 + 5 * x1 / pi + x2 - 6)^2 +
    (10 - 5 / (4 * pi)) * cos(x1) + 10
}
response <- function(x) cbind(f1(x), f2(x), f3(x))
set.seed(2025)
X <- maximinLHS(30, 2)
Y <- response(X)

# Fit a deepICMGP model with 5000 MCMC samples and set nugget to 1e-8 for deterministic case.
fit <- fit_deepICMGP(X, Y, nmcmc = 5000, true_g = 1e-8)

# Trace plots for checking convergence of lengthscale parameters.
## Inner-layer lengthscale
ts.plot(fit$theta_w)

## Outer-layer lengthscale
ts.plot(fit$theta_y)

Predict output means, variances, and within-location covariances

Description

Compute posterior predictive mean, variance and within-location output covariances at new inputs using all retained MCMC samples. Parallel computation for prediction is supported with the parallel package.

Usage

pred_deepICMGP(
  object,
  x_new,
  store_latent = FALSE,
  return_all = FALSE,
  cores = 1
)

Arguments

Details

The function samples latent values at the new inputs, then combines conditional moments of outer-layer using the laws of total expectation and covariance.

C2[i, , ] describes dependence among outputs at location i. s2_all stores conditional variances for individual MCMC samples. s2 combines their average with the variation in predicted means across samples to obtain the posterior predictive variance.

Value

The supplied deepICMGP object with the following fields provided. Let nt be the number of retained samples and q the number of outputs. The number of prediction locations is n_new = nrow(x_new).

• x_new: the n_{\mathrm{new}} \times d prediction matrix.

• mean: n_{\mathrm{new}} \times q matrix of posterior predictive means.

• s2: n_{\mathrm{new}} \times q matrix of posterior predictive variances.

• C2: n_{\mathrm{new}} \times q \times q array of within-location output covariances. C2[i, , ] is the covariance matrix at row i.

• s2_all: when return_all = TRUE, a list of nt variance matrices, each n_{\mathrm{new}} \times q.

• w_new: when store_latent = TRUE, a list of nt sampled latent matrices, each n_{\mathrm{new}} \times D.

• quad: when store_latent = TRUE, an n_{\mathrm{new}} \times nt matrix of outer-kernel quadratic terms used for active learning.

• time: Total elapsed computation time in seconds, including the time stored in object and the current prediction call.

Code provenance

This implementation adapts code from the deepgp package and extends it to deep intrinsic coregionalization models. See Sauer, Gramacy, and Higdon (2023) for the underlying deep GP methodology, and https://CRAN.R-project.org/package=deepgp for the source package.

References

Chang, C.-Y. and Sung, C.-L. (2026). Deep Intrinsic Coregionalization Multi-Output Gaussian Process Surrogate with Active Learning. International Journal for Uncertainty Quantification, 16(4), 29–52. doi:10.1615/Int.J.UncertaintyQuantification.2026062570.

Sauer, A., Gramacy, R. B., and Higdon, D. (2023). Active Learning for Deep Gaussian Process Surrogates. Technometrics, 65(1), 4–18. doi:10.1080/00401706.2021.2008505.

See Also

fit_deepICMGP, trim_deepICMGP, ALC_deepICMGP

Examples

library(lhs)

# Forrester Function: 1d input with 2 outputs
f1 <- function(x) {
  (6 * x - 2)^2 * sin(12 * x - 4)
}
f2 <- function(x) {
  0.5 * f1(x) + 10 * (x - 0.5) + 5
}
response <- function(x) cbind(f1(x), f2(x))
set.seed(2025)
X <- maximinLHS(9, 1)
Y <- response(X)

Xtest <- maximinLHS(50, 1)
Ytrue <- response(Xtest)

# Fit a deepICMGP model with 2000 MCMC samples and set nugget to 1e-8
# for deterministic case.
fit <- fit_deepICMGP(X, Y, nmcmc = 2000, true_g = 1e-8)

# Trim a fitted deepICMGP model by removing first 1500 MCMC samples and
# retaining every second sample.
trim <- trim_deepICMGP(fit, burn = 1500, thin = 2)

# Predict a trimmed deepICMGP model on Xtest.
pred <- pred_deepICMGP(trim, x_new = Xtest)

# Calculate RMSE for each output.
RMSE <- numeric(ncol(Y))
for (j in seq_len(ncol(Y))) {
  RMSE[j] <- sqrt(mean((pred$mean[, j] - Ytrue[, j])^2))
}
RMSE


# Branin Function: 2d input with 3 outputs
f1 <- function(x) {
  x1 <- x[, 1] * 15 - 5
  x2 <- x[, 2] * 15
  10 * sqrt(
    (-1.275 * (1.2 * x1 + 0.4)^2 / pi^2 +
       5 * (1.2 * x1 + 0.4) / pi + (1.2 * x2 + 0.4) - 6)^2 +
      (10 - 5 / (4 * pi)) * cos(1.2 * x1 + 0.4) + 10
  ) + 2 * (1.2 * x1 + 1.9) - 3 * (3 * (1.2 * x2 + 2.4) - 1) -
    1 - 3 * x2 + 1
}
f2 <- function(x) {
  x1 <- x[, 1] * 15 - 5
  x2 <- x[, 2] * 15
  10 * sqrt(
    (-1.275 * (x1 + 2)^2 / pi^2 + 5 * (x1 + 2) / pi + (x2 + 2) - 6)^2 +
      (10 - 5 / (4 * pi)) * cos(x1 + 2) + 10
  ) + 2 * (x1 - 0.5) - 3 * (3 * x2 - 1) - 1
}
f3 <- function(x) {
  x1 <- x[, 1] * 15 - 5
  x2 <- x[, 2] * 15
  (-1.275 * x1^2 / pi^2 + 5 * x1 / pi + x2 - 6)^2 +
    (10 - 5 / (4 * pi)) * cos(x1) + 10
}
response <- function(x) cbind(f1(x), f2(x), f3(x))
set.seed(2025)
X <- maximinLHS(30, 2)
Y <- response(X)

Xtest <- maximinLHS(100, 2)

# Fit a deepICMGP model with 5000 MCMC samples and set nugget to 1e-8
# for deterministic case.
fit <- fit_deepICMGP(X, Y, nmcmc = 5000, true_g = 1e-8)

# Predict a fitted deepICMGP model on Xtest and return sampled latent values.
pred <- pred_deepICMGP(fit, Xtest, store_latent = TRUE)
length(pred$w_new)

# For parallel prediction:
pred <- pred_deepICMGP(fit, Xtest, cores = 2)

Remove burn-in and thin stored MCMC samples

Description

Retain selected samples from a fitted deepICMGP object before prediction or active learning.

Usage

trim_deepICMGP(object, burn, thin = 1)

Arguments

Details

The number of stored samples nt must be greater than burn. After discarding burn-in, keep samples at the thinning interval. When both nt and burn are multiples of thin, the retained count is (nt - burn) / thin. For example, 1000 samples with burn = 500 and thin = 2 leave 250 samples. The actual retained count is stored in nmcmc.

Retain at least two samples for prediction, which estimates variance across samples. Invalid burn-in/thinning values and choices that retain no samples are rejected.

Value

An object of class deepICMGP with g, theta_y, theta_w, B_y, B_w, w, and ll subset to the retained samples, and nmcmc updated accordingly. Elapsed trimming time is added to time. Training data and settings are unchanged.

References

Chang, C.-Y. and Sung, C.-L. (2026). Deep Intrinsic Coregionalization Multi-Output Gaussian Process Surrogate with Active Learning. International Journal for Uncertainty Quantification, 16(4), 29–52. doi:10.1615/Int.J.UncertaintyQuantification.2026062570.

Sauer, A., Gramacy, R. B., and Higdon, D. (2023). Active Learning for Deep Gaussian Process Surrogates. Technometrics, 65(1), 4–18. doi:10.1080/00401706.2021.2008505.

See Also

fit_deepICMGP, pred_deepICMGP

Examples

library(lhs)

# Forrester Function: 1d input with 2 outputs
f1 <- function(x) {
  (6 * x - 2)^2 * sin(12 * x - 4)
}
f2 <- function(x) {
  0.5 * f1(x) + 10 * (x - 0.5) + 5
}
response <- function(x) cbind(f1(x), f2(x))
set.seed(2025)
X <- maximinLHS(9, 1)
Y <- response(X)

# Fit a deepICMGP model with 2000 MCMC samples and set nugget to 1e-8 
# for deterministic case.
fit <- fit_deepICMGP(X, Y, nmcmc = 2000, true_g = 1e-8)

# Trim a fitted deepICMGP model by removing first 1500 MCMC samples and
# retaining every second sample.
trim <- trim_deepICMGP(fit, burn = 1500, thin = 2)

# Check the number of retained MCMC samples.
trim$nmcmc