Application of second generation wavelets to blind spherical deconvolution

1.2 MB·English

by  Thomas Vareschi

Application of second generation wavelets to blind
3 spherical deconvolution
1
0
2
n T. Vareschi ∗
a
J

5
1
Abstract

]
T Weadresstheproblemofsphericaldeconvolutioninanonpara-
S metricstatisticalframework,whereboththesignalandtheoperator
.
h kernelaresubjecttoerrormeasurements.Afterapreliminarytreat-
t ment ofthe kernel,we apply a thresholding procedureto the signal
a
m in a second generation wavelet basis. Under standard assumptions
on the kernel, we study the theoritical performance of the result-
[
ing algorithm in terms of Lp losses (p 1) on Besov spaces on
≥
1 the sphere. We hereby extend the application of second generation
v
spherical wavelets to the blind deconvolution framework [16]. The
9
procedure is furthermore adaptive with regard both to the target
8
2 functionsparsityandsmoothness,andthekernelblurringeffect.We
3 endwiththestudyofaconcreteexample,putting intoevidencethe
. improvement of our procedure on the recent blockwise-SVD algo-
1
0 rithm [6].
3
1
Keywords: Blind deconvolution; blockwise SVD; spherical deconvolution;
:
v secondgeneration wavelets; nonparametricadaptiveestimation; linearin-
i
X verse problems.
r Mathematical Subject Classification: 62G05, 62G99, 65J20, 65J22.
a
1 Introduction
1.1 Statistical framework
Considerthefollowingproblem:weaimatrecoveringasignalf L2(S2).
∈
f is not observed directly, but through the action of a blurring process
∗Universit´e Denis Diderot Paris 7 and CNRS-UMR 7099, 175 rue du Chevaleret
75013 Paris, France. E-mail: [email protected]
1
modeledbyalinearoperatorK.Tothisend,weconsidertheclassicwhite
noise model, where the available information is the noisy version
g = Kf +εW˙ (1.1)
ε
of f, where W˙ is a white noise on S2 and K : L2(S2) L2(S2) is a
→
measurable operator. We further restrict the shape of K by assuming
that K is a convolution operator on L2(S2), a classic framework ([13],
[17] and[16]) enjoying convenient mathematical properties (see Part1.2).
This model is equivalently formulated in a density estimation framework,
in which one aims at recovering the density f of a random variable X
on S2 from a n-sample (θ X ,...,θ X ) of Z = θX (with the analogy
1 1 n n
ε n−1/2),whereθ isarandomelementinthegroupof withdensity
3
∼ SO
h ,andZ hasadensityf L2(S2).Inpractice,theblurringoperatorK
θ Z ∈
is seldom directly observable and is itself subject to measurement errors.
This covers the cases where either K is unknown but approximated via
preliminary inference, or K is known but always observed with noise for
experimental reasons. The result is a noisy version K , satisfying
δ
K = K +δB˙ (1.2)
δ
where B˙ is a gaussian white noise on L2( ), independent from W˙ .
3
SO
The relevance of this generic setting was adequately discussed in Efro-
movich and Kolchinskii [8] and Hoffmann and Reiß [14], and covers nu-
merous fieldsof applications. Let usmention, forexample, image process-
ing,afieldwhich covers astronomy as well aselectronic microscopy where
an image, assimilated to a function f L2([0,1]2) is observed through
∈
its convolution with the Point Spread Function of the measuring device,
which hence requires to be estimated in first instance (see [23],[1]).
For u,u′,v,v′,w,w′ L2(S2), observable quantities obtained from
∈
1.1 and 1.2 hence take the form Kf,u + εα(u) (signal) and
h i
Ku,v + δβ(v,w) (operator) where α(u) (0, u ), β(v,w)
2
h i ∼ N k k ∼
(0, v 2 w 2) and E[α(u)α(u′)] = u,u′ L2(S2), E[β(v,w)β(v′,w′)] =
N k k k k h i
v,v′ L2(S2) w,w′ L2(S2).
h i h i
As we stated, we deal with a convolution on the 2-dimensional sphere.
Namely,ifZadmitsadensityhon withrespecttotheHaarmeasure,
3
SO
then Kf has the following expression
Kf(ω)= f(g−1ω)h(g)dg (1.3)
ZSO3
where dg is the Haar measure on . That is, f is averaged on a neigh-
3
SO
bourhood of ω with weight h(g) for each rotation g−1 applied to ω. This
problem, together with the introduction of needlets, is for example well
2
illustrated by the study of ultra high energy cosmic rays (UHECR).
An UHECR is a radiation hitting the earth with very high energy, and
whose physical origin is still unknown. Yet the understanding the mech-
anisms at work in this phenomenon is fundamental. Current hypothesis
involve pulsars, hypernovaes or black holes. Robust statistic tools are
heavily required, in order to properly estimate the density shape of the
radiation, which is highly related to the physical processes at stake in its
formation. One could ask, for example, whether the density is uniformly
distributed among the sphere, indicating a cosmological cause, or if it
is the superposition of localized spikes. In the latter case, it is crucial
to determine precisely the positions of this spikes. In practice however,
observations (X ,...,X ) of such radiations are often subject to various
1 n
physical perturbations, translated through the impulse response of the
measuring device. We modelize these by a random rotation θ, which is
to say we actually observe (θ X ,...θ X ) realisations of the random
1 1 n n
variable Z = θX. The difficulty of the problem is characterized by the
spreading of h around the identity : the less localized it is, the more dif-
θ
ficult the estimation of f should be. Moreover, the law of θ is not known
in general, even if some assumptions can restrict its shape. In this case,
preliminary inference is necessary, and leads to an estimator K of K.
δ
Case of a known operator
We shall concentrate here on the case where δ = 0, and describe the
path which finally led to the introduction and use of needlets in this
setting. Spherical harmonics constitute the most natural set of functions
to expand a target function f L2(S2), and present a structure highly
∈
compatible with deconvolution problems. It prompted Healy et al. [13]
to solve the deconvolution problem with their use, hereby reaching op-
timal L2 rates of convergence (Kim and Koo [17]). Unfortunately their
performances can prove quite poor in general cases, since they lack local-
ization in the spatial domain (see [12]). More recently, spherical wavelets
were introduced (Shr¨oder and Sweldens [25] and Narcowich and Ward
[19]) and have found various applications for a direct estimation of f, in-
cluding geophysics or atmospherics sciences (see for example Freeden and
Schreiner [9] or Freeden and Michel [10]). However, these wavelets, which
rely on a spatial construction, have an infinite support in the frequency
domain, and hence are not suited for the case of spherical deconvolution,
unlike spherical harmonics. The solution to this problem was brought by
Narcowich and Ward [19], who introduced a new set of functions, called
3
needets, which preserve the frequency localization of spherical harmonics
as well as the compatibility with inverse problems, all the while remedy-
ing their lack of spatial localization. Since then, needlets became widely
used in astrophysics (Marinucci et al. [18] or Guilloux et al. [12]) or brain
shape modeling (Tournier et al. [26]). In particular, Kerkyacharian et al.
[16] reached near-minimax rates of convergence for Lp losses (1 p )
≤ ≤ ∞
in the present spherical deconvolution setting.
Resolution when K is unknown : Galerkin projection
In the case of unknown operator K, the main methods involve SVD,
WVD and Galerkin schemes (see [3],[4],[14] for example). We now give
an overview of the so called Galerkin method and present its application
to blind-deconvolution. It is based upon on a discretization of 1.1 and
1.2 through the choice of appropriate test functions. Suppose we want
to recover f from the observation g = Kf. Let X,Y L2(S2) be two
⊂
finite dimensional subset which admit the respective orthogonal bases
ϕ = (ϕ ) and Φ =(Φ ) . The Galerkin approximation f of
k k=1,...,n k k∈1,...,n n
f is the solution of the equation
Kf ,v = g,v v Y
n
h i h i ∀ ∈
Kϕ ,Φ f ,ϕ = g,Φ k′ n (1.4)
k k′ n k k′
⇔ h ih i h i ∀ ≤
k≤n
X
f is easily computable, as the equivalent solution of the finite dimen-
n
sional linear system g = K f where g is the vector whose components
n n n n
are( g,ϕ ) andK thematrixwithentries( Kϕ ,ϕ′ ) .Hence,
h ki k≤n n h k ki k,k′≤n
this method relies on the discretization of the operator K, together with
the discretization of the function f.
Galerkin projection were successfully applied to blind deconvolution to
reach optimal rates of convergence on generic Hilbert spaces ([8]) or on
Besov spaces through wavelet-thresholding technics ([14], [4]).
Its remainstohandletwo practical problems:thealgorithm mustinclude
and articulate two essential steps, namely theinversion of K and thereg-
ularization of the datas through a projection/thresholding scheme . Note
that both the signal and the operator K can be subject to regularization
(see [14],[6]).
The second practical problem remains in chosing the right functions ϕ
k
andΦ .Thischoiceshouldideally answerthedilemmatofindasetwhich
k
is both compatible with the representation of f (through the belonging
to a certain set of functions) and with the structure of K (see [14], [6]).
4
Sphericalharmonicsrespondoptimalytotheprobleminthecaseofspher-
ical deconvolution on Sobolev spaces for a L2 error, since they realize a
blockwise-SVD decomposition of K, as shown in 1.1. More importantly
here, when δ is non negative, they allow a fine treatment of K thanks to
δ
the sparse structure of the original operator K, which allowed Delattre
et al. [6] to reach optimal L2-rates of convergence for a natural class of
operators and functions. Thus, we should always seek to preserve this
property of sparsity whenever possible.
1.2 Harmonic analysis on and S2
3
SO
The next part provides preliminary tools in order to apply a blockwise
scheme to the case of spherical deconvolution. It is a quick overview of
harmonic analysis on the spaces S2 and which is mostly inspired by
3
SO
Healy et al. [13]. Let us define the Euler matrices
cosϕ sinϕ 0 cosθ 0 sinθ
−
u(ϕ) = sinϕ cosϕ 0 , a(θ)= 0 1 0
   
0 0 1 sinθ 0 cosθ
−
   
where ϕ [0,2π), θ [0,π).
∈ ∈
Every rotation g in is the product of 3 elementary rotations :
3
SO
ε= u(ϕ)a(θ)u(ψ) (1.5)
where ϕ,ψ [0,2π), θ [0,π) are the Euler angles of g. Let l N and
∈ ∈ ∈
l m,n l. We also define the rotational harmonics
− ≤ ≤
Rl (ϕ,θ,ψ) = e−i(mϕ+nψ)Pl (cos(θ)) (1.6)
mn mn
where Pl are the second type Legendre functions.
mn
The functions Rl , l N, m , n 6 l are the eigenfunctions of the
mn ∈ | | | |
Laplace-Beltrami operatoron ,associated withtheeigenvalues 2l+1.
3
SO
Therefore, the system (√2l+1Rl )l forms a complete orthonormal
mn mn
basis of L2( ). Let h L2( ). For all l 0, the projection of h on
3 3
SO ∈ SO ≥
the space of rotational harmonics with degree l is
l
hˆl Rl
mn mn
m,n=−l
X
where hˆl is the (l,m,n) Fourier coefficient of h, defined by
mn
hˆl = h(g)Rl (g)dµg (1.7)
mn mn
ZSO3
5
An analogous study is available on S2. Any point ω S2 is determined
∈
by its spherical coordinates (θ,ϕ):
ω = (sin(θ)cos(ϕ),sin(θ)sin(ϕ),cos(θ)) (1.8)
where θ [0,π) and ϕ [0,2π). Let l a positive integer, m,n two inte-
∈ ∈
gers ranking from l to l. Define the following functions, known as the
−
spherical harmonics, on S2 :
2l+1(l m)!
Yl = ( 1)m − Pl (cos(θ)) (1.9)
m − s 4π (l+m)! m
where Pl are the Legendre functions. The set (Yl) constitutes an or-
m m
thonormal basis of L2(S2). Note H the space of spherical harmonics of
l
degree l and P the orthogonal projector onto H . For every function
l l
f L2(S2),
∈
l
P f = fˆl Yl
l m m
m=−l
X
where fˆl is the (l,m) Fourier coefficient of f, defined by
m
fˆl = f(ω)Yl (ω)dω
m m
S2
Z
The term of Blockwise-SVD finds its roots in the following proposition,
which expresses the link between Fourier coefficients of h f and those of
∗
h and f. A proof is present in [13].
Proposition 1.1 (Blockwiseproperty). Leth L2( )and f L2(S2)
3
∈ SO ∈
The Fourier coefficients of h f are
∗
l l
(h[f)l = hˆl fˆl = h Yl,Yl fˆl
∗ m mn n h ∗ n mi n
n=−l n=−l
X X
Hence, if K is a convolution operator over L2(S2) and f L2(S2), and
∈
if we note, Kl thematrix KYl,Yl M (C)andfl the vec-
h n mi |m|,|n|≤l ∈ 2l+1
tor (hf,Yml i)|m|≤l, Propos(cid:0)ition 1.1 tra(cid:1)nslates into Kf l = Klfl. Hence,
turning back to the Galerkin projection of K, take ϕ = Φ = (Yl ) ,
(cid:0) (cid:1) m m,l
m l. Proposition 1.1 actually implies that the Galerkin matrix
| | ≤
KYl1 ,Yl2 is sparse, with blocks Kl on its diagonal.
h m1 m2i li≥0,|mi|≤li,i=1,2
This justifies the denomination of blockwise-SVD decomposition. In the
(cid:0) (cid:1)
sequel, if f L2(S2) and K : L2(S2) L2(S2) is a convolution opera-
∈ →
tor, we will refer indifferently to P f or fl, and to P KP or Kl. Be-
l l l
sides, due to Parseval’s formula, we also have kPlfkL2(S2) = kflkℓ2(C2l+1)
6
and kPlKPlkL2(S2)→L2(S2) = kKlkℓ2(C2l+1)→ℓ2(C2l+1). Turning back to the
original problem and reminding Proposition 1.1, we can reformulate the
equivalent problem,obtained by projecting1.1 and1.2 onevery spaceH :
l
l 0, gl =Klfl+εW˙ l (1.10)
∀ ≥ ε
l 0, Kl =Kl+δB˙ l (1.11)
∀ ≥ δ
where W˙ l is a centered gaussian vector with covariance I , and B˙ l is
2l+1
a (2l+1) (2l+1) matrix whose entries are iid variables with common
×
law (0,1).
N
As we said, spherical harmonics however show great inconvenients when
used in the estimation of a generic function f L2(S2). We turn to the
∈
presentation of far more accute functions to this end.
2 Needlets
2.1 Construction of Needlets
Needlets were introduced in Narcowich et al. [20], and used in the frame-
work of density estimation on the sphere by Kerkyacharian and Picard
[15], Baldi et al. [2] and Kerkyacharian et al. [16]. As their construction
relies on a rearrangement of spherical harmonics, they inherit the very
useful stability properties of the latter in inverse problems (as expressed
in 1.1). In addition, whereas spherical harmonics’ supports spread all
over the sphere, needlets are almost exponentially localized around their
respective centers, thus allowing a fine multi-resolution analysis and a
description of very general regularity spaces on S2.
Needlet framework
∞
As we have seen, the following decomposition holds : L2(S2) = H .
l
l=0
The orthogonal projector P on H can be written M
l l
l
P (f) = L ( x,y )f(y)dy = Yl(x)Yl(y)f(y)dy
l l h i m m
Z Z m=−l
X
where L is the Legendre polynomial of degree l, and . stands for the
l
hi
usualscalar producton R2.Finally, thefact thatP is aprojector implies
l
7
the identity
L ( x,y )L ( y,z )dz =δ L ( x,z ) (2.1)
l k l,k l
S2 h i h i h i
Z
Littlewood-Paley decomposition
Let a be a ∞(R) symetric function, compactly supported in [ 1,1],
C −
decreasing on R+, such that for all x R, 0 a(x) 1 and for all
∈ ≤ ≤
x 1/2, a(x) = 1. Define , for all x R,
| | ≤ | | ∈
x
b2(x) = a( ) a(x)
2 −
b2 is a positive function, supported in [ 2; 1/2] [1/2;2], satisfying by
− −
construction
x S
x 1, b2( ) =1 (2.2)
∀| |≥ 2j
j≥0
X
Define the kernels
l l
Λ (x,y) = b2( )L ( x,y ), and M (x,y) = b( )L ( x,y )
j 2j l h i j 2j l h i
l≥0 l≥0
X X
(2.3)
and the associated operators
J
B f = Λ (x,y)f(y)dy and A f = B f (2.4)
j j J j
S2
Z j=−1
X
with the convention B f = P f. Note that the two sums in 2.3 are
−1 0
finite since b( l ) = 0 if l / L , where L is the set of integers between
2j ∈ j j
2j−1 and 2j+1 1. It is straightforward to show that, for all f L2(S2),
− ∈
f 2 = f,ψ 2 (2.5)
k k2 h j,ηi
Xj≥0ηX∈Zj
One of the main results in Narcowich et al. [21] is that A also mimicks
J
the bestpolynomial approximation of f with respect to . for all p 1,
p
kk ≥
as expressed in the following theorem:
Theorem 2.1. For all p [1, [, if f Lp(S2), then
∈ ∞ ∈
A f f 0 as J
J p
k − k → → ∞
, with uniform convergence if f 0(S2).
∈ C
8
Space discretization
Proposition 2.2 (Quadrature formula). . Note the set of polynoms
l
P
with degree less than l on S2. For each l 0, there exists a finite set Z
l
≥
of cubature points, and non negative reals (λ ) such that
η η∈Zl
f , f(ω)dω = λ f(η)
l η
∀ ∈ P S2
Z ηX∈Zl
Since b( l ) = 0 only if 2j−1 l < 2j+1 the function z M (x,z)
2j 6 ≤ 7→ j
belongs to , and z M (x,z)M (z,y) is an element of .
2j+1−1 j j 2j+2−2
P 7→ P
For more convenience, we will note Z = . Hence, B writes
2j+2−2 j j
Z
B (f)= λ M (x,η)M (η,y)dz f(y)dy
j η j j
S2
Z (cid:0)ηX∈Zj (cid:1)
= λ M (x,η) λ M (η,y)f(y)dy
η j η j
S2
ηX∈Zjp Z p
The functions ψ = λ M (.,η), are called needlets. Furthermore, it
j,η η j
can be prooved that the cubature points η and weights λ can be chosen
η
p
so that the two following conditions are verified,
c−122j card( ) c22j and c−12−2j λ c2−2j (2.6)
j η
≤ Z ≤ ≤ ≤
with a constant c> 0
2.2 Besov spaces
Properties of needlets
By construction, needlets are well localized in frequence ( ∞, compactly
C
supported). A crucial result proved by Narcowich et al. [21] shows that
they are furthermore near exponentially localized in space.
Theorem 2.3. Let j 0, η . For all M > 0, there exists C > 0
j M
≥ ∈ Z
such that
C 2j
x S2, ψ (x) M (2.7)
∀ ∈ | j,η | ≤ (1+2jd(x,η))M
where d(x,y) = arccos( x,y ) is the geodesic distance on the sphere.
h i
To illustrate this point, we represented two needlets of level j = 2,3 on
figure 1. The following properties are all consequences of this localization
property.
9
Figure 1: ASphericalrepresentationoftwoneedlets(levelj =2,3from lefttoright)
centeredaroundthepoint(0,0,1).Thedarkenedzonescorrespondtotheregionswhere
the needlet is high.
Proposition 2.4. For all p 1 (with the convention 1/ = 0), there
≥ ∞
exists c ,C ,D > 0 such that
p p p
cp2j(21−p1) ≤ kψj,ηkp ≤ Cp2j(21−p1) (2.8)
Proposition 2.5. For all p [1,+ ], there exists a constant C such
p
∈ ∞
that for all f Lp(S2),
∈
B (f) C λ ψ (2.9)
k j kp ≤ p | η|k j,ηkp η∈Zj ℓp
Moreover, λ ψ(cid:13)(cid:0) (cid:1) f(cid:13) (2.10)
| η|k(cid:13)j,ηkp η∈Zj ℓp ≤ k(cid:13)kp
(cid:13)(cid:0) (cid:1) (cid:13)
(cid:13) (cid:13)
Construction of Besov spaces
Besov spaces on the sphere naturally generalize the usual approximation
properties of regular functions, all the while being simply characterized
with the help of needlets. A complete description, and the proofs of the
resultsclaimedinthispartcanbefoundinNarcowichetal.[20]orKerky-
acharian and Picard [15]. Let f : S2 R be a measurable function and
7→
let E (π 1) be the distance of f to with respect to . , that is
k,π k Lπ
≥ P kk
E = inf f P
k,π Lπ
P∈Pkk − k
Theorem 2.6. Let 0 < s < , 1 p < and 0 < r .Let f Lπ.
∞ ≤ ∞ ≤ ∞ ∈
10