Basis Splines¶
This chapter describes functions for the computation of smoothing basis splines (Bsplines). 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 which is subdivided into smaller intervals by defining a nondecreasing sequence of knots. Between any two adjacent knots, a Bspline of order is a polynomial of degree , and so Bsplines defined over the whole interval 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.
Bsplines of order are basis functions for spline functions of the same order, so that all possible spline functions on some interval can be expressed as linear combinations of Bsplines,
Here, the coefficients are known as control points and are often determined from a least squares fit. The are the basis splines of order , and depend on a knot vector , which is a nondecreasing sequence of knots
On each knot interval , the Bsplines are polynomials of degree . The basis splines of order are defined by
for . The common case of cubic Bsplines is given by . The above recurrence relation can be evaluated in a numerically stable way by the de Boor algorithm.
In what follows, we will drop the and simply refer to the spline functions as and for simplicity, with the understanding that each set of spline basis functions depends on the chosen knot vector . We will also occasionally drop the order and use and interchangeably.
Vector Splines¶
The spline above can be generalized to vector control points,
where the may have some length .
Knot vector¶
In practice, one will often identify locations inside the interval 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
The interior knots are in . In order to have a full basis of Bsplines, we require an additional exterior knots, and . These knots must satisfy
but can otherwise be arbitrary.
Allocation of Bsplines¶

type gsl_bspline_workspace¶
The computation of Bspline 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 Bsplines of order
k
. The number of breakpoints is given bynbreak
. Alternatively, the number of control points may be specified inncontrol
. The relationship is . Cubic Bsplines are specified by . The size of the workspace is .

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 using the nondecreasing interior knot vector
tau
of length by choosing the exterior knots to be equal to the first and last interior knot. The total knot vector iswhere is specified by the th element of
tau
.This function does not verify that the
tau
vector is nondecreasing, 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 and sets the exterior knots to the interval endpoints. The total knot vector is
with

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 Bspline on with a period . There are several ways to choose a periodic knot vector. Here, we adopt the approach of Dierckx, 1982, and require
We further define a uniform spacing between each consecutive knot pair. The knot values are given by

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 abscissa values provided in
x
, which must be in ascending order. The knot vector for interpolation must satisfy the SchoenbergWhitney conditions,This function selects the knots according to
which satisfy the SchoenbergWhitney 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 abscissa values provided in
x
, which must be in ascending order. The parameternderiv
specifies the maximum derivative order which will be interpolated.The knot vector has length , where , and is given by,
This function uses the convention that and . 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 usersupplied vector
t
which must have length and be nondecreasing.
Properties of Bsplines¶

size_t gsl_bspline_ncontrol(const gsl_bspline_workspace *w)¶
This function returns the number of Bspline control points given by .

size_t gsl_bspline_order(const gsl_bspline_workspace *w)¶
This function returns the Bspline order .

size_t gsl_bspline_nbreak(const gsl_bspline_workspace *w)¶
This function returns the number of breakpoints (interior knots) defined for the Bspline.
Evaluation of Bsplines¶

int gsl_bspline_calc(const double x, const gsl_vector *c, double *result, gsl_bspline_workspace *w)¶
This function evaluates the Bspline
at the point
x
given the control pointsc
. The result is stored inresult
. Ifx
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 Bspline
at the point
x
given the control vectors . The vector is stored in the th column ofc
. The result is stored inresult
. The parametersc
andresult
must have the same number of rows. Ifx
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 Bspline basis functions at the position
x
and stores them in the vectorB
of length , so that the th element is for . The outputistart
is guaranteed to be in . By returning only the nonzero basis functions, this function allows quantities involving linear combinations of the to be computed without unnecessary terms (such linear combinations occur, for example, when evaluating an interpolated function).
Evaluation of Bspline 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 Bspline derivative
at the point
x
given the control pointsc
. The derivative order is specified bynderiv
. The result is stored inresult
. Ifx
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 Bspline derivative
at the point
x
given the control vectors . The vector is stored in the th column ofc
. The result is stored inresult
. The parametersc
andresult
must have the same number of rows. The derivative order is specified innderiv
. Ifx
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 Bspline basis function derivatives of orders through
nderiv
(inclusive) at the positionx
and stores them in the matrixdB
. The th element ofdB
is . The last row ofdB
contains . The matrixdB
must be of size by at least . The outputistart
is guaranteed to be in . Note that function evaluations are included as the zeroth order derivatives indB
. By returning only the nonzero basis functions, this function allows quantities involving linear combinations of the and their derivatives to be computed without unnecessary terms.
Evaluation of Bspline 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
using the coefficients stored in
c
and the integration limits (a
,b
). The integral value is stored inresult
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
tob
of each basis function ,storing the results in the vector
y
, of lengthncontrol
. The algorithm uses GaussLegendre quadrature on each knot interval.
Least Squares Fitting with Bsplines¶
A common application is to fit a Bspline to a set of data points 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 and the result is a smooth spline which minimizes the cost function,
where is the Bspline defined above, is the number of control points for the spline, , and is the least squares design matrix. The parameters are optional weights which can be assigned to each data point . Because each basis spline has local support for , the least squares design matrix has only nonzero entries per row. Therefore, the normal equations matrix is symmetric and banded, with lower bandwidth . The GSL routines below use a banded Cholesky factorization to solve the normal equations system for the unknown coefficient vector,
The normal equations approach is often avoided in least squares problems, since if is badly conditioned, then 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 Bspline basis is independent of the knot vector and bounded, depending only on . Numerical evidence has shown that for reasonably small 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 Bspline to the input data (
x
,y
) which must be the same length . Optional weights may be given in thewts
vector. The banded normal equations matrix is assembled and factored with a Cholesky factorization, and then solved, storing the coefficients in the vectorc
, which has length . The cost function is output inchisq
. If the normal equations matrix is singular, the error codeGSL_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 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 functionsgsl_bspline_covariance()
andgsl_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 and right hand side vector , which are stored in
XTX
andXTy
on output. The diagonal weight matrix is specified by the input vectorwts
, which may be set toNULL
, in which case . The inputsx
,y
, andwts
represent the data and weights to be fitted in a leastsquares sense, and they all have length . While the full normal equations matrix has dimension by, it has a significant number of zeros, and so the output matrixXTX
is stored in banded symmetric format. Therefore, its dimensions are by. The vectorXTy
has length . The output matrixXTX
can be passed directly into a banded Cholesky factorization routine, such asgsl_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 and right hand side vectors , which are stored in
XTX
andXTY
on output. The diagonal weight matrix is specified by the input vectorwts
, which may be set toNULL
, in which case . The inputY
is a matrix of the right hand side vectors, of size bynrhs. The inputsx
andwts
represent the data and weights to be fitted in a leastsquares sense, and they have length . While the full normal equations matrix has dimension by, it has a significant number of zeros, and so the output matrixXTX
is stored in banded symmetric format. Therefore, its dimensions are by. The matrixXTY
has size bynrhs. The output matrixXTX
can be passed directly into a banded Cholesky factorization routine, such asgsl_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 by evaluating the spline using the coefficients
c
.

int gsl_bspline_covariance(const gsl_matrix *cholesky, gsl_matrix *cov, gsl_bspline_workspace *w)¶
This function computes the by covariance matrix using the banded symmetric Cholesky decomposition stored in
cholesky
, which has dimensions by. The output is stored incov
. The functiongsl_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),
at the point
x
using the covariance matrixcov
calculated previously withgsl_bspline_covariance()
. The derivative order is specified bynderiv
. The error is stored inerr
on output, which is defined as,where is the vector of Bspline values evaluated at , and 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 1norm) of the normal equations matrix , using its Cholesky decomposition, which must be computed previously by calling
gsl_linalg_cholesky_band_decomp()
. The reciprocal condition number estimate, defined as , is stored inrcond
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 is called periodic on if it satisfies the conditions
These conditions can be satisfied with the following constraints on the knot vector and the control points (Dierckx, 1982),
where is the period. The constraint on the control points means that there are only independent control points, and then the remaining are determined. The independent control points can be determined by solving the following least squares problem (Dierckx, 1982):
where and

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 Bspline to the input data (
x
,y
) which must be the same length . Optional weights may be given in thewts
vector. The input vectorx
must be sorted in ascending order, in order to efficiently compute the decomposition of the least squares system. The algorithm computes the 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 Bspline basis. On output, the coefficients are stored in the vectorc
, which has length . The cost function is output inchisq
.

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 factor and vector for a periodic spline fitting problem on the input data (
x
,y
) with optional weightswts
. The input vectorx
must be sorted in ascending order. On output, the by factor is stored inR
, while the first elements of the vector are stored inQTy
. The residual norm is stored inrnorm
.
Gram Matrix¶
Let be the vector whose elements are the Bspline basis functions evaluated at . Then the by outer product matrix of the Bspline basis at is defined as
This can be generalized to higher order derivatives,
Then, the by Gram matrix of the Bspline basis is defined as
where the integral is usually taken over the full support of the Bsplines, i.e. from to . However it can be computed over any interval .
The Bspline outer product and Gram matrices are symmetric and banded with lower bandwidth . These matrices are positive semidefinite, and in the case of , and are positive definite. These matrices have applications in leastsquares 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 at the location
x
, where the derivative order is specified bynderiv
. The matrix is stored inA
in symmetric banded format, and soA
has dimensions by. Since the Bspline 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 semidefinite Gram matrix , where the derivative order is specified by
nderiv
. The matrix is stored inG
in symmetric banded format, and soG
has dimensions by. The integrals are computed over the full support of the Bspline basis.The algorithm uses GaussLegendre 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 semidefinite Gram matrix , where the derivative order is specified by
nderiv
. The matrix is stored inG
in symmetric banded format, and soG
has dimensions by. The integrals are computed over the interval .The algorithm uses GaussLegendre quadrature on each knot interval. See algorithm 5.22 of L. Schumaker, “Spline Functions: Basic Theory”.
Interpolation with Bsplines¶
Given a set of data points, , we can define an interpolating spline with control points which passes through each data point,
This is a square system of equations which can be solved by inverting the by collocation matrix, whose elements are . In order for the collocation matrix to be invertible, the and spline knots must satisfy the SchoenbergWhitney conditions, which are
If these conditions are met, then the collocation matrix is banded,
with both the upper and lower bandwidths equal to . 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 abscissa values provided in
x
. The values inx
must be sorted in ascending order, although this is not checked explicitly by the function. The collocation matrix is stored inXB
on output in banded storage format suitable for use with the banded LU routines. Therefore,XB
must have dimension by.
Hermite Interpolation with Bsplines¶
Given a set of data interpolation sites, , and a set of interpolation values,
then there is a unique interpolating spline of order which will match the given function values (and their derivatives). The spline has control points and is given by,
This spline will satisfy,
Mummy (1989) provides an efficient algorithm to determine the spline coefficients .

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 and first derivatives at the interpolation sites . The inputs
x
,y
,dy
are all of the same length and contain the , , and respectively. On output, the spline coefficients are stored inc
, which must have length .This function is specifically designed for cubic Hermite splines, and so the spline order must be . 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 Bspline Basis¶
A function can be projected onto the subspace spanned by the Bsplines as follows,
If is a piecewise polynomial of degree less than , then will lie in and . The expansion coefficients can be calculated as follows,
where is the Bspline Gram Matrix. The above equation can be solved using a Cholesky factorization of .

int gsl_bspline_proj_rhs(const gsl_function *F, gsl_vector *y, gsl_bspline_workspace *w)¶
This function computes the right hand side vector for projecting a function
F
onto the Bspline basis. The output is stored iny
which must have sizencontrol
. The integrals are computed using GaussLegendre quadrature, which will produce exact results if is a piecewise polynomial of degree less than .
Greville abscissae¶
The Greville abscissae are defined to be the mean location of consecutive knots in the knot vector for each basis spline function of order . With the first and last knots in the knot vector excluded, there are Greville abscissae for any given Bspline basis. These values are often used in Bspline collocation applications and may also be called MarsdenSchoenberg points.

double gsl_bspline_greville_abscissa(const size_t i, gsl_bspline_workspace *w)¶
Returns the location of the th Greville abscissa for the given Bspline basis. For the illdefined case when , 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 . The knot vector is
where the first and last knot are repeated times, where is the spline order. The output is shown in Fig. 40.
This figure makes it clear that each Bspline has local support and spans knots. The Bsplines near the interval endpoints appear to span less knots, but note that the endpoint knots have multiplicities of .
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 Bsplines and their derivatives using breakpoints and uniform knots on the interval . All derivatives up to order are computed. The output is shown in Fig. 41.
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 Bsplines with uniform breakpoints. The data is generated from the curve with Gaussian noise added. Splines are fitted with and 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.
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 uniform breakpoints and computing splines of order .
The data and fitted models are shown in Fig. 43.
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 Bspline to data generated from the Gaussian
on the interval with noise added. Two data gaps are deliberately introduced on and , which makes the problem illposed. Fig. 44 shows that the ordinary least squares solution exhibits large errors in the data gap regions due to the illconditioned 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:
where is the Gram matrix of second derivatives of the Bsplines. Therefore, our regularized least squares problem is
where the regularization parameter represents a tradeoff between minimizing the data misfit and minimizing the spline curvature. The solution of this leastsquares problem is
The example program below solves this problem without regularization and with regularization for comparison. We also plot the confidence intervals for both solutions using the diagonal elements of the covariance matrix .
Fig. 44 shows the result of regularizing the spline over the full interval, using the command
./a.out > data.txt
Fig. 45 shows the result of regularizing the spline over the smaller interval , correcting only the second gap, using the command
./a.out 0 0.6 > data.txt
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 splines to the Runge function,
with 20 uniform breakpoints on the interval . 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, and . Note that
where is the outer product matrix of first derivatives. Therefore, our regularized least squares problem is
where we have chosen to use the same regularization parameter to regularize both ends of the spline. In the example program below, we compute the unregularized solution () and regularized solution (). 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.735857e02, 5.371758e03]
showing that the regularization has forced the first derivative smaller at both endpoints. The results are shown in Fig. 46.
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 Bspline curves to a dataset which is periodic. Both a nonperiodic and periodic spline are fitted to the data for comparison. The data is derived from the function on ,
with Gaussian noise added. The program fits order splines to the dataset.
The program additionally computes the derivatives up to order at the endpoints and :
=== Nonperiodic spline endpoint derivatives ===
deriv 0: [ 8.697939e01, 8.328553e01]
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 , with the th derivative being noncontinuous.
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); /* nonperiodic 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 nonperiodic knots on [a, b] */
gsl_bspline_init_uniform(a, b, w);
/* solve least squares problem for nonperiodic 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, "=== Nonperiodic 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 Bspline basis¶
The following example program projects the function
onto a cubic Bspline basis with uniform knots on the interval . Because the function is cubic, it can be exactly decomposed in the Bspline basis.
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 Bspline 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.
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*(k1) + 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[n1]; 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), SpringerVerlag, ISBN 0387953663.
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.
Dierckx, Curve and surface fitting with splines, Oxford University Press, 1995.
M. S. Mummy, Hermite interpolation with Bsplines, Computer Aided Geometric Design, 6, 177–179, 1989.
Further information of Greville abscissae and Bspline collocation can be found in the following paper,
Richard W. Johnson, Higher order Bspline collocation at the Greville abscissae. Applied Numerical Mathematics. vol.: 52, 2005, 63–75.
A large collection of Bspline routines is available in the PPPACK library available at http://www.netlib.org/pppack, which is also part of SLATEC.