Basis Splines¶

This chapter describes functions for the computation of smoothing basis splines (B-splines). A smoothing spline differs from an interpolating spline in that the resulting curve is not required to pass through each data point. For information about interpolating splines, see Interpolation.

The header file gsl_bspline.h contains the prototypes for the bspline functions and related declarations.

Overview¶

Basis splines are commonly used to fit smooth curves to data sets. They are defined on some interval $[a,b]$ which is sub-divided into smaller intervals by defining a non-decreasing sequence of knots. Between any two adjacent knots, a B-spline of order $k$ is a polynomial of degree $k-1$ , and so B-splines defined over the whole interval $[a,b]$ are called piecewise polynomials. In order to generate a smooth spline spanning the full interval, the different polynomial pieces are often required to match one or more derivatives at the knots.

B-splines of order $k$ are basis functions for spline functions of the same order, so that all possible spline functions on some interval $[a,b]$ can be expressed as linear combinations of B-splines,

$f(x; \mathbf{t}) = \sum_{i=1}^n c_i B_{i,k}(x; \mathbf{t}), \quad a \le x \le b$

Here, the $n$ coefficients $c_i$ are known as control points and are often determined from a least squares fit. The $B_{i,k}(x; \mathbf{t})$ are the basis splines of order $k$ , and depend on a knot vector $\mathbf{t}$ , which is a non-decreasing sequence of $n + k$ knots

$\mathbf{t} = \left\{ t_1, t_2, \dots, t_{n + k} \right\}$

On each knot interval $t_i \le x < t_{i+1}$ , the B-splines $B_{i,k}(x;\mathbf{t})$ are polynomials of degree $k-1$ . The $n$ basis splines of order $k$ are defined by

$B_{i,1}(x) &= \left\{ \begin{array}{cc} 1, & t_i \le x < t_{i+1} \\ 0, & \textrm{else} \end{array} \right. \\ B_{i,k}(x) &= {(x - t_i) \over (t_{i+k-1} - t_i)} B_{i,k-1}(x) + {(t_{i+k} - x) \over (t_{i+k} - t_{i+1})} B_{i+1,k-1}(x)$

for $i = 1, \ldots, n$ . The common case of cubic B-splines is given by $k = 4$ . The above recurrence relation can be evaluated in a numerically stable way by the de Boor algorithm.

In what follows, we will drop the $\mathbf{t}$ and simply refer to the spline functions as $f(x)$ and $B_{i,k}(x)$ for simplicity, with the understanding that each set of spline basis functions depends on the chosen knot vector $\mathbf{t}$ . We will also occasionally drop the order $k$ and use $B_i(x)$ and $B_{i,k}(x)$ interchangeably.

Vector Splines¶

The spline above can be generalized to vector control points,

$\mathbf{f}(x) = \sum_{i=1}^n \mathbf{c}_i B_{i,k}(x)$

where the $\mathbf{c}_i$ may have some length $m$ .

Knot vector¶

In practice, one will often identify locations inside the interval $[a,b]$ where it makes sense to end one polynomial piece and begin another. These locations are called interior knots or breakpoints. It may be desirable to increase the density of knots in regions with rapid changes, and decrease the density when the underlying data is changing slowly. For some applications it may be simply best to distribute the knots with equal spacing on the interval, which are called uniform knots.

The full knot vector will look like

$\mathbf{t} = \{ \underbrace{t_1, \dots, t_{k-1}}_{k-1}, \underbrace{t_k, \dots, t_{n+1},}_{n-k+2 \textrm{ interior knots}} \underbrace{t_{n+2}, \dots, t_{n+k}}_{k-1} \}$

The $n-k+2$ interior knots are in $[a,b]$ . In order to have a full basis of B-splines, we require an additional $2(k-1)$ exterior knots, $\{ t_1,\dots,t_{k-1}\}$ and $\{ t_{n+2},\dots,t_{n+k}\}$ . These knots must satisfy

$t_1 & \le \cdots \le t_{k-1} = a \\ b &= t_{n+2} \le \cdots \le t_{n+k}$

but can otherwise be arbitrary.

Allocation of B-splines¶

type gsl_bspline_workspace¶: The computation of B-spline functions requires a preallocated workspace.

gsl_bspline_workspace *gsl_bspline_alloc(const size_t k, const size_t nbreak)¶
gsl_bspline_workspace *gsl_bspline_alloc_ncontrol(const size_t k, const size_t ncontrol)¶: These functions allocate a workspace for computing B-splines of order k. The number of breakpoints is given by nbreak. Alternatively, the number of control points $n$ may be specified in ncontrol. The relationship is $\textrm{nbreak} = n - k + 2$ . Cubic B-splines are specified by $k = 4$ . The size of the workspace is $O(2k^2 + 5k + \textrm{nbreak})$ .

void gsl_bspline_free(gsl_bspline_workspace *w)¶: This function frees the memory associated with the workspace w.

Initializing the knots vector¶

int gsl_bspline_init_augment(const gsl_vector *tau, gsl_bspline_workspace *w)¶

This function initializes the knot vector $\mathbf{t}$ using the non-decreasing interior knot vector tau of length $n - k + 2$ by choosing the exterior knots to be equal to the first and last interior knot. The total knot vector is

$\mathbf{t} = \{ \underbrace{\tau_1, \dots, \tau_1}_{k-1 \textrm{ copies}}, \tau_1, \tau_2, \dots, \tau_{n-k+2}, \underbrace{\tau_{n-k+2}, \dots, \tau_{n-k+2}}_{k-1 \textrm{ copies}} \}$

where $\tau_i$ is specified by the $i$ -th element of tau.

This function does not verify that the tau vector is non-decreasing, so the user must ensure this.

int gsl_bspline_init_uniform(const double a, const double b, gsl_bspline_workspace *w)¶

This function constructs a knot vector using uniformly spaced interior knots on $[a,b]$ and sets the exterior knots to the interval endpoints. The total knot vector is

$\mathbf{t} = \{ \underbrace{a, \dots, a}_{k-1 \textrm{ copies}}, \tau_1, \tau_2, \dots, \tau_{n-k+2}, \underbrace{b, \dots, b}_{k-1 \textrm{ copies}} \}$

with

$\tau_i = a + \frac{i-1}{n-k+1} (b-a), \qquad i = 1, \cdots, n-k+2$

int gsl_bspline_init_periodic(const double a, const double b, gsl_bspline_workspace *w)¶

This function constructs a knot vector suitable for defining a periodic B-spline on $[a,b]$ with a period $P = b - a$ . There are several ways to choose a periodic knot vector. Here, we adopt the approach of Dierckx, 1982, and require

$\begin{array}{ll} t_i = t_{i + n - k + 1} - P, & i = 1, \dots, k-1 \\ t_{n-k+1+i} = t_i + P, & i = k+1, \dots, 2k-1 \end{array}$

We further define a uniform spacing between each consecutive knot pair. The knot values are given by

$t_i = \frac{i-k}{n-k+1} (b-a) + a, \qquad i = 1, \dots, n+k$

int gsl_bspline_init_interp(const gsl_vector *x, gsl_bspline_workspace *w)¶

This function calculates a knot vector suitable for interpolating data at the $n$ abscissa values provided in x, which must be in ascending order. The knot vector for interpolation must satisfy the Schoenberg-Whitney conditions,

$t_i < x_i < t_{i+k}, \quad i = 1, \dots, n$

This function selects the knots according to

$t_i = \left\{ \begin{array}{ll} x_1 & i=1,\dots,k \\ \frac{1}{k-1} \sum_{j=i-k+1}^{i-1} x_j & i=k+1,\dots,n \\ x_n & i=n+1,\dots,n+k \end{array} \right.$

which satisfy the Schoenberg-Whitney conditions, although this may not be the optimal knot vector for a given interpolation problem.

int gsl_bspline_init_hermite(const size_t nderiv, const gsl_vector *x, gsl_bspline_workspace *w)¶

This function calculates a knot vector suitable for performing Hermite interpolation data at the $n$ abscissa values provided in x, which must be in ascending order. The parameter nderiv specifies the maximum derivative order which will be interpolated.

The knot vector has length $(m+1)(n+2)$ , where $m = nderiv$ , and is given by,

$\mathbf{t} = \{ \underbrace{x_0, \dots, x_0}_{m+1 \textrm{ copies}}, \underbrace{x_1, \dots, x_1}_{m+1}, \dots, \underbrace{x_n, \dots, x_n}_{m+1}, \underbrace{x_{n+1}, \dots, x_{n+1}}_{m+1} \}$

This function uses the convention that $x_0 = x_1$ and $x_{n+1} = x_n$ . See Mummy (1989) for more details.

int gsl_bspline_init(const gsl_vector *t, gsl_bspline_workspace *w)¶: This function sets the knot vector equal to the user-supplied vector t which must have length $n+k$ and be non-decreasing.

Properties of B-splines¶

size_t gsl_bspline_ncontrol(const gsl_bspline_workspace *w)¶: This function returns the number of B-spline control points given by $n = \textrm{nbreak} + k - 2$ .

size_t gsl_bspline_order(const gsl_bspline_workspace *w)¶: This function returns the B-spline order $k$ .

size_t gsl_bspline_nbreak(const gsl_bspline_workspace *w)¶: This function returns the number of breakpoints (interior knots) defined for the B-spline.

Evaluation of B-splines¶

int gsl_bspline_calc(const double x, const gsl_vector *c, double *result, gsl_bspline_workspace *w)¶

This function evaluates the B-spline

$f(x) = \sum_{i=1}^n c_i B_i(x)$

at the point x given the $n$ control points c. The result $f(x)$ is stored in result. If x lies outside the knot interval, a Taylor series approximation about the closest knot endpoint is used to extrapolate the spline.

int gsl_bspline_vector_calc(const double x, const gsl_matrix *c, gsl_vector *result, gsl_bspline_workspace *w)¶

This function evaluates the vector valued B-spline

$\mathbf{f}(x) = \sum_{i=1}^n \mathbf{c}_i B_i(x)$

at the point x given the $n$ control vectors $\mathbf{c}_i$ . The vector $\mathbf{c}_i$ is stored in the $i$ -th column of c. The result $\mathbf{f}(x)$ is stored in result. The parameters c and result must have the same number of rows. If x lies outside the knot interval, a Taylor series approximation about the closest knot endpoint is used to extrapolate the spline.

int gsl_bspline_basis(const double x, gsl_vector *B, size_t *istart, gsl_bspline_workspace *w)¶: This function evaluates all potentially nonzero B-spline basis functions at the position x and stores them in the vector B of length $k$ , so that the $i$ -th element is $B_{(istart+i)}(x)$ for $i=0,1,\dots,k-1$ . The output istart is guaranteed to be in $[0, n - k]$ . By returning only the nonzero basis functions, this function allows quantities involving linear combinations of the $B_i(x)$ to be computed without unnecessary terms (such linear combinations occur, for example, when evaluating an interpolated function).

Evaluation of B-spline derivatives¶

int gsl_bspline_calc_deriv(const double x, const gsl_vector *c, const size_t nderiv, double *result, gsl_bspline_workspace *w)¶

This function evaluates the B-spline derivative

$\frac{d^j}{dx^j} f(x) = \sum_{i=1}^n c_i \frac{d^j}{dx^j} B_i(x)$

at the point x given the $n$ control points c. The derivative order $j$ is specified by nderiv. The result $d^j f(x) / dx^j$ is stored in result. If x lies outside the knot interval, a Taylor series approximation about the closest knot endpoint is used to extrapolate the spline.

int gsl_bspline_vector_calc_deriv(const double x, const gsl_matrix *c, const size_t nderiv, gsl_vector *result, gsl_bspline_workspace *w)¶

This function evaluates the vector valued B-spline derivative

$\frac{d^j}{dx^j} \mathbf{f}(x) = \sum_{i=1}^n \mathbf{c}_i \frac{d^j}{dx^j} B_i(x)$

at the point x given the $n$ control vectors $\mathbf{c}_i$ . The vector $\mathbf{c}_i$ is stored in the $i$ -th column of c. The result $d^j \mathbf{f}(x) / dx^j$ is stored in result. The parameters c and result must have the same number of rows. The derivative order $j$ is specified in nderiv. If x lies outside the knot interval, a Taylor series approximation about the closest knot endpoint is used to extrapolate the spline.

int gsl_bspline_basis_deriv(const double x, const size_t nderiv, gsl_matrix *dB, size_t *istart, gsl_bspline_workspace *w)¶: This function evaluates all potentially nonzero B-spline basis function derivatives of orders $0$ through nderiv (inclusive) at the position x and stores them in the matrix dB. The $(i,j)$ -th element of dB is $d^jB_{(istart+i)}(x)/dx^j$ . The last row of dB contains $d^jB_{istart+k-1}(x)/dx^j$ . The matrix dB must be of size $k$ by at least $nderiv + 1$ . The output istart is guaranteed to be in $[0, n - k]$ . Note that function evaluations are included as the zeroth order derivatives in dB. By returning only the nonzero basis functions, this function allows quantities involving linear combinations of the $B_i(x)$ and their derivatives to be computed without unnecessary terms.

Evaluation of B-spline integrals¶

int gsl_bspline_calc_integ(const double a, const double b, const gsl_vector *c, double *result, gsl_bspline_workspace *w)¶

This function computes the integral

$\int_a^b f(x) dx = \sum_{i=1}^n c_i \int_a^b B_i(x) dx$

using the coefficients stored in c and the integration limits (a, b). The integral value is stored in result on output.

int gsl_bspline_basis_integ(const double a, const double b, gsl_vector *y, gsl_bspline_workspace *w)¶

This function computes the integral from a to b of each basis function $B_i(x)$ ,

$y_i = \int_a^b B_i(x) dx$

storing the results in the vector y, of length ncontrol. The algorithm uses Gauss-Legendre quadrature on each knot interval.

Least Squares Fitting with B-splines¶

A common application is to fit a B-spline to a set of $m$ data points $(x_i,y_i)$ without requiring the spline to pass through each point, but rather the spline fits the data in a least square sense. In this case, the number of control points can be much less than $m$ and the result is a smooth spline which minimizes the cost function,

$\chi^2 = \sum_{i=1}^m w_i \left( y_i - f(x_i) \right)^2 = || y - X c ||_W^2$

where $f(x) = \sum_{j=1}^n c_j B_{j,k}(x)$ is the B-spline defined above, $n$ is the number of control points for the spline, $W = \textrm{diag}(w)$ , and $X_{ij} = B_{j,k}(x_i)$ is the least squares design matrix. The parameters $w_i$ are optional weights which can be assigned to each data point $(x_i,y_i)$ . Because each basis spline $B_{j,k}(x)$ has local support $(B_{j,k}(x) = 0$ for $x \notin \left[ t_j,\dots,t_{j+k} \right])$ , the least squares design matrix $X$ has only $k$ nonzero entries per row. Therefore, the normal equations matrix $X^T W X$ is symmetric and banded, with lower bandwidth $k - 1$ . The GSL routines below use a banded Cholesky factorization to solve the normal equations system for the unknown coefficient vector,

$c = \left( X^T W X \right)^{-1} X^T W y$

The normal equations approach is often avoided in least squares problems, since if $X$ is badly conditioned, then $X^T W X$ is much worse conditioned, which can lead to loss of accuracy in the solution vector. However, a theorem of de Boor [de Boor (2001), XI(4)] states that the condition number of the B-spline basis is independent of the knot vector $\mathbf{t}$ and bounded, depending only on $k$ . Numerical evidence has shown that for reasonably small $k$ the normal equations matrix is well conditioned, and a banded Cholesky factorization can be used safely.

int gsl_bspline_lssolve(const gsl_vector *x, const gsl_vector *y, gsl_vector *c, double *chisq, gsl_bspline_workspace *w)¶

int gsl_bspline_wlssolve(const gsl_vector *x, const gsl_vector *y, const gsl_vector *wts, gsl_vector *c, double *chisq, gsl_bspline_workspace *w)¶

These functions fit a B-spline to the input data (x, y) which must be the same length $m$ . Optional weights may be given in the wts vector. The banded normal equations matrix is assembled and factored with a Cholesky factorization, and then solved, storing the coefficients in the vector c, which has length $n$ . The cost function $\chi^2$ is output in chisq. If the normal equations matrix is singular, the error code GSL_EDOM is returned. This could occur, for example, if there is a large gap in the data vectors so that no data points constrain a particular $B_{i,k}(x)$ basis function.

On output, w->XTX contains the Cholesky factor of the symmetric banded normal equations matrix, and this can be passed directly to the functions gsl_bspline_covariance() and gsl_bspline_rcond().

int gsl_bspline_lsnormal(const gsl_vector *x, const gsl_vector *y, const gsl_vector *wts, gsl_vector *XTy, gsl_matrix *XTX, gsl_bspline_workspace *w)¶

This function builds the normal equations matrix $X^T W X$ and right hand side vector $X^T W y$ , which are stored in XTX and XTy on output. The diagonal weight matrix $W$ is specified by the input vector wts, which may be set to NULL, in which case $W = I$ . The inputs x, y, and wts represent the data and weights to be fitted in a least-squares sense, and they all have length $m$ . While the full normal equations matrix has dimension $n$ -by- $n$ , it has a significant number of zeros, and so the output matrix XTX is stored in banded symmetric format. Therefore, its dimensions are $n$ -by- $k$ . The vector XTy has length $n$ . The output matrix XTX can be passed directly into a banded Cholesky factorization routine, such as gsl_linalg_cholesky_band_decomp().

This function is normally called by gsl_bspline_wlssolve() but is available for users with more specialized applications.

int gsl_bspline_lsnormalm(const gsl_vector *x, const gsl_matrix *Y, const gsl_vector *wts, gsl_matrix *XTY, gsl_matrix *XTX, gsl_bspline_workspace *w)¶: This function builds the normal equations matrix $X^T W X$ and right hand side vectors $X^T W Y$ , which are stored in XTX and XTY on output. The diagonal weight matrix $W$ is specified by the input vector wts, which may be set to NULL, in which case $W = I$ . The input Y is a matrix of the right hand side vectors, of size $m$ -by-nrhs. The inputs x and wts represent the data and weights to be fitted in a least-squares sense, and they have length $m$ . While the full normal equations matrix has dimension $n$ -by- $n$ , it has a significant number of zeros, and so the output matrix XTX is stored in banded symmetric format. Therefore, its dimensions are $n$ -by- $k$ . The matrix XTY has size $n$ -by-nrhs. The output matrix XTX can be passed directly into a banded Cholesky factorization routine, such as gsl_linalg_cholesky_band_decomp().

int gsl_bspline_residuals(const gsl_vector *x, const gsl_vector *y, const gsl_vector *c, gsl_vector *r, gsl_bspline_workspace *w)¶: This function computes the residual vector $r_i = y_i - f(x_i)$ by evaluating the spline $f(x)$ using the coefficients c.

int gsl_bspline_covariance(const gsl_matrix *cholesky, gsl_matrix *cov, gsl_bspline_workspace *w)¶: This function computes the $n$ -by- $n$ covariance matrix $\left( X^T W X \right)^{-1}$ using the banded symmetric Cholesky decomposition stored in cholesky, which has dimensions $n$ -by- $k$ . The output is stored in cov. The function gsl_linalg_cholesky_band_decomp() must be called on the normal equations matrix to compute the Cholesky factor prior to calling this function.

int gsl_bspline_err(const double x, const size_t nderiv, const gsl_matrix *cov, double *err, gsl_bspline_workspace *w)¶

This function computes the standard deviation error of a spline (or its derivative),

$f^{(q)}(x) = \sum_{j=1}^n c_j B^{(q)}_j(x)$

at the point x using the covariance matrix cov calculated previously with gsl_bspline_covariance(). The derivative order is specified by nderiv. The error is stored in err on output, which is defined as,

$\delta f = \sqrt{\mathbf{B}^{(q)}(x)^T C \mathbf{B}^{(q)}(x)}$

where $\mathbf{B}^{(q)}(x) = \left( B_1^{(q)}(x), \dots, B_n^{(q)}(x) \right)$ is the vector of B-spline values evaluated at $x$ , and $C$ is the covariance matrix.

int gsl_bspline_rcond(const gsl_matrix *XTX, double *rcond, gsl_bspline_workspace *w)¶: This function estimates the reciprocal condition number (using the 1-norm) of the normal equations matrix $X^T W X$ , using its Cholesky decomposition, which must be computed previously by calling gsl_linalg_cholesky_band_decomp(). The reciprocal condition number estimate, defined as $1 / (||X^T W X||_1 \cdot ||(X^T W X)^{-1}||_1)$ , is stored in rcond on output.

Periodic Spline Fitting¶

Some applications require fitting a spline to a dataset which has an underlying periodicity associated with it. In these cases it is desirable to construct a spline which exhibits the same periodicity. A spline $f(x)$ is called periodic on $[a,b]$ if it satisfies the conditions

$f^{(j)}(a) = f^{(j)}(b), \qquad j = 0, \dots, k - 2$

These conditions can be satisfied with the following constraints on the knot vector and the control points (Dierckx, 1982),

$t_i &= t_{i + n - k + 1} - P, \qquad i = 1, \dots, k-1 \\ t_{n-k+1+i} &= t_i + P, \qquad i = k+1, \dots, 2k-1 \\ c_i &= c_{n - k + 1 + i}, \qquad i = 1, \dots, k - 1$

where $P = b - a$ is the period. The constraint on the control points $c_i$ means that there are only $n-k+1$ independent control points, and then the remaining $k-1$ are determined. The independent control points $c^{*} = \{ c_1, \dots, c_{n-k+1} \}$ can be determined by solving the following least squares problem (Dierckx, 1982):

$\min_{c^{*}} \left|\left| y - X^{*} c^{*} \right|\right|_W^2$

where $X^{*}_{ij} = B^{*}_{j,k}(x_i)$ and

$B^{*}_{j,k}(x) = \left\{ \begin{array}{ll} B_{j,k}(x) + B_{n-k+j+1,k}(x), & j = 1, \dots, k-1 \\ B_{j,k}(x), & j = k, \dots, n-k+1 \end{array} \right.$

int gsl_bspline_plssolve(const gsl_vector *x, const gsl_vector *y, gsl_vector *c, double *chisq, gsl_bspline_workspace *w)¶
int gsl_bspline_pwlssolve(const gsl_vector *x, const gsl_vector *y, const gsl_vector *wts, gsl_vector *c, double *chisq, gsl_bspline_workspace *w)¶: These functions fit a periodic B-spline to the input data (x, y) which must be the same length $m$ . Optional weights may be given in the wts vector. The input vector x must be sorted in ascending order, in order to efficiently compute the $QR$ decomposition of the least squares system. The algorithm computes the $QR$ decomposition of the least squares matrix, one row at a time using Givens transformations, taking advantage of sparse structure due to the local properties of the B-spline basis. On output, the coefficients are stored in the vector c, which has length $n$ . The cost function $\chi^2 = || y - X c ||_W^2$ is output in chisq.

int gsl_bspline_plsqr(const gsl_vector *x, const gsl_vector *y, const gsl_vector *wts, gsl_matrix *R, gsl_vector *QTy, double *rnorm, gsl_bspline_workspace *w)¶: This function calculates the $R$ factor and $Q^T y$ vector for a periodic spline fitting problem on the input data (x, y) with optional weights wts. The input vector x must be sorted in ascending order. On output, the $(n-k+1)$ -by- $(n-k+1)$ $R$ factor is stored in R, while the first $n-k+1$ elements of the vector $Q^T y$ are stored in QTy. The residual norm $||y - X^{*} c^{*}||_W$ is stored in rnorm.

Gram Matrix¶

Let $\mathbf{B}(x) = (B_1(x), B_2(x), \dots, B_n(x))^T$ be the $n$ -vector whose elements are the B-spline basis functions evaluated at $x$ . Then the $n$ -by- $n$ outer product matrix of the B-spline basis at $x$ is defined as

$A(x) = \mathbf{B}(x) \otimes \mathbf{B}(x) = \mathbf{B}(x) \mathbf{B}(x)^T, \quad A_{ij}(x) = B_i(x) B_j(x)$

This can be generalized to higher order derivatives,

$A_{ij}^{(q)}(x) = \left( \frac{d^q}{dx^q} B_i(x) \right) \left( \frac{d^q}{dx^q} B_j(x) \right)$

Then, the $n$ -by- $n$ Gram matrix of the B-spline basis is defined as

$G_{ij}^{(q)} = \int_a^b A_{ij}^{(q)}(x) dx = \int_a^b \left( \frac{d^q}{dx^q} B_i(x) \right) \left( \frac{d^q}{dx^q} B_j(x) \right) dx$

where the integral is usually taken over the full support of the B-splines, i.e. from $\min(\mathbf{t})$ to $\max(\mathbf{t})$ . However it can be computed over any interval $[a,b]$ .

The B-spline outer product and Gram matrices are symmetric and banded with lower bandwidth $k - 1$ . These matrices are positive semi-definite, and in the case of $q = 0$ , $A^{(0)}$ and $G^{(0)}$ are positive definite. These matrices have applications in least-squares fitting.

int gsl_bspline_oprod(const size_t nderiv, const double x, gsl_matrix *A, gsl_bspline_workspace *w)¶: This function calculates the outer product matrix $A^{(q)}(x)$ at the location x, where the derivative order $q$ is specified by nderiv. The matrix is stored in A in symmetric banded format, and so A has dimensions $n$ -by- $k$ . Since the B-spline functions have local support, many elements of the outer product matrix could be zero.

int gsl_bspline_gram(const size_t nderiv, gsl_matrix *G, gsl_bspline_workspace *w)¶

This function calculates the positive semi-definite Gram matrix $G_{ij}^{(q)}$ , where the derivative order $q$ is specified by nderiv. The matrix is stored in G in symmetric banded format, and so G has dimensions $n$ -by- $k$ . The integrals are computed over the full support of the B-spline basis.

The algorithm uses Gauss-Legendre quadrature on each knot interval. See algorithm 5.22 of L. Schumaker, “Spline Functions: Basic Theory”.

int gsl_bspline_gram_interval(const double a, const double b, const size_t nderiv, gsl_matrix *G, gsl_bspline_workspace *w)¶

This function calculates the positive semi-definite Gram matrix $G_{ij}^{(q)}$ , where the derivative order $q$ is specified by nderiv. The matrix is stored in G in symmetric banded format, and so G has dimensions $n$ -by- $k$ . The integrals are computed over the interval $[a,b]$ .

The algorithm uses Gauss-Legendre quadrature on each knot interval. See algorithm 5.22 of L. Schumaker, “Spline Functions: Basic Theory”.

Interpolation with B-splines¶

Given a set of $n$ data points, $(x_i,y_i)$ , we can define an interpolating spline with $n$ control points which passes through each data point,

$f(x_i) = \sum_{j=1}^n c_j B_j(x_i) = y_i, \quad i = 1,\dots,n$

This is a square system of equations which can be solved by inverting the $n$ -by- $n$ collocation matrix, whose elements are $B_j(x_i)$ . In order for the collocation matrix to be invertible, the $x_i$ and spline knots must satisfy the Schoenberg-Whitney conditions, which are

$t_i < x_i < t_{i+k}, \quad i = 1,\dots,n$

If these conditions are met, then the collocation matrix is banded, with both the upper and lower bandwidths equal to $k-1$ . This property leads to savings in storage and computation when solving interpolation problems. To construct a knot vector which satisfies these conditions, see the function gsl_bspline_init_interp().

int gsl_bspline_col_interp(const gsl_vector *x, gsl_matrix *XB, gsl_bspline_workspace *w)¶: This function constructs the banded collocation matrix to interpolate the $n$ abscissa values provided in x. The values in x must be sorted in ascending order, although this is not checked explicitly by the function. The collocation matrix is stored in XB on output in banded storage format suitable for use with the banded LU routines. Therefore, XB must have dimension $n$ -by- $3(k-1) + 1$ .

Hermite Interpolation with B-splines¶

Given a set of $n$ data interpolation sites, $x_1 < x_2 < \dots < x_n$ , and a set of interpolation values,

$f_i^{(j)}, \qquad i=1,\dots,n, j=0,1,\dots,m$

then there is a unique interpolating spline of order $k = 2m+2$ which will match the given function values (and their derivatives). The spline has $(m+1)n$ control points and is given by,

$f(x) = \sum_{j=1}^{(m+1)n} c_j B_j(x)$

This spline will satisfy,

$f^{(j)}(x_i) = f_i^{(j)}, \qquad i=1,\dots,n, j=0,1,\dots,m$

Mummy (1989) provides an efficient algorithm to determine the spline coefficients $c_j$ .

int gsl_bspline_interp_chermite(const gsl_vector *x, const gsl_vector *y, const gsl_vector *dy, gsl_vector *c, const gsl_bspline_workspace *w)¶

This function will calculate the coefficients of an interpolating cubic Hermite spline which match the provided function values $y_i$ and first derivatives $y'_i$ at the interpolation sites $x_i$ . The inputs x, y, dy are all of the same length $n$ and contain the $x_i$ , $y_i$ , and $y'_i$ respectively. On output, the spline coefficients are stored in c, which must have length $2n$ .

This function is specifically designed for cubic Hermite splines, and so the spline order must be $k = 4$ . The function gsl_bspline_init_hermite() must be called to initialize the knot vector. The algorithm used is provided by Mummy (1989).

Projection onto the B-spline Basis¶

A function $f(x)$ can be projected onto the subspace $V = \textrm{span}(B_i)$ spanned by the B-splines as follows,

$f(x) \approx \textrm{proj}_V(f) = \sum_j c_j B_j(x)$

If $f(x)$ is a piecewise polynomial of degree less than $k$ , then $f(x)$ will lie in $V$ and $f(x) = \textrm{proj}_V(f)$ . The expansion coefficients $c_j$ can be calculated as follows,

$f(x) B_i(x) &= \sum_j c_j B_j(x) B_i(x) \\ \underbrace{\int f(x) B_i(x) dx}_{y_i} &= \sum_j c_j \underbrace{\int B_i(x) B_j(x) dx}_{G_{ij}} \\ G c &= y$

where $G$ is the B-spline Gram Matrix. The above equation can be solved using a Cholesky factorization of $G$ .

int gsl_bspline_proj_rhs(const gsl_function *F, gsl_vector *y, gsl_bspline_workspace *w)¶: This function computes the right hand side vector $y_i = \int f(x) B_i(x) dx$ for projecting a function F onto the B-spline basis. The output is stored in y which must have size ncontrol. The integrals are computed using Gauss-Legendre quadrature, which will produce exact results if $f(x)$ is a piecewise polynomial of degree less than $k$ .

Greville abscissae¶

The Greville abscissae are defined to be the mean location of $k-1$ consecutive knots in the knot vector for each basis spline function of order $k$ . With the first and last knots in the knot vector excluded, there are $n$ Greville abscissae for any given B-spline basis. These values are often used in B-spline collocation applications and may also be called Marsden-Schoenberg points.

double gsl_bspline_greville_abscissa(const size_t i, gsl_bspline_workspace *w)¶: Returns the location of the $i$ -th Greville abscissa for the given B-spline basis. For the ill-defined case when $k = 1$ , the implementation chooses to return breakpoint interval midpoints.

Examples¶

Example 1: Basis Splines and Uniform Knots¶

The following program computes and outputs linear, quadratic, cubic, and quartic basis splines using uniform knots on the interval $[0,1]$ . The knot vector is

$\mathbf{t} = \left( \underbrace{0, \dots, 0}_{k}, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, \underbrace{1, \dots, 1}_{k} \right)$

where the first and last knot are repeated $k$ times, where $k$ is the spline order. The output is shown in Fig. 40.

Fig. 40 Linear, quadratic, cubic, and quartic basis splines defined on the interval $[0,1]$ with uniform knots. Black circles denote the location of the knots.¶

This figure makes it clear that each B-spline has local support and spans $k + 1$ knots. The B-splines near the interval endpoints appear to span less knots, but note that the endpoint knots have multiplicities of $k$ .

The program is given below.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>

void
print_basis(FILE *fp_knots, FILE *fp_spline, const size_t order, const size_t nbreak)
{
  gsl_bspline_workspace *w = gsl_bspline_alloc(order, nbreak);
  const size_t p = gsl_bspline_ncontrol(w);
  const size_t n = 300;
  const double a = 0.0;
  const double b = 1.0;
  const double dx = (b - a) / (n - 1.0);
  gsl_vector *B = gsl_vector_alloc(p);
  size_t i, j;

  /* use uniform breakpoints on [a, b] */
  gsl_bspline_init_uniform(a, b, w);

  gsl_vector_fprintf(fp_knots, w->knots, "%f");

  for (i = 0; i < n; ++i)
    {
      double xi = i * dx;

      gsl_bspline_eval_basis(xi, B, w);

      fprintf(fp_spline, "%f ", xi);

      for (j = 0; j < p; ++j)
        fprintf(fp_spline, "%f ", gsl_vector_get(B, j));

      fprintf(fp_spline, "\n");
    }

  fprintf(fp_knots, "\n\n");
  fprintf(fp_spline, "\n\n");

  gsl_vector_free(B);
  gsl_bspline_free(w);
}

int
main (void)
{
  const size_t nbreak = 11; /* number of breakpoints */
  FILE *fp_knots = fopen("bspline1_knots.txt", "w");
  FILE *fp_spline = fopen("bspline1.txt", "w");

  print_basis(fp_knots, fp_spline, 2, nbreak); /* linear splines */
  print_basis(fp_knots, fp_spline, 3, nbreak); /* quadratic splines */
  print_basis(fp_knots, fp_spline, 4, nbreak); /* cubic splines */
  print_basis(fp_knots, fp_spline, 5, nbreak); /* quartic splines */

  fclose(fp_knots);
  fclose(fp_spline);

  return 0;
}

Example 2: Derivatives of Basis Splines¶

The following program computes and outputs cubic B-splines and their derivatives using $6$ breakpoints and uniform knots on the interval $[0,1]$ . All derivatives up to order $3$ are computed. The output is shown in Fig. 41.

Fig. 41 Cubic B-splines and derivatives up to order 3. Black circles denote the location of the knots.¶

The program is given below.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_bspline.h>

int
main (void)
{
  const size_t nbreak = 6;
  const size_t spline_order = 4;
  gsl_bspline_workspace *w = gsl_bspline_alloc(spline_order, nbreak);
  const size_t p = gsl_bspline_ncontrol(w);
  const size_t n = 300;
  const double a = 0.0;
  const double b = 1.0;
  const double dx = (b - a) / (n - 1.0);
  gsl_matrix *dB = gsl_matrix_alloc(p, spline_order);
  size_t i, j, k;

  /* uniform breakpoints on [a, b] */
  gsl_bspline_init_uniform(a, b, w);

  /* output knot vector */
  gsl_vector_fprintf(stdout, w->knots, "%f");
  printf("\n\n");

  for (i = 0; i < spline_order; ++i)
    {
      for (j = 0; j < n; ++j)
        {
          double xj = j * dx;

          gsl_bspline_eval_deriv_basis(xj, i, dB, w);

          printf("%f ", xj);

          for (k = 0; k < p; ++k)
            printf("%f ", gsl_matrix_get(dB, k, i));

          printf("\n");
        }

      printf("\n\n");
    }

  gsl_matrix_free(dB);
  gsl_bspline_free(w);

  return 0;
}

Example 3: Least squares and breakpoints¶

The following program computes a linear least squares fit to data using cubic B-splines with uniform breakpoints. The data is generated from the curve $y(x) = \cos{x} \exp{\left( -x/10 \right)}$ with Gaussian noise added. Splines are fitted with $10$ and $40$ breakpoints for comparison.

The program output is shown below:

$ ./a.out > bspline_lsbreak.txt
40 breakpoints: chisq/dof = 1.008999e+00
10 breakpoints: chisq/dof = 1.014761e+00

The data and fitted models are shown in Fig. 42.

Fig. 42 Cubic B-spline least squares fit¶

The program is given below.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>

int
main (void)
{
  const size_t n = 500;     /* number of data points to fit */
  const size_t k = 4;       /* spline order */
  const double a = 0.0;     /* data interval [a,b] */
  const double b = 15.0;
  const double sigma = 0.2; /* noise */
  gsl_bspline_workspace *work1 = gsl_bspline_alloc(k, 40);   /* 40 breakpoints */
  gsl_bspline_workspace *work2 = gsl_bspline_alloc(k, 10);   /* 10 breakpoints */
  const size_t p1 = gsl_bspline_ncontrol(work1);             /* number of control points */
  const size_t p2 = gsl_bspline_ncontrol(work2);             /* number of control points */
  const size_t dof1 = n - p1;                                /* degrees of freedom */
  const size_t dof2 = n - p2;                                /* degrees of freedom */
  gsl_vector *c1 = gsl_vector_alloc(p1);
  gsl_vector *c2 = gsl_vector_alloc(p2);
  gsl_vector *x = gsl_vector_alloc(n);
  gsl_vector *y = gsl_vector_alloc(n);
  gsl_vector *w = gsl_vector_alloc(n);
  size_t i;
  gsl_rng *r;
  double chisq1, chisq2;

  gsl_rng_env_setup();
  r = gsl_rng_alloc(gsl_rng_default);

  /* this is the data to be fitted */
  for (i = 0; i < n; ++i)
    {
      double xi = (b - a) / (n - 1.0) * i + a;
      double yi = cos(xi) * exp(-0.1 * xi);
      double dyi = gsl_ran_gaussian(r, sigma);

      yi += dyi;

      gsl_vector_set(x, i, xi);
      gsl_vector_set(y, i, yi);
      gsl_vector_set(w, i, 1.0 / (sigma * sigma));

      printf("%f %f\n", xi, yi);
    }

  /* use uniform breakpoints on [a, b] */
  gsl_bspline_init_uniform(a, b, work1);
  gsl_bspline_init_uniform(a, b, work2);

  /* solve least squares problem */
  gsl_bspline_wlssolve(x, y, w, c1, &chisq1, work1);
  gsl_bspline_wlssolve(x, y, w, c2, &chisq2, work2);

  fprintf(stderr, "40 breakpoints: chisq/dof = %e\n", chisq1 / dof1);
  fprintf(stderr, "10 breakpoints: chisq/dof = %e\n", chisq2 / dof2);

  printf("\n\n");

  /* output the spline curves */
  {
    double xi;

    for (xi = a; xi <= b; xi += 0.1)
      {
        double result1, result2;

        gsl_bspline_calc(xi, c1, &result1, work1);
        gsl_bspline_calc(xi, c2, &result2, work2);
        printf("%f %f %f\n", xi, result1, result2);
      }
  }

  gsl_rng_free(r);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(w);
  gsl_vector_free(c1);
  gsl_vector_free(c2);
  gsl_bspline_free(work1);
  gsl_bspline_free(work2);

  return 0;
}

Example 4: Least squares and spline order¶

The following program computes least squares fits to the same dataset in the previous example, using $10$ uniform breakpoints and computing splines of order $1,2,3,4,5$ .

The data and fitted models are shown in Fig. 43.

Fig. 43 B-spline least squares fit with varying spline orders¶

The program is given below.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>

int
main (void)
{
  const size_t n = 500;              /* number of data points to fit */
  const double a = 0.0;              /* data interval [a,b] */
  const double b = 15.0;
  const size_t max_spline_order = 5; /* maximum spline order */
  const size_t nbreak = 10;          /* number of breakpoints */
  const double sigma = 0.2;          /* noise */
  gsl_bspline_workspace **work = malloc(max_spline_order * sizeof(gsl_bspline_workspace *));
  gsl_vector **c = malloc(max_spline_order * sizeof(gsl_vector *));
  gsl_vector *x = gsl_vector_alloc(n);
  gsl_vector *y = gsl_vector_alloc(n);
  gsl_vector *w = gsl_vector_alloc(n);
  size_t i;
  gsl_rng *r;
  double chisq;

  gsl_rng_env_setup();
  r = gsl_rng_alloc(gsl_rng_default);

  /* this is the data to be fitted */
  for (i = 0; i < n; ++i)
    {
      double xi = (b - a) / (n - 1.0) * i + a;
      double yi = cos(xi) * exp(-0.1 * xi);
      double dyi = gsl_ran_gaussian(r, sigma);

      yi += dyi;

      gsl_vector_set(x, i, xi);
      gsl_vector_set(y, i, yi);
      gsl_vector_set(w, i, 1.0 / (sigma * sigma));

      printf("%f %f\n", xi, yi);
    }

  printf("\n\n");

  for (i = 0; i < max_spline_order; ++i)
    {
      /* allocate workspace for this spline order */
      work[i] = gsl_bspline_alloc(i + 1, nbreak);
      c[i] = gsl_vector_alloc(gsl_bspline_ncontrol(work[i]));

      /* use uniform breakpoints on [a, b] */
      gsl_bspline_init_uniform(a, b, work[i]);

      /* solve least squares problem */
      gsl_bspline_wlssolve(x, y, w, c[i], &chisq, work[i]);
    }

  /* output the spline curves */
  {
    double xi;

    for (xi = a; xi <= b; xi += 0.1)
      {
        printf("%f ", xi);
        for (i = 0; i < max_spline_order; ++i)
          {
            double result;
            gsl_bspline_calc(xi, c[i], &result, work[i]);
            printf("%f ", result);
          }
        printf("\n");
      }
  }

  for (i = 0; i < max_spline_order; ++i)
    {
      gsl_vector_free(c[i]);
      gsl_bspline_free(work[i]);
    }

  gsl_rng_free(r);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(w);
  free(work);
  free(c);

  return 0;
}

Example 5: Least squares, regularization and the Gram matrix¶

The following example program fits a cubic B-spline to data generated from the Gaussian

$g(x) = e^{-x^2}$

on the interval $[-1.5,1.5]$ with noise added. Two data gaps are deliberately introduced on $[-1.1,-0.7]$ and $[0.1,0.55]$ , which makes the problem ill-posed. Fig. 44 shows that the ordinary least squares solution exhibits large errors in the data gap regions due to the ill-conditioned least squares matrix.

One way to correct the situation is to minimize the data misfit but also minimize the curvature of the spline averaged over the full interval. We can define the average curvature using the second derivative of the spline as follows:

$\textrm{curvature} &\approx \int \left| \frac{d^2}{dx^2} f(x) \right|^2 dx \\ &= \int \left| \frac{d^2}{dx^2} \sum_{i=1}^n c_i B_i(x) \right|^2 dx \\ &= \int \left| \sum_{i=1}^n c_i B''_i(x) \right|^2 dx \\ &= \int \left| c^T B''(x) \right|^2 dx \\ &= \int c^T B''(x) B''^T(x) c dx \\ &= c^T \left( \int B''(x) B''^T(x) dx \right) c \\ &= c^T G^{(2)} c$

where $G^{(2)}$ is the Gram matrix of second derivatives of the B-splines. Therefore, our regularized least squares problem is

$\min_{c} || y - X c ||_W^2 + \lambda^2 c^T G^{(2)} c$

where the regularization parameter $\lambda^2$ represents a tradeoff between minimizing the data misfit and minimizing the spline curvature. The solution of this least-squares problem is

$c = \left( X^T W X + \lambda^2 G^{(2)} \right)^{-1} X^T W y$

The example program below solves this problem without regularization $(\lambda^2 = 0)$ and with regularization $(\lambda^2 = 0.1)$ for comparison. We also plot the $1 \sigma$ confidence intervals for both solutions using the diagonal elements of the covariance matrix $\textrm{Cov} = \left( X^T W X + \lambda^2 G^{(2)} \right)^{-1}$ .

Fig. 44 shows the result of regularizing the spline over the full interval, using the command

./a.out > data.txt

Fig. 44 B-spline least squares fit with exact (green), unregularized (red with $1 \sigma$ confidence intervals), and fully regularized (blue with $1 \sigma$ confidence intervals) solutions.¶

Fig. 45 shows the result of regularizing the spline over the smaller interval $[0,0.6]$ , correcting only the second gap, using the command

./a.out 0 0.6 > data.txt

Fig. 45 B-spline least squares fit with exact (green), unregularized (red with $1 \sigma$ confidence intervals), and partially regularized (blue with $1 \sigma$ confidence intervals) solutions.¶

The program is given below.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int
main (int argc, char * argv[])
{
  const size_t n = 500;                                     /* number of data points to fit */
  const size_t k = 4;                                       /* spline order */
  const double a = -1.5;                                    /* data interval [a,b] */
  const double b = 1.5;
  const double sigma = 0.02;                                /* noise */
  const double lambda_sq = 0.1;                             /* regularization parameter */
  gsl_bspline_workspace *work = gsl_bspline_alloc(k, 25);   /* 25 breakpoints */
  const size_t p = gsl_bspline_ncontrol(work);              /* number of control points */
  gsl_vector *c = gsl_vector_alloc(p);                      /* control points for unregularized model */
  gsl_vector *creg = gsl_vector_alloc(p);                   /* control points for regularized model */
  gsl_vector *x = gsl_vector_alloc(n);                      /* x data */
  gsl_vector *y = gsl_vector_alloc(n);                      /* y data */
  gsl_vector *w = gsl_vector_alloc(n);                      /* data weights */
  gsl_matrix *XTX = gsl_matrix_alloc(p, k);                 /* normal equations matrix */
  gsl_vector *XTy = gsl_vector_alloc(p);                    /* normal equations rhs */
  gsl_matrix *cov = gsl_matrix_alloc(p, p);                 /* covariance matrix */
  gsl_matrix *cov_reg = gsl_matrix_alloc(p, p);             /* covariance matrix for regularized model */
  gsl_matrix *G = gsl_matrix_alloc(p, k);                   /* regularization matrix */
  double reg_a = -1.0;                                      /* regularization interval [reg_a,reg_b] */
  double reg_b = -1.0;
  int reg_interval = 0;                                     /* apply regularization on smaller interval */
  size_t i;
  gsl_rng *r;

  if (argc > 2)
    {
      reg_a = atof(argv[1]);
      reg_b = atof(argv[2]);
      reg_interval = 1;
    }

  gsl_rng_env_setup();
  r = gsl_rng_alloc(gsl_rng_default);

  /* this is the data to be fitted */
  i = 0;
  while (i < n)
    {
      double ui = gsl_rng_uniform(r);
      double xi = (b - a) * ui + a;
      double yi, dyi;

      /* data gaps between [-1.1,-0.7] and [0.1,0.5] */
      if ((xi >= -1.1 && xi <= -0.7) || (xi >= 0.1 && xi <= 0.55))
        continue;

      dyi = gsl_ran_gaussian(r, sigma);
      yi = exp(-xi*xi) + dyi;

      gsl_vector_set(x, i, xi);
      gsl_vector_set(y, i, yi);
      gsl_vector_set(w, i, 1.0 / (sigma * sigma));

      printf("%f %f\n", xi, yi);

      ++i;
    }

  printf("\n\n");

  /* use uniform breakpoints on [a, b] */
  gsl_bspline_init_uniform(a, b, work);

  gsl_bspline_lsnormal(x, y, w, XTy, XTX, work); /* form normal equations matrix */
  gsl_linalg_cholesky_band_decomp(XTX);          /* banded Cholesky decomposition */
  gsl_linalg_cholesky_band_solve(XTX, XTy, c);   /* solve for unregularized solution */
  gsl_bspline_covariance(XTX, cov, work);        /* compute covariance matrix */

  /* compute regularization matrix */
  if (reg_interval)
    gsl_bspline_gram_interval(reg_a, reg_b, 2, G, work);
  else
    gsl_bspline_gram(2, G, work);

  /* multiply by lambda^2 */
  gsl_matrix_scale(G, lambda_sq);

  gsl_bspline_lsnormal(x, y, w, XTy, XTX, work);  /* form normal equations matrix */
  gsl_matrix_add(XTX, G);                         /* add regularization term */
  gsl_linalg_cholesky_band_decomp(XTX);           /* banded Cholesky decomposition */
  gsl_linalg_cholesky_band_solve(XTX, XTy, creg); /* solve for regularized solution */
  gsl_bspline_covariance(XTX, cov_reg, work);     /* compute covariance matrix */

  /* output the spline curves */
  {
    double xi;

    for (xi = a; xi <= b; xi += 0.01)
      {
        double result_unreg, result_reg;
        double err_unreg, err_reg;

        /* compute unregularized spline value and error at xi */
        gsl_bspline_calc(xi, c, &result_unreg, work);
        gsl_bspline_err(xi, 0, cov, &err_unreg, work);

        /* compute regularized spline value and error at xi */
        gsl_bspline_calc(xi, creg, &result_reg, work);
        gsl_bspline_err(xi, 0, cov_reg, &err_reg, work);

        printf("%f %e %e %e %e %e\n",
               xi,
               exp(-xi*xi),
               result_unreg,
               err_unreg,
               result_reg,
               err_reg);
      }
  }

  gsl_rng_free(r);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(w);
  gsl_vector_free(c);
  gsl_vector_free(creg);
  gsl_bspline_free(work);
  gsl_matrix_free(XTX);
  gsl_vector_free(XTy);
  gsl_matrix_free(cov);
  gsl_matrix_free(cov_reg);
  gsl_matrix_free(G);

  return 0;
}

Example 6: Least squares and endpoint regularization¶

Data fitting with splines can often produce undesired effects at the spline endpoints, where they are less constrained by data. As an example, we fit order $k = 10$ splines to the Runge function,

$g(x) = \frac{1}{1 + 25 x^2}$

with 20 uniform breakpoints on the interval $[a,b] = [-1,1]$ . Fig. 46 shows that the fitted spline (red) exhibits a sharp upward feature at the left endpoint. This is due to the lack of data constraining this part of the spline, and would likely result in errors if the user wanted to extrapolate the spline outside the fit interval. One method to correct this is to apply regularization to the spline at the endpoints by minimizing one or more of its derivatives. In this example, we minimize the first derivative of the spline at both endpoints, $|f'(a)|^2$ and $|f'(b)|^2$ . Note that

$\left| \frac{d}{dx} f(x) \right|^2 &= \left| \frac{d}{dx} \sum_{i=1}^n c_i B_i(x) \right|^2 \\ &= \left| c^T B'(x) \right|^2 \\ &= c^T B'(x) B'^T(x) c \\ &= c^T A^{(1)}(x) c$

where $A^{(1)}$ is the outer product matrix of first derivatives. Therefore, our regularized least squares problem is

$\min_{c} || y - X c ||_W^2 + \lambda^2 c^T A^{(1)}(a) c + \lambda^2 c^T A^{(1)}(b) c$

where we have chosen to use the same regularization parameter $\lambda^2$ to regularize both ends of the spline. In the example program below, we compute the unregularized solution ( $\lambda^2 = 0$ ) and regularized solution ( $\lambda^2 = 10$ ). The program outputs the first derivative of the spline at both endpoints:

unregularized endpoint deriv 1: [ -1.081170e+01,  -2.963725e+00]
  regularized endpoint deriv 1: [ -2.735857e-02,  -5.371758e-03]

showing that the regularization has forced the first derivative smaller at both endpoints. The results are shown in Fig. 46.

Fig. 46 Runge function B-spline least squares fit with exact (green), unregularized (red) and regularized (blue) solutions. The zoomed inset views show the spline endpoints with $1 \sigma$ confidence intervals.¶

The program is given below.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

double
f(double x)
{
  return (1.0 / (1.0 + 25.0*x*x));
}

int
main (void)
{
  const size_t n = 500;                                     /* number of data points to fit */
  const size_t k = 10;                                      /* spline order */
  const double a = -1.0;                                    /* data interval [a,b] */
  const double b = 1.0;
  const double sigma = 0.03;                                /* noise */
  const double lambda_sq = 10.0;                            /* regularization parameter */
  const size_t nderiv = 1;                                  /* derivative order to regularize at endpoints */
  gsl_bspline_workspace *work = gsl_bspline_alloc(k, 20);   /* 20 breakpoints */
  const size_t p = gsl_bspline_ncontrol(work);              /* number of control points */
  gsl_vector *c = gsl_vector_alloc(p);                      /* control points for unregularized model */
  gsl_vector *creg = gsl_vector_alloc(p);                   /* control points for regularized model */
  gsl_vector *x = gsl_vector_alloc(n);                      /* x data */
  gsl_vector *y = gsl_vector_alloc(n);                      /* y data */
  gsl_vector *w = gsl_vector_alloc(n);                      /* data weights */
  gsl_matrix *XTX = gsl_matrix_alloc(p, k);                 /* normal equations matrix */
  gsl_vector *XTy = gsl_vector_alloc(p);                    /* normal equations rhs */
  gsl_matrix *cov = gsl_matrix_alloc(p, p);                 /* covariance matrix */
  gsl_matrix *cov_reg = gsl_matrix_alloc(p, p);             /* covariance matrix for regularized model */
  gsl_matrix *A1 = gsl_matrix_alloc(p, k);                  /* left regularization matrix */
  gsl_matrix *A2 = gsl_matrix_alloc(p, k);                  /* right regularization matrix */
  size_t i;
  gsl_rng *r;

  gsl_rng_env_setup();
  r = gsl_rng_alloc(gsl_rng_default);

  /* this is the data to be fitted */
  i = 0;
  while (i < n)
    {
      double ui = gsl_rng_uniform(r);
      double xi = (b - a) * ui + a;
      double yi, dyi;

      dyi = gsl_ran_gaussian(r, sigma);
      yi = f(xi) + dyi;

      gsl_vector_set(x, i, xi);
      gsl_vector_set(y, i, yi);
      gsl_vector_set(w, i, 1.0 / (sigma * sigma));

      printf("%f %f\n", xi, yi);

      ++i;
    }

  printf("\n\n");

  /* use uniform breakpoints on [a, b] */
  gsl_bspline_init_uniform(a, b, work);

  gsl_bspline_lsnormal(x, y, w, XTy, XTX, work); /* form normal equations matrix */
  gsl_linalg_cholesky_band_decomp(XTX);          /* banded Cholesky decomposition */
  gsl_linalg_cholesky_band_solve(XTX, XTy, c);   /* solve for unregularized solution */
  gsl_bspline_covariance(XTX, cov, work);        /* compute covariance matrix */

  /* compute regularization matrices */
  gsl_bspline_oprod(nderiv, a, A1, work);
  gsl_bspline_oprod(nderiv, b, A2, work);

  /* multiply by lambda^2 */
  gsl_matrix_scale(A1, lambda_sq);
  gsl_matrix_scale(A2, lambda_sq);

  gsl_bspline_lsnormal(x, y, w, XTy, XTX, work);  /* form normal equations matrix */
  gsl_matrix_add(XTX, A1);                        /* add regularization terms */
  gsl_matrix_add(XTX, A2);
  gsl_linalg_cholesky_band_decomp(XTX);           /* banded Cholesky decomposition */
  gsl_linalg_cholesky_band_solve(XTX, XTy, creg); /* solve for regularized solution */
  gsl_bspline_covariance(XTX, cov_reg, work);     /* compute covariance matrix */

  /* output the spline curves */
  {
    double xi;

    for (xi = a; xi <= b; xi += 0.001)
      {
        double result_unreg, result_reg;
        double err_unreg, err_reg;

        /* compute unregularized spline value and error at xi */
        gsl_bspline_calc(xi, c, &result_unreg, work);
        gsl_bspline_err(xi, 0, cov, &err_unreg, work);

        /* compute regularized spline value and error at xi */
        gsl_bspline_calc(xi, creg, &result_reg, work);
        gsl_bspline_err(xi, 0, cov_reg, &err_reg, work);

        printf("%f %e %e %e %e %e\n",
               xi,
               f(xi),
               result_unreg,
               err_unreg,
               result_reg,
               err_reg);
      }
  }

  {
    double result0, result1;

    gsl_bspline_calc_deriv(a, c, nderiv, &result0, work);
    gsl_bspline_calc_deriv(b, c, nderiv, &result1, work);
    fprintf(stderr, "unregularized endpoint deriv %zu: [%14.6e, %14.6e]\n", nderiv, result0, result1);

    gsl_bspline_calc_deriv(a, creg, nderiv, &result0, work);
    gsl_bspline_calc_deriv(b, creg, nderiv, &result1, work);
    fprintf(stderr, "  regularized endpoint deriv %zu: [%14.6e, %14.6e]\n", nderiv, result0, result1);
  }

  gsl_rng_free(r);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(w);
  gsl_vector_free(c);
  gsl_vector_free(creg);
  gsl_bspline_free(work);
  gsl_matrix_free(XTX);
  gsl_vector_free(XTy);
  gsl_matrix_free(cov);
  gsl_matrix_free(cov_reg);
  gsl_matrix_free(A1);
  gsl_matrix_free(A2);

  return 0;
}

Example 7: Least squares and periodic splines¶

The following example program fits B-spline curves to a dataset which is periodic. Both a non-periodic and periodic spline are fitted to the data for comparison. The data is derived from the function on $[0,2\pi]$ ,

$f(x) = \sin{x} - \cos{2x}$

with Gaussian noise added. The program fits order $k=6$ splines to the dataset.

Fig. 47 Data (black) with non-periodic spline fit (red) and periodic spline fit (green).¶

The program additionally computes the derivatives up to order $5$ at the endpoints $x=0$ and $x=2\pi$ :

=== Non-periodic spline endpoint derivatives ===
deriv 0: [ -8.697939e-01,  -8.328553e-01]
deriv 1: [ -2.423132e+00,   4.277439e+00]
deriv 2: [  3.904362e+01,   3.710592e+01]
deriv 3: [ -2.096142e+02,   2.036125e+02]
deriv 4: [  7.240121e+02,   7.384888e+02]
deriv 5: [ -1.230036e+03,   1.294264e+03]
=== Periodic spline endpoint derivatives ===
deriv 0: [ -1.022798e+00,  -1.022798e+00]
deriv 1: [  1.041910e+00,   1.041910e+00]
deriv 2: [  4.466384e+00,   4.466384e+00]
deriv 3: [ -1.356654e+00,  -1.356654e+00]
deriv 4: [ -2.751544e+01,  -2.751544e+01]
deriv 5: [  5.055419e+01,   2.674261e+01]

We can see that the periodic fit matches derivative values up to order $4$ , with the $5$ -th derivative being non-continuous.

The program is given below.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>

int
main (void)
{
  const size_t n = 500;              /* number of data points to fit */
  const double a = 0.0;              /* data interval [a,b] */
  const double b = 2.0 * M_PI;
  const size_t spline_order = 6;     /* spline order */
  const size_t ncontrol = 15;        /* number of control points */
  const double sigma = 0.2;          /* noise */
  gsl_bspline_workspace *w = gsl_bspline_alloc_ncontrol(spline_order, ncontrol);
  gsl_bspline_workspace *wper = gsl_bspline_alloc_ncontrol(spline_order, ncontrol);
  gsl_vector *c = gsl_vector_alloc(ncontrol);    /* non-periodic coefficients */
  gsl_vector *cper = gsl_vector_alloc(ncontrol); /* periodic coefficients */
  gsl_vector *x = gsl_vector_alloc(n);
  gsl_vector *y = gsl_vector_alloc(n);
  gsl_vector *wts = gsl_vector_alloc(n);
  size_t i;
  gsl_rng *r;
  double chisq;

  gsl_rng_env_setup();
  r = gsl_rng_alloc(gsl_rng_default);

  /* this is the data to be fitted */
  for (i = 0; i < n; ++i)
    {
      double xi = (b - a) / (n - 1.0) * i + a;
      double yi = sin(xi) - cos(2.0 * xi);
      double dyi = gsl_ran_gaussian(r, sigma);

      yi += dyi;

      gsl_vector_set(x, i, xi);
      gsl_vector_set(y, i, yi);
      gsl_vector_set(wts, i, 1.0 / (sigma * sigma));

      printf("%f %f\n", xi, yi);
    }

  printf("\n\n");

  /* use uniform non-periodic knots on [a, b] */
  gsl_bspline_init_uniform(a, b, w);

  /* solve least squares problem for non-periodic spline */
  gsl_bspline_wlssolve(x, y, wts, c, &chisq, w);

  /* use periodic knots on [a, b] */
  gsl_bspline_init_periodic(a, b, wper);

  /* solve least squares problem for periodic spline */
  gsl_bspline_pwlssolve(x, y, wts, cper, &chisq, wper);

  /* output the spline curves */
  {
    double xi;

    for (xi = a; xi <= b; xi += 0.01)
      {
        double result, result_per;
        gsl_bspline_calc(xi, c, &result, w);
        gsl_bspline_calc(xi, cper, &result_per, wper);
        printf("%f %f %f\n", xi, result, result_per);
      }
  }

  fprintf(stderr, "=== Non-periodic spline endpoint derivatives ===\n");

  for (i = 0; i < spline_order; ++i)
    {
      double result0, result1;

      gsl_bspline_calc_deriv(a, c, i, &result0, w);
      gsl_bspline_calc_deriv(b, c, i, &result1, w);

      fprintf(stderr, "deriv %zu: [%14.6e, %14.6e]\n", i, result0, result1);
    }

  fprintf(stderr, "=== Periodic spline endpoint derivatives ===\n");

  for (i = 0; i < spline_order; ++i)
    {
      double result0, result1;

      gsl_bspline_calc_deriv(a, cper, i, &result0, wper);
      gsl_bspline_calc_deriv(b, cper, i, &result1, wper);

      fprintf(stderr, "deriv %zu: [%14.6e, %14.6e]\n", i, result0, result1);
    }

  gsl_vector_free(c);
  gsl_vector_free(cper);
  gsl_bspline_free(w);
  gsl_bspline_free(wper);
  gsl_rng_free(r);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(wts);

  return 0;
}

Example 8: Projecting onto the B-spline basis¶

The following example program projects the function

$f(x) = 3 x^3 - 2 x^2 - 7 x$

onto a cubic B-spline basis with uniform knots on the interval $[-2,2]$ . Because the function is cubic, it can be exactly decomposed in the B-spline basis.

Fig. 48 Exact function (red) and spline projection (green). Small offset deliberately added to green curve.¶

The program is given below.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>

/* function to be projected */
double
f(double x, void * params)
{
  (void) params;
  return (x * (-7.0 + x*(-2.0 + 3.0*x)));
}

int
main (void)
{
  const size_t k = 4;                                       /* spline order */
  const double a = -2.0;                                    /* spline interval [a,b] */
  const double b = 2.0;
  gsl_bspline_workspace *work = gsl_bspline_alloc(k, 10);   /* 10 breakpoints */
  const size_t n = gsl_bspline_ncontrol(work);              /* number of control points */
  gsl_vector *c = gsl_vector_alloc(n);                      /* control points for spline */
  gsl_vector *y = gsl_vector_alloc(n);                      /* rhs vector */
  gsl_matrix *G = gsl_matrix_alloc(n, k);                   /* Gram matrix */
  gsl_function F;

  F.function = f;
  F.params = NULL;

  /* use uniform breakpoints on [a, b] */
  gsl_bspline_init_uniform(a, b, work);

  /* compute Gram matrix */
  gsl_bspline_gram(0, G, work);

  /* construct rhs vector */
  gsl_bspline_proj_rhs(&F, y, work);

  /* solve system */
  gsl_linalg_cholesky_band_decomp(G);
  gsl_linalg_cholesky_band_solve(G, y, c);

  /* output the result */
  {
    double x;

    for (x = a; x <= b; x += 0.01)
      {
        double result;
        gsl_bspline_calc(x, c, &result, work);
        printf("%f %f %f\n", x, GSL_FN_EVAL(&F, x), result);
      }
  }

  gsl_vector_free(y);
  gsl_vector_free(c);
  gsl_matrix_free(G);
  gsl_bspline_free(work);

  return 0;
}

Example 9: Interpolation¶

The following example program computes an interpolating cubic B-spline for a set of data. It constructs the collocation matrix for the dataset and factors it with a banded LU decomposition. The results are shown in Fig. 49.

Fig. 49 Data points are shown in black and interpolated spline in light blue.¶

The program is given below.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>

int
main (void)
{
  const size_t n = 9;                  /* number of data to interpolate */
  const double x_data[] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
  const double y_data[] = { 3.0, 2.9, 2.5, 1.0, 0.9, 0.8, 0.5, 0.2, 0.1 };
  gsl_vector_const_view xv = gsl_vector_const_view_array(x_data, n);
  gsl_vector_const_view yv = gsl_vector_const_view_array(y_data, n);
  const size_t k = 4;                                 /* spline order */
  gsl_vector *c = gsl_vector_alloc(n);                /* control points for spline */
  gsl_bspline_workspace *work = gsl_bspline_alloc_ncontrol(k, n);
  gsl_matrix * XB = gsl_matrix_alloc(n, 3*(k-1) + 1); /* banded collocation matrix */
  gsl_vector_uint * ipiv = gsl_vector_uint_alloc(n);
  double t;
  size_t i;

  /* initialize knots for interpolation */
  gsl_bspline_init_interp(&xv.vector, work);

  /* compute collocation matrix for interpolation */
  gsl_bspline_col_interp(&xv.vector, XB, work);

  /* solve linear system with banded LU */
  gsl_linalg_LU_band_decomp(n, k - 1, k - 1, XB, ipiv);
  gsl_linalg_LU_band_solve(k - 1, k - 1, XB, ipiv, &yv.vector, c);

  /* output the data */
  for (i = 0; i < n; ++i)
    {
      double xi = x_data[i];
      double yi = y_data[i];
      printf("%f %f\n", xi, yi);
    }

  printf("\n\n");

  /* output the spline */
  for (t = x_data[0]; t <= x_data[n-1]; t += 0.005)
    {
      double result;
      gsl_bspline_calc(t, c, &result, work);
      printf("%f %f\n", t, result);
    }

  gsl_vector_free(c);
  gsl_bspline_free(work);
  gsl_vector_uint_free(ipiv);
  gsl_matrix_free(XB);

  return 0;
}

References and Further Reading¶

Further information on the algorithms described in this section can be found in the following references,

C. de Boor, A Practical Guide to Splines (2001), Springer-Verlag, ISBN 0-387-95366-3.
L. L. Schumaker, Spline Functions: Basic Theory, 3rd edition, Cambridge University Press, 2007.
P. Dierckx, Algorithms for smoothing data with periodic and parametric splines, Computer Graphics and Image Processing, 20, 1982.
1. Dierckx, Curve and surface fitting with splines, Oxford University Press, 1995.
M. S. Mummy, Hermite interpolation with B-splines, Computer Aided Geometric Design, 6, 177–179, 1989.

Further information of Greville abscissae and B-spline collocation can be found in the following paper,

Richard W. Johnson, Higher order B-spline collocation at the Greville abscissae. Applied Numerical Mathematics. vol.: 52, 2005, 63–75.

A large collection of B-spline routines is available in the PPPACK library available at http://www.netlib.org/pppack, which is also part of SLATEC.