activegp {activegp} R Documentation

## 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')

```

[Package activegp version 1.0.6 Index]