@cindex basis splines, B-splines
@cindex splines, basis

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 datapoint.  @xref{Interpolation}, for information
about interpolating splines.

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

@menu
* Overview of B-splines::
* Initializing the B-splines solver::
* Constructing the knots vector::
* Evaluation of B-spline basis functions::
* Evaluation of B-spline basis function derivatives::
* Working with the Greville abscissae::
* Example programs for B-splines::
* B-Spline References and Further Reading::
@end menu

@node Overview of B-splines
@section Overview
@cindex basis splines, overview

B-splines are commonly used as basis functions to fit smoothing
curves to large data sets. To do this, the abscissa axis is
broken up into some number of intervals, where the endpoints
of each interval are called @dfn{breakpoints}. These breakpoints
are then converted to @dfn{knots} by imposing various continuity
and smoothness conditions at each interface. Given a nondecreasing
knot vector
@c{$t = \{t_0, t_1, \dots, t_{n+k-1}\}$}
@math{t = @{t_0, t_1, @dots{}, t_@{n+k-1@}@}},
the @math{n} basis splines of order @math{k} are defined by
@tex
\beforedisplay
$$
B_{i,1}(x) = \left\{\matrix{1, & t_i \le x < t_{i+1}\cr
                            0, & else}\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)
$$
\afterdisplay
@end tex
@ifinfo

@example
B_(i,1)(x) = (1, t_i <= x < t_(i+1)
             (0, else
B_(i,k)(x) = [(x - t_i)/(t_(i+k-1) - t_i)] B_(i,k-1)(x)
              + [(t_(i+k) - x)/(t_(i+k) - t_(i+1))] B_(i+1,k-1)(x)
@end example

@end ifinfo
@noindent
for @math{i = 0, @dots{}, n-1}. The common case of cubic B-splines
is given by @math{k = 4}. The above recurrence relation can be
evaluated in a numerically stable way by the de Boor algorithm.

If we define appropriate knots on an interval @math{[a,b]} then
the B-spline basis functions form a complete set on that interval.
Therefore we can expand a smoothing function as
@tex
\beforedisplay
$$
f(x) = \sum_{i=0}^{n-1} c_i B_{i,k}(x)
$$
\afterdisplay
@end tex
@ifinfo

@example
f(x) = \sum_i c_i B_(i,k)(x)
@end example

@end ifinfo
@noindent
given enough @math{(x_j, f(x_j))} data pairs. The coefficients
@math{c_i} can be readily obtained from a least-squares fit.

@node Initializing the B-splines solver
@section Initializing the B-splines solver
@cindex basis splines, initializing

The computation of B-spline functions requires a preallocated
workspace of type @code{gsl_bspline_workspace}.

@deftypefun {gsl_bspline_workspace *} gsl_bspline_alloc (const size_t @var{k}, const size_t @var{nbreak})
@tindex gsl_bspline_workspace
This function allocates a workspace for computing B-splines of order
@var{k}. The number of breakpoints is given by @var{nbreak}. This
leads to @math{n = nbreak + k - 2} basis functions. Cubic B-splines
are specified by @math{k = 4}. The size of the workspace is
@math{O(2k^2 + 5k + nbreak)}.
@end deftypefun

@deftypefun void gsl_bspline_free (gsl_bspline_workspace * @var{w})
This function frees the memory associated with the workspace @var{w}.
@end deftypefun

@node Constructing the knots vector
@section Constructing the knots vector
@cindex knots, basis splines

@deftypefun int gsl_bspline_knots (const gsl_vector * @var{breakpts}, gsl_bspline_workspace * @var{w})
This function computes the knots associated with the given breakpoints
and stores them internally in @code{w->knots}.
@end deftypefun

@deftypefun int gsl_bspline_knots_uniform (const double @var{a}, const double @var{b}, gsl_bspline_workspace * @var{w})
This function assumes uniformly spaced breakpoints on @math{[a,b]}
and constructs the corresponding knot vector using the previously
specified @var{nbreak} parameter. The knots are stored in
@code{w->knots}.
@end deftypefun

@node Evaluation of B-spline basis functions
@section Evaluation of B-splines
@cindex basis splines, evaluation

@deftypefun int gsl_bspline_eval (const double @var{x}, gsl_vector * @var{B}, gsl_bspline_workspace * @var{w})
This function evaluates all B-spline basis functions at the position
@var{x} and stores them in the vector @var{B}, so that the @math{i}-th element
is @math{B_i(x)}. The vector @var{B} must be of length
@math{n = nbreak + k - 2}.  This value may also be obtained by calling
@code{gsl_bspline_ncoeffs}.
Computing all the basis functions at once is more efficient than
computing them individually, due to the nature of the defining
recurrence relation.
@end deftypefun

@deftypefun int gsl_bspline_eval_nonzero (const double @var{x}, gsl_vector * @var{Bk}, size_t * @var{istart}, size_t * @var{iend}, gsl_bspline_workspace * @var{w})
This function evaluates all potentially nonzero B-spline basis
functions at the position @var{x} and stores them in the vector @var{Bk}, so
that the @math{i}-th element is @c{$B_{(istart+i)}(x)$}
@math{B_(istart+i)(x)}.
The last element of @var{Bk} is @c{$B_{iend}(x)$}
@math{B_(iend)(x)}.  The vector @var{Bk} must be
of length @math{k}.  By returning only the nonzero basis functions,
this function
allows quantities involving linear combinations of the @math{B_i(x)}
to be computed without unnecessary terms
(such linear combinations occur, for example,
when evaluating an interpolated function).
@end deftypefun

@deftypefun size_t gsl_bspline_ncoeffs (gsl_bspline_workspace * @var{w})
This function returns the number of B-spline coefficients given by
@math{n = nbreak + k - 2}.
@end deftypefun

@node Evaluation of B-spline basis function derivatives
@section Evaluation of B-spline derivatives
@cindex basis splines, derivatives

@deftypefun int gsl_bspline_deriv_eval (const double @var{x}, const size_t @var{nderiv}, gsl_matrix * @var{dB}, gsl_bspline_workspace * @var{w})
This function evaluates all B-spline basis function derivatives of orders
@math{0} through @math{nderiv} (inclusive) at the position @var{x}
and stores them in the matrix @var{dB}.  The @math{(i,j)}-th element of @var{dB}
is @math{d^jB_i(x)/dx^j}.  The matrix @var{dB} must be
of size @math{n = nbreak + k - 2} by @math{nderiv + 1}.
The value @math{n} may also be obtained
by calling @code{gsl_bspline_ncoeffs}.  Note that function evaluations
are included as the zeroth order derivatives in @var{dB}.
Computing all the basis function derivatives at once is more efficient
than computing them individually, due to the nature of the defining
recurrence relation.
@end deftypefun

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

@node Working with the Greville abscissae
@section Working with the Greville abscissae
@cindex basis splines, Greville abscissae
@cindex basis splines, Marsden-Schoenberg points

The Greville abscissae are defined to be the mean location of @math{k-1}
consecutive knots in the knot vector for each basis spline function of order
@math{k}.  With the first and last knots in the @code{gsl_bspline_workspace}
knot vector excluded, there are @code{gsl_bspline_ncoeffs} 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.

@deftypefun double gsl_bspline_greville_abscissa (size_t @var{i}, gsl_bspline_workspace *@var{w});
Returns the location of the @math{i}-th Greville abscissa for the given
B-spline basis.  For the ill-defined case when @math{k=1}, the implementation
chooses to return breakpoint interval midpoints.
@end deftypefun

@c See https://savannah.gnu.org/bugs/index.php?34361
@c @deftypefun int gsl_bspline_knots_greville (const gsl_vector * @var{abscissae}, gsl_bspline_workspace * @var{w}, double * @var{abserr});
@c Given target Greville abscissae values in @var{abscissae} and a workspace
@c @var{w} where @code{abscissae->size == gsl_bspline_ncoeffs(w)}, this functions
@c computes and stores the knots required for the workspace to best approximate
@c the target abscissae.  The approximation is optimal in that the first and last
@c values in @var{abscissae} are preserved exactly while the 2-norm of the error
@c in any other abscissae is minimized.  If not-@code{NULL}, the sum of the
@c absolute approximation errors over each abscissa is returned in @var{abserr}.
@c
@c The workspace order must satisfy @math{k > 1} and @var{abscissae} should be
@c monotonically increasing.  Beware that when @code{w->nbreak} is small relative
@c to @code{w->k} the best approximation may still be of poor quality for
@c non-uniformly spaced @var{abscissae}.  This function has memory and runtime
@c overhead that scales like a QR-based linear least squares solution on a
@c @code{(abscissae->size - 2)} by @code{(w->nbreak - 2)} problem.
@c @end deftypefun

@node Example programs for B-splines
@section Examples
@cindex basis splines, examples

The following program computes a linear least squares fit to data using
cubic B-spline basis functions with uniform breakpoints. The data is
generated from the curve @math{y(x) = \cos{(x)} \exp{(-x/10)}} on
the interval @math{[0, 15]} with Gaussian noise added.

@example
@verbatiminclude examples/bspline.c
@end example

The output can be plotted with @sc{gnu} @code{graph}.

@example
$ ./a.out > bspline.txt
chisq/dof = 1.118217e+00, Rsq = 0.989771
$ graph -T ps -X x -Y y -x 0 15 -y -1 1.3 < bspline.txt > bspline.ps
@end example

@iftex
@sp 1
@center @image{bspline,3.4in}
@end iftex

@node B-Spline References and Further Reading
@section B-Spline References and Further Reading

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

@itemize @w{}
@item
C. de Boor, @cite{A Practical Guide to Splines} (1978), Springer-Verlag,
ISBN 0-387-90356-9.
@end itemize

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

@itemize @w{}
@item
Richard W. Johnson, Higher order B-spline collocation at the Greville
abscissae.  @cite{Applied Numerical Mathematics}. vol.@: 52, 2005, 63--75.
@end itemize

@noindent
A large collection of B-spline routines is available in the
@sc{pppack} library available at @uref{http://www.netlib.org/pppack},
which is also part of @sc{slatec}.