D \(\mathit{B}\)-spline basis

Given that the solution of the optimization problem is a natural cubic spline with \(n-2\) interior knots, we can represent it in terms of \(\mathit{B}\)-spline basis functions.

We can write \(s(\mathbf{x}) = \sum_1^{n+2}\gamma_d\mathit{B}_d(\mathbf{x})\), where \(\gamma_j\) are coefficients and the \(\mathit{B}_d\) are the cubic \(\mathit{B}\)-spline basis functions. We define the \(n \times (n+2)\) matrix \(\mathbf{B}\) and the \((n+2) \times (n+2)\) matrix \(\mathbf{\Omega}\) by \[\mathit{B}_{id} = \mathit{B}_d(x_i)\] and \[\Omega_{ii'} = \int \mathit{B}_i''(\mathbf{x}) \mathit{B}_{i'}''(\mathbf{x}) dx\]

The optimization criterion \[\sum_{i=1}^n||y_i - s(x_i)||^2 + \lambda \int_a^b{s''(t)}^2dt,\] can be rewrite as:

\[(\mathbf{y} - \mathbf{B}\boldsymbol{\gamma})^T(\mathbf{y} - \mathbf{B}\boldsymbol{\gamma}) + \lambda\boldsymbol{\gamma}^T\boldsymbol{\Omega}\boldsymbol{\gamma} \label{eq:optiBspline}\]

We set the derivative of to 0 with respect to \(\boldsymbol{\gamma}\) to get the solution:

\[\frac{\partial [(\mathbf{y} - \mathbf{B}\boldsymbol{\gamma})^T(\mathbf{y} - \mathbf{B}\boldsymbol{\gamma}) + \lambda\boldsymbol{\gamma}^T\mathbf{\Omega}\boldsymbol{\gamma}]}{\partial \boldsymbol{\gamma}} = 0\]

\[\frac{\partial [\mathbf{y}^T\mathbf{y} - 2 \mathbf{B}^T\mathbf{y}\boldsymbol{\gamma} + \mathbf{B}^T\boldsymbol{\gamma}^T \mathbf{B}\boldsymbol{\gamma} + \lambda\boldsymbol{\gamma}^T\boldsymbol{\Omega}\boldsymbol{\gamma}]}{\partial \boldsymbol{\gamma}} = 0\]

\[-2\mathbf{B}^T\mathbf{y} + 2\mathbf{B}^T\mathbf{B} + 2\lambda\boldsymbol{\Omega}\boldsymbol{\gamma} = 0\]

\[(\mathbf{B}^T\mathbf{B} + \lambda\boldsymbol{\Omega})\boldsymbol{\gamma} = \mathbf{B}^T\mathbf{y}\]

\[\hat{\boldsymbol{\gamma}} = (\mathbf{B}^T\mathbf{B} + \lambda\boldsymbol{\Omega})^{-1}\mathbf{B}^T\mathbf{y}\]

For computational purpose, it can be shown (Hastie and Tibshirani 1990) that the \(\mathit{B}\)-spline basis function of size \(n \times (\mathit{K} + 4)\) can be expressed as a \(n \times (\mathit{K} + 2)\) basis matrix \(\mathbf{N}\) for the natural cubic splines with the same interior and boundary knots at the extreme of \(\mathbf{x}\).

The solution vector \(\hat{s}\) can be write as

\[\hat{\mathbf{S}} = \mathbf{N}\hat{\mathbf{B}eta} = \mathbf{N}(\mathbf{N}^T\mathbf{N} + \lambda\mathbf{\Omega})^{-1}\mathbf{N}^T\mathbf{y} = (\mathbf{I} + \lambda\mathbf{K})^{-1}\mathbf{y}\]

where \(\mathbf{K} = \mathbf{N}^{-T}\mathbf{\Omega}\mathbf{N}^{-1}\) and \(\hat{\mathbf{B}eta}\) the transformed version of \(\hat{\boldsymbol{\gamma}}\) corresponding to the change in basis. In terms of the candidate fitted vector \(\mathbf{f}\) and \(\mathbf{K}\), the cubic smoothing spline \(\hat{f}\) minimizes \[(\mathbf{y}-\mathbf{s})^T(\mathbf{y}-\mathbf{s}) + \lambda\mathbf{s}^T\mathbf{K}\mathbf{s}\] over all vectors \(\mathbf{s}\).

To compute all of this efficiently, the natural spline basis functions should be chosen so that \(\mathbf{N}\) and (4) are band limited which thereby allow the fitted values to be computed in \(O(n)\) calculations (Eubank 1999). Specific ways to obtain such band limited structures are given in (Reinsch 1967). In our case, where the natural spline basis is a cubic spline basis, we can use the piecewise polynomial representation for the estimator describe in (Boor 1975) to show that

\[\hat{\mathbf{S}} = \mathbf{y} - \lambda\mathbf{Q}(\lambda\mathbf{Q}^T\mathbf{Q}+ \Delta)^{-1}\mathbf{Q}^T\mathbf{y},\]

where \(\mathbf{Q}^T\) is an \((n-2) \times n\) tridiagonal matrix with \(\textit{i}\)th row \[(\underbrace{0, \dots, 0}_{i-1}, \frac{1}{t_{i+1}- t_i}, -\frac{1}{t_{i+2}- t_{i+1}} - \frac{1}{t_{i+1}- t_i}, \frac{1}{t_{i+2}- t_{i+1}}, \underbrace{0,\dots,0}_{n-i-2})\] and \(\Delta\) is symmetric, \((n-2) \times (n-2)\), tridiagonal matrix having first and last rows \((t_2 - t_1, t_3-t_2, \underbrace{0,\dots,0}_{n-4})\) and \((\underbrace{0,\dots,0}_{n-4}, t_{n-1}-t_{n-2}, t_n - t_{n_1}),\) and with \(i^{th}\) row \[(\underbrace{0,\dots,0}_{i-2}, t_{i+1} - t_i, 2(t_{i+2}-t_i), t_{i+2}-t_{i+1}, \underbrace{0,\dots,0}_{n-i-3})\]

for \(i = 2, \dots, n-3\).

The fitted values for the cubic smoothing splines can therefore be obtained in \(O(n)\) operations by first solving the 5-banded system

\[(\lambda\mathbf{Q}^T\mathbf{Q} + \Delta)\gamma = \mathbf{Q}^T\mathbf{y}\] and then using \(\hat{\mathbf{S}} = \mathbf{y} - \lambda\mathbf{Q}\boldsymbol{\gamma}\).

References

Boor, Carl de. 1975. A Practical Guide to Splines. Springer. Springer. http://www.springer.com/gb/book/9780387953663.

Eubank, Randall L. 1999. Nonparametric Regression and Spline Smoothing, Second Edition. CRC Press.

Hastie, T. J., and R. J. Tibshirani. 1990. Generalized Additive Models. CRC Press.

Reinsch, Christian H. 1967. “Smoothing by Spline Functions.” Numer. Math. 10 (3): 177–83. http://dx.doi.org/10.1007/BF02162161.