Note: This page has been translated by MathWorks. Click here to see

To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

To view all translated materials including this page, select Country from the country navigator on the bottom of this 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

$$E(f)={\displaystyle \sum _{j=1}^{n}w(j)|(y(:,j)-f(x(j))){|}^{2}}$$

with the default choice for the weights `w`

making *E*(*f*)
the composite trapezoidal rule approximation to $$\int}_{\mathrm{min}(x)}^{\mathrm{max}(x)}|y-f{|}^{2$$, and |*z*|^{2} denoting
the sum of squares of the entries of *z*.

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

$$F({D}^{m}f)={\displaystyle \underset{\mathrm{min}(x)}{\overset{\mathrm{max}(x)}{\int}}\lambda (t)|{D}^{m}}f(t)|{}^{2}dt$$

where *D ^{m}f* denotes
the

`m`

th derivative of `m`

is `2`

, the
default value for the roughness measure weight When `tol`

is nonnegative, then the spline *f* is
determined as the unique minimizer of the expression ρ*E*(*f*)
+ *F*(*D ^{m}f*),
with the smoothing
parameter ρ (optionally returned) so chosen that

`tol`

. Hence, when `m`

is `2`

,
then, after conversion to ppform, the result should be the same (up
to round-off) as obtained by csaps(`tol`

is zero, then the “natural”
or variational spline interpolant
of order 2`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.

Reinsch's approach [1] 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.