# 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 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\left(f\right)=\sum _{j=1}^{n}w\left(j\right)|\left(y\left(:,j\right)-f\left(x\left(j\right)\right)\right){|}^{2}$

with the default choice for the weights `w` making E(f) the composite trapezoidal rule approximation to ${\int }_{\mathrm{min}\left(x\right)}^{\mathrm{max}\left(x\right)}|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\left({D}^{m}f\right)=\underset{\mathrm{min}\left(x\right)}{\overset{\mathrm{max}\left(x\right)}{\int }}\lambda \left(t\right)|{D}^{m}f\left(t\right)|dt$

where Dmf denotes the `m`th 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`.

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

collapse all

### Algorithms

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

## References

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