Spatially adaptive splines for statistical linear inverse problems.
This paper introduces a new nonparametric estimator based on penalized regression splines for linear operator equations when the data are noisy. A local roughness penalty that relies on local support properties of B-splines is introduced in order to deal with spatial heterogeneity of the function to be estimated. This estimator is shown to be consistent under weak conditions on the asymptotic behaviour of the singular values of the linear operator. Furthermore, in the usual non-parametric settings, it is shown to attain optimal rates of convergence. Then its good performances are confirmed by means of a simulation study.
INTRODUCTION
Statistical linear inverse problems consist of indirect noisy observations
of a parameter (a function generally) of interest. Such problems occur
in many areas of science such as genetics with DNA sequences (Mendel-
sohn and Rice 1982), optics and astronomy with image restoration (Craig
and Brown 1986), biology and natural sciences (Tikhonov and Goncharsky
1987). Then, the data are a linear transform of an original signal f corrupted
by noise, so that we have:
where K is some known compact linear operator defined on a separable
Hilbert space H (supposed in the following to be L 2 [0; 1]; the space of
square integrable functions defined on [0,1]) and ffl i is a white noise with
unknown variance oe 2 : These problems are also called ill-posed problems
because the operator K is compact and consequently equation (1) can not
be inverted directly since K \Gamma1 is not a bounded operator. The reader is
referred to Tikhonov and Arsenin (1977) for a seminal book on ill-posed
operator equations and O'Sullivan (1986) for a review of the statistical
perspective on ill-posed problems. In the following, we will restrict ourself
to integral equations with kernel k(s;
which include deconvolution
There is a vast literature in numerical analysis (Hansen 1998, Neumaier
1998 and references therein) and in statistics dealing with inverse problems
(e.g Wahba 1977, Mendelsohn and Rice 1982, Nychka and Cox 1989,
Abramovich and Silverman 1998 among others). Actually, since model (1)
can not be inverted directly, even if the data are not corrupted by noise,
one has to regularize the estimator by adding a constraint in the estimation
procedure (Tikhonov and Arsenin 1977). The regularization can be linear
and is generally based on a windowed singular value decomposition (SVD)
of K: Indeed, since K is compact, it admits the following decomposition
are orthonormal bases of H
and the singular values are sorted by decreasing order,  1
Most estimators of f proposed in the literature are based either on a finite
rank, depending on the sample size, approximation of K achieved by truncating
the basis expansions obtained by means of the SVD or by adding
a regularization (smoothing) parameter to the eigenvalues that makes K
The sequence of filtering coefficients f j controls the regularity of the so-
for the Tikhonov method and f I [jk] for the truncated
SVD method (see Hansen 1998, for an exhaustive review of these
methods). The rate of decay of the singular values indicates the degree of
ill-posedness of the problem: the more the singular values decrease rapidly
the more ill-posed the problem is.
Nevertheless, estimators based on the SVD have two main defaults. On
the one hand, the basis functions depend explicitly on the operator K and
not on the function of interest. For instance, it is well known that Fourier
basis are the singular functions of K for deconvolution problems and that
they can not provide a parsimonious approximation of the true function
if it is smooth in some regions and rapidly oscillates in other regions. On
the other hand, the usual regularization procedures do not allow to deal
with spatial heterogeneity of the function to be recovered. Several authors
have proposed spatially adaptive estimators based on wavelet decomposition
(Donoho 1995, Abramovich and Silverman 1998) that attain minimax
rates of convergence for particular operators K such as homogeneous oper-
ators. Our approach is quite different and relies on spline fitting with local
roughness penalties.
Until now, spatially adaptive splines were computed by means of knots
selection procedures that require sophisticated algorithms (Friedman 1991,
Stone et al. 1997, Denison et al. 1998). This paper does not address the
topic of knots selection and the estimator proposed below is a penalized
regression splines whose original idea traces back to O'Sullivan (1986) and
Ruppert and Carroll (1999). Actually, Ruppert and Carroll's method con-
sists in penalizing the jumps of the function at the interior knots, each
being controlled by a smoothing parameter, in order to manage both the
highly variable part and the smooth part of the estimator. In this arti-
cle, a similar approach is proposed. Using the fact that B-spline functions
have local supports and that the derivative of a B-spline of order q is the
combination of two B-splines of order are able to define local
measures of the squared norm of a given order derivative of the function
of interest. Thus the curvature of the estimator can be controlled locally
by means of smoothing parameters associated to these local measures of
roughness. Some asymptotic properties of the estimator are given. These
local penalties are controlled by local smoothing parameters whose values
must be chosen very carefully in practical situations in order to get accurate
estimates. The generalized cross validation (GCV) criterion is widely used
for nonparametric regression and generally allows to select "good" values
of the smoothing parameter (Green and Silverman, 1994). Unfortunately,
GCV seems to fail to select effective smoothing parameter values in the
framework of adaptive splines for inverse problems by giving too often un-
dersmoothed estimates. Further investigation is needed to cope with this
important practical topic but that is beyond the scope of this paper. Nevertheless
a small Monte Carlo experiment has been performed to show the
potential of this new approach.
The organization of the paper is as follows. In section 2, the spatial
adaptive regression splines estimator is defined. In section 3, upper bounds
for the rates of convergence are given. The particular case where
(the usual nonparametric framework) is also tackled and the spatially
adaptive estimator is shown to attain optimal rates of convergence. Then,
in section 4, a simulation study compares the behavior of this estimator
to the penalized regression splines proposed by O'Sullivan (1986). Finally,
section 5 gathers the proofs.
Splus programs for carrying out the estimation are available on request.
2. SPATIALLY ADAPTIVE SPLINE ESTIMATES
The estimator proposed below is based on spline functions. Let's now
briefly recall the definition and some known properties of these functions.
Suppose that q and k are integers and let S qk be the space of spline functions
defined on [0; 1]; of order q (q  2); with k equispaced interior knots. The
set S qk is then the set of functions s defined as :
ffl s is a polynomial of degree on each interval
ffl s is continuously differentiable on [0; 1].
The space S qk is known to be of dimension q +k and one can derive a basis
by means of normalized B-splines fB q
1978 or Dierckx, 1993). These functions are non negative and have local
support:
where
Furthermore, a remarkable property of B-splines is that the derivative of a
B-spline of order q can be expressed as a linear combination of two B-splines
of order precisely, if
kj is the jth normalized B-spline of S (q\Gamma1)k and B qk is the vector
of all the B-splines of S qk : Let's define by D qk the weighted differentiation
which gives the coordinates in S (q\Gamma1)k of the
derivative of a function of S qk :
Then, by iterating this process, one can easily obtain the coordinates of a
given order derivative of a function of S qk by applying the (k+q \Gammam)\Theta(k+q)
matrix \Delta (m) defined as follows:
We consider a penalized least squares regression estimator with penalty
proportional to the weighted squared norm of a given order m (m
derivative of the functional coefficient, the effect of which being to give
preference for a certain local degree of smoothness. Using the local support
properties of B-splines, this adaptive roughness penalty is controlled by
local positive smoothing parameters ae that may
take spatial heterogeneity into account. Our penalized B-splines estimate
of f is thus defined as
qk
' is a solution of the following minimization problem
min
'2R q+kn
j is the jth element of ' (m) and k:k denotes the usual L 2 [0; 1] norm.
Let An be the n \Theta (q
C qk the (k +q) \Theta (k +q) matrix whose generic element is the inner product
between two B-splines:
Z 1B q
ki (t)B q
Let's define
I ae C (q\Gammam)k I ae \Delta (m)
where I ae is the diagonal matrix with diagonal elements Then, the solution
" ' of the minimization problem (10) is given by:
n;aen
where Y is the vector of R n with elements
Remark 2. 1. If ae then the estimator defined
by (9) is the same as the estimator proposed by O'Sullivan (1986):
min
'2R q+kn
Furthermore, in the usual nonparametric settings (i.e
kind of penalized regression splines have already been used for different
purposes. Kelly and Rice (1991) have used them to nonparametrically
estimate dose-response curves and Cardot and Diack (1998) have demonstrated
they could attain optimal rates of convergence. Besse et al. (1997)
have performed the principal components analysis of unbalanced longitudinal
data and Cardot (2000) has studied the asymptotic convergence of
the principal components analysis of sampled noisy functional data.
Remark 2. 2. The local penalty defined in (10) may be viewed as a
kind of discrete version of the continuous penalty defined as follows:
the local roughness being continously controlled by function ae(t):
3. ASYMPTOTIC RESULTS
Our convergence results are derived according to the semi-norm induced
by the linear operator K (named K-norm in the
and the empirical norm
Then we have kfk belongs to the null space of K and thus
can not be estimated. This norm allows us to measure the distance of the
estimate from the recoverable part of function f and has been considered
by Wahba (1977). Let's define
It is easily
seen that K(m) is the space of polynomial functions defined on [0; 1] with
degree less than m:
To ensure the existence and the convergence of the estimator we need
the following assumptions on the regularity of f; on the repartition of the
design points, the moments of the noise and on the operator
(H.2) The ffl i 's are independent and distributed as ffl where
Let's denote by Fn the empirical distribution of the design se-
quence, ng ae [0; 1] and suppose it converges to a design
measure F that has a continuous, bounded, and strictly positive density h
on [0; 1]: Furthermore, let's suppose that there exists a sequence fdng of
positive numbers tending to zero such that
sup
(H.4) The kernel k(s; t) belongs to L 2 ([0; 1] \Theta [0; 1]) and, for fixed s,
the function t 7! k(s; t) is a continous function whose derivative belongs to
It exists C ? 0 such that 8g 2 K(m); kKgk  C kgk :
In other words, assumption (H.5) means that the null space of K should
not contain a (non null) polynomial whose degree is less than m: This
condition is rather weak when dealing with deconvolution problems but
excludes some operator equations such as differentiation. By assumption
(H.3), the norm of L 2 ([0; 1]; dF (t)) is equivalent to the L 2 ([0; 1]; dt) norm
with respect to the Lebesgue measure. Assumptions (H.5) and (H.3) ensure
the invertibility of G n;ae and hence the unicity of "
f provided that n is
sufficiently large. Finally assumption (H.4) is a technical assumption that
ensures a certain amount of regularity for operator K but that can be
relaxed for particular operator equations. More precisely, it implies that K
is a Hilbert-Schmidt operator, i.e
is the sequence
of singular values of K:
Let's define
0: We can state
now the two main theorems of this article:
Theorem 3.1. Suppose that n tends to infinity and k=o(n);
then under hypotheses (H.1)-(H.5) we have:
f
K;n
O
The best upper bound is
f
K;n
O
\Gamma2p
It is obtained when
There is
no strong assumption on the decay of
ae as n goes to infinity and actually
ae
can be as small as we want. However, it is well known that in practical situations
a too small value of
ae leads to very bad estimates having undesirable
oscillations. Thus the empirical K-norm should not be considered as an
effective criterion to evaluate the asymptotic performance of this estimator.
Theorem 3.2. Suppose that n tends to infinity and k=o(n);
one has the following upper
f
O
ae 6
Remark 3. 1. Upper bounds for the empirical and the K-norm are
different and surprinsigly, that difference is entirely caused by the bias
term whereas one should expect it would be the result of variance. The
bounds obtained in the K-norm error depend directly on how accurately the
empirical measure Fn of the design points approximates the true measure F:
Furthermore, a larger amount of regularization is needed for the estimator
to be convergent. For instance, if the sequence dn decreases at the usual
rate one choose
ae  n \Gamma(p+m)=(2p+1) as before then the
f is not consistent since ae
ae and then d 2
goes to infinity.
If we choose
ae  ae  n \Gamma(p+m)=(4p+3m) and k  n 1=(4p+3m) , then the
asymptotic error is
f
O
Remark 3. 2. This bound may not be optimal for particular operator
equations since the demonstration relies on general arguments without
assuming any particular decay of the singular values of K (excepted the
implicit conditions imposed by H.4). Thus it must be interpreted as an
upper bound for the rates of convergence: under assumptions (H.4) and
on operator K; the rate of convergence is at least the one given in
(15).
Remark 3. 3. The consistency of the estimator (eq. 12) proposed by
O'Sullivan (1986) is a direct consequence of Theorem 3.2. Upper bounds
for the rates of convergence are those obtained in (15).
3. 4. We have supposed that the interior knots were equispaced
but Theorem 3.2 remains true provided that the distance between two
successive knots satisfies the asymptotic condition :
Remark 3. 5. In the usual nonparametric framework,
the estimator "
qk
' is defined as follows :
min
'2R q+kn
Writing
I ae C (q\Gammam)k I ae \Delta (m) ;
defined as in (11). The demonstration
of the convergence of this estimator is an immediate consequence
of the convergence of (9) since it only remains to study the asymptotic behaviour
of Cn : This has already been done by Agarwall and Studden (1980)
who have shown that if
we get under (H1), (H2) and (H3):
f
O
and the usual optimal rates of convergence are attained if k  n 1=(2p+1)
and
Note that there are no conditions on ae since
Cn is a well conditioned matrix.
4. A SIMULATION STUDY
In this section, a small Monte Carlo experiment has been performed in
order to compare the behaviour of the two estimators defined in section 2.
We have simulated ns = 100 samples, each being composed of
noisy measurements of the convoluted function at equidistant design points
in [0; 1]:
Z 1k
ds
with
\Gamma0:5
\Gamma0:5
\Gamma0:5
\Gamma0:5
x
x
Y
-226f
adapt pens
(b)
FIG. 1. (a) Convoluted function Kf and its noisy observation. (b) True function
f , adaptive spline and O'Sullivan's penalized spline estimates with median fit.
and
The noise ffl has gaussian distribution with standard deviation 0:2 so that
the signal-to-noise-ratio is 1 to 8. The integral equation (18) is practically
approximated by means of a quadrature rule. Function f is drawn in Fig.
1, it is flat in some regions and oscillates in others.
We need to choose the smoothing parameter values to compute the es-
timates. These tuning parameters which control the regularity of the estimators
are numerous: the number of knots, the order q of the splines,
the order m of derivation involved in the roughness penalty and the vector
smoothing parameters. Fortunately, all these
parameters have not the same importance to control the behaviour of the
estimators. Indeed, it appears in the usual nonparametric settings that
the most crucial parameters are the elements of ae which are regularization
parameters. The number of knots and their locations are of minor
importance (Eilers and Marx 1996, Besse et al. 1997), provided they are
numerous enough to capture the variability of the true function f: The
number of derivatives used in (10) controls the roughness penalty is rather
important since two different values of m lead to two different estimators.
It may have mechanical interpretation and its value can be chosen by the
practitioners. Here, it was fixed to and the order of the splines to
as it is the case in a lot of applications in the literature. We consider
a set equispaced knots in [0; 1] to build the estimator and thus we
deal with 44\Theta44 square matrices. Nevertheless, the number of smoothing
parameters remains very large: ae 2 R 42
To face this problem, we used the
method proposed by Ruppert and Carrol (1999) which consists to select
a subset of N k ; N k ! k; smoothing parameters ae
including the "edges" ae
The criterion (GCV,
AIC, .) used to select the smoothing parameter values is then optimized
according to this subset of variables, the values of the other smoothing
parameters being determined by linear interpolation: if
\Gammaae
(j
This subset of smoothing parameters
may be chosen a priori if one has some a priori knowledge of the spatial
variability of the true function but in the following we will consider N
equispaced "quantile" smoothing parameters.
14 CARDOT0.2Pens Adapt
MSE
FIG. 2. Boxplot of the mean square error of estimation for each estimate. The
median error is 0.23 for the penalized spline and 0.08 for the adaptive spline.
We first consider a generalized cross validation criterion in order to
choose the values of ae
because it is computationally fast, widely used as
an automatic procedure and has been proved to be efficient in many statistical
settings (Green and Silverman, 1994). Unfortunately, it seems to fail
here (if there is more than one smoothing parameter) and systematically
gives too small smoothing parameter values that lead to undersmoothed
estimates. Actually, we think it would be better to consider a penalized
version of the GCV that takes into account the number of smoothing parameters
and our future work will go in that direction.
Thus, we have defined the exact empirical risk
in order to evaluate the accuracy of estimate "
Smoothing parameter
values are chosen by minimizing the above risk so that we compare the best
attainable penalized splines defined in (12) and adaptive splines estimators
from samples y
Boxplots of this empirical risk are drawn in Fig. 2 and show that the use
of local penalties may lead to substantial improvements of the estimate:
the median error is 0.08 for the adaptative spline whereas it is 0.23 for
the penalized spline. The penalized spline estimates whose curvature is
only controlled by one parameter can not manage both the flat regions and
the oscillatories regions of function f: That's why undesirable oscillations
of the penalized spline estimate appear in the intervals [0; 0:2] and [0:7; 1]
whereas the use of local smoothing parameters allows to cope efficiently
this problem (see Fig. 1).
5. PROOFS
Let's decompose the mean square error into a squared bias and a variance
term according to the x norm which is successively fK; ng and fKg:
f
x
f
x
f
x
and let's study each term separately. Technical Lemmas are gathered at
the end of the section.
5.1. Bias term
5.1.1. Empirical Bias
Define by
qk
n;aen A 0
thermore, it is easy to show that
' k is the solution of the minimization
problem
min
'2R q+kn
Criterion (20) can be written equivalently with the empirical K-norm :
min
K;n
From Theorem XII.1 in De Boor (1978) and regularity assumption (H.1),
there exists
qk such that
sup
where constant C does not depend on k:
Furthermore we have:
K;n
'Z
'Z
Ck \Gamma2p 1
because lim n (1=n)
ds dF (t) ! +1 by
assumptions (H.3) and (H.4).
On the other hand, since
f k is the solution of (20), we have:
K;n
K;n
From Lemma 5.2 and k' s
Finally, the empirical bias is bounded as follows:
f
K;n
5.1.2. K-norm bias
Let's denote by f
ae the solution of the following optimization
problem:
min
ae is the solution of (23), using the continuity of K; one gets with
Lemma 5.2:
where function s is defined in (21).
Writing now
it remains to study the last term of the right side of this inequality to
complete the proof. Let's denote by K k the (q
elements
I ae C (q\Gammam)k I ae
We have
qk
with
n;ae
n IE(Y): It is easy to check that
classical interpolation
theory that
Pursuing the calculus begun in (25) and appealing to Lemmas 5.1 and
5.3, one gets
ae
n;ae
ae
n;ae
O
O(d 2
ae 4
O
O(d 2
O
ae 6
that completes the proof.
5.2. Variance
5.2.1. Empirical variance
Let's denote by K n;k the (q+k)\Theta(q+k) matrix with elements
It is exactly the matrix n \Gamma1 A 0
Before embarking on the calculus, let's
notice that I q+k
I ae C (q\Gammam)k I ae \Delta (m) is a nonnegative
matrix and hence
tr
n;ae K n;k
Furthermore, the largest eigenvalues of G \Gamma1
n;ae K n;k is less than one and thus
for any (q nonnegative matrix A; one has tr(G \Gamma1
n;ae K n;k
et al. 1998, Lemma 6.5). Thus, under (H.2), the empirical
variance term is bounded as follows :
f
K;n
n;aen
tr
n;ae K n;k
tr
n;ae K n;k
5.2.2. K-norm variance
Using the same decomposition as in (27), one obtains readily:
f
n;aen
tr
On the other hand, one can easily check that
n;ae
n;ae K n;k
O
O
O
by Lemma 5.1.
Using equations (28), (29) and the condition finally gets
f
O
5.3. Technical Lemmas
Lemma 5.1. Under assumptions (H.3) and (H.4) one has
An
Proof. Let g be a function of S qk ; then
By integration by
parts followed by the Holder inequality and invoking assumption (H.3), we
An
K;n
Z 1jKg(t)j jDKg(t)jdt
Z 1@k
@t
(s; ds:
One the other hand, from lemma 5.3 we have kKgk
and, with similar arguments, since by assumption (H.4) DK
is a bounded operator, we also get kD Kgk
Matching previous remarks, one finally obtains the desired result :
An
Let's denote by N It is the null
space of
Lemma 5.2. There are two positive constants c 1 and c 2 such that:
I ae C (q\Gammam)k I ae \Delta (m) u  c 1 k
I ae C (q\Gammam)k I ae \Delta (m) u:
Proof.
The Grammian matrix C (q\Gammam)k is positive and from Agarwall and Studden
(1980) it exists two positive constants c 3 and c 4 such that:
Then, it easy to check that the matrix I ae C (q\Gammam)k I ae is positive, its largest
eigenvalue is proportional to  ae 2 k \Gamma1 and its smallest eigenvalue is proportional
to ae remains to study
to complete the proof.
Let's begin with the first point.
is defined in (7). Furthermore,
writing this weighted difference as follows:
'
and remembering that \Delta (m) is obtained by iterations of the differentiation
process, one gets
that completes the proof of the first point.
Let's suppose now that u 2 N
since Furthermore, the sum of (u
can be expressed in a matrix way:
SPATIALLY ADAPTIVE SPLINES 21
where matrix L is a kind of discretized Laplacian matrix defined as follows
. :::
Its eigenvalues are 2 [1 \Gamma cos
bill, 1969). The null space of L is spanned by the constant vector and the
smallest non null eigenvalue is proportional to  2 (k +q) \Gamma2 . Hence, we have
Since the smallest eigenvalue of I ae C (q\Gammam)k I ae is proportional to ae
gets the desired result for 1: The proof is complete by iterating these
calculus for
Lemma 5.3.
ae
n;ae
Proof.
Let's denote by K   the adjoint operator of K; then, for ' 2 R q+k
kK   Kk
By inequality (31), one gets kC qk
bounded since K is continuous and the first point is now complete.
22 CARDOT
Let's recall that G
I ae C (q\Gammam)k I ae \Delta (m) and decompose any
function
From Lemma 5.2, one gets:
I ae C (q\Gammam)k I ae \Delta (m) ' g2  c ae
Then, under (H.5), one has
min
and thus the smallest eigenvalue of G ae satisfies  min Writing
now G n;ae \Gamma G
we get with Lemma 5.1:
Consequently the smallest eigenvalue of G n;ae satisfies
and then  min (G n;ae )  ae



--R

Wavelet decomposition approaches to statistical inverse problems.
Asymptotic integrated mean square error using least squares and bias minimizing splines.
Simultaneousnonparametricregressions of unbalanced longitudinal data.
Convergence en moyenne quadratique de l'estimateur de la r'egression par splines hybrides.
Nonparametric estimation of smoothed principal components analysis of sampled noisy functions.
Inverse Problems in Astronomy.
A Practical Guide to Splines.
Automatic Bayesian curve fitting.
Curve and Surface Fitting with Splines.
Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition
Flexible Smoothing with B-splines and Penalties (with discussion)
Multivariate adaptive regression splines (with discussion).
Introduction to Matrices With Applications in Statistics.
Nonparametric Regression and Generalized Linear Models.

Monotone smoothing with application to dose-response curves and the assessment of synergism
Deconvolution of Microfluometric Histograms with B-splines
Solving Ill-Conditioned and Singular Linear Systems: a Tutorial on Regularization
Convergence rates for regularized solutions of integral equations from discrete noisy data.
A Statistical Perspective on Ill-Posed Inverse Problems

Polynomial splines and their tensor product in extended linear modeling.
Solutions of Ill-posed problems

Practical Approximate Solutions to Linear Operator Equations when the Data are Noisy.
Local Asymptotics for Regression Splines and Confidence Regions.
--TR
Curve and surface fitting with splines
Simultaneous non-parametric regressions of unbalanced longitudinal data
Solving Ill-Conditioned and Singular Linear Systems
Rank-deficient and discrete ill-posed problems

--CTR
Herv Cardot , Pacal Sarda, Estimation in generalized linear models for functional data via penalized likelihood, Journal of Multivariate Analysis, v.92 n.1, p.24-41, January 2005
deconvolution;convergence;linear inverse problems;regularization;local roughness penalties;spatially adaptive estimators;regression splines;integral equations