Usage Notes¶
The majority of the functions in this chapter produce cubic piecewise
polynomial or general spline functions that either interpolate or
approximate given data or support the evaluation and integration of these
functions. Two major subdivisions of functions are provided. The cubic
spline functions begin with the prefix “cubSpline” and use the
piecewise polynomial representation described below. The spline functions
begin with the prefix “spline” and use the B-spline representation
described below. Most of the spline functions are based on routines in the
book by de Boor (1978).
We provide a few general purpose routines for general least-squares fit to data and a routine that produces an interpolant to two-dimensional scattered data.
Piecewise Polynomials¶
A univariate piecewise polynomial (function) p is specified by giving its breakpoint sequence ξ∈ℜ, the order k (degree k − 1) of its polynomial pieces, and the k × (n - 1) matrix c of its local polynomial coefficients. In terms of this information, the piecewise polynomial (ppoly) function is given by
The breakpoint sequence ξ is assumed to be strictly increasing, and we extend the ppoly function to the entire real axis by extrapolation from the first and last intervals. This representation is redundant when the ppoly function is known to be smooth. For example, if p is known to be continuous, then we can compute \(c_{1,i+1}\) from the \(c_{ii}\) as follows:
For smooth ppoly, we prefer to use the nonredundant representation in terms of the “basis” or B-splines, at least when such a function is first to be determined.
Splines and B-Splines¶
B-splines provide a particularly convenient and suitable basis for a given class of smooth ppoly functions. Such a class is specified by giving its breakpoint sequence, its order k, and the required smoothness across each of the interior breakpoints. The corresponding B-spline basis is specified by giving its knot sequence \(t \epsilon \Re^M\). The specification rule is as follows: If the class is to have all derivatives up to and including the j-th derivative continuous across the interior breakpoint \(\xi_i\), then the number \(\xi_i\) should occur \(k − j − 1\) times in the knot sequence. Assuming that \(\xi_1\) and \(\xi_n\) are the endpoints of the interval of interest, choose the first k knots equal to \(\xi_1\) and the last k knots equal to \(ξ_n\). This can be done because the B-splines are defined to be right continuous near \(ξ_1\) and left continuous near \(\xi_n\).
When the above construction is completed, a knot sequence t of length M is generated, and there are \(m := M − k\) B-splines of order k, for example \(B_0, ..., B_{m-1}\), spanning the ppoly functions on the interval with the indicated smoothness. That is, each ppoly function in this class has a unique representation
as a linear combination of B-splines. A B-spline is a particularly compact ppoly function. \(B_i\) is a nonnegative function that is nonzero only on the interval \([t_i, t_{i+k}]\). More precisely, the support of the i-th B-spline is \([t_i, t_{i+k}]\). No ppoly function in the same class (other than the zero function) has smaller support (i.e., vanishes on more intervals) than a B-spline. This makes B-splines particularly attractive basis functions since the influence of any particular B-spline coefficient extends only over a few intervals. When it is necessary to emphasize the dependence of the B-spline on its parameters, we will use the notation \(B_{i,k,t}\) to denote the i-th B-spline of order k for the knot sequence t.
Cubic Splines¶
Cubic splines are smooth (i.e., \(C^0\), \(C^1\) or \(C^2\)), fourth-order ppoly functions. For historical and other reasons, cubic splines are the most heavily used ppoly functions. Therefore, we provide special functions for their construction and evaluation. These routines use the ppoly representation as described above for general ppoly functions (with \(k = 4\)).
We provide three cubic spline interpolation functions:
cubSplineInterpECnd,
cubSplineInterpShape, and
cubSplineTcb. The function
cubSplineInterpECnd allows the user to specify various endpoint
conditions (such as the value of the first or second derivative at the right
and left points). The natural cubic spline, for example, can be obtained
using this function by setting the second derivative to zero at both
endpoints. The function cubSplineInterpShape is designed so that the
shape of the curve matches the shape of the data. In particular, one option
of this function preserves the convexity of the data while the default
attempts to minimize oscillations. The function cubSplineTcb allows the
user to specify tension, continuity and bias parameters at each data point.
It is possible that the cubic spline interpolation functions will produce unsatisfactory results. For example, the interpolant may not have the shape required by the user, or the data may be noisy and require a least-squares fit. The interpolation function splineInterp is more flexible, as it allows you to choose the knots and order of the spline interpolant. We encourage the user to use this routine and exploit the flexibility provided.
Tensor Product Splines¶
The simplest method of obtaining multivariate interpolation and approximation functions is to take univariate methods and form a multivariate method via tensor products. In the case of two-dimensional spline interpolation, the derivation proceeds as follows. Let \(t_x\) be a knot sequence for splines of order \(k_x\), and \(t_v\) be a knot sequence for splines of order \(k_v\). Let \(N_x\) + \(k_x\) be the length of \(t_x\), and \(N_v\) + \(k_x\) be the length of \(t_v\). Then, the tensor-product spline has the following form.
Given two sets of points
and
for which the corresponding univariate interpolation problem can be solved, the tensor-product interpolation problem finds the coefficients \(c_{nm}\) so that
This problem can be solved efficiently by repeatedly solving univariate interpolation problems as described in de Boor (1978, p. 347). Three-dimensional interpolation can be handled in an analogous manner. This chapter provides functions that compute the two-dimensional, tensor-product spline coefficients given two-dimensional interpolation data (spline2dInterp) and that compute the two-dimensional, tensor-product spline coefficients for a tensor-product, least-squares problem (spline2dLeastSquares). In addition, we provide evaluation, differentiation, and integration functions for the two-dimensional, tensor-product spline functions. The relevant functions are spline2dValue and spline2dIntegral.
Scattered Data Interpolation¶
The PyIMSL Math Library provides one function, scattered2dInterp, that returns values of an interpolant to scattered data in the plane. This function is based on work by Akima (1978), which uses \(C^1\) piecewise quintics on a triangular mesh.
Multi-dimensional Interpolation¶
splineNdInterp computes a piecewise polynomial interpolant, of up to 15-th degree, to a function of up to 7 variables, defined on a multi-dimensional grid.
Least Squares¶
The PyIMSL Math Library includes functions for smoothing noisy data. The function userFcnLeastSquares computes regressions with user-supplied functions. The function splineLeastSquares computes a least-squares fit using splines with fixed knots or variable knots. These functions produce cubic spline, least-squares fit by default. Optional arguments allow the user to choose the order and the knot sequence. PyIMSL Math Library also includes a tensor-product spline regression function (spline2dLeastSquares), mentioned above. The function radialScatteredFit computes an approximation to scattered data in \(\Re^N\) using radial-basis functions.
In addition to the functions listed above, several functions in Chapter 10, “Statistics and Random Number Generation”, provide for polynomial regression and general linear regression.
Smoothing by Cubic Splines¶
One ‘‘smoothing spline’’ function is provided. The default action of cubSplineSmooth estimates a smoothing parameter by cross-validation and then returns the cubic spline that smooths the data. If the user wishes to supply a smoothing parameter, then this function returns the appropriate cubic spline.