discretized laplacian smoothing by fourier iviethods - Semantic Scholar
DISCRETIZED LAPLACIAN SMOOTHING BY FOURIER IVIETHODS
by
Finbarr O'Sullivan
TECHNICAL REPORT No. 162 June 1989
Department of Statistics, GN-22 University of Washington Seattle, Wasbington9819S USA
Abstract An approach to multi-dimensional smoothing is introduced which is based on a penalized with a modified The choice of penalty stanaard multi-dim~nsiQnal.Lapladan SQ11aI'e errQrcharacteristics at-least on the interior of hyper-rectangular~omains, which has wide restorationproblems,computar tions are carried out using fast Fourier transforms. Non-linear smoothing is accomplished by iterative application of the linear smoothing technique. The iterative procedure can be accelerated by means of preconditioned conjugate gradients. The methods are implemented in one- and two-dimensional settings. Adaptive choice of the amount of smoothing is based on approximate cross-validation type scores. Some illustrations are given relating to scatter-plot smoothing, estimation of a logistic regression surface and density estimation. The asymptotic mean square error characteristics of the linear smoother are derived and the
smoothing splines even when these boundary conditions are not satisfied.
Discretized Laplacian Smoothing by Fourier Methods Finbarr O'Sullivan Dept. of Biostatistics and Statistics University of Washington Seattle, WA 98195 March 30, 1989
1
Introduction Laplacian Smoothing and Discretization
1.1
Consider a region of il c Rd. Suppose we have data (Xi, zi, Wi; i = 1,2,... where Xi E il, EZi = U(Xi) and Var(zi) ex: wil. A class of multi-dimensional smoothers are defined as minimizers of the penalized likelihoods of the form
(1) where X > 0 and L is a linear differential operator, see Cox [4]. The parameter ,\ > 0 controls the smoothness of the estimator. Particular choices for L lead to thin-plate and Laplacian smoothing splines which have been studied by Meinguet [9], Wahba and Wendelberger, [15] and others. For example, in l-dimension with Lu = u" we obtain cubic smoothing splines. In d-dimensions we let L = ~ - the Laplacian in Rd d
Lu(x)
= ~u(x) = ~
02 u(x ) OX]
(2)
A discretized version of the multi-dimensional smoother with Laplacian penalty is defined as follows: For some parameter h (which might be a vector) let il h be a finite dimensional grid defined on il. Suppose we have data Zh,Wh specified on ilh where EZh = Uh and Wh gives the weights associated with Zh (Zh and Wh could be obtained by projecting the observed potentially scattered data onto the grid ilh). We wish to estimate Uh. Let ~h be a matrix of order 4im(ilh) x di.m(rt h) which yieldsadiscretizedapproximation to the Laplacian, i.e, Uh. is a discretized versiQn of a function U defined on il then
(3) Let Bh = ~~~h' An s'th order Discreiized Laplacian smoother is defined as the minimizer of the penalized likelihood criterion
(4) Here Wh is regarded as a diagonal matrix. We note that the smooth solves the estimating equation (5) any real For ease of notation we IS
8
but for smootnmz purposes we are only interested in the 8 from
8
> O.
1.2
Outline
dnUellSHJUS greater than one on ( are O'Sullivan [11) and Wahba Wendelberger The goal of paper is to show a wen chosen discretization can simplifications in computations yet without comprimising performance cnarachenstics. The paper only considers the case where the region of interest is a d~dimensional rectangle. This is a situation having immediate application to a range of image processing problems. We begin in section 2 by describing a particular choice for Doh whose spectral decomposition is amenable to Fourier analysis. This greatly simplifies the computations. Extensions to non-constant weights is achieved by an iterative scheme whose convergence properties are analyzed. Conjugate gradient methods may be used to accelerate convergence. Adaptive .choice of the smoothing parameter is done using approximate cross-validation assessments. Non-linear .SJ.fioothingby iterativere-weighting is considered in section 3. Various illustrations .and timing information are presentedin section 4. Computing times are dominated by the fast Fourier transform algorithm. one-dimensional Monte-Carlo experiments shows that there may be be little or no difference between a cross-validated discretized Laplacian smoother and a corresponding cross-validated spline smoother. Some asymptotic mean square error performance characteristics are worked out in section 4. The asymptotic analysis uses only very basic tools. In the case that the target function is locally constant at the boundary of the region of interest the asymptotic characteristics are seen to match those of multidimensional smoothing splines.
2
2
Computational Methods amenaote to
Fourier In case the weiehts the computation trivial. non-constant weizhts an iterative algorithm convergence of is it is techniques may be used for acceleration.
2.1
Discretized Laplacian and Spectral Analysis
The region of interest n is a d-dimensional rectangle which for notational convenience we take to be [O,l]d. A lattice is defined on il as follows. Let h = (hI, h 2 , ••• , h d ) where hj = 1/ nj for j = 1,2, ... , d. ilh consists of the set of d-vectors Xhi = (XiI" •• ,Xij' ••• , Xid) with = .hj (6) for 1 :::; i :::;
»s. and 1:::; j
:::; d. There area total of
(7) grid points in this lattice. Let D j be nj x nj-matrix which yields a discrete approximation to {;z. Thus if v = (v(Xd,V(X2),.' .v(x n j ) ) is a nj-vector which is a discretized version of a function v(x) then D j has the property that
for i
= 1,2, ... , nj.
Given Dil a discrete approximation for 6.h is the N x N-matrix d
6.h
= I: 11 ® 12 ® ... Dj."
® t,
(9)
j=l
where Ij choices for D j are possible corresponding to Periodic, Dirichlet boundary conditions for example, see Buzbee et. al. [3]. A modified choice related to the Neuman discretization is -1 1 Dj=
1
0
1 1 -2
1
1
( 1 -2 1 1 -2 1 1 -1
0 3
np. as 1 -+ co. Proof:
iterative scheme may be exnressed as M u(l+1) = b + Du(l)
where M = (>.B h + fiJI], D stage 1, so from (24)
= -[Wh -
(24)
and b = WhZh. Let £(1)
= Uh.\ -
u(l)
be the error at
(25) for 1 = 1,2, .... Such a method is convergent whenever the largest eigenvalue of NI- 1 D is less than 1. The analysis is straightforward. Clearly y'Bf.y 2: 0 so provided ib :» 0 y'M y
=
>,y'.8hY
>
fiJy'y
+ fiJy'Y (26)
This implies y'M- 2y ::; fiJ-2 y'y for all N-vectors y. From here
II M-IDy 11 2 II y 11 2
y'DM-2Dy y'y
w- 2 y' D 2 y
11
= [1 + ~~II](P+280:)/4Iuhlll then
'2) 1('" ,1..2) L...
"'( ~A~~)2 1 28 = ( L... 1 + ,U hll [1-1- ~2lu»J(280:)/2 (/>11 II I
2
=
14
II
'Yll
(63)
(64)
v
Since b is positive sup;
(bi~:~2) 2
is finite for 0 ::S
0:.
v
b depends on h through
>"hv
but by using the bounds in (60) we can find a constant M independent of h such that supv
(br~:~2) 2 ::S M.
Thus for p ::S t ::S 2 with
0:
= (65)
Note that if D.U is continuous and u satisfies the Neuman boundary conditions then for Ihl sufficiently small with t = 2
(66) 5.2.2
r
Bounding the Expected Squared Error
is orthogonal so Var£h
= (T2I.
Thus the expected squared error is
Ihl[1
(67)
Using (60)
I: Ihl[1 + >"~vlP/2[1 + >..>..~~t2(T2 v
(68) where a = c~ and b cis. The summation may be analyzed by slightly generalizing the mtegral approximation argument in Craven and Wahba Let t Vj >..1/4s(Vj - 1) for Vi 1,2, ... ,nj + 1 and j 1,2, ... ,d. Also >..d/4s. Now for o ::S >.. ::S >"0
=
=
+
=
< u
+
- tvJ
)..~/ ", The summation over integral over the positive orthant
is a Riemann sum approximation to the d-dimensional
1/
(70)
Transforming to polar co-ordinates the integral is bounded by
1
00
Jyf
o
[c + ar 4 ]P/ 2 d-l [ dr 1 + br 4 S ]2 r
(71)
where At[ is a positive constant. (The r d - 1 term comes from the Jacobian of the transformation to polar co-ordinates.) Note the latter integral is convergent provided p/2s < 2 - d/4s. Thus wecan find a constant M (independent of h) for which
E as )..
-+
II U k:A -
Eu k:A
11 2kp-< MN- 1 )..-(p/2s+d/4s)
(72)
O. From here we have the following result:
Theorem 2 There is a constant M E
II Uk
- Uk>.
> 0 such that for p < min(t,2s - d/2)
lI~p::; M·
()..(t-p)/s
II Uk
lI~t +N- 1 )..-
(P/ 2S+ d/ 4S»)
(73)
An optimal upper bound on the asymptotic mean square error is obtained by choosing ).. so that the squared bias and variance are on the same order. If Llu is continuous and u satisfies fhe Neuman boundary conditions then this calculation gives (for t = 2 and p < min(t,2s d/2)) E II Uk - Uk>. lI~p::; M . N- 2 (t - p )/ (4s+ d) (74) which occurs when X = O(N-4s/(2t+d»). We conjecture that this upper bound also holds for t E (0,2) provided p < min(t,28 - d/2). These results parallel those obtained for multidimensional smoothing splines by Cox [4].
6
Acknowledgement
I am gratefull to Professor Gene Golub for a some valuable discussions related to this work.
References Akaike H. Information theory an extension of maximum likelihood principle. In Petrov B.N. and Caski F _, editors, Second Iniernation Symposium on Information Theory, pages 267-281, Budapest, 1973. Akademiai Kiado. [2J Buzbee B.L., Door F.vV., George J .A. and Golub G.H. The direct solutions of the discrete poisson equation on irregular regions. SL4iV! J. Numer. Anal., 8:722-736, 1971. [3J Buzbee B.L., Golub G.H. and Nielson C.W. On direct methods for solving poissons equation. SIAlY! J. Numer. Ana!., 7:627-656, 1970. [4J Cox D.D. Multivariate smoothing spline functions. SIAM J. Numer. Anal., 21:789-813, 1984.
[5) Craven P. and Wahba G. Smoothing noisy data with spline functions: estimating the correct degree of smoothing by the method of generalized cross-validation. Numer. Math., 31:377-403, 1979. [6) Golub G.H. and Van Loan C.F. Matrix Computations. North Oxford Academic, New York, 1983. [7] Huber P.J. Robust smoothing. In Launer R. L. and Wilkinson G.H., editors, Robustness in Statistics. Academic Press, New York, 1979. [8) Mardia K.V., Kent J.T. and Bibby J.M. London, 1979.
Multivariate Analysis. Academic Press,
[9] Meinguet J. Multivariate interpolation at arbitrary knots made simple. J. Appl. Math. and Physics, 30:292-304, 1979. [10) O'Sullivan F., Yandell B. and Raynor W. J. Automatic smoothing of regression functions in generalized linear models. J. Amer. Statist. Assoc., 81:96-104, 1986.
[11) O'Sullivan F. An iterative approac.h to 2-dimensional Laplacian smoothing with application to image restoration. J. Amer. Statist. Assoc., 1989(to appear). [12) O'Sullivan F. Fast computation of fully automated log-density and log-hazard estimators. Siam J. Sci. Statist. Comp., 9:363-379, 1988. [13) Silverman B. W. On the estimation of a probability density function by the maximum perrauseu li..-.';:a.LllVVU metnoa. Ann. 10:795-810, pacxage of fortran subprograms National
v :.-
v.:>
.,
0
(]l
oi
o
o
-1 I
-""
o
I
o (]l
9 o
o
1.0
-l .....
c
co
I
I
-""
-""
a
a
I
I
a
a
01
01
F
x
OJ OJ
a a
(j)
a a
a
01
01
.....
..... a
0.0
0.02
0.04
0.06
..... I
a
I
a 01
r vJ
x
x
a
a
a a
a 01
.
(j)
a.
--1 F
,.:::>
.....
C
CD
by
Finbarr O'Sullivan
TECHNICAL REPORT No. 162 June 1989
Department of Statistics, GN-22 University of Washington Seattle, Wasbington9819S USA
Abstract An approach to multi-dimensional smoothing is introduced which is based on a penalized with a modified The choice of penalty stanaard multi-dim~nsiQnal.Lapladan SQ11aI'e errQrcharacteristics at-least on the interior of hyper-rectangular~omains, which has wide restorationproblems,computar tions are carried out using fast Fourier transforms. Non-linear smoothing is accomplished by iterative application of the linear smoothing technique. The iterative procedure can be accelerated by means of preconditioned conjugate gradients. The methods are implemented in one- and two-dimensional settings. Adaptive choice of the amount of smoothing is based on approximate cross-validation type scores. Some illustrations are given relating to scatter-plot smoothing, estimation of a logistic regression surface and density estimation. The asymptotic mean square error characteristics of the linear smoother are derived and the
smoothing splines even when these boundary conditions are not satisfied.
Discretized Laplacian Smoothing by Fourier Methods Finbarr O'Sullivan Dept. of Biostatistics and Statistics University of Washington Seattle, WA 98195 March 30, 1989
1
Introduction Laplacian Smoothing and Discretization
1.1
Consider a region of il c Rd. Suppose we have data (Xi, zi, Wi; i = 1,2,... where Xi E il, EZi = U(Xi) and Var(zi) ex: wil. A class of multi-dimensional smoothers are defined as minimizers of the penalized likelihoods of the form
(1) where X > 0 and L is a linear differential operator, see Cox [4]. The parameter ,\ > 0 controls the smoothness of the estimator. Particular choices for L lead to thin-plate and Laplacian smoothing splines which have been studied by Meinguet [9], Wahba and Wendelberger, [15] and others. For example, in l-dimension with Lu = u" we obtain cubic smoothing splines. In d-dimensions we let L = ~ - the Laplacian in Rd d
Lu(x)
= ~u(x) = ~
02 u(x ) OX]
(2)
A discretized version of the multi-dimensional smoother with Laplacian penalty is defined as follows: For some parameter h (which might be a vector) let il h be a finite dimensional grid defined on il. Suppose we have data Zh,Wh specified on ilh where EZh = Uh and Wh gives the weights associated with Zh (Zh and Wh could be obtained by projecting the observed potentially scattered data onto the grid ilh). We wish to estimate Uh. Let ~h be a matrix of order 4im(ilh) x di.m(rt h) which yieldsadiscretizedapproximation to the Laplacian, i.e, Uh. is a discretized versiQn of a function U defined on il then
(3) Let Bh = ~~~h' An s'th order Discreiized Laplacian smoother is defined as the minimizer of the penalized likelihood criterion
(4) Here Wh is regarded as a diagonal matrix. We note that the smooth solves the estimating equation (5) any real For ease of notation we IS
8
but for smootnmz purposes we are only interested in the 8 from
8
> O.
1.2
Outline
dnUellSHJUS greater than one on ( are O'Sullivan [11) and Wahba Wendelberger The goal of paper is to show a wen chosen discretization can simplifications in computations yet without comprimising performance cnarachenstics. The paper only considers the case where the region of interest is a d~dimensional rectangle. This is a situation having immediate application to a range of image processing problems. We begin in section 2 by describing a particular choice for Doh whose spectral decomposition is amenable to Fourier analysis. This greatly simplifies the computations. Extensions to non-constant weights is achieved by an iterative scheme whose convergence properties are analyzed. Conjugate gradient methods may be used to accelerate convergence. Adaptive .choice of the smoothing parameter is done using approximate cross-validation assessments. Non-linear .SJ.fioothingby iterativere-weighting is considered in section 3. Various illustrations .and timing information are presentedin section 4. Computing times are dominated by the fast Fourier transform algorithm. one-dimensional Monte-Carlo experiments shows that there may be be little or no difference between a cross-validated discretized Laplacian smoother and a corresponding cross-validated spline smoother. Some asymptotic mean square error performance characteristics are worked out in section 4. The asymptotic analysis uses only very basic tools. In the case that the target function is locally constant at the boundary of the region of interest the asymptotic characteristics are seen to match those of multidimensional smoothing splines.
2
2
Computational Methods amenaote to
Fourier In case the weiehts the computation trivial. non-constant weizhts an iterative algorithm convergence of is it is techniques may be used for acceleration.
2.1
Discretized Laplacian and Spectral Analysis
The region of interest n is a d-dimensional rectangle which for notational convenience we take to be [O,l]d. A lattice is defined on il as follows. Let h = (hI, h 2 , ••• , h d ) where hj = 1/ nj for j = 1,2, ... , d. ilh consists of the set of d-vectors Xhi = (XiI" •• ,Xij' ••• , Xid) with = .hj (6) for 1 :::; i :::;
»s. and 1:::; j
:::; d. There area total of
(7) grid points in this lattice. Let D j be nj x nj-matrix which yields a discrete approximation to {;z. Thus if v = (v(Xd,V(X2),.' .v(x n j ) ) is a nj-vector which is a discretized version of a function v(x) then D j has the property that
for i
= 1,2, ... , nj.
Given Dil a discrete approximation for 6.h is the N x N-matrix d
6.h
= I: 11 ® 12 ® ... Dj."
® t,
(9)
j=l
where Ij choices for D j are possible corresponding to Periodic, Dirichlet boundary conditions for example, see Buzbee et. al. [3]. A modified choice related to the Neuman discretization is -1 1 Dj=
1
0
1 1 -2
1
1
( 1 -2 1 1 -2 1 1 -1
0 3
np. as 1 -+ co. Proof:
iterative scheme may be exnressed as M u(l+1) = b + Du(l)
where M = (>.B h + fiJI], D stage 1, so from (24)
= -[Wh -
(24)
and b = WhZh. Let £(1)
= Uh.\ -
u(l)
be the error at
(25) for 1 = 1,2, .... Such a method is convergent whenever the largest eigenvalue of NI- 1 D is less than 1. The analysis is straightforward. Clearly y'Bf.y 2: 0 so provided ib :» 0 y'M y
=
>,y'.8hY
>
fiJy'y
+ fiJy'Y (26)
This implies y'M- 2y ::; fiJ-2 y'y for all N-vectors y. From here
II M-IDy 11 2 II y 11 2
y'DM-2Dy y'y
w- 2 y' D 2 y
11
= [1 + ~~II](P+280:)/4Iuhlll then
'2) 1('" ,1..2) L...
"'( ~A~~)2 1 28 = ( L... 1 + ,U hll [1-1- ~2lu»J(280:)/2 (/>11 II I
2
=
14
II
'Yll
(63)
(64)
v
Since b is positive sup;
(bi~:~2) 2
is finite for 0 ::S
0:.
v
b depends on h through
>"hv
but by using the bounds in (60) we can find a constant M independent of h such that supv
(br~:~2) 2 ::S M.
Thus for p ::S t ::S 2 with
0:
= (65)
Note that if D.U is continuous and u satisfies the Neuman boundary conditions then for Ihl sufficiently small with t = 2
(66) 5.2.2
r
Bounding the Expected Squared Error
is orthogonal so Var£h
= (T2I.
Thus the expected squared error is
Ihl[1
(67)
Using (60)
I: Ihl[1 + >"~vlP/2[1 + >..>..~~t2(T2 v
(68) where a = c~ and b cis. The summation may be analyzed by slightly generalizing the mtegral approximation argument in Craven and Wahba Let t Vj >..1/4s(Vj - 1) for Vi 1,2, ... ,nj + 1 and j 1,2, ... ,d. Also >..d/4s. Now for o ::S >.. ::S >"0
=
=
+
=
< u
+
- tvJ
)..~/ ", The summation over integral over the positive orthant
is a Riemann sum approximation to the d-dimensional
1/
(70)
Transforming to polar co-ordinates the integral is bounded by
1
00
Jyf
o
[c + ar 4 ]P/ 2 d-l [ dr 1 + br 4 S ]2 r
(71)
where At[ is a positive constant. (The r d - 1 term comes from the Jacobian of the transformation to polar co-ordinates.) Note the latter integral is convergent provided p/2s < 2 - d/4s. Thus wecan find a constant M (independent of h) for which
E as )..
-+
II U k:A -
Eu k:A
11 2kp-< MN- 1 )..-(p/2s+d/4s)
(72)
O. From here we have the following result:
Theorem 2 There is a constant M E
II Uk
- Uk>.
> 0 such that for p < min(t,2s - d/2)
lI~p::; M·
()..(t-p)/s
II Uk
lI~t +N- 1 )..-
(P/ 2S+ d/ 4S»)
(73)
An optimal upper bound on the asymptotic mean square error is obtained by choosing ).. so that the squared bias and variance are on the same order. If Llu is continuous and u satisfies fhe Neuman boundary conditions then this calculation gives (for t = 2 and p < min(t,2s d/2)) E II Uk - Uk>. lI~p::; M . N- 2 (t - p )/ (4s+ d) (74) which occurs when X = O(N-4s/(2t+d»). We conjecture that this upper bound also holds for t E (0,2) provided p < min(t,28 - d/2). These results parallel those obtained for multidimensional smoothing splines by Cox [4].
6
Acknowledgement
I am gratefull to Professor Gene Golub for a some valuable discussions related to this work.
References Akaike H. Information theory an extension of maximum likelihood principle. In Petrov B.N. and Caski F _, editors, Second Iniernation Symposium on Information Theory, pages 267-281, Budapest, 1973. Akademiai Kiado. [2J Buzbee B.L., Door F.vV., George J .A. and Golub G.H. The direct solutions of the discrete poisson equation on irregular regions. SL4iV! J. Numer. Anal., 8:722-736, 1971. [3J Buzbee B.L., Golub G.H. and Nielson C.W. On direct methods for solving poissons equation. SIAlY! J. Numer. Ana!., 7:627-656, 1970. [4J Cox D.D. Multivariate smoothing spline functions. SIAM J. Numer. Anal., 21:789-813, 1984.
[5) Craven P. and Wahba G. Smoothing noisy data with spline functions: estimating the correct degree of smoothing by the method of generalized cross-validation. Numer. Math., 31:377-403, 1979. [6) Golub G.H. and Van Loan C.F. Matrix Computations. North Oxford Academic, New York, 1983. [7] Huber P.J. Robust smoothing. In Launer R. L. and Wilkinson G.H., editors, Robustness in Statistics. Academic Press, New York, 1979. [8) Mardia K.V., Kent J.T. and Bibby J.M. London, 1979.
Multivariate Analysis. Academic Press,
[9] Meinguet J. Multivariate interpolation at arbitrary knots made simple. J. Appl. Math. and Physics, 30:292-304, 1979. [10) O'Sullivan F., Yandell B. and Raynor W. J. Automatic smoothing of regression functions in generalized linear models. J. Amer. Statist. Assoc., 81:96-104, 1986.
[11) O'Sullivan F. An iterative approac.h to 2-dimensional Laplacian smoothing with application to image restoration. J. Amer. Statist. Assoc., 1989(to appear). [12) O'Sullivan F. Fast computation of fully automated log-density and log-hazard estimators. Siam J. Sci. Statist. Comp., 9:363-379, 1988. [13) Silverman B. W. On the estimation of a probability density function by the maximum perrauseu li..-.';:a.LllVVU metnoa. Ann. 10:795-810, pacxage of fortran subprograms National
v :.-
v.:>
.,
0
(]l
oi
o
o
-1 I
-""
o
I
o (]l
9 o
o
1.0
-l .....
c
co
I
I
-""
-""
a
a
I
I
a
a
01
01
F
x
OJ OJ
a a
(j)
a a
a
01
01
.....
..... a
0.0
0.02
0.04
0.06
..... I
a
I
a 01
r vJ
x
x
a
a
a a
a 01
.
(j)
a.
--1 F
,.:::>
.....
C
CD
