Smoothing

Marie Mørch

Theory

Aim: Implement a spline smoother using LOOCV for selecting the parameter \(\lambda\). If we minimize \[ RSS(s)=\sum_{i=1}^n(y_i-s(x_i))^2+\lambda\int s''(t)^2dt \] we get a cubic spline with knots in \(x_i\). A cubic spline is piecewise polynomial of degree 3 on the form \[ f=\sum_i \beta_i \phi_i \] where \(\phi_i\) are basis functions.

Linear Algebra

In vector notation we have \[ \hat{s} = \Phi \hat{\beta} \] where \(\phi_{ij}=\phi_j(x_i)\) and \[ RSS(s,\lambda) = (Y-\Phi \beta )^T(Y-\Phi \beta) + \lambda \beta^T \Omega \beta \] where \(\Omega_{ij}= \int \phi_i''(t)\phi_j''(t)dt\). This is minimized by \[ \hat{\beta}=(\Phi^T \Phi + \lambda \Omega)^{-1} \Phi^T Y \] and the resulting smoother is \[ \hat{s}=\Phi (\Phi^T \Phi + \lambda \Omega)^{-1} \Phi^T Y \]

LOOCV

The spline smoother is a linear smoother (\(\hat{s}=SY\) ) so we can use the formula \[ CVRSS = \sum_i (Y_i- \hat{s}_i^{(-i)})^2 = \sum_i \left(\frac{Y_i-\hat{s}_i}{1- S_{ii}}\right)^2 \]

First implementation

Phi <- function(x){
  knot <- c(rep(min(x),3), x, rep(max(x),3))
  splineDesign(knots=knot,x)
}

Omega <- function(x){
  bsplinepen(create.bspline.basis(range(x), breaks = sort(x)))
}

Smatrix <- function(phi, omega, l){
  phi%*%solve(t(phi)%*%phi+l*omega)%*%t(phi)
}

cvrss <- function(y, l, phi, omega){
  S<- Smatrix(phi, omega, l)
  shat <- S%*%y
  sum((y-shat)^2/(1-diag(S))^2)
}

lambda.opt <- function(y, phi, omega){
  optim(1,cvrss, phi=phi, omega=omega, y=y, method = "BFGS")$par
}

First implementation

smoothspline <- function(x,y){
  phi <- Phi(x)
  omega <- Omega(x)
  lambda <- lambda.opt(y,phi, omega)
  S <- Smatrix(phi,omega, lambda)
  shat <- S%*%y
  mySpline <- structure(list(x=x, s=shat, y = y, lambda =lambda, Omega=omega, Phi=phi) ,class = "mySpline")
  mySpline
}

Methods

summary.mySpline <- function(x){
  cat(paste("Optimal lambda:", x$lambda, "\n \n"))
  tmp<-data.frame( 
    list(xval = x$x,
         yval = x$y,
         smoothval = x$s))
  colnames(tmp) <- c("x values", "y values", "Smoothed values")
  tmp
}

print.mySpline <- function(x){
  cat(paste("Call: smoothspline \n \n"))
  cat("Smoothed values: \n")
  cat(x$s)
}

Methods

plot.mySpline <- function(x, ...){
  p<-ggplot() + geom_point(aes(x=x$x, y=x$y)) + geom_line( aes(x=x$x, y=x$s))
  p
}

lambdaplot <- function(x) UseMethod("lambdaplot")
lambdaplot.mySpline <- function(sm){
  lambdaopt <-sm$lambda
  yopt <- cvrss(sm$y, sm$lambda, sm$Phi, sm$Omega)
  lambdas <- seq(from = 0.001, to= lambdaopt*20, length.out = 300)
  cvrssval <-  sapply(lambdas, function(lambdas) cvrss(sm$y, lambdas,sm$Phi, sm$Omega))
  p<-ggplot()+ geom_line(aes(x=lambdas, y=cvrssval)) 
  p+ labs(title = paste("Optimal lambda", round(lambdaopt,2))) + geom_point(aes(x=lambdaopt, y=yopt), col="red") +xlab("Lambdas") + ylab("CVRSS") 
}

Test

Test

Comparison with the R-smoother:

Simulation

x_sim1 <- seq(-100,100,3)
y_sim1 <- x_sim1^2 + rnorm(length(x_sim1), mean = 0, sd = 3000) 

x_sim2 <- seq(-6,6,0.1) 
y_sim2 <- cos(x_sim2) + rnorm(length(x_sim2), mean = 0, sd = 1)

Profiling

http://rpubs.com/mamo3007/325276

Bottleneck:

Computation of the matrix product

Optimizing

spline <- function(x,y, omega, phi, phiTphi){
  structure(list(x=x, shat=NULL, y = y, lambda = NULL, Omega = omega, Phi = phi, phiTphi = phiTphi, knots=NULL) ,class = "mySpline")
}

newSmatrix <- function(phi, phiTphi, omega, l){
  phi%*%solve(phiTphi+l*omega)%*%t(phi)
}

newcvrss <- function(y, l, phi,phiTphi, omega){
  S<- newSmatrix(phi, phiTphi, omega, l)
  shat <- S%*%y
  sum((y-shat)^2/(1-diag(S))^2)
}

newoptlambda <- function(y, phi, phiTphi, omega){
  optim(1,newcvrss, phi=phi, phiTphi= phiTphi, omega=omega, y=y, method = "BFGS")$par
}

Optimizing

smoothsplineOpt <- function(x,y){
  omega <- Omega(x)
  phi <- Phi(x)
  phiTphi <- crossprod(phi)
  mySpline <- spline(x,y, omega, phi, phiTphi)
  lambda <- newoptlambda(mySpline$y,mySpline$Phi,mySpline$phiTphi, mySpline$Omega)
  S <- newSmatrix(mySpline$Phi,mySpline$phiTphi,mySpline$Omega, lambda)
  shat <- S%*%y
  mySpline$shat <- shat
  mySpline$lambda <- lambda
  mySpline
}

Benchmarking

library(microbenchmark)
microbenchmark(smoothspline(radius,temperature), smoothsplineOpt(radius,temperature))
## Unit: milliseconds
##                                  expr      min       lq      mean   median
##     smoothspline(radius, temperature) 80.41388 87.05506 100.05890 92.40779
##  smoothsplineOpt(radius, temperature) 68.66599 74.19768  83.04732 79.03629
##         uq      max neval
##  101.07965 274.4574   100
##   85.78186 236.5452   100

Sparsity

Plot of the \(\Phi\):

Plot of the \(\Omega\):

Sparsity

sparsePhi <- function(x){
  knot <- c(rep(min(x),3), x, rep(max(x),3))
  splineDesign(knots=knot,x, sparse=TRUE)
}

sparseOmega <- function(x){
  Matrix(bsplinepen(create.bspline.basis(range(x), breaks = sort(x))), sparse = TRUE)
}

sparseSmatrix <- function(phi, phiTphi, omega, l){
  Matrix(phi%*%solve(phiTphi+l*omega)%*%t(phi), sparse = TRUE)
}

sparsecvrss <- function(y, l, phi,phiTphi, omega){
  S<- sparseSmatrix(phi, phiTphi, omega, l)
  shat <- S%*%y
  sum((y-shat)^2/(1-diag(S))^2)
}

sparseoptlambda <- function(y, phi, phiTphi, omega){
  optim(1,sparsecvrss, phi=phi, phiTphi= phiTphi, omega=omega, y=y, method = "BFGS")$par
}

sparsesmoothspline <- function(x,y){
  omega <- sparseOmega(x)
  phi <- sparsePhi(x)
  phiTphi <- crossprod(phi)
  mySpline <- spline(x,y, omega, phi, phiTphi)
  lambda <- sparseoptlambda(mySpline$y,mySpline$Phi,mySpline$phiTphi, mySpline$Omega)
  S <- sparseSmatrix(mySpline$Phi,mySpline$phiTphi,mySpline$Omega, lambda)
  shat <- as.numeric(S%*%y)
  mySpline$shat <- shat
  mySpline$lambda <- lambda
  mySpline
}

Benchmarking

microbenchmark(sparsesmoothspline(radius,temperature), smoothsplineOpt(radius,temperature) )
## Unit: milliseconds
##                                     expr       min        lq     mean
##  sparsesmoothspline(radius, temperature) 130.74860 146.97782 171.3157
##     smoothsplineOpt(radius, temperature)  76.92858  87.00962 102.8447
##    median        uq      max neval
##  162.3281 176.55588 512.9753   100
##   92.5502  99.50594 428.8546   100
microbenchmark(sparsesmoothspline(xsim,ysim), smoothsplineOpt(xsim,ysim), smoothspline(xsim,ysim), times = 10L)
## Unit: seconds
##                            expr      min       lq     mean   median
##  sparsesmoothspline(xsim, ysim) 2.159418 2.288117 2.569798 2.534678
##     smoothsplineOpt(xsim, ysim) 4.819525 5.107443 5.381594 5.224918
##        smoothspline(xsim, ysim) 6.217133 6.362084 6.782426 6.656636
##        uq      max neval
##  2.882165 3.136082    10
##  5.484642 6.513259    10
##  6.736122 7.949948    10

Number of knots

smoothsplineOpt <- function(x,y, nknots=NULL){
  if(is.null(nknots)){
    knot<-c(rep(min(x),3), x, rep(max(x),3))
    phi<- splineDesign(knots=knot,x)
    omega <- bsplinepen(create.bspline.basis(range(x),breaks = x))
  }else{
    index <- round(seq.int(1, length(x), length.out = nknots-6))
    knot <- c(rep(min(x[index]),3), x[index], rep(max(x[index]),3))
    phi<- splineDesign(knots=knot,x)
    omega <- bsplinepen(create.bspline.basis(range(x),breaks = x[index]))
  }
  phiTphi <- crossprod(phi)
  mySpline <- spline(x,y, omega, phi, phiTphi)
  lambda <- newoptlambda(mySpline$y,mySpline$Phi,mySpline$phiTphi, mySpline$Omega)
  S <- newSmatrix(mySpline$Phi,mySpline$phiTphi,mySpline$Omega, lambda)
  shat <- S%*%y
  mySpline$shat <- shat
  mySpline$lambda <- lambda
  mySpline$knots <- knot
  mySpline
}

Number of knots and running time

microbenchmark(smoothsplineOpt(radius,temperature, nknots=10), smoothsplineOpt(radius,temperature, nknots=20), smoothsplineOpt(radius,temperature, nknots=40), smoothsplineOpt(radius,temperature, nknots=60), smoothsplineOpt(radius,temperature))
## Unit: milliseconds
##                                               expr       min        lq
##  smoothsplineOpt(radius, temperature, nknots = 10)  5.924302  6.663621
##  smoothsplineOpt(radius, temperature, nknots = 20)  9.815475 11.029920
##  smoothsplineOpt(radius, temperature, nknots = 40) 21.987687 24.885275
##  smoothsplineOpt(radius, temperature, nknots = 60) 49.606955 56.096676
##               smoothsplineOpt(radius, temperature) 70.562831 77.327953
##       mean    median        uq       max neval
##   7.677046  7.084069  9.082674  11.41140   100
##  12.383532 11.833853 13.482996  19.79536   100
##  33.047422 26.799630 29.276666 328.31105   100
##  62.593037 59.547345 62.897350 208.56993   100
##  80.610207 79.634486 83.202838  99.43354   100

Number of knots and precision

plot((s80$shat- s20$shat)^2,col="purple", type="l",xlab=("radius"),ylab=("temperature"))
lines((s80$shat- s40$shat)^2, col="blue")
lines((s80$shat- s60$shat)^2, col="red")

plot((s80$shat- s40$shat)^2,col="blue", type="l",xlab=("radius"),ylab=("temperature"))
lines((s80$shat- s60$shat)^2, col="red")

Rsm<-smooth.spline(radius,temperature)
Rsm$fit$nk
## [1] 59