C_GP: Active Subspace Matrix closed form expression for a GP.

Description

Computes the integral over the input domain of the outer product of the gradients of a Gaussian process. The corresponding matrix is the C matrix central in active subspace methodology.

Usage

C_GP(
  modelX,
  y,
  measure = "lebesgue",
  xm = NULL,
  xv = NULL,
  S = NULL,
  verbose = TRUE
)

Value

a const_C object with elements

model: GP model provided or estimated;
mat: C matrix estimated;
Wij: list of W matrices, of size number of variables;
ct: covariance type (1 for "Gaussian", 2 for "Matern3_2", 3 for "Matern5_2").

Arguments

modelX: This may be either 1) a homGP or hetGP GP model, see hetGP-package containing, e.g., a vector of thetas, type of covariance ct, an inverse covariance matrix Ki, a design matrix X0, and response vector Z0. 2) A matrix of design locations, in which case a vector of responses must be given as the y argument, and this function will fit a default model for you.
y: A vector of responses corresponding to the design matrix; may be ommited if a GP fit is provided in the modelX argument.
measure: One of c("lebesgue", "gaussian", "trunc_gaussian", "sample", "discrete"), indiciating the probability distribution with respect to which the input points are drawn in the definition of the active subspace. "lebesgue" uses the Lebesgue or Uniform measure over the unit hypercube [0,1]^d. "gaussian" uses a Gaussian or Normal distribution, in which case xm and xv should be specified. "trunc_gaussian" gives a truncated Gaussian or Normal distribution over the unit hypercube [0,1]^d, in which case xm and xv should be specified. "sample" gives the Sample or Empirical measure (dirac deltas located at each design point), which is equivalent to calculating the average expected gradient outer product at the design points. "discrete" gives a measure which puts equal weight at points in the input space specified via the S parameter, which should be a matrix with one row for each atom of the measure.
xm: If measure is "gaussian" or "trunc_gaussian", gives the mean vector.
xv: If measure is "gaussian" or "trunc_gaussian", gives the marginal variance vector. The covariance matrix is assumed to be diagonal.
S: If measure is "discrete", gives the locations of the measure's atoms. S is a matrix, each row of which gives an atom.
verbose: Should we print progress?

References

N. Wycoff, M. Binois, S. Wild (2019+), Sequential Learning of Active Subspaces, preprint.

P. Constantine (2015), Active Subspaces, Philadelphia, PA: SIAM.

Examples

Run this code

################################################################################
### Active subspace of a Gaussian process
################################################################################
# \donttest{ 
library(hetGP); library(lhs)
set.seed(42)

nvar <- 2
n <- 100

# theta gives the subspace direction
f <- function(x, theta, nugget = 1e-3){
  if(is.null(dim(x))) x <- matrix(x, 1)
  xact <- cos(theta) * x[,1] - sin(theta) * x[,2]
  return(hetGP::f1d(xact) + rnorm(n = nrow(x), sd = rep(nugget, nrow(x))))
}

theta_dir <- pi/6
act_dir <- c(cos(theta_dir), -sin(theta_dir))

# Create design of experiments and initial GP model
design <- X <- matrix(signif(maximinLHS(n, nvar), 2), ncol = nvar)
response <- Y <- apply(design, 1, f, theta = theta_dir)
model <- mleHomGP(design, response, known = list(beta0 = 0))

C_hat <- C_GP(model)

# Subspace distance to true subspace:
print(subspace_dist(C_hat, matrix(act_dir, nrow = nvar), r = 1))
plot(design %*% eigen(C_hat$mat)$vectors[,1], response, 
  main = "Projection along estimated active direction")
plot(design %*% eigen(C_hat$mat)$vectors[,2], response, 
  main = "Projection along estimated inactive direction")
  
# For other plots:
# par(mfrow = c(1, 3)) # uncomment to have all plots together
plot(C_hat)
# par(mfrow = c(1, 1)) # restore graphical window

# }