activegp {activegp}

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