Choice of wavelet smoothness, primary resolution and threshold in wavelet shrinkage.
This article introduces a fast cross-validation algorithm that performs wavelet shrinkage on data sets of arbitrary size and irregular design and also simultaneously selects good values of the primary resolution and number of vanishing moments.We demonstrate the utility of our method by suggesting alternative estimates of the conditional mean of the well-known Ethanol data set. Our alternative estimates outperform the Kovac-Silverman method with a global variance estimate by 25% because of the careful selection of number of vanishing moments and primary resolution. Our alternative estimates are simpler than, and competitive with, results based on the Kovac-Silverman algorithm equipped with a local variance estimate.We include a detailed simulation study that illustrates how our cross-validation method successfully picks good values of the primary resolution and number of vanishing moments for unknown functions based on Walsh functions (to test the response to changing primary resolution) and piecewise polynomials with zero or one derivative (to test the response to function smoothness).
Introduction
Wavelet shrinkage is a technique for estimating curves in the presence of noise which is
appealing because it is nearly minimax for a wide range of functions, computationally
practical and spatially adaptive (see the seminal work of Donoho et al. (1995)). This paper
assumes a familiarity with wavelets and wavelet shrinkage up to the level of Nason and
Department of Mathematics, University of Bristol, University Walk, BRISTOL, BS8 1TW, UK
Silverman (1994). We also rely heavily on developments in Kovac and Silverman (2000).
For a recent survey of research in the area see Vidakovic (1999) or Abramovich, Bailey and
Sapatinas (2000).
Most early wavelet shrinkage techniques relied on Mallat's (1989) pyramid algorithm
for computing the discrete wavelet transform (DWT) which in its standard form requires
data to be equally spaced and contain 2 J values. For this limited data situation wavelet
shrinkage works by taking the DWT, then thresholding or shrinking the coefficients, and
then taking the inverse transformation. A great deal of research effort has been expended
on methods to choose the threshold value of wavelet shrinkage, again see Vidakovic (1999)
for an excellent overview. However, the threshold, albeit important, is but one parameter
involved in wavelet shrinkage. For successful shrinkage the following criteria need also to
be chosen well:
the primary resolution Thresholding of coefficients is applied to coefficients whose resolution
level is equal to or finer than the primary resolution. The primary resolution
parameter is similar to the usual bandwidth parameter in linear smoothing methods.
For wavelet shrinkage choice of the primary resolution was first investigated by Hall
and Patil (1995). Hall and Nason (1997) suggest that actually choosing the primary
resolution on a continuous scale may be advantageous. However, even if the primary
resolution is to be chosen on a discrete scale, as in standard wavelet shrinkage, it is
critical to use good values (just as the bandwidth is critical in linear smoothing).
the analysing wavelet. Very little detailed attention has been paid to the problem of which
wavelet should be used in wavelet shrinkage. The Daubechies' (1988) series of compactly
supported wavelets provide a family of mother wavelets of varying smooth-
ness, V , where V is the number of vanishing moments and V 2 C V (R) where
0:2 and C  (R) is the space of  times continuously differentiable functions on
R. We only consider the Daubechies' extremal phase wavelets in this article ranging
from the discontinuous Haar wavelet to the smooth
denote the quadrature mirror filter associated with Daubechies' extremal phase
wavelet of order V and N hV be the length of this filter.
The type of wavelet transform is also important. For example, the translation-invariant
transform of Coifman and Donoho (1995) often gives better results than using the
DWT.
how the threshold is applied Early work described in Donoho et al. (1995) considered
two methods for applying the threshold to wavelet coefficients, the "keep-or-kill"
hard thresholding (similar to model selection) and "shrink-or-kill" soft threshold-
ing. Other techniques have been suggested such as the firm thresholding of Gao
and Bruce (1997) and more recently Bayesian wavelet shrinkage which also has a
thresholding interpretation, see Chipman et al. (1997) or Abramovich et al. (1998).
This article proposes a fast-update cross-validation method which aims to find good combinations
of the threshold, number of vanishing moments (wavelet smoothness), V , and
primary resolution parameters. Often there is no unique "optimal" combination as the
parameters interact and sometimes quite different combinations give similar results. The
cross-validation method described below could be extended to incorporate different choices
of threshold application and type of wavelet transform but it is not clear that such choices
could be implemented in a fast-update algorithm. Cross-validation for threshold selection
in wavelet shrinkage was proposed, especially for functions sparsely represented by
wavelets, by Nason (1996).
Generalized cross-validation for wavelet shrinkage was proposed by Jansen, Malfait
and Vial (1997) and other cross-validation techniques were proposed by Wang (1996) and
Weyrich and Warhola (1998) but of course cross-validation as a general technique has been
around for a very long time; see Stone (1974) for further details. Recently, papers such
as Hall and Turlach (1997), Sardy et al. (1999) and Kovac and Silverman (2000) have
adapted the wavelet shrinkage methodology to data sets with arbitrary size and irregular
design. See also Antoniadis, Gr-egoire and Vial (1997) and Antoniadis and Pham (1998)
for work on fast linear wavelet methods for random designs. We should also mention that
recent algorithms based on the lifting transform show great promise for curve and surface
estimation for irregular data (see for example, Daubechies et al. (1999) for an excellent
review).
Our cross-validatory method is a development of Kovac and Silverman (2000) and as
such works with data sets of irregular design and arbitrary size but ours additionally gives
useful information on which wavelet and primary resolution to use as well as choosing the
threshold value.
1.1 Wavelet shrinkage and the Kovac-Silverman algorithm
First we establish some notation and describe the data model specified by Kovac and Silverman
(2000). We suppose that f(x) for x 2 (0; 1) is the function that we wish to estimate
and that we observe n data points g(x i ) according to the model
i=0 is i.i.d. noise with mean zero and variance  2 and fx
i=0 are not
necessarily equally spaced.
Kovac and Silverman (2000) propose choosing a new equally spaced grid t
on (0; 1) where and interpolate the observed data onto the new
grid. They propose choosing t
such that ng. Throughout their article they linearly interpolate the
original data to new values, y k , on the grid by
(2)
where although they admit that higher order interpolants
or other reweighting schemes might also be of some use. Writing the original and interpolated
data as vectors the linear transform
described by (2) can be written in matrix form by
where the interpolation matrix R depends on t and x. Each row of R always contains either
one or two non-zero entries which always sum to one. Interpolation to a grid is a useful
technique but certainly not new see, for example, Jones and Lotwick (1983) or Silverman
(1986). Kovac and Silverman (2000) then apply wavelet shrinkage to the interpolated
data, y, which first involves taking the DWT by
is the N  N orthogonal matrix associated with Daubechies' extremal phase
wavelet with V vanishing moments (in practice Mallat's fast algorithm is used but the
matrix multiplication representation is mathematically convenient).
Given model (1) above and in particular the i.i.d. assumption on the noise the covariance
of the interpolated data is given by
Kovac and Silverman (2000) exploit the fact that for the linear interpolation scheme described
above  Y is actually a band matrix. After applying the DWT to y Kovac and
Silverman (2000) show that the variances of the individual wavelet coefficients can be
computed exactly, up to knowledge of  2 which has to be estimated from the data, using
a fast algorithm of computational order O(b2 J ) where Y is the
bandwidth of the matrix  Y . A useful consequence of Kovac and Silverman's work is that
the variance of all the wavelet coefficients can be computed with no more effort than computing
the wavelet coefficients themselves, which is O(N hV 2 J ) and also fast. Vannucci
and Corradi (2000) also present a fast algorithm to compute the variance-covariance matrix
of the wavelet coefficients and link it to the two-dimensional DWT.
Kovac and Silverman (2000) show how knowledge of the wavelet coefficient variances
permits extension of the universal and SURE thresholds of Donoho and Johnstone (1994;
1995) to their interpolated data situation. We are also interested in choosing the threshold,
but by cross-validation, and also simultaneously selecting good values for the number of
vanishing moments, V , and primary resolution, p. The next section shows that it is possible
to apply full leave-one-out cross-validation to the Kovac-Silverman set-up and still retain a
fast algorithm.
cross-validation
Cross-validation is a well-established technique for assessing model prediction error and,
in our situation, selecting good choices of the threshold, number of vanishing moments,
and primary resolution. In the following sections we describe how to obtain a leave-one-
out estimate of the prediction error. That is, our wavelet shrinkage estimator, ^
f t;V;p (x) with
threshold t, number of vanishing moments, V , and primary resolution p aims to minimize
the mean integrated squared error (MISE)
Z 1n
dx:
With the methodology below we could easily choose another form of loss function. In par-
ticular, with wavelet shrinkage, we might be interested in doing better near known discontinuities
or inhomogeneities for example. However, using MISE for now, we can estimate
M by
t;V;p is an estimate of f constructed from all the data except the ith point. To
find good values of (t; V; p) we
M . However, we do not believe that the first
minimizer we come across is in any sense "optimal". Unlike, say, the cross-validation
score developed in Nason (1996) the score ^
M has multiple minima and as many as possible
should be investigated further. The next few sections describe the construction and efficient
computation of the leave-one-out predictor ^
using all the
original data points apart from The key to the efficiency is that removing the ith
original data point only changes grid points and thus only wavelet coefficients local to x i .
2.1 Leave-one-out and interpolation
Removing the ith original data point only has a very local effect on the interpolated data
points because with linear interpolation only those grid points that lie in an interval with
the ith point as one of the end points are affected. For these points the interpolated point t k
is either to the left or the right of x i . If removed (or
removed) then new grid points to the left of x 1 are updated to take the value g 1 (or those to
the right of xn 2 take the value g n 2 ).
Assume that x i 1  t k  x i . We can compute the value of the updated interpolated
points
y k from the old ones y k by the simple formula
where
For each removed point, L i only has to be computed once and only
local to the ith
original point have to be recomputed. If N is chosen well then only a few y k need to be
updated. We record the indices of those t k whose y k value has been updated and pass them
onto the next stage.
2.2 Updating the wavelet transform
The previous section tells which of the y changed. The Mallat
algorithm is a recursive algorithm which takes the fy k g N 1
k=0 as input and
computes coarser versions (called father wavelet coefficients) and detail coefficients between
successive levels of coarse coefficients. For the father wavelet coefficients the formula
for computing a coarser approximation to the data c j 1 from a finer approximation
c j is given by
The finest scale approximation of c J
initializes the algorithm and then coarser ap-
proximations, c J are generated using successive applications of (3).
The mother wavelet coefficients, d lost when moving from a finer
scale c j to a coarser scale c j 1 and are computed by a similar formula to (3) except that
the smoothing filter h V is replaced by a "detail extraction" filter g V , we refer to Nason and
Silverman (1994) or Vidakovic (1999) for further details.
The key point for efficiency is that changing a single y k only affects those wavelet
coefficients which are derived from the y k and only those coefficients local to t k will be
changed. In summary, changing a y k changes very few wavelet coefficients.
More specifically, if the single father wavelet coefficient c j
k is changed then only the
coarser father wavelet coefficients c j 1
' where
l k N hV
need to be recomputed (here dxe is the smallest integer greater than or equal to x, and bxc
is the largest integer less than or equal to x). Recall that the DWT is recursive starting with
the fy k g N 1
as the input so formula (4) shows which coefficients need to be recomputed
at a coarser resolution level and then supplies the indices of those changed recursively to
the same routine for the next coarsest level.
Further efficiency gains can be achieved by noting the range of the changed c j
coefficients
and only recomputing those coarser c j 1 that are involved. For example, if c j
k has
changed for k min  k  k max then we only need to recompute c j 1
' for
l k min N hV
in other words 1(k max k min +N hV +1) coarser coefficients at resolution level j 1 need to
be computed from k max k min +1 coefficients at level j. Since each coefficient computation
takes O(N hV ) the recursive update of wavelet coefficients is effectively O(N hV J) and
hence extremely fast (i.e. effectively O(1) with respect to N and n). Note that a similar
algorithm can be developed for the DWT inversion.
2.3 Updating the wavelet coefficient variance factors
Given the covariance matrix  Y of y the wavelet coefficients' covariance matrix is given
by W V  Y W T
. Removing the ith point alters the covariance matrix  RR T because
the N  n matrix R changes to a N  (n 1) matrix
R. Let r k denote the kth row of R.
has a non-zero entry at position ig. Then rows not in R i in both R
and
R will be the same except that those in
R will be one entry missing where there was a
zero in R. However, rows k 2 R i will be different in R and
R. Therefore, the difference
R
R T will be zero apart from a cross-shaped region where any
row or column of D i in R i can, in general, be non-zero. However, since both  Y and
are band matrices most of the entries in the cross-shaped region in D i , except those close
to the main diagonal (less than b away from the main diagonal) will be zero.
Summarizing we can compute the covariance matrix
Y for the interpolated data from
using
where D i is almost all zero apart from some of the rows/columns near the main diagonal
and in R i . To compute the wavelet coefficient variance for the new interpolated data
we only need to consider W
since we already know W V  Y W T
from the Kovac-
Silverman algorithm. Computation of W V D i W T
can be easily performed using the updating
wavelet transform described in the previous section since multiplication by W or W T
is simply an application of the DWT. Since most of the entries in rows or columns of D i
are zero the updating algorithm can be executed first with a zero transform, and then with
those non-zero entries in each row/column of D i . Again application of such a transform is
extremely fast: O(N hV J).
2.4 Thresholding, inversion and optimisation
Using the above information about which wavelet coefficients' variance change one can
identify those coefficient positions where the quantity
d jk =^ jk has changed (where
jk are the updated variances of the wavelet coefficients,
jk are the variance
factors as in Kovac and Silverman (2000), ^  2 is some estimate of  2 and ^
d jk are the up-dated
empirical coefficients.) In wavelet shrinkage an estimate of  is typically computed
by using a robust estimator such as the median absolute deviation (MAD) of the wavelet
coefficients at the finest level, divided by 0.6745. The estimate ^
can itself be updated
quickly by keeping track of which coefficients in the finest resolution level change.
For thresholding after point removal we only need note which  jk have changed. If all
of those  jk that changed were previously thresholded, and if after point i removal they
are subsequently all smaller in absolute size than the threshold then the estimate does not
change. This means that inversion does not have to be performed and the prediction error
simply taken from the non-removed point estimate. However, if any  jk have changed their
status since last time the estimate ^
has to be recomputed using the efficient inverse
algorithm described in section 2.2.
For optimisation we have found that a grid search algorithm works extremely effectively
for finding minimal values of ^
M . We also have used a golden section search method but

Table

1: Table showing values of ^
M (1000 to 3 s.f.) with t  3:11 fixed for various
values of the primary resolution, p, and number of vanishing moments, V .
Number of vanishing moments, V
this tends to get stuck in one of the multiple minima. Another strategy that we adopt is
to condition on the universal threshold t
log N and then optimise ^
to
find good values of (V; p). Then using the good values of (V; p) we optimise over the
threshold t. This strategy is effective in practice because the universal threshold makes a
useful starting value for the optimiser as its value is independent of (V; p).
3 Example: the ethanol data
Before we describe a simulation study we present an applied example in detail. The well-studied
ethanol data from Brinkman (1981) has been analysis by Cleveland et al. (1993)
and Hastie (1993) but more importantly for our purposes by Kovac and Silverman (2000).
The data consist of measurements from an experiment where ethanol was burned
in a single cylinder engine. The concentration of the total amount of nitric oxide and
nitrogen dioxide in engine exhaust, normalized by work done by the engine is related to
the "equivalence ratio", a measure of the richness of the air ethanol mixture. Note that the
range of the x-axis or "Equivalence ratio" variable is (0:535; 1:232) so this was linearly
shifted to (0; 1). The ethanol data are plotted in the top left-hand corner of Figure 1.
For the purposes of this example we fixed the threshold value t to be equal to the
universal threshold value of Donoho and Johnstone (1994). The default Kovac-Silverman
method chooses equally-spaced grid points to interpolate the original data. Thus
the universal threshold computed was
2 log N  3:11. With this threshold fixed

Table

shows our computed value of ^
M for all values of the primary resolution ranging
from 0 to 6 and for numbers of vanishing moments (smoothness) ranging from 1 to 10 from
the Daubechies "extremal phase" series. From the table it is clear that the lowest value of
M occurs for
or 4 might also be of interest (one could continue to higher values of ^
Next, we conditioned on the three different pairs of (V; p) and optimized over the
Equivalence Ratio
Equivalence ratio
Equivalence ratio
Equivalence ratio

Figure

1: Top left: ethanol data with NOx emission versus Equivalence ratio.
Clockwise from top right: estimates ^
f t;V;p (x) for (t; V; p) equal to (3:11; 8; 2), (3:13; 10; 5)
and (2:88; 2; 3).
threshold value to minimize ^
M . In the cases (8; 2), (10; 5) and (2; 3) the minimizing
thresholds were 3:11, 3:13 and 2:88 all of which are not actually too far from the universal
threshold. The estimated curves ^
f t;V;p (x) are shown in Figure 1. The top-right plot in

Figure

shows our "best" estimate for the underlying curve.
Kovac and Silverman (2000) use Daubechies' ``extremal phase'' wavelet with
vanishing moments and a primary resolution of 3. Referring to Table 1 again one can
see that in terms of minimizing M this combination is only ranked 31st out of all the
combinations tried. Therefore we obtain an approximate 25% improvement on Kovac
and Silverman by using our best combination. The improvement using our search method
is also demonstrated by comparing the plot for our best estimate and their two best in
their bottom row in Figure 1 of Kovac and Silverman (2000). Both of their plots use a
combination of (5; 3) with the universal threshold. Their bottom-left uses a global
estimate of variance, their bottom-right uses a local estimate of variance by noticing that
the variance of the NOx variable decreases with Equivalence ratio. Their use of the
more complex local variance estimate is motivated by the fact that their estimate with the
global variance contains a small spike at about 0.8 (like the one in our bottom-left plot of

Figure

1). However, after examining our top-right plot of Figure 1 which tracks the double
bump and does not have a spike at 0.8 we claim that actually Kovac and Silverman (2000)
need not go to the trouble of forming a local variance estimate but merely need change the
number of vanishing moments V to 8 and their primary resolution to 2. We repeated the
Kovac-Silverman (2000) analysis with the new number of vanishing moments and primary
resolution and the resulting estimate is significantly better than theirs and looks more like
our top-right plot in Figure 1 (note it is not exactly the same since they use soft thresholding
and choose t to be exactly the universal threshold, whereas we use hard thresholding and
optimize the value of t).
In no way are we trying to denigrate Kovac and Silverman (2000). Indeed, this paper is
based on their extremely useful methodology. However, we have used the above example
to stress that choice of number of vanishing moments and primary resolution is extremely
important, probably as important as choice of threshold but considerably neglected by much
of the literature and available software.
4 Simulation study
We performed several simulations that show that once the number of vanishing moments
and primary resolution are correctly selected that our cross-validation method produces
broadly similar results to GlobalSURE type thresholding (and also universal thresholding,
although the goal of universal thresholding is not MISE minimisation). The GlobalSURE
thresholding method was introduced by Nason (1996) and is a single threshold version
of the level-dependent technique based on Stein's unbiased risk estimation for wavelet
shrinkage introduced by Donoho and Johnstone (1995).
However, as the ethanol example demonstrates in the previous section it is important
to show that methodology can adapt to features such as the smoothness of the underlying
function or the scale of the features (which require particular choices of numbers
of vanishing moments or primary resolution). The first simulation ("adapting to primary
scale") concentrates on choice of primary resolution, the next three simulations ("adapting
to smoothness") concentrate on choice of number of vanishing moments.
The "adapting to smoothness" simulations fix the threshold to be the universal threshold
using sample sizes of and 500 and for each sample size perform 10
simulations to find the best combinations of number of vanishing moments and primary
resolution. In each case we use Gaussian noise with zero mean and the variance,  2 , is
specified in each section below. The "adapting to smoothness" simulations demonstrate
how the primary resolution is most influenced by discontinuities in lowest-order derivatives
and that the number of vanishing moments chosen by the cross-validation algorithm
is influenced by the underlying smoothness of the true function. However, these results are
not hard and fast and occasionally the cross-validation technique gets it wrong.
4.1 Adapting to primary scale
The underlying true function for this simulation is the Walsh function W (p; x) defined on
which is a piecewise constant function taking only the values 0 or 1 starting at 0
and then switching to 1 and then back to 0 and so on. The number of switches in the interval
is parametrised by p and the distance between each switch is 1=p. A convenient way to
think about Walsh functions is as a sine wave that has been "blocked"! The parameter p is
more formally known as the sequency number of the Walsh function and it is akin to the
frequency parameter of a sine wave, but not exactly the same as the Walsh function is not
always periodic on (0; 1). See Stoffer (1991) for further information on Walsh functions
and their applications in statistics.

Table

2 demonstrates that the selected primary resolution increases with the fineness of
the true Walsh function although there appears to be quite some variability in the selected
primary resolution values at 2. However, notice that for (e.g.) 2the width of
the Walsh peaks is 3=64 and so a primary resolution matching this of 4 or 5 (corresponding
to nearest widths of 8/256 or 4/256 of the Haar wavelet at this resolution level) might
have been expected. However, our algorithm chooses primary resolution of 7, and indeed
the other primary resolutions also "over-estimate" in this way. This effect is presumably
because of the addition of noise which causes the procedure to be conservative and use
finer scale wavelets. Conceptually, the best wavelet basis for representing this set of Walsh
functions should be the Haar basis. Table 3 shows that our cross-validation method nearly
always selects Haar to be the best basis in this situation.
The ten simulations in Tables 2 and 3 are based on sample sizes of
variance of  . The deliberately large signal to noise ratio in this simulation
is to verify that in low-noise situations the cross-validation procedure chooses reasonable
values of the primary resolution and number of vanishing moments. Clearly as the signal
to noise ratio decreases our procedure will choose the "best" values far less often. We leave
to further work the interesting questions: for which values of the signal to noise ratio does
it become very hard to select good parameters and whether competitors such as SURE can
do a better job.

Table

2: Best primary resolution for Walsh function W (p; x) at sequency number 2Simulation Number

Table

3: Best number of vanishing moments for Walsh function W (p; x) at sequency number
2Simulation Number
4.2 Adapting to smoothness: no derivatives
The underlying true function for this example is
x for x 2 [0; 1
This function is continuous on [0; 1] but is not differentiable everywhere (the first derivative
has a discontinuity at 1). The variance of the noise for these simulations  . The
best primary resolutions and numbers of vanishing moments for each simulation/sample
size combination are shown in Tables 4 and 5. Over all sample sizes the modal number
of vanishing moments is 3 and the associated approximate wavelet smoothness is
0:6. However, at smaller sample sizes the wavelet with 7 vanishing moments
is selected 4 times out of 10, as many times as the wavelet. The modal primary
resolution appears to be 4 but chosen.

Table

4: Best primary resolution for function with no derivatives at different sample sizes.

Table

5: Best number of vanishing moments for function with no derivatives at different
sample sizes.
Simulation Number

Table

primary resolution for function with one derivative at different sample sizes.
Simulation Number
4.3 Adapting to smoothness: one derivative
The underlying true function for this example is
This function is continuous on [0; 1], its first derivative is continuous on [0; 1] but the first
derivative is not differentiable everywhere (the second derivative has discontinuities at 1and 3). The variance of the added noise for these simulations was  . The best
primary resolutions and numbers of vanishing moments for each simulation/sample size
combination are shown in Tables 6 and 7. Over all sample sizes the modal number of vanishing
moments is 9 and the associated approximate wavelet smoothness is approximately
1:8. However, at smaller sample sizes the wavelet with 10 vanishing
moments is selected 5 times out of 10. It is interesting to note that the smoothness of the
wavelets selected for the example with one derivative is approximately twice that in the example
with no derivatives, which is what one would expect. The modal primary resolution

Table

7: Best number of vanishing moments for function with one derivative at different

Table

8: Best primary resolution for function with mixed derivatives at different sample
sizes.
Simulation Number

Table

Best number of vanishing moments for function with mixed derivatives at different
sample sizes.
Simulation Number
is 2.
4.4 Adapting to smoothness: mixed derivatives
The underlying true function for this example mixes the two functions, f 0 (x) and f 1 (x)
from the previous two examples.
f mixed
This function is continuous on [0; 1] but the first derivative has a discontinuity at 1and
the second derivative has a discontinuities at 5and 7. The variance of the added noise
for these simulations was  . The best primary resolutions and numbers of
vanishing moments for each simulation/sample size combination are shown in Tables 8
and 9. The modal primary resolution is 4 for the smaller sample sizes (agreeing with the
primary resolution for the "no" derivative case in section 4.2) but 5 for the
size. The primary resolution is most greatly influenced by the lowest-order derivative,
as one might expect from the work of Hall and Patil (1995). The number of vanishing
moments at the sample size is 3 for four simulations and around 6/7 for most of
the others. At lower sample sizes the number of vanishing moments is generally larger than
the ones selected for the f 0 (x) example, but not as high as for the f 1 (x) example which
indicates that maybe some sort of compromise is being made.
5 Conclusions and further work
In this article we have introduced a fast cross-validation method that performs wavelet
shrinkage on data sets of irregular design and arbitrary size and also selects good values of
the number of vanishing moments and primary resolution.
Our cross-validation method has been shown to work well on the Ethanol data set
and on simulated data where the scale (primary resolution) and smoothness (vanishing mo-
ments) of the underlying true function can be controlled. Further work could be performed
to investigate the conditions under which our method would break down both in terms of
diminishing signal to noise ratio and in non-Gaussian and correlated noise situations. Our
method could easily be extended to use level-dependent thresholds which would be of use
with correlated data. It would also be interesting to see how well a MISE estimator such as
GlobalSURE would perform in place of the cross-validation estimate.
Herrick (1999) uses cross-validation with the Kovac-Silverman (2000) algorithm in the
two-dimensional case, although not using the fast version described above and as such
the implementation is slow. It remains to be seen whether a fast version is plausible: the
fast point insertion/deletion techniques of Green and Sibson (1978) for the Voronoi dia-
gram/Delaunay triangulation used for data interpolation would certainly be of value.

Acknowledgments

The author would like to particularly thank Arne Kovac for helpful discussions and access
to his code. He would like to thank Bernard Silverman, David Herrick, others in the Bristol
Statistics Group and the audience of the wavelets contributed papers session at the 1999 ISI
Helsinki meeting for helpful conversations and suggestions. The author would also like to
thank the referees and Editors for supplying extremely helpful comments and suggestions.



--R

Wavelet analysis and its statistical applications.

Random design wavelet curve smoothing.
Statistics and Probability Letters

Computational Statistics and Data Analysis

Adaptive Bayesian Wavelet Shrink- age
Models.

Orthonormal bases of compactly supported wavelets.
Wavelets on irregular point sets.
spatial adaptation by wavelet shrinkage.
Adapting to unknown smoothness via wavelet shrinkage.
Wavelet shrinkage: asymptopia?
WaveShrink with firm shrinkage.
Computing Dirichlet tesselations in the plane.
On choosing a non-integer resolution level when using wavelet methods
for mean integrated squared error of nonlinear wavelet-based density estimators

Generalized Additive Models.
Wavelet methods for curve and surface estimation.
Generalised cross validation for wavelet thresholding.
On the errors involved in computing the empirical characteristic function.
Extending the scope of wavelet regression methods by coefficient-dependent thresholding
A theory for multiresolution signal decomposition: the wavelet rep- resentation
Wavelet shrinkage using cross-validation
The discrete wavelet transform in S.
Wavelet shrinkage for unequally spaced data.
Density estimation.



Statistical modeling by Wavelets.
Function estimation via wavelet shrinkage for long-memory data
Wavelet shrinkage and generalized cross validation for image denoising.
--TR

--CTR
Matthew A. Nunes , Marina I. Knight , Guy P. Nason, Adaptive lifting for nonparametric regression, Statistics and Computing, v.16 n.2, p.143-159, June      2006
ethanol data;kovac-silverman algorithm;fast wavelet shrinkage updates;cross-validation