This is taken from the jupyter notebook by Dino Sejdinovic, Department of Statistics, Oxford. Slides and other materials available at: http://www.stats.ox.ac.uk/~sejdinov/atsml19/
Some of the topics in that course we cover here also. The level is roughly the same.
Let's generate a function that is nonlinear in $X$.
fstar <- function(x){
return(x-1.2*x^2-0.8*x^3+0.6*cos(2*pi*(x^3)*3)) # play with the coefficient (frequency of the cosine) and the magnitudes
}
#
sgma <- 0.1 # play with noise level
n <- 100 # play with the sample size - how does this impact the choice of lambda and bandwidth?
x <- runif(n)
y <- fstar(x)+sgma*rnorm(n)
plot(x,y)
gridx=0:200/200
fgridx <- fstar(gridx)
lines(gridx,fgridx,col='red')
We will use the 'kernlab' library and fast CV package to fit kernel ridge regression with a Gaussian kernel.
Important The package CVST/kernlab we use below uses the convention$ $$k(x,y)=e^(−\sigma∥x−x′∥^2)$$ for the Gaussian RBF kernel.
In the lecture notes we used the $\sigma$ in the denominator so you will see the opposite effect of the bandwidth here. Be careful to check this for packages you use!
So, here larger values of $\sigma$ correspond to rougher kernels. First let us keep $\sigma$ fixed and vary the regularisation parameter $\lambda$ . The blue lines are fitted regression functions.
Note how very small $\lambda$ leads to overfitting (not enough regularisation), while very large $\lambda$ leads to underfitting (too much regularisation).
library(CVST) # Why do we need fast CV? Extensive grid to search for both regularization and bandwidth.
# First let's just run a fit for a subset of lambdas
krr <- constructKRRLearner()
dat <- constructData(x,y)
dat_tst <- constructData(gridx,0)
par(mfrow=c(3,3),oma = c(5,4,0,0) + 0.1,mar = c(0,0,1,1) + 0.1)
lambdas= 10^(seq(-8,0,by=1)) # you can change the density of the lambda grid here.
for(lambda in lambdas){
param <- list(kernel="rbfdot", sigma=50, lambda=lambda)
krr.model <- krr$learn(dat,param)
pred <- krr$predict(krr.model,dat_tst)
plot( x, y, xaxt='n', yaxt='n', main=paste('lambda =',signif(lambda,digits=3)) )
lines(gridx,fgridx,col='red')
lines(gridx,pred,col='blue')
}
Let's fix the regularization and now explore the bandwidth (or inverse bandwidth since the $\sigma$ appears in the numerator in the kernel expressionabove). That is, very small $\sigma$ leads to to underfitting (kernel too smooth), while very large $\sigma$ leads to overfitting (kernel too rough).
par(mfrow=c(3,3),oma = c(5,4,0,0) + 0.1,mar = c(0,0,1,1) + 0.1)
sigmas=10^(seq(1,9,by=1)/3) # change the span or grid density here by varying the step size or the range in seq()
#
for(sigma in sigmas){
param <- list(kernel="rbfdot", sigma=sigma, lambda=0.01)
krr.model <- krr$learn(dat,param)
pred <- krr$predict(krr.model,dat_tst)
plot(x,y,xaxt='n',yaxt='n',main=paste('sigma =',signif(sigma,digits=3)))
lines(gridx,fgridx,col='red')
lines(gridx,pred,col='blue')
}
Let's use cross-validation to select the optimal $\lambda$ and $\sigma$.
#
lambas= 10^(seq(-8,0,by=.5)) # you can change the density of the lambda grid here.
sigmas=10^(seq(1,5,by=.1)/3) # change the span or grid density here by varying the step size or the range in seq()
#
params <- constructParams(kernel="rbfdot", sigma=sigmas, lambda=lambdas) # this is what's going to cost you with dense grids...
opt <- CV(dat, krr, params, fold=10, verbose=FALSE)
param <- list(kernel="rbfdot", sigma=opt[[1]]$sigma, lambda=opt[[1]]$lambda)
#lambda=1*10^-2)
#lambda=opt[[1]]$lambda)
#lambda=0 - no regularization - what happens?
#
krr.model <- krr$learn(dat,param)
pred <- krr$predict(krr.model,dat_tst)
plot(x,y,xaxt='n',yaxt='n',
main=paste("selected values: sigma=",signif(param$sigma,digits=3),", lambda=",signif(param$lambda,digits=3)))
lines(gridx,fgridx,col='red')
lines(gridx,pred,col='blue')
To try at home
Also investigate the stability of the CV selection as a function of e.g. sample size. What if you use no regularization? What happens?