Learn more about Spline Toolbox

spaps - Smoothing spline

Syntax

sp = spaps(x,y,tol) [sp,values] = spaps(x,y,tol) [sp,values,rho] = spaps(x,y,tol) [...] = spaps(x,y,tol,arg1,arg2,...) [...] = spaps({x1,...,xr},y,tol,...)

Description

sp = spaps(x,y,tol) returns the B-form of the smoothest function that lies within the given tolerance tol of the given data points (x(j), y(:,j)), j=1:length(x). The data values y(:,j) may be scalars, vectors, matrices, even ND-arrays. Data points with the same data site are replaced by their weighted average, with its weight the sum of the corresponding weights, and the tolerance tol is reduced accordingly.

[sp,values] = spaps(x,y,tol) also returns the smoothed values, i.e., values is the same as fnval(sp,x).

Here, the distance of the function from the given data is measured by

with the default choice for the weights w making the composite trapezoidal rule approximation to , and denoting the sum of squares of the entries of .

Further, smoothest means that the following roughness measure is minimized:

where denotes the mth derivative of . The default value for m is 2, the default value for the roughness measure weight is the constant 1, and this makes a cubic smoothing spline.

When tol is nonnegative, then the spline is determined as the unique minimizer of the expression , with the smoothing parameter (optionally returned) so chosen that equals tol. Hence, when m is 2, then, after conversion to ppform, the result should be the same (up to roundoff) as obtained by . Further, when tol is zero, then the ``natural" or variational spline interpolant of order is returned. For large enough tol, the least-squares approximation to the data by polynomials of degree <m is returned.

When tol is negative, then ρ is taken to be -tol.

The default value for the weight function in the roughness measure is the constant function 1. But you may choose it to be, more generally, a piecewise constant function, with breaks only at the data sites. Assuming the vector x to be strictly increasing, you specify such a piecewise constant by inputting tol as a vector of the same size as x. In that case, tol(i) is taken as the constant value of on the interval (x(i-1) .. x(i)), i=2:length(x), while tol(1) continues to be used as the specified tolerance.

[sp,values,rho] = spaps(x,y,tol) also returns the actual value of used as the third output argument.

[...] = spaps(x,y,tol,arg1,arg2,...) lets you specify the weight vector w and/or the integer m, by supplying them as an argi. For this, w must be a nonnegative vector of the same size as x; m must be 1 (for a piecewise linear smoothing spline), or 2 (for the default cubic smoothing spline), or 3 (for a quintic smoothing spline).

If the resulting smoothing spline, sp, is to be evaluated outside its basic interval, it should be replaced by fnxtr(sp,m) to ensure that its m-th derivative is zero outside that interval.

[...] = spaps({x1,...,xr},y,tol,...) returns the B-form of an r-variate tensor-product smoothing spline that is roughly within the specified tolerance to the given gridded data. (For scattered data, use tpaps.) Now y is expected to supply the corresponding gridded values, with size(y) equal to [length(x1),...,length(xr)] in case the function is scalar-valued, and equal to [d,length(x1),...,length(xr)] in case the function is d-valued. Further, tol must be a cell array with r entries, with tol{i} the tolerance used during the i-th step when a univariate (but vector-valued) smoothing spline in the i-th variable is being constructed. The optional input for m must be an r-vector (with entries from the set {1,2,3}), and the optional input for w must be a cell array of length r, with w{i} either empty (to indicate that the default choice is wanted) or else a positive vector of the same length as xi.

Examples

The statements

w = ones(size(x)); w([1 end]) = 100;
sp = spaps(x,y, 1.e-2, w, 3);

give a quintic smoothing spline approximation to the given data that close to interpolates the first and last datum, while being within about 1.e-2 of the rest.

x = -2:.2:2; y=-1:.25:1; [xx,yy] = ndgrid(x,y); rand('seed',39);
z = exp(-(xx.^2+yy.^2)) + (rand(size(xx))-.5)/30;
sp = spaps({x,y},z,8/(60^2));  fnplt(sp), axis off

produces the figure below, showing a smooth approximant to noisy data from a smooth bivariate function. Note the use of ndgrid here; use of meshgrid would have led to an error.

Algorithm

Reinsch's approach References is used (including his clever way of choosing the equation for the optimal smoothing parameter in such a way that a good initial guess is available and Newton's method is guaranteed to converge and to converge fast).