fit_deepICMGP: Fit a two-layer multi-output Gaussian process

In deepICMGP: Deep Intrinsic Coregionalization Gaussian Processes for Multiple Outputs

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

x	Numeric input matrix with n rows and d columns. A numeric vector is converted to a one-column matrix. Each column is automatically scaled to [0,1] using its training minimum and range.
y	Numeric response matrix with n rows and q columns, in the same row order as x. Each column represents one output and is automatically centered and divided by its training sample standard deviation. For a single output, supply a one-column matrix, i.e., matrix(y, ncol = 1).
nmcmc	Length of the MCMC chain. Must be an integer greater than 1. Default is 5000 MCMC runs.
D	Positive integer giving the number of latent variables. Defaults to max(ncol(x), ncol(y)).
trace	Logical. Print progress for every 500 iterations when TRUE.
true_g	A positive scalar specifying a fixed value for the outer-layer nugget g. Default is NULL, in which case g is sampled by MCMC with eps as its lower bound. A small positive value approximates a noise-free response model, e.g., true_g = 1e-8.
settings	List of proposal and prior settings. Default is NULL, which uses the defaults described in Details. One may provide alternative values for any settings for customizing.
eps	A small positive constant for numerical stability. The default is sqrt(.Machine$double.eps), which is approximately 1e-8.

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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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)

deepICMGP documentation built on Nov. 3, 2026, 5:08 p.m.