This is machine translation

Translated by Microsoft
Mouse over text to see original. Click the button below to return to the English verison of the page.


Smoothing spline


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,...)


sp = spaps(x,y,tol) returns the B-form of the smoothest function f 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 f from the given data is measured by


with the default choice for the weights w making E(f) the composite trapezoidal rule approximation to min(x)max(x)|yf|2, and |z|2 denoting the sum of squares of the entries of z.

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


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

When tol is nonnegative, then the spline f is determined as the unique minimizer of the expression ρE(f) + F(Dmf), with the smoothing parameter ρ (optionally returned) so chosen that E(f) equals tol. Hence, when m is 2, then, after conversion to ppform, the result should be the same (up to roundoff) as obtained by csaps(x,y,ρ/(ρ + 1)). Further, when tol is zero, then the "natural" or variational spline interpolant of order 2m 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.


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); rng(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.

More About

collapse all


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).


[1] C. Reinsch, "Smoothing by spline functions", Numer. Math. 10 (1967), 177–183.

See Also

| | |

Was this topic helpful?