GCV spline smoothing in Matlab

This GitHub source code repository provides a Matlab interface to Woltring’s classic generalized, cross-validatory (GCV) spline smoothing and differentiation code ¹. First, Woltring’s original Fortran implementation ² has been converted to C using Bell Laboratories’ f2c converter ³. Second, two Matlab MEX wrappers have been implemented making the resulting code available to a wider range of users.

Below you will find a brief description of the theory of spline smoothing and a documentation of the Matlab interface provided here. For further information, please refer to the original publications ^{1 4 5}.

• Theoretical minimum
• Matlab interface
• Matlab example
• Licenses
• References

Theoretical minimum

Let $t_i$ ($i=1,\ldots,n\ge2m$) be a set of strictly increasing, not necessarily equidistant abscissa values with corresponding ordinates $x_i$ and positive weight factors $w_i$. A natural spline $s_p(t)$ is a piecewise polynomial function on knot positions $t_i$ which minimizes the criterion function $C_p$ for suitably selected regularization parameter $p\ge0$

Through the regularization parameter $p$, a trade-off can be effectuated between the smoothness of the spline (its $m$th derivative) and the goodness-of-fit to the given data. For $p=0$, an exactly interpolating spline is obtained. In the limiting case $p\rightarrow\infty$, the spline becomes an $m$th order polynomial which fits the data in a weighted least-squares sense. Woltring ¹ implemented ways to determine the regularization parameter $p$ by generalized cross-validation or from a given predicted mean squared error.

Using the basis of B-splines ⁵, $B_i(t)$, the spline can be written as the linear combination $s_p(t)=\sum_{i=1}^nc_{p,i}B_i(t)$ with spline coefficients $c_{p,i}$. The resulting polynomials have order $\le2m$ between the knots, and order $m$ outside the knot range. The polynomials are continuous at the knots up to and including the $(2m-2)$th derivative with vanishing $m$th and higher derivatives at the terminal knots $t_1$ and $t_n$.

Key characteristics of some frequently used spline representations

Back to top ↑

Matlab interface

For the construction and evaluation of splines, Woltring’s original Fortran code provided two separate functions, gcvspl and splder. Access to these functions is made available here by the use of two Matlab MEX wrappers. In order to use these wrappers in your Matlab installation, you first (and only once) need to compile them. The compilation requires the presence of a C/C++ compiler software on your computer, for example, the GNU Compiler Collection gcc or similar. If you do not know if this is the case, please refer to the Matlab documentation on how to install, setup, and test a proper MEX building environment before continuing (ask someone from your IT-section if necessary).

Once a C compiler is available on your computer, start Matlab, change to the path where the files of this repository are located and run the commands listed below (advanced users can use the Makefile provided, instead):

mex -v gcvsplmex.c gcvspl.c
mex -v spldermex.c splder.c

Matlab commands to compile the two MEX wrappers

The resulting binaries gcvsplmex.mexa64 and spldermex.mex64 (perhaps with different filename extensions on non-Unix platforms) are Matlab executables which can be used as is. That is, you can use the functions gcvsplmex and spldermex in Matlab as direct replacements of their original Fortran implementations. However, this repository also provides (yet another) set of wrappers with additional error handling. The use of these wrappers is recommended. Below you will find short documentation for each of them.

Spline construction with gcvspl

The gcvspl function computes a natural B-spline using the generalized cross-validation and mean-squared prediction error criteria of Craven & Wahba ⁴. The model assumes uncorrelated, additive noise and essentially smooth, underlying functions. The independent coordinates may be spaced non-equidistantly.

% compute spline coefficients
%
% c = gcvspl( x, t, m, v, w = ones )
% 
% INPUT
% x : ordinate values to be smoothed (numeric matrix [K, N])
% t : corresonding abscissa values (numeric row [1, N])
% m : half order of the resulting spline (numeric scalar)
% v : expected variance of the ordinates (numeric scalar)
% w : weight factors of the ordinates (numeric row [1, N])
%
% OUTPUT
% c : spline coefficients (numeric matrix [K, N])
%
% REMARKS
% - input abscissa values t must be sorted
% - input half orders m = 1, 2, 3, 4 correspond to the output of a linear, cubic, quintic, heptic spline (etc.)
% - input of negative variance v (any value) switches to generalized cross-validation
% - input weight factors w default to TODO
%
% REQUIREMENTS
% - the binary gcvsplmex must be accessible from the Matlab search path

Type help gcvspl to see this output (or similar) in your Matlab installation

Spline evaluation with splder

The splder function computes the values of a B-spline at selected evaluation points. Woltring’s underlying code is an adoption of Lyche et al. ⁶.

% compute spline values and derivatives
%
% x = splder( c, t, m, s, n )
%
% INPUT
% c : spline coefficients (numeric matrix [K, N])
% t : original abscissa values (numeric row [1, N])
% m : original spline half order (numeric scalar)
% s : evaluation points (numeric row [1, M])
% n : order of derivative (numeric scalar)
%
% OUTPUT
% x : spline values or derivatives (numeric matrix [K, M])
%
% REMARKS
% - input variables t (abscissa values) must be identical to those used for spline creation
% - input half order m must be identical to that used for spline creation
% - orders of derivate n = 0, 1, 2 correspond to position, velocity, acceleration values (etc.)
%
% REQUIREMENTS
% - the binary spldermex must be accessible from Matlab'c search path

Type help splder to see this output (or similar) in your Matlab installation

Zero crossings with splzer

The splzer function computes the location (and sign) of zero crossings of a B-spline.

% compute spline zeros
%
% [z, s] = splzer( t, c, m, n, xmin = min( t ), xmax = max( t ), ofs = 0 )
%
% INPUT
% t : independent variables (numeric row [1, N])
% c : spline coefficients (numeric row [1, N])
% m : spline half order (numeric scalar)
% n : order of derivative (numeric scalar)
% xmin, xmax : search interval (numeric scalar)
% ofs : spline offset (numeric scalar)
%
% OUTPUT
% z : zero locations (numeric row [1, Z])
% s : zero signs (numeric row [1, Z])
%
% REMARKS
% - array of independent variables t must be sorted
% - spline half orders m = 1, 2, 3, 4 correspond to linear, cubic, quintic, heptic splines (etc.)
% - physically, orders of derivate n = 0, 1, 2 correspond to position, velocity, acceleration values (etc.)
% - this function requires the presence of the Curve Fitting Toolbox
%
% REQUIREMENTS
% - Matlab'c Curve Fitting Toolbox must be installed on your computer

Type help splzer to see this output (or similar) in your Matlab installation

Back to top ↑

Matlab example

As an example of usage, consider the following few lines of code, which compute the (fifth-order) quintic spline values and their first two derivatives (velocity and acceleration) of a given set of data points (x,y):

c = gcvspl( x, y, 3 ); % compute spline coefficients

y0 = splder( x, c, 3, x, 0 ); % compute spline values
y1 = splder( x, c, 3, x, 1 ); % compute first derivatives (velocity)
y2 = splder( x, c, 3, x, 2 ); % compute second derivatives (acceleration)

z = splzer( x, c, 3, 1 ); % compute zero crossings of the first derivative

A simple example on how to use the splines interface

Note that by specification of more (or less) than x evaluation points in the fourth argument to the function splder, re-sampling of the original data is possible.

Back to top ↑

Licenses

For archiving and documentation purposes, Woltring’s original Fortran file (gcvspl.f) is stored within this repository, as well as the converted C files (gcvspl.c, gcvspl.h). Please adhere to their original copyright (unrestricted non-commercial use) ². Other source code files published in the repository are released to the public domain. If you intend to cite this software in an academic context, please make use the following metadata:

• Kuberski, Stephan R. (2023). GCV spline smoothing in Matlab [Software]. github.com/srkuberski/gcvspl

Back to top ↑

References