MathSource Reviews

Splines

Edited by Matthew M. Thomas

At a symposium on the approximation of functions held at General   Motors in 1964, one Philip J. Davis is said to have remarked: "Spline approximation contains the delicious paradox of Prokofieff's Classical Symphony: it seems as though it might have been written several centuries ago, but of course it could not have been" [Schoenberg 1973, p. v]. As this journal enters its eighth volume, it is well nigh time that the topic of splines received its due in this space. Given the importance of splines to the field of approximation problems (of both data-fitting and operator-equation flavors), this particular review might itself have been written several volumes ago, but of course it could not have covered recent MathSource spline-related items had it appeared in the premiere volume. As it now appears, this review covers items submitted to MathSource as recently as mid-1998. It unfolds in threefold manner: It begins   with a nod to the spline literature, useful to any reader whose appetite for   splines is whetted by the herein; it proceeds with the promised examination of MathSource spline-related items, emphasizing item highlights; it concludes with an assessment of what had been examined, acknowledging the broad scope of the examined items.

Background

Spline function approximation is a topic worthy of entire textbooks, and indeed such textbooks have been written. Here, then, is the threepenny synopsis: A spline function approximation of a complex mathematical function entails selection of piecewise-continuous at-least-once-differentiable functions in adjacent subintervals that span the independent-variable range of the function to be approximated. The cubic polynomial is a popular and useful choice for a spline function, as it offers superior representation of approximated function curvature. Spline functions should be continuous at points where adjacent approximations meet, and one or more of their derivatives should also be equal at these same points. In general spline parlance, "knot" is the term given to those points of adjacent approximation confluence and nth derivative equality, while "node" denotes abscissae at which ordinates of the complex function to be approximated are specified. Cubic polynomial coefficient determination is usually achieved through matrix inversion, in which the matrix elements are functions of subinterval width (knot separation distance), and in which the matrix itself is often tridiagonal. Tridiagonal matrix inversion routines are thus  
commonly found among spline packages.

Use of the term "splines" in a function approximation context began with a suggestion from L. H. Thomas of the Ohio State University to I. J.  
Schoenberg of the Ballistic Research Laboratories -- Aberdeen Proving Ground in Maryland. Used for drawing smooth curves, a spline is a thin elastic bar, held in place at various points by weights ("dogs" or "rats") which  
force the spline to pass through the prescribed points [Schoenberg 1946a, p.  
67]. A later work illustrates a mechanical spline [Schumaker 1981, p. 6], though "duck" supplants "dog" in that description; this later work also notes the importance of the mechanical spline to naval architecture. A stylish essay labeled "Why splines?" [Shikin and Plis 1995, pp. 1-4] muses upon that importance. As Schoenberg's 1946 work notes, spline function approximation found an important early and practical use at Aberdeen, in relating projectile drag coefficient to its velocity. In taking spline  
curves from polynomial to analytic approximations, Schoenberg turns to that seen regularly in this space as of late -- conductive heat transfer theory: Indeed, Schoenberg makes reference to the work of Horatio Scott Carslaw, and his 1946 efforts [Schoenberg 1946a, 1946b] sandwich a note on error function erf(x) by John Conrad Jaeger (who, with Carslaw, would write the authoritative analytical text on conductive heat transfer).

Other important works on spline function approximation laid the groundwork for extensive post-1960 publication in the field. An oft-cited one-paragraph abstract [Curry and Schoenberg 1947] relates spline and Pólya distribution functions; a subsequent paper [Schoenberg and Whitney 1953] furthers that relationship. The fundamental theorem of algebra has also been shown to hold for "monosplines" derived from spline functions [Schoenberg 1958]. In a retrospective monograph [Schoenberg 1973] that he describes   as an homage to Leonhard Euler, Schoenberg 1) credits Euler with contributing to the central parts of spline theory, 2) ties together many of the aforementioned works, and 3) notes that winning World War II emphasized spline performance over spline mathematical underpinnings: "In 1944-45 applicable methods were needed, and there was no time to develop their mathematical properties" [p. v]. Urgent goals tend to have that   effect. Some worthwhile latter-day texts [de Boor 1978, Schumaker 1981, Burden and Faires 1993] either have been or will be cited in this review.

What inherent treatment does Mathematica give to splines? The first Mathematica book reference to "splines" occurs in Section 1.9.4 - Redrawing and Combining Plots. Section 1.9.4 states that the Plot[] function uses actual points and avoids "smoothed or splined curves, which can give misleading results in mathematical graphics". Nevertheless, Mathematica works with splines through its Interpolation[] and InterpolatingFunction[] functions, which operate in conjunction with the NumericalMath`SplineFit` package. Note that the default value of Interpolation[] option InterpolationOrder is 3, denoting a cubic polynomial as the default approximating function for Interpolation[]. What Mathematica itself omits, MathSource helps to provide. The following section examines what is provided, with details on four items in particular.

Splines in MathSource

If we use our browser to gain access to MathSource via the World-Wide Web (URL www.mathsource.com), we can search both new MathSource (3.0) and old MathSource (2.2 and earlier) for documents containing the term "splines". Such a search reveals five items in the new domain, and the same number in the old. We shall briefly examine four of the five items in the new domain; bear in mind, though, that three of these five are also among the five   spline-related items found in the old domain. Interestingly enough, the two old-domain   items absent from the new domain were contributed by Wolfram Research insiders (specifically, John M. Novak and the late Jerry B. Keiper). Novak's "Splines" (/Enhancements/Graphics/2D/0200-619), dated July 1992, comprises a 5-kB package, a 4-kB documentation text file, and a 481-byte sample input text file; the package is used in the George Beck contribution to be examined. Keiper's "Tridiagonal Matrix Routines" (/Enhancements/LinearAlgebra/0200-721) comprises a 2-kB package, a 2-kB documentation text file, and a 324-byte sample input text file. The new-domain entry excluded from our examination but present in the old domain is "Splines and Tridiagonal Matrix Solutions" (/Enhancements/Graphics/2D/0202-307) by Silvio Levy, founding editor of The Mathematica Journal: This December 1989 item comprises 18-kB Splines.m and 1-kB   Tridiagonal.m packages. While these packages get no further mention herein, they are nevertheless very useful and effective. The spline-inclined Mathematica user is well advised to visit the Novak,  Keiper, and Levy items in MathSource first, and then examine what we shall now examine.

George Beck

We begin with an item in which splines are but part of the whole. "High School Notebooks" (/Applications/Education/Secondary/0207-649), by George Beck, comprises fourteen notebooks, only one of which -- the 58-kB   splines.nb file -- involves splines. It is one of the three items also found in the old MathSource domain. Beck writes for high school teachers, and his aim is to set   forth examples of how these teachers might craft notebooks for use with their students. The notebooks tackle topics such as long multiplication and division of positive integers (cleverly coded but better suited, one would hope, for the earlier grades), tweaking slope a and y-intercept b to shift line y = a x + b, similarly tweaking the coefficients in parabola y = a   x^2 + b x + c, scrutinizing Pi, nesting regular polygons and circles, solving quadratic and linear equations, differentiation (also cleverly coded   and very well suited for calculus neophytes), and complex numbers. The notebooks are surely worthy of the attention of Beck's target audience.

Splines, for the most part, are beyond the scope of the high schooler; Beck, accordingly, sidesteps spline history, formal spline usage, and all other formal discussion associated with splines. His splines.nb notebook, entitled "Drawing Letters with Splines", shows how one might use the Novak "Splines" package to create letters of the alphabet, in the style of cursive writing. This notebook contains no original code: Novak's package does  
all of the heavy lifting. Splines.nb loads this package as such:

Needs["Graphics`Spline`"]

The notebook draws the letters V, W, and R. The Plot[] function is used to illustrate mathematical functions which resemble letters V and W; presumably, one is to bear witness to these graphs, obtain knot- and node-coordinates from them, and use Spline[] to connect the resulting dots. No ready-made mathematical function resembles the letter R; one must obtain the following coordinates intuitively:

r = {{0, 0}, {0, 1}, {0, 2}, {1, 1.5}, {0, 1}, {1.5, 0}};

Having obtained these coordinates, let us both plot and connect the  
resulting dots:

Show[Graphics[{PointSize[0.03], Point /@ r}], Graphics[Spline[r, Cubic]]];

This result may not be what 3Com® engineers had   in mind when they created the Graffiti® text input system for their Palm Computing® connected organizers. Nevertheless, the result does strongly resemble the letter R,  and the notebook on the whole provides a non-threatening introduction to   splines that the mathematically inclined high school student might later find useful.

Joseph M. Herrmann

"A Natural Cubic Spline coefficient command and a Clamped Cubic Spline coefficient command for Mathematica and examples of its use" (/Applications/Mathematics/0207-447) is the lengthy title of Joseph M. Herrmann's 1994 contribution to MathSource. This item is without counterpart in the old MathSource domain. It comprises two notebooks: 26-kB cspline.nb, and 29-kB nspline.nb. These notebooks create, respectively, the clamped and natural cubic spline coefficients for a set of 2-D data points, based upon a published algorithm [see Burden and Faires 1993, pp. 133-139]. Given   function f defined on [a, b], and a set of nodes that fall on and between a and b, a cubic spline interpolant S for f is a function that meets a certain set of conditions. Among these are two possible boundary conditions: If the second derivatives of S with respect to the nodes at a and b are both zero, the boundary is a natural boundary; if the first derivative of S with respect to the node at a equals the first derivative of f with respect to that   node (and the same equality also holds at node b), the boundary is a clamped boundary. The natural boundary condition leads to a natural cubic spline interpolant; the clamped boundary condition leads to a clamped cubic spline interpolant. Clamped boundary conditions generally lead to more accurate   approximations of f by S, but they require the values of the first derivatives (or an accurate approximation thereof) at the endpoints as inputs.

The structure of cspline.nb is virtually parallel to that of nspline.nb. Both notebooks make reference to Burden et al (an earlier version of [Burden and Faires 1993]), warn of an absence of input data error checking, provide commands CSplineCoef[S,db,de,a,b,c,d,n,t] and NSplineCoef[S,a,b,c,d,n,t] along with well-written explanatory text, and apply the commands to the same set of input data. That application creates a set of clamped/natural cubic spline functions, which are first plotted individually,  then plotted as one, and finally plotted atop the original set of input data.  Below, in order, are the set of input data, and the final plots from the clamped and natural interpolants, respectively:

B = {{46, 40}, {49, 50}, {51, 55}, {52, 63}, {54, 72}, {56, 70}, {57, 77}, 
    {58, 73}, {59, 90}, {60, 93}, {61, 96}, {62, 88}, {63, 99}, {64, 110}, 
    {66, 113}, {67, 120}, {68, 127}, {71, 137}, {72, 132}};

The interpolant behavior becomes interesting for x > 68; more about that in the Assessment.

Edward Neuman

"Spline Functions" (/Enhancements/Numerical/0209-562) is the July 1998 contribution to MathSource by Edward Neuman, whose affiliation is with the Southern Illinois University at Carbondale Department of Mathematics. It too is without counterpart in the old MathSource domain; in fact, it is new enough to warrant entry (through 1998, at least) in the What's New directory of MathSource. This item comprises the 13-kB SplineFunctions.m package and the 395-kB SplineFunctionsTest.nb notebook. The package contains the Usage and Warning statements that one would hope to find in a well-written Mathematica package, but few comments within the functional code. The notebook contains an introduction, a short paragraph on B-splines and their applications (and four related examples), a similar paragraph on convex interpolating splines (and a related example), a similar paragraph on quadratic histosplines (and a related example), and a similar paragraph on box splines on s+1-directional meshes (and the by-now-expected related example). The notebook refers interested readers to a textbook [de Boor   1978] for more information on B-splines; that text, in turn, credits the aforementioned abstract [Curry and Schoenberg 1947] with originally introducing B-splines.

As one might expect from a notebook occupying hundreds of kilobytes of memory, SplineFunctionsTest.nb is replete with text, equations, and explanatory graphics. To illustrate but a portion of this notebook, we begin at the beginning. The first example in the notebook generates a family of quadratic or parabolic B-splines. To replicate this family, we first load SplineFunctions.m:

Needs["SplineFunctions`"]

We next define spline degree (k = 2 for the quadratic case) and the knot sequence t:

k = 2; t = {0, 1, 1, 3, 4, 6, 6, 6}; tp = Partition[t, k + 2, 1]; 
  le = Length[t] - k - 1;

We then make use of SplineFunctions.m functions bpol[] and plotpdf[] to generate the resulting parabolic B-splines, in a manner piecewise defined:

pprep = (bpol[1, k, #1, x] & ) /@ tp; (plotpdf[#1, x] & ) /@ pprep;

Those familiar with the literature on splines [de Boor 1978, p. 112, Figure IX.2] will recognize the above sequence of plots; Neuman, to his credit, gives credit where credit is due. Carl de Boor uses this sequence to make an important point about the connection between knot multiplicity and spline smoothness: A preponderance of the former leads to a deficiency in the latter.

Susumu Sakakibara

"Spline Wavelet Analysis" (/Applications/Mathematics/0206-873) is the May 1994 contribution to MathSource by Susumu Sakakibara, whose affiliation was with Iwaki Meisei University at the time of the contribution. As is the Beck contribution, so is the Sakakibara work found in the old MathSource domain. This MathSource item comprises a 10-kB SplineWavelet.m package, a 1.2-MB SplineWaveletManual.nb notebook, and a 6320-number data file entitled exp136b.03 that occupies 74 kB of memory and is used by the notebook. The package contains not only a Usage statement for each newly-defined function and Warning statements for situations gone awry, but also explanatory comments within the functional portion of the code. The package is   based upon algorithms from an introductory wavelet text [Chui 1992], but since two chapters of that text deal explicitly with cardinal splines (polynomial spline functions with equidistant simple knots), it is evident that   said text has its basis within an earlier work on splines [indeed, Chui 1988].

There are nine sections (in both a rhetorical and a Mathematica-notebook sense) within SplineWavelet.m: theory; initialization; the structure of the data; the scaling function; wavelet; interpolation; decomposition; data smoothing; references. The "theory" section defines scaling functions and wavelets, based upon cardinal B-splines as   presented by [Chui 1992]. The "initialization" section loads the SplineWavelet.m package. The "structure of the data" section explains that data is dealt with in the form of nested lists, wherein the first list element is a data sequence index, and the second list element is the data sequence itself. The "scaling function" section provides background material on the scaling function, as well as examples of SplineWavelet.m function use in   illustrating said scaling function. The "wavelet" section does much of the same,   but for the wavelet in lieu of the scaling function. The "interpolation" section determines and uses a scaling function in a data interpolation example,   while the "decomposition" section does the same with wavelets in a data decomposition example involving a noisy data set -- both sections feature comparisons of newly-created function output with original data. Below is an example from the decomposition section, in which the "sum" of the three plots below (twice-smoothed, noisy extract from once-smoothed, noisy extract from original data) yields the original data:

The "data smoothing" section applies the SplineWavelet.m functions illustrated in the previous sections to the exp136b.03 data set, and   not only plots but also uses ListPlay[] to play the results. And, true to its name, the "references" section cites relevant references.

Assessment

Except for the Beck contribution, whose involvement with splines is trivial, one could write a complete review on not only the three other items examined above, but also the three aforementioned spline-related items (from Novak, Keiper, and Levy) that were not examined above. Of those other examined items, that from Joseph M. Herrmann shows the least mathematical rigor, but compensates through offering the highest degree of practicality: By using Herrmann's notebooks in conjunction with the text by Burden and Faires, one comes by a good working knowledge of spline function approximation. One unfortunate aspect of these notebooks is that, given how they are stored in MathSource, the user must evaluate the input presented in them to see and   understand the resulting output. Requiring input evaluation is a fine teaching strategy if the user is reading the notebook with Mathematica. If, however, not Mathematica but MathReader is in use, the user cannot evaluate the input, and is thus unable to see (let alone understand) the resulting output. The MathReader-hidden output is unfortunate, given the care and thoughtfulness behind the Herrmann notebooks. Also: In the examples from those notebooks, it is evident that the clamped interpolant behaves more erratically than its natural counterpart for x > 68 (even though the cited textbook ascribes better fitting to clamped interpolants). Kudos to Herrmann for providing two answers that raise a question in the mind of the curious student.

The package from Edward Neuman suffers only in that it tends to lack (* in-code comments *) to complement its ample Usage and Warning statements ... the main flaw in an otherwise well-presented package. Neuman's notebook serves as a very good tool for demonstrating the capabilities of the   package, with functions from the package applied to a number of example problems, and with ample text to describe each example. References to Neuman's own published work and to work cited elsewhere in this review punctuate that text, leaving little doubt that SplineFunctionsTest.nb is a serious notebook by a serious mathematician working with spline function approximation. The example parabolic B-spline family shown above complements its counterpart in de Boor very well, and that is a favorable outcome, since Neuman's   intent was to reproduce that family and nothing more. de Boor, however,   illustrated that family by superimposing each of the five B-spline functions upon the same graph, with a separate y-axis for each spline but with the same x-axis for all five splines. The result was to prove de Boor's point: Knot multiplicity (two at x = 1; three at x = 6) leads to spline discontinuity. Reproducing de Boor's exact plot format would be a challenge -- perhaps more so for Mathematica than for Neuman.

The package from Susumu Sakakibara is very well documented, with not only Usage and Warning statements but also the (* in-code comments *) that help explain the code. As Neuman's notebook serves to cast Neuman's package in good light, so does Sakakibara's notebook serve for Sakakibara's package. Whereas Hermann's work with splines is oriented toward the practical, and Neuman's work with splines is oriented toward the mathematical, Sakakibara's work is wavelet-oriented, with splines as a framework but with wavelets as the primary focus. Accordingly, one might initially steer clear of   Sakakibara if pure spline function approximation is the aim. Nevertheless, Sakakibara does an extensive job of bringing Charles K. Chui's wavelet-based algorithms to bear within Mathematica, using results both visible and audible to convey output to the user. Neuman offers many more references in SplineFunctionsTest.nb than Sakakibara provides in SplineWaveletManual.nb. To his credit, however, Sakakibara illustrates the results of tackling a "realistic" data set,   containing 6320 numbers in 74 kB of memory, and requiring 10 MB for the Mathematica kernel to do the processing. Coding popular algorithms in Mathematica, using the code to solve real-world problems, and using Mathematica to present the solutions (in both traditional and innovative formats) is a definite formula for creating value. Sakakibara follows this formula; others should do likewise.

References

BURDEN, RICHARD L. and FAIRES, J. DOUGLAS. Numerical Analysis, 5th ed. PWS-Kent Publishing Company, Boston (1993).

CHUI, CHARLES K. Multivariate Splines. Society for Industrial and Applied Mathematics [Conference Board of Mathematical Sciences - National Science Foundation (CBMS - NSF) Regional Conference Series in Applied Mathematics, Volume 54], Philadelphia   (1988).

CHUI, CHARLES K. An Introduction to Wavelets. Academic Press, Boston (1992).

CURRY, H. B. and SCHOENBERG, I. J. "On Spline Distributions and Their Limits: The Pólya Distribution Function", abstract 380t, Bulletin of the American Mathematical Society, 53(11), 1114 (1947). 

DE BOOR, CARL. A Practical Guide to Splines. Springer-Verlag, New York (1978).

SCHOENBERG, I. J. "Contributions to the Problem of Approximation of Equidistant Data by Analytic Functions. Part A - On the Problem of Smoothing or Graduation. A First Class of   Analytic Approximation Formulae", Quarterly of Applied Mathematics, 4(1), 45-99 (1946a).

SCHOENBERG, I. J. "Contributions to the   Problem of Approximation of Equidistant Data by Analytic Functions. Part B - On the Problem of Osculatory Interpolation. A Second Class of Analytic   Approximation Formulae", Quarterly of Applied Mathematics, 4(2), 112-141 (1946b).

SCHOENBERG, I. J. Cardinal Spline Interpolation. Society for Industrial and Applied Mathematics [Conference Board of Mathematical Sciences - National Science Foundation (CBMS - NSF) Regional Conference Series in Applied Mathematics, Volume 12], Philadelphia   (1973).

SCHOENBERG, I. J. and WHITNEY, ANNE. "On Pólya Frequency Functions. III. The Positivity of Translation Determinants with an Application to the Interpolation Problem by Spline Curves", Transactions of the American Mathematical Society, 74(2), 246-259 (1953).

SCHOENBERG, I. J. "Spline Functions, Convex Curves, and Mechanical Quadrature", Bulletin of the American Mathematical Society, 64(6), 352-357 (1958).

SCHUMAKER, LARRY L. Spline Functions: Basic Theory. John Wiley and Sons, New York (1981).

SHIKIN, EUGENE V. and PLIS, ALEXANDER I. Handbook on Splines for the User. CRC Press, Boca Raton, Florida (1995).

About The Author

Matthew M. Thomas completed his doctorate in chemical engineering at Washington University in St. Louis, in December 1995. He and   compatriots Thad Podgajny and Steve Buckner are the co-inventors responsible for U.S. Patent 5,712,613, issued on 27 January 1998; the patent features Mathematica in a product that reduces electromagnetic wave scatter.

thomas@wuche2.wustl.edu