Kernel Ridge Regression demo¶

This first part taken from the jupyter notebook by Dino Sejdinovic, Department of Statistics, Oxford.

Let's generate a function that is nonlinear in $X$.

In [1]:
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')
No description has been provided for this image

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

In [3]:
install.packages("CVST")
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')
}

Installing package into 'C:/Users/jornsten/AppData/Local/R/win-library/4.3'
(as 'lib' is unspecified)

also installing the dependency 'kernlab'



package 'kernlab' successfully unpacked and MD5 sums checked
package 'CVST' successfully unpacked and MD5 sums checked

The downloaded binary packages are in
	C:\Users\jornsten\AppData\Local\Temp\RtmpWaZEMW\downloaded_packages

Loading required package: kernlab

Loading required package: Matrix
No description has been provided for this image

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

In [4]:
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')
}
No description has been provided for this image

Let's use cross-validation to select the optimal $\lambda$ and $\sigma$.

In [5]:
# 
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)
In [6]:
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')
No description has been provided for this image

To try at home

  • different n: what if you don't have much data?
  • different frequency components in the function - how "rough" is the true function
  • different density grids

Also investigate the stability of the CV selection as a function of e.g. sample size. What if you use no regularization? What happens?

In case you want to try this from scratch, here is some simple code for that.

In [6]:
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
trn <- 100 # play with the sample size - how does this impact the choice of lambda and bandwidth?
x <- runif(trn)
y <- fstar(x)+sgma*rnorm(trn)
plot(x,y)

ten <- 2000 # test data
xt <- runif(ten)
yt <- fstar(xt)+sgma*rnorm(ten)


gridx=0:200/200
fgridx <- fstar(gridx)
lines(gridx,fgridx,col='red')
No description has been provided for this image
In [21]:
# KRR from scratch
sigm<-46.4 # try different values here
lam<-0.00001
itr<-seq(1,trn)
alldata <- c(x,xt)
ally <- c(y,yt)
K <- exp(-((as.matrix(dist(alldata)))^2)*sigm/2)
dim(K)
alpha <- solve(K[itr,itr]+lam*diag(rep(1,length(itr))))%*%y
fpred <- K[-itr,itr]%*%alpha

plot(x,y)
lines(gridx,fgridx,col='red')
lines(sort(xt),fpred[sort.list(xt)],col="blue")
  1. 2100
  2. 2100
No description has been provided for this image
In [17]:
# CV on training data to select sigma and lambda
lambas= 10^(seq(-8,0,by=.25)) # you can change the density of the lambda grid here.
sigmas=10^(seq(1,5,by=.05)/3)

CVerr<-matrix(0,length(lambas)*length(sigmas),10) # 10 fold
folds <- sample(rep(seq(1,10),ceiling(trn/10))[1:trn],trn)
ct<-0
param<-matrix(0,length(lambas)*length(sigmas),2)
for (ll in (1:length(lambas))) {
    for (ss in (1:length(sigmas))) {
        ct<-ct+1
        param[ct,]<-c(sigmas[ss],lambas[ll])
        for (zz in (1:10)) {
          alldata<-x
          K <- exp(-((as.matrix(dist(alldata)))^2)*sigmas[ss]/2)
          alpha <- solve(K[folds!=zz,folds!=zz]+lambas[ll]*diag(rep(1,length(folds[folds!=zz]))))%*%y[folds!=zz]
          fpred <- K[folds==zz,folds!=zz]%*%alpha
          CVerr[ct,zz]<-sum((y[folds==zz]-fpred)^2)
    }
  }
}
In [18]:
plot(param[,1],apply(CVerr,1,mean),xlab="sigmas",ylab="CVerror",col=as.factor(param[,2]))
plot(param[,2],apply(CVerr,1,mean),xlab="lambdas",ylab="CVerror",col=as.factor(param[,1]))
puse <- param[which.min(apply(CVerr,1,mean)),]
alldata <- c(x,xt)
K <- exp(-((as.matrix(dist(alldata)))^2)*puse[1]/2)
alpha <- solve(K[itr,itr]+puse[2]*diag(rep(1,length(itr))))%*%y
fpred <- K[-itr,itr]%*%alpha

plot(x,y)
lines(gridx,fgridx,col='red')
lines(sort(xt),fpred[sort.list(xt)],col="blue")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

In the next demo I will now illustrate how you can use gradient descent to obtain KRR and Kernel logistic regression models. This is a lead in to gradient boosting.

The simple L2-loss with ridge penalty $\alpha´K\alpha$ results in the negative gradient $Ky-K(K+\lambda I)\alpha$. I initialize with $\alpha=0$ everywhere which acts as a regularizer in itself (if training slowly).

Kernel logistic regression looks almost like the KRR gradient because the gradient of the sigmoid cancels out in the chain rule.