B-spline derivatives

[Jan-Ole Koslik]

Hi, 

I would like to implement a model with a B-spline, but the basis function evaluations are parameter dependent. Hence to make this efficient and stable, one would usually evaluate the design matrix once on a fine grid once before model fitting and interpolate inside the likelihood function for new data.

In theory, this is particularly easy for B-splines as the derivative of f(x) = B(x) %*% coef

is just f`(x) = B`(x) %*% coef, where B` is the design matrix of one degree lower. Obtaining these matrices is also easy through splines2 for example. 

I have the code below and it seems to work, but the obvious question is what to do about higher order derivatives. Again, conceptually they are easy. The derivative in x is just a B-spline one order lower again and the second derivative in the coefficients is already zero.

library(splines2)
library(RTMB)

xfine <- seq(0, 10, length = 1000)
k <- 8 # number of basis functions
knots <- seq(0, 10, length = k-2)
knots <- knots[2:(length(knots)-1)]

# B-spline design matrix
B_fine <- bSpline(xfine, knots = knots, Boundary.knots = c(0,10), degree = 3, intercept = TRUE)
# First derivative of the B-spline design matrix
dB_fine <- bSpline(xfine, knots = knots, Boundary.knots = c(0,10), degree = 3, intercept = TRUE,
                   derivs = 1)

# function that evaluates f(x, par) = B(x) %*% par using linear interpolation
f <- function(xpar) {
  n <- length(xpar) - k
  x <- xpar[1:n]
  par <- xpar[(n + 1):(n + k)]
 
  B_interp <- sapply(1:ncol(B_fine), function(j) approx(xfine, B_fine[, j], xout = x, rule = 2)$y)
  as.vector(B_interp %*% par)
}

# function that evaluates f'(x, par) = d / d c(x, par) f(x, par) using linear interpolation
df <- function(xpar, y, dy) {
  n <- length(xpar) - k
  x <- xpar[1:n]
  par <- xpar[(n + 1):(n + k)]
 
  B_interp <- sapply(1:k, function(j) approx(xfine, B_fine[, j], xout = x, rule = 2)$y)
  Bprime_interp <- sapply(1:ncol(dB_fine), function(j) approx(xfine, dB_fine[, j], xout = x, rule = 2)$y)
 
  # derivative in x
  grad_x <- dy * rowSums(Bprime_interp * matrix(par, nrow = n, ncol = k, byrow = TRUE))
  # derivative in par
  grad_par <- colSums(B_interp * matrix(dy, nrow = n, ncol = k, byrow = FALSE))
 
  c(grad_x, grad_par)
}

# Create the ADjoint function using the wrappers
f_spline <- ADjoint(f, df, name = "f_spline")

x <- 1:9
xpar <- c(x, par)

F_spline <- MakeTape(f_spline, xpar)
F_spline(xpar)
F_spline$jacobian(xpar)

[Phillip Vetter]

Aren't you already there?

Higher order derivatives automatically work provided that df is composed by functions that RTMB already knows how to differentiate.

F_spline$jacfun()$jacobian(xpar)?

with the only trouble being the use of stats::approx which RTMB doesn't know? So you should perhaps interpolate with RTMB::interpol1Dfun.

[Kasper Kristensen]

To follow up on Phillip's answer, if you use RTMB::interpol1Dfun you might not have to worry about derivatives at all. I guess you could just create an interpolating function for each basis function:

li <- split(B_fine, col(B_fine))
B_list <- lapply(li, interpol1Dfun, xlim=c(0,10))
B <- function(x) do.call("cbind", lapply(B_list, function(f) f(x)))
F <- MakeTape(B, 1:9) ## Works

Now you can calculate B(x) %*% par without writing adjoint code.

[Jan-Ole Koslik]

Thanks to both of you. 

I already tried the interpol1Dfun approach and it was very unstable and did not seem to work. Hence in this case I would like to use the fact that B-spline derivatives are "easy" to have the correct interpolation.

I worked on it and figured out the only implementing a function for the (flattened) basis matrix is probably easier and it seems to work perfectly:

# second derivative with derivative

ddb <- ADjoint(
  function(x) {
    ddB_interp <- sapply(1:ncol(ddB_fine), function(j) approx(xfine, ddB_fine[, j], xout = x, rule = 2)$y)
    as.numeric(ddB_interp)
  },
  function(x, y, dy) {
    dddB_interp <- sapply(1:ncol(dddB_fine), function(j) approx(xfine, dddB_fine[, j], xout = x, rule = 2)$y)
    dy_mat <- matrix(dy, ncol = k)
    rowSums(dy_mat * dddB_interp)
  },
  name = "ddb"
)

# derivative with derivative
db <- ADjoint(
  function(x) {
    dB_interp <- sapply(1:ncol(dB_fine), function(j) approx(xfine, dB_fine[, j], xout = x, rule = 2)$y)
    as.numeric(dB_interp)
  },
  function(x, y, dy) {
    ddB_interp <- matrix(ddb(x), ncol = k)
    dy_mat <- matrix(dy, ncol = k)
    rowSums(dy_mat * ddB_interp)
  },
  name = "db"
)

# B-spline with derivative
b <- ADjoint(
  function(x) {

    B_interp <- sapply(1:ncol(B_fine), function(j) approx(xfine, B_fine[, j], xout = x, rule = 2)$y)

    as.numeric(B_interp)
  },
  function(x, y, dy) {
    dB_interp <- matrix(db(x), ncol = k)
    dy_mat <- matrix(dy, ncol = k)
    rowSums(dy_mat * dB_interp)
  },
  name = "b"
)

# matrix format
B <- function(x) matrix(b(x), ncol = k)

par <- sin(1:k)
x <- 1:9

f <- function(xpar) {
  x <- xpar[1:(length(xpar) - k)]
  par <- xpar[(length(xpar) - k + 1):length(xpar)]
  B(x) %*% par
}
xpar <- c(x, par)

F <- MakeTape(f, xpar)
FF <- F$jacfun(sparse = TRUE)
FF(xpar)
FFF <- F$jacfun(sparse = TRUE)
FFF(xpar)

B_fine, dB_fine, ddB_fine, and dddB_fine are pre-evaluated basis matrices.

The only issue I have with the approach above is that it doesn't allow for sparsity, which I also need. Do you think that's easy to incorporate?

[Kasper Kristensen]

You need two changes to get sparsity. First, prevent that B(x[i]) links to all x by evaluating one point at a time:

B <- function(x) t(sapply(x, b))

However, that will generate length(x) calls to the b function, which might be expensive in R...

Second, prevent that B %*% par links to all B by temporarily disabling atomic matrix multiply:

TapeConfig(matmul="plain")
B(x) %*% par
TapeConfig(matmul="atomic") ## restore

[Phillip Vetter]

Just out of curiosity Kasper, is that TapeConfig(atomic="disable") and TapeConfig(atomic="enable")? I can't find anything on matmul.

[Phillip Vetter]

Just to answer my own question: Yes it is.

https://github.com/kaskr/RTMB/blob/09a691a29df1bd8c1aaf0e0db695f4255ef6a50e/RTMB/R/advector.R#L559

    ## Set atomic
    atomic <- match.arg(atomic)
    atomic <- c("NA"="NA", "disable"="plain", "enable"="atomic")[atomic]
    args$matmul <- atomic