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
)

Arguments

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").

References

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

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

See Also

print.const_C, plot.const_C

Examples

################################################################################
### Active subspace of a Gaussian process
################################################################################
 
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

 

CI on Eigenvalues via Monte Carlo/GP

Description

CI on Eigenvalues via Monte Carlo/GP

Usage

C_GP_ci(model, B = 100)

Arguments

Value

A list with elements ci giving 95

Examples

################################################################################
## Example of uncertainty quantification on C estimate
################################################################################
library(hetGP); library(lhs)
set.seed(42)

nvar <- 2
n <- 20
nits <- 20

# theta gives the subspace direction
f <- function(x, theta, nugget = 1e-6){
  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))

res <- C_GP_ci(model)

plot(c(1, 2), log(c(mean(res$eigen_draws[,1]), mean(res$eigen_draws[,2]))),
  ylim = range(log(res$eigen_draws)), ylab = "Eigenvalue", xlab = "Index")
  segments(1, log(res$ci[1,1]), 1, log(res$ci[2,1]))
  segments(2, log(res$ci[1,2]), 2, log(res$ci[2,2]))

Equivalent of C_GP using RcppArmadillo

Description

Equivalent of C_GP using RcppArmadillo

Usage

C_GP_cpp(design, response, theta, Ki, ct)

Arguments

Value

The active subspace matrix C.

Active subspace for second order linear model

Description

Active subspace for second order linear model

Usage

C_Q(design, response)

Arguments

Value

A matrix corresponding to the active subspace C matrix.

Examples

set.seed(42) 
A <- matrix(c(1, -1, 0, -1, 2, -1.5, 0, -1.5, 4), nrow = 3, byrow = TRUE)
b <- c(1, 4, 9)

# Quadratic function
ftest <- function(x, sd = 1e-6){
   if(is.null(dim(x))) x <- matrix(x, nrow = 1)
   return(3 + drop(diag(x %*% A %*% t(x)) + x %*% b) + 
     rnorm(nrow(x), sd = sd))
}

ntrain <- 10000
design <- 2 * matrix(runif(ntrain * 3), ntrain) - 1
response <- ftest(design)

C_hat <- C_Q(design, response)

plot(design %*% eigen(C_hat)$vectors[,1], response)

# Test 
gfun <- function(x){2 * A %*% t(x) + matrix(b, nrow = nrow(A), ncol = nrow(x))}
grads <- gfun(design)
C_MC <- tcrossprod(grads)/ntrain
C_true <- 4/3 * A %*% A + tcrossprod(b)
subspace_dist(eigen(C_MC)$vectors[,1:2], eigen(C_true)$vectors[,1:2]) 

Expected variance of trace of C

Description

Expected variance of trace of C

Usage

C_tr(C, xnew, grad = FALSE)

Arguments

Value

A real number giving the expected variance of the trace of C given the current design.

References

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

Examples

################################################################################
### Variance of trace criterion landscape
################################################################################
    library(hetGP)
    set.seed(42)
    nvar <- 2
    n <- 20

    # theta gives the subspace direction
    f <- function(x, theta = pi/6, nugget = 1e-6){
     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))))
    }

    design <- matrix(signif(runif(nvar*n), 2), ncol = nvar)
    response <- apply(design, 1, f)
    model <- mleHomGP(design, response, lower = rep(1e-4, nvar),
                      upper = rep(0.5,nvar), known = list(g = 1e-4))
                      
    C_hat <- C_GP(model)

    ngrid <- 101
    xgrid <- seq(0, 1,, ngrid)
    Xgrid <- as.matrix(expand.grid(xgrid, xgrid))
    filled.contour(matrix(f(Xgrid), ngrid))

    Ctr_grid <- apply(Xgrid, 1, C_tr, C = C_hat)
    filled.contour(matrix(Ctr_grid, ngrid), color.palette = terrain.colors,
                   plot.axes = {axis(1); axis(2); points(design, pch = 20)})

Element-wise Cn+1 variance

Description

Element-wise Cn+1 variance

Usage

C_var(C, xnew, grad = FALSE)

Arguments

Value

A real number giving the expected elementwise variance of C given the current design.

References

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

Examples

################################################################################
### Norm of the variance of C criterion landscape
################################################################################
library(hetGP)
set.seed(42)
nvar <- 2
n <- 20

# theta gives the subspace direction
f <- function(x, theta = pi/6, nugget = 1e-6){
 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))))
}

design <- matrix(signif(runif(nvar*n), 2), ncol = nvar)
response <- apply(design, 1, f)
model <- mleHomGP(design, response, lower = rep(1e-4, nvar),
                  upper = rep(0.5,nvar), known = list(g = 1e-4))
                  
C_hat <- C_GP(model)

ngrid <- 51
xgrid <- seq(0, 1,, ngrid)
Xgrid <- as.matrix(expand.grid(xgrid, xgrid))
filled.contour(matrix(f(Xgrid), ngrid))

cvar_crit <- function(C, xnew){
 return(sqrt(sum(C_var(C, xnew)^2)))
}

Cvar_grid <- apply(Xgrid, 1, cvar_crit, C = C_hat)
filled.contour(matrix(Cvar_grid, ngrid), color.palette = terrain.colors,
               plot.axes = {axis(1); axis(2); points(design, pch = 20)})

Alternative Variance of Update

Description

Defined as E[(C - E[C])^2], where A^2 = AA (not elementwise multiplication).

Usage

C_var2(C, xnew, grad = FALSE)

Arguments

Value

A real number giving the expected variance of C defined via matrix multiplication given the current design.

References

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

Examples

################################################################################
### Norm of the variance of C criterion landscape
################################################################################
 
library(hetGP)
set.seed(42)
nvar <- 2
n <- 20

# theta gives the subspace direction
f <- function(x, theta = pi/6, nugget = 1e-6){
 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))))
}

design <- matrix(signif(runif(nvar*n), 2), ncol = nvar)
response <- apply(design, 1, f)
model <- mleHomGP(design, response, lower = rep(1e-4, nvar),
                  upper = rep(0.5,nvar), known = list(g = 1e-4))
                  
C_hat <- C_GP(model)

ngrid <- 51
xgrid <- seq(0, 1,, ngrid)
Xgrid <- as.matrix(expand.grid(xgrid, xgrid))
filled.contour(matrix(f(Xgrid), ngrid))

cvar_crit <- function(C, xnew){
 return(sqrt(sum(C_var(C, xnew)^2)))
}

Cvar_grid <- apply(Xgrid, 1, cvar_crit, C = C_hat)
filled.contour(matrix(Cvar_grid, ngrid), color.palette = terrain.colors,
               plot.axes = {axis(1); axis(2); points(design, pch = 20)})

Active Subspace Prewarping

Description

Computes a matrix square root of C = Lt

Usage

Lt_GP(..., C)

Arguments

Value

The matrix Lt which can be used for sensitivity prewarping, i.e. by computing Xw = X

Covariance of kernel computations

Description

Computes Int(kappa_i(X, design) . kappa_j(design, X)). This function is preferred for initialization

Usage

W_kappa_ij(design, theta, i1, i2, ct)

Arguments

Value

The matrix representing the result of the integration.

Covariance of kernel computations

Description

Computes Int(kappa_i(X, design1) . kappa_j(design2, X)). This function is preferred for initialization

Usage

W_kappa_ij2(design1, design2, theta, i1, i2, ct)

Arguments

Value

matrix of size nrow(design1) x nrow(design2)

Covariance of kernel computations

Description

Computes Int(kappa_i(X, design) . kappa_j(design, X)). This function is preferred for updates

Usage

W_kappa_ij_up(W, design, theta, i1, i2, start, ct)

Arguments

Value

W is modified in-place.

Covariance of kernel computations

Description

Creates a submatrix of the same tensor as W_kappa_ij, but this time l,k give the indices of the observations, not variables.

Usage

W_kappa_lk(design, theta, i1, i2, ct)

Arguments

Value

The matrix representing the result of the integration.

Package activegp

Description

Active subspace estimation with Gaussian processes

Details

The primary function for analysis of the active subspace given some set of function evaluations is C_GP.

C_var, C_var2, and C_tr give three possible acquisition functions for sequential design. Either C_var or C_var2 is recommended, see Wycoff et al for details and the example below for usage.

Author(s)

Nathan Wycoff, Mickael Binois

References

N. Wycoff, M. Binois, S. Wild (2019+), Sequential Learning of Active Subspaces, preprint.
P. Constantine (2015) Active Subspaces: Emerging Ideas for Dimension Reduction in Parameter Studies, SIAM Spotlights

Examples

################################################################################
### Sequential learning of active subspaces
################################################################################
    library(hetGP); library(lhs)
    set.seed(42)

    nvar <- 2
    n <- 20
    nits <- 20

    # theta gives the subspace direction
    f <- function(x, theta, nugget = 1e-6){
      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))))
      return(100*erf((xact + 0.5)*5) + 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, lower = rep(1e-4, nvar),
                      upper = rep(0.5,nvar), known = list(g = 1e-6, beta0 = 0))

    C_hat <- C_GP(model)

    ngrid <- 51
    xgrid <- seq(0, 1,, ngrid)
    Xgrid <- as.matrix(expand.grid(xgrid, xgrid))
    filled.contour(matrix(f(Xgrid, theta = theta_dir), ngrid))
    ssd <- rep(NA, nits)

    # Main loop
    for(nit in 1:nits) {
      cat(nit)
      cat(" ")
      
      af <- function(x, C) C_var(C, x, grad = FALSE)
      af_gr <- function(x, C) C_var(C, x, grad = TRUE)
      Ctr_grid <- apply(Xgrid, 1, af, C = C_hat) # CVAR
      
      # Best candidate point
      opt_cand <- matrix(Xgrid[which.max(Ctr_grid),], 1)
      
      # Refine with gradient based optimization
      opt <-  optim(opt_cand, af, af_gr, method = 'L-BFGS-B', lower = rep(0, nvar), C = C_hat,
                    upper = rep(1, nvar), hessian = TRUE, 
                    control = list(fnscale=-1, trace = 0, maxit = 10))
      
      # Criterion surface with best initial point and corresponding local optimum
      filled.contour(matrix(Ctr_grid, ngrid), color.palette = terrain.colors,
                     plot.axes = {axis(1); axis(2); points(X, pch = 20);
                                  points(opt_cand, pch = 20, col = 'blue');
                                  points(opt$par, pch = 20, col = 'red')})
      
      X <- rbind(X, opt$par)
      Ynew <- f(opt$par, theta = theta_dir)
      Y <- c(Y, Ynew)
      model <- update(model, Xnew = opt$par, Znew = Ynew)
      
      ## periodically restart model fit
      if(nit %% 5 == 0){
        mod2 <- mleHomGP(X = list(X0 = model$X0, Z0 = model$Z0, mult = model$mult), Z = model$Z,
                         known = model$used_args$known, lower = model$used_args$lower,
                         upper = model$used_args$upper)
        if(mod2$ll > model$ll) model <- mod2
      }
      C_hat <- C_GP(model)
      # Compute subspace distance
      ssd[nit] <- subspace_dist(C_hat, matrix(act_dir, nrow = nvar), r = 1)
    }
    plot(ssd, type = 'b')

Extract Matrix

Description

Given a const_C object, extracts the actual matrix itself.

Usage

## S3 method for class 'const_C'
as.matrix(x, ...)

Arguments

Value

The mat entry of C, a matrix.

Get Integerodifferential Quantities

Description

Returns an element related to an integral of a derivative of the desired kernel. See references for details.

Usage

d1(X, x, sigma, type)

Arguments

Value

The scalar value given the integrated derivative.

References

Mickael Binois, Robert B. Gramacy, Michael Ludkovski (2017), Practical Heteroscedastic Gaussian Process Modeling for Large Simulation Experiments, Journal of Computational and Graphical Statistics

Rectangular Domain -> Unbounded Domain

Description

Given an m dimensional function whose inputs live in bounded intervals [a1, b1], ..., [am, bm], return a wrapped version of the function whose inputs live in R^m. Transformed using the logit function.

Usage

domain_to_R(f, domain)

Arguments

Value

A function wrapping f.

Change a function's inputs to live in [-1, 1]

Description

Given an m dimensional function whose inputs live in bounded intervals [a1, b1], ..., [am, bm], return a wrapped version of the function whose inputs live in [-1, 1], ..., [-1, 1].

Usage

domain_to_unit(f, domain)

Arguments

Value

A function wrapping f.

Quantities for Acquisition Functions

Description

Create a single element of the BETA/GAMMA matrix. Used to compute acquisition functions and their gradients.

Usage

get_betagamma(C, xnew, grad = FALSE)

Arguments

Value

If grad == FALSE, A numeric vector of length 2, whose first element of beta_ij and the second gamma_ij. If grad == TRUE, a list with 3 numeric vector elements, the first giving the gradient for beta_ij, and the second for gamma_ij, and the third is the same vector as would have been returned if grad was FALSE: simply the values of beta and gamma.

Covariance of kernel computations

Description

Computes gradient of Int(kappa_i(X, design1) . kappa_j(design2, X)) with respect to the first argument.

Usage

grad_W_kappa_ij2(design1, design2, theta, i1, i2, ct)

Arguments

Value

matrix of size nrow(design1) x nrow(design2)

Covariance of kernel computations

Description

Computes gradient of Int(kappa_i(X, design1) . kappa_j(design2, X)) with respect to the second argument.

Usage

grad_W_kappa_ij2_w2(design1, design2, theta, i1, i2, ct)

Arguments

Value

matrix of size nrow(design1) x nrow(design2)

Estimate the Active Subspace of a Cheap Function using Gradients

Description

Looks between [-1, 1]

Usage

grad_est_subspace(f, r, m, M = NULL, scale = FALSE)

Arguments

Value

A list with sub, the active subspace, sv, the singular values (all m of them), fs, which gives function values, gs, function grads, and X, which gives sampled locations.

Hessian of the log-likelihood with respect to lengthscales hyperparameters Works for homGP and hetGP models from the hetGP package for now.

Description

Hessian of the log-likelihood with respect to lengthscales hyperparameters Works for homGP and hetGP models from the hetGP package for now.

Usage

logLikHessian(model)

Arguments

Value

A matrix giving the Hessian of the GP loglikelihood.

f:[-1, 1] -> R Becomes f:[0,1] -> R

Description

f:[-1, 1] -> R Becomes f:[0,1] -> R

Usage

n11_2_01(f)

Arguments

Value

The same function with domain shifted.

Plot const_C objectc

Description

Plot const_C objectc

Usage

## S3 method for class 'const_C'
plot(x, output = c("all", "matrix", "logvals", "projfn"), ...)

Arguments

Print const_C objects

Description

Print const_C objects

Usage

## S3 method for class 'const_C'
print(x, ...)

Arguments

Covariance of kernel computations

Description

Computes Int(kappa_i(X, design) . kappa_j(design, X)). This function is preferred for initialization

Usage

quick_C(measure, design, Ki, Kir, theta, xm, xv, ct, verbose)

Arguments

Value

The matrix representing the result of the integration.

Get the distance between subspaces defined as the ranges of A and B

Description

Get the distance between subspaces defined as the ranges of A and B

Usage

subspace_dist(A, B, r)

Arguments

Value

A nonnegative scalar giving the cosine of the first principle angle between the two subspaces.

C update with new observations

Description

Update Constantine's C with new point(s) for a GP

Usage

## S3 method for class 'const_C'
update(object, Xnew, Znew, ...)

Arguments

Value

The updated const_C object originally provided.

See Also

C_GP to generate const_C objects from mleHomGP objects; update_C2 for an update using faster expressions.

Examples

 
################################################################################
### Active subspace of a Gaussian process
################################################################################
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)

print(C_hat)
print(subspace_dist(C_hat, matrix(act_dir, nrow = nvar), r = 1))

# New designs
Xnew <- matrix(runif(2), 1)
Znew <- f(Xnew, theta_dir)

C_new <- update(C_hat, Xnew, Znew)
print(C_new)
subspace_dist(C_new, matrix(act_dir, nrow = nvar), r = 1)

Update Constantine's C, using update formula

Description

Update Constantine's C, using update formula

Usage

update_C2(C, xnew, ynew)

Arguments

Value

Updated C matrix, a const_C object.

References

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