Ex. derive natural spline bases from B-spline bases

[szcf-weiya]

The question can be found in p36 of Hastie, T., & Tibshirani, R. (1990). Generalized additive models (Vol. 43).

[szcf-weiya]

[szcf-weiya]

illustration with R

let the number of interior knots K = 3, and the total number of knots is 3 + 2 = 5,

> s = smooth.spline(rnorm(5))
> Aknots = s$fit$knot
> Aknots
[1] 0.00 0.00 0.00 0.00 0.25 0.50 0.75 1.00 1.00 1.00 1.00

bs()

calculate the B-spline basis matrix

> B.bs = bs(Aknots[4:8], knots = Aknots[5:7], intercept = T)
> B.bs
     1    2         3         4         5    6 7
[1,] 1 0.00 0.0000000 0.0000000 0.0000000 0.00 0
[2,] 0 0.25 0.5833333 0.1666667 0.0000000 0.00 0
[3,] 0 0.00 0.1666667 0.6666667 0.1666667 0.00 0
[4,] 0 0.00 0.0000000 0.1666667 0.5833333 0.25 0
[5,] 0 0.00 0.0000000 0.0000000 0.0000000 0.00 1
attr(,"degree")
[1] 3
attr(,"knots")
[1] 0.25 0.50 0.75
attr(,"Boundary.knots")
[1] 0 1
attr(,"intercept")
[1] TRUE
attr(,"class")
[1] "bs"     "basis"  "matrix"

which is equivalent to

> B.sD = splineDesign(Aknots, Aknots[4:8])
> B.sD
     [,1] [,2]      [,3]      [,4]      [,5] [,6] [,7]
[1,]    1 0.00 0.0000000 0.0000000 0.0000000 0.00    0
[2,]    0 0.25 0.5833333 0.1666667 0.0000000 0.00    0
[3,]    0 0.00 0.1666667 0.6666667 0.1666667 0.00    0
[4,]    0 0.00 0.0000000 0.1666667 0.5833333 0.25    0
[5,]    0 0.00 0.0000000 0.0000000 0.0000000 0.00    1

ns()

calculate the natural spline basis via ns()

> B.ns = ns(Aknots[4:8], knots = Aknots[5:7], intercept = T)
> B.ns
              1         2           3         4           5
[1,] -0.2672612 0.0000000 -0.21428571 0.6428571 -0.42857143
[2,]  0.5910913 0.1666667 -0.06059511 0.1817853 -0.12119021
[3,]  0.1589087 0.6666667  0.14854170 0.0543749 -0.03624993
[4,]  0.0000000 0.1666667  0.59523810 0.2142857  0.02380952
[5,]  0.0000000 0.0000000 -0.14285714 0.4285714  0.71428571
attr(,"degree")
[1] 3
attr(,"knots")
[1] 0.25 0.50 0.75
attr(,"Boundary.knots")
[1] 0 1
attr(,"intercept")
[1] TRUE
attr(,"class")
[1] "ns"     "basis"  "matrix"

We can derive it from B-spline basis as follows. Firstly, calculate the matrix C

> const = splineDesign(Aknots, c(0, 1), ord = 4L, derivs = c(2L, 2L))
> const
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]   96 -144   48    0    0    0    0
[2,]    0    0    0    0   48 -144   96

then perform QR decomposition,

> qr.const <- qr(t(c))
> qr.const
$qr
             [,1]         [,2]
[1,] -179.5995546    0.0000000
[2,]   -0.8017837 -179.5995546
[3,]    0.2672612    0.0000000
[4,]    0.0000000    0.0000000
[5,]    0.0000000    0.2672612
[6,]    0.0000000   -0.8017837
[7,]    0.0000000    0.5345225

$rank
[1] 2

$qraux
[1] 1.534522 1.000000

$pivot
[1] 1 2

attr(,"class")
[1] "qr"
> qr.Q(qr.const, complete = T)
           [,1]       [,2]       [,3] [,4]        [,5]      [,6]        [,7]
[1,] -0.5345225  0.0000000 -0.2672612    0 -0.21428571 0.6428571 -0.42857143
[2,]  0.8017837  0.0000000  0.1396433    0 -0.15529755 0.4658927 -0.31059511
[3,] -0.2672612  0.0000000  0.9534522    0 -0.03732123 0.1119637 -0.07464246
[4,]  0.0000000  0.0000000  0.0000000    1  0.00000000 0.0000000  0.00000000
[5,]  0.0000000 -0.2672612  0.0000000    0  0.92857143 0.2142857 -0.14285714
[6,]  0.0000000  0.8017837  0.0000000    0  0.21428571 0.3571429  0.42857143
[7,]  0.0000000 -0.5345225  0.0000000    0 -0.14285714 0.4285714  0.71428571
> qr.R(qr.const, complete = T)
          [,1]      [,2]
[1,] -179.5996    0.0000
[2,]    0.0000 -179.5996
[3,]    0.0000    0.0000
[4,]    0.0000    0.0000
[5,]    0.0000    0.0000
[6,]    0.0000    0.0000
[7,]    0.0000    0.0000

note that the above QR decomposition is for Rectangular matrix (https://en.wikipedia.org/wiki/QR_decomposition#Rectangular_matrix) where Q2 can be the basis of the null space of C (see also: Null-space of a rectangular dense matrix), then the matrix N can be calculated as

> B.ns.scratch = B.bs %*% qr.Q(qr.const, complete = T)[,-(1:2)]
> B.ns.scratch
           [,1]      [,2]        [,3]      [,4]        [,5]
[1,] -0.2672612 0.0000000 -0.21428571 0.6428571 -0.42857143
[2,]  0.5910913 0.1666667 -0.06059511 0.1817853 -0.12119021
[3,]  0.1589087 0.6666667  0.14854170 0.0543749 -0.03624993
[4,]  0.0000000 0.1666667  0.59523810 0.2142857  0.02380952
[5,]  0.0000000 0.0000000 -0.14285714 0.4285714  0.71428571

which also can be calculated as

> t(qr.qty(qr.const, t(B.bs)))[, -(1:2)]
              3         4           5         6           7
[1,] -0.2672612 0.0000000 -0.21428571 0.6428571 -0.42857143
[2,]  0.5910913 0.1666667 -0.06059511 0.1817853 -0.12119021
[3,]  0.1589087 0.6666667  0.14854170 0.0543749 -0.03624993
[4,]  0.0000000 0.1666667  0.59523810 0.2142857  0.02380952
[5,]  0.0000000 0.0000000 -0.14285714 0.4285714  0.71428571

we can compare the differences,

> sum((B.ns.scratch - B.ns)^2)
[1] 1.651485e-32

Actually, it can be checked that the source code of ns is exactly following the above way except for usingqr.qty.

> ns
function (x, df = NULL, knots = NULL, intercept = FALSE, Boundary.knots = range(x)) 
{
...
    else basis <- splineDesign(Aknots, x, ord = 4L)
    const <- splineDesign(Aknots, Boundary.knots, ord = 4L, derivs = c(2L, 
        2L))
    if (!intercept) {
        const <- const[, -1, drop = FALSE]
        basis <- basis[, -1, drop = FALSE]
    }
    qr.const <- qr(t(const))
    basis <- as.matrix((t(qr.qty(qr.const, t(basis))))[, -(1L:2L), 
        drop = FALSE])
...
<bytecode: 0x555e95afe8d0>
<environment: namespace:splines>