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.