Overview
The goal of outerbase is to make the production of
near-interpolators easy, stable, and scalable. It is based on a
C++ backend and interfaced with R via
Rcpp. Under the hood, it leverages unique, custom linear
algebra using RcppArmadillo (Armadillo) and omp.
The overall structure is designed to be interacted with in an
object-oriented manner using Rcpp modules.
There are ways to interact with outerbase for those who
are uncomfortable with (or just actively detest) object oriented
programming.
To begin, load the package.
Note that if you built this package from source, make sure you use a
compiler that can process omp commands to access the entire
speed benefits.
Simple prediction
To understand how to get started with using the package, we predict
using data from a function with an eight dimensional input commonly
known as the Borehole function.
We begin by generating 1000 points using a test function
?obtest_borehole8d built into the package.
sampsize = 400
d = 8
x = matrix(runif(sampsize*d),ncol=d) #uniform samples
y = obtest_borehole8d(x) + 0.5*rnorm(sampsize)Our goal will be to design a predictor for y given
x that is a near-interpolator.
The simplest way to interact with this package is using
?obfit (fitting outerbase). The function requires
x, y and two other objects.
The value of numb is the number of basis functions you
want to use. The choice of numb is still under research,
but generally you want it to be large as tolerable. Play around!
The underlying concepts of this approach come from Gaussian
processes. Thus the core building block of predictors will be
covariances functions. The choice of the covariances needed for
obfit is a list of strings corresponding to each column in
x. Type listcov() to discover what covariance
functions are are currently deployed.
listcov()
#> [1] "mat25pow" "mat25" "mat25ang"If you are curious about any of them, type,
e.g. ?covf_mat25pow.
Note obfit has some checks in place to prevent serious
damage. They are not foolproof.
obmodel = obfit(x, y, covnames=rep("elephant",8))
#> Error in .checkcov(covnames[k], x[, k]):
#> covariances must be from listcov()
obmodel = obfit(x, y[1:200], covnames=rep("mat25pow",5))
#> Error in obfit(x, y[1:200], covnames = rep("mat25pow", 5)):
#> x and y dims do not align
obmodel = obfit(x[1:2,], y[1:2], covnames=rep("mat25pow",8))
#> Error in obfit(x[1:2, ], y[1:2], covnames = rep("mat25pow", 8)):
#> dimension larger than sample size has not been tested
obmodel = obfit(x, y, numb = 2, covnames=rep("mat25pow",8))
#> Error in obfit(x, y, numb = 2, covnames = rep("mat25pow", 8)):
#> number of basis functions should be less than twice the dimension
obmodel = obfit(100*x, y, covnames=rep("mat25pow",8))
#> Error in .checkcov(covnames[k], x[, k]):
#> x ranges exceed limits of covariance functions
#> the limits are between 0 and 1
#> try rescaling
obmodel = obfit(0.001*x, y, covnames=rep("mat25pow",8))
#> Error in .checkcov(covnames[k], x[, k]):
#> x are too small for ranges
#> the limits are between 0 and 1
#> try rescalingBelow is one correct deployment, where mat25pow is used
for all dimensions. This should take a bit to run, but it should be
around a second on most modern computers.
ptm = proc.time()
obmodel = obfit(x, y, numb=300, covnames=rep("mat25pow",8),
verbose = 3)
#> doing partial optimization
#> max number of cg steps set to 100
#>
#> ########started BFGS#######
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 0 213.968 NA NA 0.1
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 1 -304.83 -518.742 -136.219 0.4
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 2 -809.594 -504.693 -5779.22 0.109596
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 3 -989.343 -179.73 -8.95883 0.273429
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 4 -1091.09 -101.699 -5478.9 0.155643
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 5 -1168.83 -77.7316 -144.603 0.187464
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 6 -1220.9 -52.069 -26.577 0.221628
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 7 -1271.16 -50.2562 -43.5984 0.257669
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 8 -1302.52 -31.3541 -27.7063 0.295091
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 9 -1323.82 -21.2922 -18.7852 0.333396
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 10 -1338.56 -14.7402 -17.5528 0.372104
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 11 -1349.4 -10.8448 -15.595 0.41077
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 12 -1354.93 -5.52821 -6.03393 0.448992
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 13 -1358.64 -3.70348 -4.06504 0.486424
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 14 -1360.97 -2.3326 -3.25464 0.522773
#> num iter: 15 obj start: 213.9682 obj end: -1362.3
#> final learning rate: 0.5578043
#> approx lower bound (not achieved): -1365.812
#> #########finished BFGS########
#>
#> doing optimization 1
#> max number of cg steps set to 312
#>
#> ########started BFGS#######
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 0 -1371.48 NA NA 0.278902
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 1 -1374.76 -3.27751 -0.537417 0.278902
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 2 -1381.21 -6.44684 -5.66844 0.316889
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 3 -1384.73 -3.51968 -3.28774 0.355481
#> num iter: 4 obj start: -1371.485 obj end: -1386.449
#> final learning rate: 0.3942162
#> approx lower bound (not achieved): -1391.91
#> #########finished BFGS########
#>
#> doing optimization 2
#> max number of cg steps set to 315
#>
#> ########started BFGS#######
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 0 -1386.39 NA NA 0.197108
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 1 -1386.74 -0.349188 -0.170362 0.197108
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 2 -1387.06 -0.316986 -0.191856 0.231864
#> iter.no obj.value wolfe.cond.1 wolfe.cond.2 learning.rate
#> 3 -1387.3 -0.23807 -0.194395 0.268356
#> num iter: 4 obj start: -1386.393 obj end: -1387.446
#> final learning rate: 0.3060832
#> approx lower bound (not achieved): -1388.044
#> #########finished BFGS########
print((proc.time() - ptm)[3])
#> elapsed
#> 3.952Note that the package is made using custom parallelization at the
linear-algebra level. The package relies on omp for
parallelization, so if the package was not compiled with that in place
there will be no benefits. The default call of obfit grabs
all available threads, which is ideal for desktops/laptops. It might be
less ideal for large clusters where the CPU might be shared.
We can adjust the number of threads manually. Below we reduce ourselves to a single thread, which should slow things down.
ptm = proc.time()
obmodel = obfit(x, y, numb=300, covnames=rep("mat25pow",8),
nthreads=1) #optional input
print((proc.time() - ptm)[3])
#> elapsed
#> 5.71We can then predict using ?obpred. While it is not an
exact interpolator, it is close.
predtr = obpred(obmodel, x)
rmsetr = sqrt(mean((y-predtr$mean)^2))
plot(predtr$mean, y,
main=paste("training \n RMSE = ", round(rmsetr,3)),
xlab="prediction", ylab = "actual")
Since we generated this data, we can show that outerbase
can reasonably predict ground truth, meaning overfitting is not an
issue.
ytrue = obtest_borehole8d(x)
rmsetr = sqrt(mean((ytrue-predtr$mean)^2))
plot(predtr$mean, ytrue,
main=paste("oracle \n RMSE = ", round(rmsetr,3)),
xlab="prediction", ylab="actual")
1000 test points generated the same way as our original
data can also serve as a verification process.
xtest = matrix(runif(1000*d),ncol=d) #prediction points
ytest = obtest_borehole8d(xtest) + 0.5*rnorm(1000)The predictions at these new points are also quite good. Not quite as good as the residuals on the test set, but we are extrapolating here.
predtest = obpred(obmodel, xtest)
rmsetst = sqrt(mean((ytest-predtest$mean)^2))
plot(predtest$mean, ytest,
main=paste("testing \n RMSE = ", round(rmsetst,3)),
xlab="prediction", ylab="actual")
This package also produces variances on the predictions which we can use to test reasonableness. The fact that the second histogram looks like a standard Normal is promising that the predictions are reasonable.
hist((ytest-predtest$mean),
main="testing \n residuals", xlab="residuals")
hist((ytest-predtest$mean)/sqrt(predtest$var),
main="testing \n standarized residuals",
xlab="standarized residuals")