38. Basis Splines

This chapter describes functions for the computation of smoothing basis splines (B-splines). The header file ‘gsl_bspline.h’ contains prototypes for the bspline functions and related declarations.

38.1 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 breakpoints. These breakpoints are then converted to knots by imposing various continuity and smoothness conditions at each interface. Given a nondecreasing knot vector , the basis splines of order are defined by for . The common case of cubic B-splines is given by . 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 then the B-spline basis functions form a complete set on that interval. Therefore we can expand a smoothing function as given enough data pairs. The can be readily obtained from a least-squares fit.

38.2 Initializing the B-splines solver

Function: gsl_bspline_workspace * gsl_bspline_alloc (const size_t k, const size_t nbreak)

This function allocates a workspace for computing B-splines of order k. The number of breakpoints is given by nbreak. This leads to basis functions. Cubic B-splines are specified by . The size of the workspace is .

Function: void gsl_bspline_free (gsl_bspline_workspace * w)

This function frees the memory associated with the workspace w.

38.3 Constructing the knots vector

Function: int gsl_bspline_knots (const gsl_vector * breakpts, gsl_bspline_workspace * w)

This function computes the knots associated with the given breakpoints and stores them internally in w->knots.

Function: int gsl_bspline_knots_uniform (const double a, const double b, gsl_bspline_workspace * w)

This function assumes uniformly spaced breakpoints on and constructs the corresponding knot vector using the previously specified nbreak parameter. The knots are stored in w->knots.

38.4 Evaluation of B-splines

Function: int gsl_bspline_eval (const double x, gsl_vector * B, gsl_bspline_workspace * w)

This function evaluates all B-spline basis functions at the position x and stores them in B, so that the th element of B is . B must be of length . This value is also stored in w->n. It is far more efficient to compute all of the basis functions at once than to compute them individually, due to the nature of the defining recurrence relation.

38.5 Example programs for B-splines

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 on with gaussian noise added.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* number of data points to fit */
#define N        200

/* number of fit coefficients */
#define NCOEFFS  8

/* nbreak = ncoeffs + 2 - k = ncoeffs - 2 since k = 4 */
#define NBREAK   (NCOEFFS - 2)

int
main (void)
{
  const size_t n = N;
  const size_t ncoeffs = NCOEFFS;
  const size_t nbreak = NBREAK;
  size_t i, j;
  gsl_bspline_workspace *bw;
  gsl_vector *B;
  double dy;
  gsl_rng *r;
  gsl_vector *c, *w;
  gsl_vector *x, *y;
  gsl_matrix *X, *cov;
  gsl_multifit_linear_workspace *mw;
  double chisq;

  gsl_rng_env_setup();
  r = gsl_rng_alloc(gsl_rng_default);

  /* allocate a cubic bspline workspace (k = 4) */
  bw = gsl_bspline_alloc(4, nbreak);
  B = gsl_vector_alloc(ncoeffs);

  x = gsl_vector_alloc(n);
  y = gsl_vector_alloc(n);
  X = gsl_matrix_alloc(n, ncoeffs);
  c = gsl_vector_alloc(ncoeffs);
  w = gsl_vector_alloc(n);
  cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
  mw = gsl_multifit_linear_alloc(n, ncoeffs);

  printf("#m=0,S=0\n");
  /* this is the data to be fitted */
  for (i = 0; i < n; ++i)
    {
      double sigma;
      double xi = (15.0/(N-1)) * i;
      double yi = cos(xi) * exp(-0.1 * xi);

      sigma = 0.1;
      dy = gsl_ran_gaussian(r, sigma);
      yi += dy;

      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 [0, 15] */
  gsl_bspline_knots_uniform(0.0, 15.0, bw);

  /* construct the fit matrix X */
  for (i = 0; i < n; ++i)
    {
      double xi = gsl_vector_get(x, i);

      /* compute B_j(xi) for all j */
      gsl_bspline_eval(xi, B, bw);

      /* fill in row i of X */
      for (j = 0; j < ncoeffs; ++j)
        {
          double Bj = gsl_vector_get(B, j);
          gsl_matrix_set(X, i, j, Bj);
        }
    }

  /* do the fit */
  gsl_multifit_wlinear(X, w, y, c, cov, &chisq, mw);

  /* output the smoothed curve */
  {
    double xi, yi, yerr;

    printf("#m=1,S=0\n");
    for (xi = 0.0; xi < 15.0; xi += 0.1)
      {
        gsl_bspline_eval(xi, B, bw);
        gsl_multifit_linear_est(B, c, cov, &yi, &yerr);
        printf("%f %f\n", xi, yi);
      }
  }

  gsl_rng_free(r);
  gsl_bspline_free(bw);
  gsl_vector_free(B);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_matrix_free(X);
  gsl_vector_free(c);
  gsl_vector_free(w);
  gsl_matrix_free(cov);
  gsl_multifit_linear_free(mw);

  return 0;
} /* main() */

The output can be plotted with GNU graph.

$ ./a.out > bspline.dat
$ graph -T ps -X x -Y y -x 0 15 -y -1 1.3 < bspline.dat > bspline.ps

38.6 References and Further Reading

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

• C. de Boor, A Practical Guide to Splines (1978), Springer-Verlag, ISBN 0-387-90356-9.

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

This document was generated by Autobuild on September, 26 2007 using texi2html 1.78.