whittaker2 {albatross} | R Documentation |
Implementation notes for Whittaker smoothing and interpolation of surfaces
Description
Smooth (Eilers 2003) or estimate the baseline (Eilers 2004) of a surface measured on an arbitrary grid by minimising a sum of penalties. Combined difference orders and two different methods of preventing negative values in the output are supported.
This is not a public interface. Subject to change without further notice. Please do not call from outside albatross.
Usage
whittaker2(x, y, z, lambda, d, p, logscale, nonneg)
diffmat(x, y, d)
vandermonde(x0)
Arguments
x |
Grid values along the rows of |
y |
Grid values along the columns of |
z |
Matrix containing the surface values to smooth or |
lambda |
A vector of smoothness penalties, one for every difference order.
Must be of the same length as |
d |
|
p |
If not missing, use the asymmetric penalty method
(Eilers 2004) to estimate the baseline by penalising the
differences with weight |
logscale |
If not Such transformation prevents the resulting values from getting lower
than A typical value would be |
nonneg |
If not |
x0 |
A vector specifying the grid where the function to be differentiated is measured. Must be sorted. |
Details
Finite difference approximation
How to differentiate a function tabulated on a fixed, potentially nonuniform grid before you even know its values? Use its Taylor series.
First derivative is special because it's possible to use central
differences and get a second-order accurate result, even on a
non-uniform grid, by carefully choosing the points where the
derivative is calculated. Let x + \frac{h}{2}
and
x - \frac{h}{2}
be a pair of adjacent points from the
grid. Here's the Taylor series expansion for f
around x
,
with the Lagrange form of the reminder:
f \left(x + \frac{h}{2}\right) =
f(x) + \frac{h}{2} f'(x) + \frac{h^2}{8} f''(x) +
\frac{h^3}{48} f'''(\zeta)
f \left(x - \frac{h}{2}\right) =
f(x) - \frac{h}{2} f'(x) + \frac{h^2}{8} f''(x) -
\frac{h^3}{48} f'''(\eta)
f \left(x + \frac{h}{2}\right) -
f \left(x - \frac{h}{2}\right) = h f'(x) +
\frac{h^3}{48} \left(f'''(\zeta) + f'''(\eta)\right)
f'(x) = \frac{
f\left(x + \frac{h}{2}\right) - f\left(x - \frac{h}{2}\right)
}{h} - \frac{h^2}{48} \left(f'''(\zeta) + f'''(\eta)\right)
|\delta f'(x)| \le \max_{
\xi \in [x - \frac{h}{2}, x + \frac{h}{2}]
} \frac{h^2}{24} f'''(\xi)
Suppose the three grid points
\xi_1 = x_1 - \frac{h_1}{2}
,
\xi_2 = x_1 + \frac{h_1}{2} = x_2 - \frac{h_2}{2}
,
\xi_3 = x_2 + \frac{h_2}{2}
are adjacent on the grid, and we
know the f'
estimates in x_1
and
x_2
:
On the one hand, Taylor series
expansion for f'(x)
around
\frac{x_1 + x_2}{2}
gives:
f'(x_2) = f'\left(\frac{x_1 + x_2}{2}\right)
+ f''\left(\frac{x_1 + x_2}{2}\right)\frac{x_2 - x_1}{2}
+ f'''\left(\frac{x_1 + x_2}{2}\right)\frac{(x_2 - x_1)^2}{8}
+ f''''(\zeta)\frac{(x_2 - x_1)^3}{48}
f'(x_1) = f'\left(\frac{x_1 + x_2}{2}\right)
- f''\left(\frac{x_1 + x_2}{2}\right)\frac{x_2 - x_1}{2}
+ f'''\left(\frac{x_1 + x_2}{2}\right)\frac{(x_2 - x_1)^2}{8}
- f''''(\eta)\frac{(x_2 - x_1)^3}{48}
f''\left(\frac{x_1 + x_2}{2}\right) =
\frac{f'(x_2) - f'(x_1)}{x_2 - x_1} -
\frac{(x_2 - x_1)^2}{48}(f''''(\zeta) + f''''(\eta))
|\delta f''(x)| \le \max_{
\xi \in [x_1, x_2]
} \frac{(x_2 - x_1)^2}{24} f''''(\xi)
On the other hand, if we substitute the estimations of f'(x)
from above, we get:
f''\left(\frac{x_1 + x_2}{2}\right) = \frac{
h_1 f(\xi_3) - h_1 f(\xi_2)
- h_2 f(\xi_2) + h_2 f(\xi_1)
}{h_1 h_2 (x_2 - x_1)} - \frac{
h_2^2 f'''(\zeta_2)
- h_1^2 f'''(\zeta_1)
}{24(x_2 - x_1)}
- \frac{(x_2 - x_1)^2}{24} f''''(\eta)
This is why we can't just keep using central differences and get
n+1
th order accurate results.
What are the general methods of finding the coefficients for the
\mathbf D
matrix? Start with a system of Taylor
series expansions for every grid point:
f(x_i) = \sum_{k=0}^{n-1} f^{(k)}(x) \frac{(x_i - x)^k}{k!}
+ f^{(n)}(\xi) \frac{(x_i - x)^{n}}{n!}
\; \forall i = 1 \dots n
We can solve this system for coefficients
c_i
giving the desired l
-th
derivative estimate with highest accuracy p
possible
(LeVeque 2007):
\sum_{i = 1}^n c_i f(x_i) = f^{(l)}(x) + o(h^p)
Substituting the approximations for
f(x_i)
into the equation, we get the
following condition for the multiplier in front of each
f^{(k)}(x)
:
\frac{1}{k!} \sum_{i = 1}^n c_i (x_i - x)^k = \mathbf{1}_{k = l}
\; \forall k = 0 \dots n-1
In the matrix form, this becomes a Vandermonde system:
V_{k,i} = \frac{(x_i - x)^k}{k!}
b_k = \mathbf{1}_{k = l}
\mathbf c = \mathbf{V}^{-1} \mathbf b
Unfortunately, this system becomes ill-conditioned for “large”
numbers of points. (Experiment shows noticeably growing
c_i
even for third derivative from 10
points and no solution for 32
points on a uniform grid.)
Fornberg (Fornberg 1988) suggests a more numerically stable
procedure, but it still breaks down for 1000
points.
It is important to note that the performance of the method depends
on the matrix \mathbf D
being sparse. While the methods
described above could give more accurate results, they do so at the cost of
providing nonzero weights for a lot of points, and the weights get larger
as the number of points increases. Therefore, with the knowledge that
difference orders above 3
are used very rarely and the interest in
simplicity and performance, we'll minimise the number of coefficients and
their values by solving the Vandermonde system for the minimally accurate
derivative estimations, taking exactly k + 1
points for k
-th
derivative.
What is the error of such estimation? Substituting the Lagrange form
of the remainder into
\mathbf{c}^\top f(\mathbf x)
, we get:
\sum_{i = 1}^n c_i f(x_i) = f^{(n-1)}(x) +
\sum_{i = 1}^n c_i f^{(n)}(\xi_i) \frac{(x_i - x)^n}{n!},
\; \xi_i \in [ x_i, x ]
Our choice of x
(middle point for odd n
, average of
middle points for even n
) lets us shave off one term from the
sum above for odd n
and get second order accurate results for
n = 2
, but other than that, the method is n
-th order
accurate.
Whittaker smoothing
Whittaker smoothing works by minimising a sum of penalties
(Eilers 2003). Interpolation can be achieved by setting weights
\mathbf w
to 0
for the missing points.
\min_{\mathbf{\hat z}} \:
(\mathbf{\hat z} - \mathbf{z})^\top
\mathrm{diag}(\mathbf w)
(\mathbf{\hat z} - \mathbf{z})
+ \lambda | \mathbf D \mathbf{\hat z} |^2
By writing down the derivatives over
\mathbf{\hat z}
and equating them to
0
, we get the normal equation:
(
\mathrm{diag}(\mathbf w) +
\lambda \mathbf{D}^\top \mathbf{D}
) \mathbf{\hat z} = \mathrm{diag}(\mathbf w) \mathbf z
The equation is then solved using solve
.
Given a one-dimensional penalty matrix
\mathbf{D}_d
of order d
obtained by solving a Vandermonde system for every successive group
of d+1
points, we can adapt it for every applicable group of
points from a fluorescence excitation-emission matrix unfolded into
a vector \mathbf z = \mathrm{vec}\, \mathbf F
by
taking Kronecker products with unit matrices:
\mathbf D = \begin{pmatrix}
\mathbf I \otimes \mathbf{D}_\mathrm{em} \\
\mathbf{D}_\mathrm{ex} \otimes \mathbf I
\end{pmatrix}
Penalties of different orders are concatenated by rows in a similar
manner (which is equivalent to adding the corresponding
\mathbf{D}^\top\mathbf D
matrices together in the normal equation).
It has been shown in (Eilers and Goeman 2004) that a combination of
first- and second-order penalty (
2 \lambda \mathbf{D}_1 + \lambda^2 \mathbf{D}_2
)
results in a non-negative impulse response, but the resulting peak
shape may be sub-optimal.
Value
whittaker2 |
A matrix of the same shape as |
diffmat |
A difference matrix |
vandermonde |
A vector |
References
Eilers PHC (2003). “A Perfect Smoother.” Analytical Chemistry, 75(14), 3631-3636. doi:10.1021/ac034173t.
Eilers PHC (2004). “Parametric Time Warping.” Analytical Chemistry, 76(2), 404-411. doi:10.1021/ac034800e.
Eilers PHC, Goeman JJ (2004). “Enhancing scatterplots with smoothed densities.” Bioinformatics, 20(5), 623-628. doi:10.1093/bioinformatics/btg454.
Fornberg B (1988). “Generation of finite difference formulas on arbitrarily spaced grids.” Mathematics of Computation, 51(184), 699-706. doi:10.1090/S0025-5718-1988-0935077-0.
LeVeque RJ (2007). “Finite Difference Approximations.” In Finite Difference Methods for Ordinary and Partial Differential Equations, chapter 1, 3-11. Society for Industrial and Applied Mathematics (SIAM). https://faculty.washington.edu/rjl/fdmbook/.
See Also
Examples
data(feems)
z <- feemscatter(feems$a, rep(25, 4), 'omit')
str(albatross:::whittaker2(
attr(z, 'emission'), attr(z, 'excitation'), z,
c(1, 1e-3), 1:2, logscale = NA, nonneg = 1
))