Computes a least-squares spline approximation.

The type Imsl_d_spline function is imsl_d_spline_least_squares.

A pointer to the structure that represents the spline fit. If a fit cannot be computed, then NULL is returned. To release this space, use imsl_free.

Let’s make the identifications

For convenience, we assume that the sequence x is increasing, although the function does not require this.

By default, k = 4, and the knot sequence we select equally distributes the knots through the distinct xi’s. In particular, the m + k knots will be generated in [x1, xn] with k knots stacked at each of the extreme values. The interior knots will be equally spaced in the interval.

Once knots t and weights w are determined (and assuming that the option IMSL_OPTIMIZE is not chosen), then the function computes the spline least-squares fit to the data by minimizing over the linear coefficients aj

where the Bj, j = 1, …, m are a (B-spline) basis for the spline subspace.

The optional argument IMSL_ORDER allows the user to choose the order of the spline fit. The optional argument IMSL_KNOTS allows user specification of knots. The function imsl_f_spline_least_squares is based on the routine L2APPR by de Boor (1978, p. 255).

If the option IMSL_OPTIMIZE is chosen, then the procedure attempts to find the best placement of knots that will minimize the least-squares error to the given data by a spline of order k with m coefficients. For this problem to make sense, it is necessary that m > k. We then attempt to find the minimum of the functional

The technique employed here uses the fact that for a fixed knot sequence t the minimization in a is a linear least-squares problem that can be easily solved. Thus, we can think of our objective function F as a function of just t by setting

A Gauss-Seidel (cyclic coordinate) method is then used to reduce the value of the new objective function G. In addition to this local method, there is a global heuristic built into the algorithm that will be useful if the data arise from a smooth function. This heuristic is based on the routine NEWNOT of de Boor (1978, pp. 184 and 258−261).

The initial guess, tg, for the knot sequence is either provided by the user or is the default. This guess must be a valid knot sequence for splines of order k with

with tg nondecreasing, and

In regard to execution speed, this function can be several orders of magnitude slower than a simple least-squares fit.

The return value for this function is a pointer of type Imsl_f_spline. The calling program must receive this in a pointer Imsl_f_spline *sp. This structure contains all the information to determine the spline (stored in B-spline form) that is computed by this function. For example, the following code sequence evaluates this spline at x and returns the value in y.

In the figure below, two cubic splines are fit to

Both splines are cubics with the same spline_space_dim = 8. The first spline is computed with the default settings, while the second spline is computed by optimizing the knot locations using the keyword IMSL_OPTIMIZE.

This example fits data generated from a trigonometric polynomial

where ɛ is a random uniform deviate over the range [-1, 1]. The data are obtained by evaluating this function at 90 equally spaced points on the interval [0, 6]. This data is fitted with a cubic spline with 12 degrees of freedom (eight equally spaced interior knots). The error at 10 equally spaced points is printed out.

This example continues with the first example in which we fit data generated from the trigonometric polynomial

where ɛ is random uniform deviate over the range [−1, 1]. The data is obtained by evaluating this function at 90 equally spaced points on the interval [0, 6]. This data was fitted with a cubic spline with 12 degrees of freedom (in this case, the default gives us eight equally spaced interior knots) and the error sum of squares was printed. In this example, the knot locations are optimized and the error sum of squares is printed. Then, the error at 10 equally spaced points is printed.