- Email: [email protected]
Sampling theorems on bounded domains
Journal of Computational and Applied Mathematics 221 (2008) 376–385 www.elsevier.com/locate/cam
Sampling theorems on bounded domains Massimo Fornasier a , Laura Gori b,∗ a Program in Applied and Computational Mathematics, Princeton University, Fine Hall, Washington Road, 08544-1000 Princeton NJ, USA b Dipartimento di Metodi e Modelli Matematici per le Scienze Applicate, Universit`a di Roma “La Sapienza”, Via A. Scarpa 16,
I-00161 Roma, Italy Received 17 January 2007; received in revised form 5 June 2007
Abstract This paper concerns with iterative schemes for the perfect reconstruction from nonuniform sampling of functions belonging to multiresolution spaces on bounded manifolds. Since the iterations converge uniformly, we can produce the corresponding iterative integration schemes that allow to recover the integral of functions belonging to multiresolution spaces from nonuniform sampling. We present also an error analysis and, in particular, we estimate the L 2 -error we produce in recovering smooth functions in H s , but not necessarily in any multiresolution space. The error analysis extends to integrals. Our results hold regardless of the dimension of the domain and for a variety of multiresolution spaces constructed from certain refinable bases formed by the so-called GP-functions. c 2007 Elsevier B.V. All rights reserved.
MSC: 65D05; 65D15; 65D32; 65T60; 94A20 Keywords: Nonuniform sampling; Quasi-interpolation; Multiresolution analysis; Iterative integration formulas
1. Introduction The reconstruction of a signal from sparse and nonuniform sampling data is a well-known problem [3] and mainly addressed to guessing or learning some relatively small missing part of a function by using the information of the relevant known part. Besides statistical learning, as refinements of the Shannon sampling theorem, mathematical methods and numerical algorithms based on deterministic techniques have been developed to compute missing parts of (band-limited) signals from few and sparse known sampling information. We refer in particular to the classical work of Feichtinger et al. [1,9– 11] and Gr¨ochenig et al. [2,14]. These methods and algorithms are essentially based on an iterative quasi-interpolation of the function defined on the Euclidean space by means of suitable series expansions of irregularly shifted (translated) basic functions or, in its discrete version [10], by discrete finite series expansion of complex exponentials. In this paper we want to present an axiomatic approach to the sampling problem for functions defined on compact subsets Ω of Rd . Because of the boundedness nature of the domain the results are nicely suited for applications to ∗ Corresponding author.
E-mail addresses: [email protected] (M. Fornasier), [email protected] (L. Gori). c 2007 Elsevier B.V. All rights reserved. 0377-0427/$ - see front matter doi:10.1016/j.cam.2007.10.037
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
377
signals in concrete situations. The reconstruction from nonuniform sampling of functions defined on intervals has also been considered in [10,14]. In [10] the functions are modeled as complex trigonometric polynomials and in [14] as restrictions of functions defined on R as linear combinations of shifted compactly supported functions. None of these two approaches can be really extended to the reconstruction of functions defined on manifolds, e.g., on the sphere. The approach that we present in this paper can be indeed generalized to compact manifolds with or without boundary and it is computationally as efficient as the methods in [10,14]. In the following, the functions are assumed to belong to level spaces of a multiresolution analysis directly constructed on the domain (or manifold) [4,7,8], and not as a restriction of functions defined on larger spaces. Moreover, we can associate wavelet bases to these multiresolution analyses to characterize function spaces of smooth functions, e.g., Sobolev spaces, and we derive certain error estimates which would not be available otherwise. The paper is organized as follows: Section 2 describes the abstract setting and the general requirements of the multiresolution spaces we use to model functions, and the corresponding sampling theorem for the reconstruction of functions from nonuniform sampling. Section 3 is devoted to the numerical implementation of the iterative reconstruction formulas derived in Section 2. We show also how the iterative formula can be used to implement iterative numerical integration schemes from nonuniform sampling. We derive error estimators for smooth functions belonging to certain Sobolev spaces in Section 4. The construction on compact manifolds of multiresolution analyses with the properties that allow the application of the sampling theorems is addressed in Section 5. We will show that we can construct a rather large variety of multiresolution spaces from certain refinable functions on the real line chosen in the class described in [13]. 2. Elementary axiomatic sampling theorems In the following Ω is a compact subset or a compact manifold in Rd . Let us assume that (MRA1) is given a multiresolution analysis of finite-dimensional spaces V = {V j } j∈N such that for all p ∈ (1, ∞) χΩ ∈ V0 (Ω ) ⊂ V1 (Ω ) ⊂ . . . V j (Ω ) ⊂ L p (Ω ),
[
V j (Ω )
L p (Ω )
= L p (Ω );
j
(MRA2) each space has a suitable basis of nonnegative continuously differentiable functions Φ j := {φ j,0 , . . . , φ j,N j } ⊂ C 1 (Ω , R+ ). Moreover, we assume that diam(supp(φ j,k )) ≤ C2− j , #{k : supp(φ j,k ) ∩ supp(φ j,k 0 ) 6= ∅} ≤ m for some m ∈ N, uniformly with respect to j, and k∇φ j,k k∞ ≤ C2 j (d/2+1) , for all k, and uniformly with respect to j; (MRA3) P j : C(Ω ) → V j (Ω ) is a bounded linear projector, i.e., P j f = f for all f ∈ V j (Ω ), for all j ∈ N. Since V j (Ω ) ⊂ C(Ω ) ⊂ L q (Ω ) for all q ∈ [1, ∞], we assume that P j are bounded on L q (Ω ) for the full range q ∈ [1, ∞]. S T M of sampling nodes in Ω is called ∆-dense if Ω = ( M B (x )) Definition 2.1. A set X = {x` }`=0 Ω , for `=0 δ ` 0 < δ := ∆−1 , where Bδ (x) is the ball centered at x of radius δ. We have the following immediate result. M ⊂ L ∞ (Ω ) with Lemma 2.2. For a given ∆-dense set X of sampling nodes, there is a partition of unity Ψ = {ψ` }`=0 the following properties
(1) 0 ≤ ψ` ≤ 1 for all ` = 0, . . . , M; (2) supp(ψ P M ` ) ⊂ Bδ (x` ); (3) `=0 ψ` ≡ 1. Associated to a ∆-dense set X of sampling nodes we can define the following quasi-interpolation operator given by Q Ψ ,X f :=
M X
f (x` )ψ` ,
for all f ∈ C(Ω ).
`=0
Let us formulate now the main result of this paper.
(2.1)
378
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
Theorem 2.3. If V = {V j } j∈N is a multiresolution analysis in L p (Ω ) for p ∈ (1, ∞) with properties (MRA1)–(MRA3), then for any j ∈ N there exists ∆ j > 0 large enough such that for any ∆ j -dense M
j sampling set X j := X (∆ j ) any function f j ∈ V j (Ω ) can be exactly recovered from its samples { f j (x j,` )}`=1 by the following iterative algorithm
(n+1)
fj
(n)
(n)
= P j Q Ψ j ,X j ( f j − f j ) + f j ,
n ≥ 1,
(0)
fj
= P j Q Ψ j ,X j f j .
(2.2)
In fact we have (n)
f j := lim f j ,
(2.3)
n→∞
where the convergence is uniform on Ω , and then it is valid in L q (Ω ) for all q ∈ [1, ∞]. For f ∈ C(Ω ) but f 6∈ V j (Ω ) (0) for any j, the iterations (2.2) applied with the initial value f j = P j Q Ψ j ,X j f will converge anyway to a function (∞)
fj
(∞)
with the property P j (Q Ψ j ,X j f j
− Q Ψ j ,X j f ) = 0.
The iterative scheme (2.2) is analogous to that appearing in [1]. Here it is adapted to our abstract setting for the reconstruction of functions defined on bounded domains. Before working out the proof of Theorem 2.3 we need to show the following technical lemma. Lemma 2.4. Under the notations and the assumptions of Theorem 2.3, for any f ∈ V j (Ω ), j ∈ N, and any δ > 0 we define the oscillation function by oscδ ( f )(x) :=
sup y∈Bδ (x)∩Ω
| f (x) − f (y)|,
x ∈ Ω.
(2.4)
Then, we have koscδ ( f )k∞ ≤ C j δ2 j (d/2+1) |Ω |k f k∞ .
(2.5)
Here and in the following |Ω | denotes the measure of the set Ω . PN j Proof. Let us assume f ∈ V j (Ω ). Therefore, f = k=0 ak ( f )φ j,k ∈ C(Ω ) ∩ L 2 (Ω ) and, if y ∈ Bδ (x), then | f (x) − f (y)| ≤
Nj X
|ak ( f )||φ j,k (x) − φ j,k (y)| ≤ k∇φ j,k k∞ |x − y|
k=0
≤ Cδ2 j (d/2+1)
Nj X
Nj X
|ak ( f )|
k=0
|ak ( f )| ≤ C(1 + N j )1/2 δ2 j (d/2+1)
Nj X
1/2 |ak ( f )|2
k=0
k=0
≤ C B(1 + N j )1/2 δ2 j (d/2+1) k f k2 ≤ C j δ2 j (d/2+1) |Ω |k f k∞ . PN j The estimate ( k=0 |ak ( f )|2 )1/2 ≤ Bk f k2 for some B > 0 is valid because Φ j is a basis for V j (Ω ). Here we have denoted C j = C B(1 + N j )1/2 . Hence (2.5) holds. Now we have all the necessary ingredients to prove Theorem 2.3. Proof. We claim that for ∆ j > 0 large enough and for any ∆ j -dense sampling set X j := X (∆ j ) we have k(I − P j Q Ψ j ,X j ) f j k∞ ≤ ηk f j k∞ ,
for all f j ∈ V j (Ω ),
(2.6)
for some 0 < η < 1. This would imply that P j Q Ψ j ,X j is a boundedly invertible operator on V j (Ω ) with inverse P n (P j Q Ψ j ,X j )−1 = ∞ n=0 (I − P j Q Ψ j ,X j ) , so that f j = (P j Q Ψ j ,X j )−1 P j Q Ψ j ,X j f j =
∞ X n=0
(I − P j Q Ψ j ,X j )n P j Q Ψ j ,X j f j ,
(2.7)
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
379
with uniform convergence. Of course, it is immediate to see that (2.2) and (2.3) are just a reformulation of (2.7). Let us show (2.6):
k(I − P j Q Ψ j ,X j ) f j k∞ = kP j (I − Q Ψ j ,X j ) f j k∞ ≤ C 00j (I − Q Ψ j ,X j ) f j ∞
Mj
Mj Mj X X X
00 00
= Cj fj ψ` − f j (x j,` )ψ` ≤ C j | f j (·) − f j (x j,` )|ψ`
`=0
`=0 `=0 ∞ ∞
Mj
X
0 ≤ C 00j |oscδ j f j (x j,` )|ψ`
≤ C j koscδ j f j k∞ .
`=0
∞
RESTORE[n, fjs ] → fcj :
INTEGRATE[n, fjs ] → I j :
fc :=
int := ω Tj M
j PQc (fjs ); j PQs (fjs );
fs := i = 0; While i ≤ n do i := i + 1 j fc := fc + PQc (fjs − fs ); j
fs := fs + PQs (fjs − fs ); enddo fcj := fc .
j f s; ψ φ˜ j j PQs (fjs );
fs := i = 0; While i ≤ n do i := i + 1 j int := int + ω Tj M ˜ (fjs − fs ); j
ψφ
fs := fs + PQs (fjs − fs ); enddo I j := int.
By applying Lemma 2.4, we finally have k(I − P j Q Ψ j ,X j ) f j k∞ ≤ C j δ j 2 j (d/2+1) |Ω |k f j k∞ . Therefore, if ∆ j = δ −1 if f ∈ C(Ω ), but not necessarily f ∈ V j (Ω ) for j > 0 is large enough, then we immediately have (2.6). Moreover, P n any j, then P j Q Ψ j ,X j f ∈ V j (Ω ) and P j Q Ψ j ,X j f = P j Q Ψ j ,X j ( ∞ n=0 (I − P j Q Ψ j ,X j ) P j Q Ψ j ,X j f ). Let us denote P (∞) (∞) n − Q Ψ j ,X j f ) = 0. fj := ( ∞ n=0 (I − P j Q Ψ j ,X j ) P j Q Ψ j ,X j f ) ∈ V j (Ω ). Then we have P j (Q Ψ j ,X j f j 3. Numerical implementation First of all it is important to discuss how to construct a possible projector P j . A natural choice is given by PN j Nj 2 ˜ ˜ ˜ Pj f = k=0 h f, φ j,k iφ j,k , where Φ j := {φ j,k }k= j ⊂ L (Ω ) is a biorthogonal dual basis for Φ j . In particular P j would result as a projection from L 2 (Ω ) onto V j (Ω ). We will assume that it is also bounded from C(Ω ) into V j (Ω ) both being endowed with the sup-norm. Therefore the problem of constructing a good projector is shifted to the construction of suitable dual bases. We will discuss this issue in Section 5. Now that we have a recipe how to construct P j , we should discuss how to implement numerically the general PM j scheme (2.2). Recalling from (2.1) that Q Ψ j ,X j f = for all f ∈ C(Ω ), we immediately `=0 f (x j,` )ψ j,` , PN j PM j PN j PM j ˜ have P j Q Ψ ,X f = h f (x j,` )ψ j,` , φ j,k iφ j,k = ( f (x j,` )hψ j,` , φ˜ j,k i)φ j,k . This suggests j
k=0
j
`=0
k=0
`=0
j to define the notations fsj := [ f j (x j,0 ), . . . , f j (x j,M j )]T , M ˜ := (hψ j,` , φ˜ j,k i)k=0,...,N j ;`=0,...,M j , 8sj := ψφ (φ j,k (x j,` ))`=0,...,M j ;k=0,...,N j , and 8cj := (φ j,k (τ i))i=i∈Zd ,τ i∈Ω ;k=0,...,N j for τ > 0 small, and we can write j
PQs (fs ) := 8sj M
j s f , ψ φ˜ j
j
PQc (fs ) := 8cj M
j s f . ψ φ˜ j
Since the iterative formula (2.2) in Theorem 2.3 converges uniformly to f j we have Z Z (n+1) f j (x)dx = lim fj (x)dx. Ω
n→∞ Ω
(3.1)
(3.2)
380
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
Fig. 1. The original surface belonging to V0 ([0, 3]2 ) = span{B0,2 (x + 2 − i)B0,2 (y + 2 − j) : i, j = 0, . . . , 4} to be reconstructed from random samples is shown on the left. On the right the result of the reconstruction is shown after 6000 iterations. Table 1 Application of the iterative integration Iter.
Approx. integral
Error. (×10−6 )
0 50 100 300 2000
1.07874 0.97298 0.953824 0.947611 0.947643
131 095.0 25 335.3 6 178.9 34.0 1.3
On the other hand Z Z Z (n) (n+1) fj (x)dx = P j Q Ψ j ,X j ( f j − f j )(x)dx + Ω
=
Ω
Mj Nj X X k=0 `=0
=
j ω Tj M ˜ (fsj ψφ
(n) (f j )s
(n)
( f (x j,` ) − f j (x j,` ))hψ j,` , φ˜ j,k i (n) − (f j )s ) +
(n) [ f j (x j,0 ), . . . ,
Z
Z Ω
(n)
Ω
f j (x)dx
Z φ j,k (x)dx +
Ω
(n)
f j (x)dx
(n)
Ω
f j (x)dx,
(n) f j (x j,M j )]T ,
R R and ω j = [ Ω φ j,0 (x)dx, . . . , Ω φ j,N j (x)dx]T . where := If f j ∈ V j (Ω ) is the function that we want to reconstruct from its vector of samples fjs = [ f j (x j,0 ), . . . , f 0 (x j,M j )]T , we can use the procedure RESTORE to approximate the vector fcj = [ f j (τ i) : i ∈ Zd , τ i ∈ Ω ]T of its sampling on a τ −1 -dense regular grid. The convergence of the discrete scheme is of course ensured by the uniform and pointwise convergence of the original iterative algorithm (2.2). See Fig. 1 for an example of the reconstruction of a surface. Moreover, we can formulate a corresponding iterative integration procedure INTEGRATE. In Table 1 we report the numerical results of the iterative integration of the surface of Fig. 1. Proposition 3.1. Under the assumptions of Theorem 2.3, the procedure INTEGRATE has the following properties: (i) It is linear, i.e., INTEGRATE[n, λfs + µgs ] = λ · INTEGRATE[n, fs ] + µ · INTEGRATE[n, gs ], for all λ, µ ∈ R. R (ii) for all f j ∈ V j (Ω ) it is Ω f j (x)dx = limn→∞ INTEGRATE[n, fsj ]. Proof. To show (i), it is sufficient to observe that all the operations executed in the procedure INTEGRATE are linear. The proof of (ii) is a straightforward application of (3.2). 4. Error analysis (∞)
In this section we address the estimate of the error k f − f j k2 , under the assumption that f has some regularity. In particular we will assume that f ∈ H s (Ω ), for s ≥ 0, where H s (Ω ) denotes the Sobolev space of order s. Moreover,
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
381
we also estimate the error that we produce by approximating the integral of f by INTEGRATE[∞, f s ]. The key tool is the characterization of H s (Ω ) by suitable (biorthogonal) wavelets: (W1) there exists a wavelet biorthogonal basis {Ψ j,k } j≥0,k∈J j associated to the MRA {V j (Ω )} j≥0 , so that any function f ∈ L 2 (Ω ) can be written as X X f = h f, φ˜ J,k iφ J,k + h f, Ψ˜ j,k iΨ j,k = PJ ( f ) + f J⊥ ; (4.1) j≥J,k∈J j
k=0,...,N J
(W2) f ∈ H s (Ω ), s ≥ 0, if and only if 1/2
c1 k f k H s (Ω ) ≤ kPJ f k2 +
X
22s j |h f, Ψ˜ j,k i|2
≤ c2 k f k H s (Ω ) ,
(4.2)
j≥J,k∈J j
for c1 , c2 > 0P independent of f . This characterization of H s (Ω ) and a straightforward computation ensure in particular that ( j≥J,k∈J j |h f, Ψ˜ j,k i|2 )1/2 ≤ c2 2−s J k f k H s (Ω ) ; dj
dj 2 (W3) kΨ R j,k k∞ ≤ CΨ 2 and #J j ≤ CJ 2 ; (W4) Ω Ψ j,k (x)dx = 0 for all j ≥ 0 and k ∈ J j .
Proposition 4.1. Let f ∈ C(Ω ) ∩ H s (Ω ) for s ≥ 0, but not necessarily f ∈ V j for any j ≥ 0. Under the assumptions P (∞) n of Theorem 2.3 we have that f j := S j f = ( ∞ n=0 (I − P j Q Ψ j ,X j ) P j Q Ψ j ,X j f ) ∈ V j (Ω ) has the property (∞)
k f − fj
k2 ≤ CkS j k L 2 →L 2 2−s j k f k H s (Ω ) .
(4.3) (∞)
Proof. With similar arguments as in the proofs of Lemma 2.4 and Theorem 2.3, we show that the map S j : f → f j is bounded with the L 2 -norm. Moreover this map coincides with the identity on V j . This together with the properties W1 and W2 imply that (∞)
k f − fj
k2 = k f − S j ( f − P j f + P j f )k2 = kS j ( f − P j f ) − ( f − S j (P j f ))k2 = kS j ( f − P j f ) − ( f − P j f )k2 ≤ (1 + kS j k L 2 →L 2 )k f − P j f k2 ≤ C 0 (1 + kS j k L 2 →L 2 )2−s j k f k H s (Ω ) .
Lemma 4.2. Let f ∈ C(Ω ). Under the assumptions of Theorem 2.3 we have that ∞ ∞ INTEGRATE[∞, f s ] ≤ |Ω |kP j k L →L k f k∞ , 1−η
(4.4)
where 0 < η < 1 is as in formula (2.6). Proof. By Theorem 2.3 and by (3.2), we have Z ! ∞ X s n INTEGRATE[∞, f ] = (I − P j Q Ψ j ,X j ) P j Q Ψ j ,X j f (x) dx Ω n=0
∞
X
|Ω |
kP j Q Ψ j ,X j f k∞ ≤ |Ω | (I − P j Q Ψ j ,X j )n P j Q Ψ j ,X j f ≤
n=0 1−η ∞
≤
|Ω |kP j k L ∞ →L ∞ kQ Ψ j ,X j k L ∞ →L ∞ 1−η
We conclude by observing that kQ Ψ j ,X j k L ∞ →L ∞ = 1. Then we have the following error estimate result.
k f k∞ .
382
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
Proposition 4.3. Assume 3d 2 < s. If the iterative integration formula is executed on a ∆ J -dense set X J , for J ≥ 0 and for ∆ J > 0 large enough, then there exists a constant C > 0 such that, for all f ∈ H s (Ω ), we have Z ( 3d 2 −s)(J +1) ∞ ∞ INTEGRATE[∞, f s ] − ≤ C |Ω |kPJ k L →L k f k H s 2 f (x)dx . 3d 1−η Ω 1 − 2 2 −s Proof. First of all observe that s > d2 . By the Sobolev embedding theorem we have f ∈ C(Ω ) and it makes sense to consider its sampling. Proposition 3.1(i) ensures the linearity of INTEGRATE and by (W1) Z X INTEGRATE[∞, f s ] − ˜ j,k i9 s f (x)dx = INTEGRATE[∞, PJ (f)s ] + INTEGRATE ∞, hf, 9 j,k Ω j≥J,k ! Z Z X − PJ ( f )(x)dx + h f, Ψ˜ j,k i Ψ j,k (x)dx . Ω Ω j≥J,k In the last equality we could exchange the integral with the sum because the sum converges absolutely due to d P P P d (W2) and (W3): j≥J,k |h f, Ψ˜ j,k i||Ψ j,k (x)| ≤ CΨ j≥J,k 2 j |h f, Ψ˜ j,k i| ≤ Ck f k H s j≥J,k 2( 2 −s) j ≤ C 0 k f k H s R P ( 3d 2 −s) j < ∞. Moreover, INTEGRATE is exact on V J . By Proposition 3.1(ii) and j≥J 2 Ω Ψ j,k (x)dx = 0 we have Z X s INTEGRATE[∞, f s ] − ˜ j,k i9 ] . f (x)dx = INTEGRATE[∞, hf, 9 j,k Ω j≥J,k By Lemma 4.2, the latter equation, and by using again (W2) and (W3), we obtain
Z
X |Ω |kPJ k L ∞ →L ∞
s INTEGRATE[∞, f ] − f (x)dx ≤ h f, Ψ˜ j,k iΨ j,k
1 − η Ω j≥J,k ≤ C0
|Ω |kPJ k L ∞ →L ∞ k f kH s 1−η
X
2( 2 −s) j = C 0 3d
j≥J
∞ ( 3d −s)(J +1) 2
2 |Ω |kPJ k L ∞ →L ∞ k f kH s . 3d 1−η 1 − 2 2 −s
5. Construction of MRA on domains and manifolds In this section we want briefly to recall the construction of multiresolution analyses on manifolds with properties MRA1–MRA3 for which the axiomatic sampling theorems illustrated in the previous sections can be applied. It is also ensured that corresponding wavelet bases with properties W1–W3 are also available. 5.1. Composite multiresolution analyses In [4,7,8] the construction of spline-type multiresolution analyses on rather general manifolds has been proposed. Such construction is based on the decomposition of the manifold Ω into submanifolds Ωi smoothly parametrized by M Ω , Ω = κ (), i = 1, . . . , M, where Ω ∩ Ω = ∅ for i 6= j, and κ : Rd → Rd 0 , := (0, 1)d , i.e., Ω = ∪i=1 i i i i j i ˜ Zi 0 d ≤ d are smooth regular functions. The idea is to define biorthogonal multiresolution analyses V () and V˜ Z i () j
j
on adjusted with suitable complementary boundary conditions Z i and Z˜ i [7], for each i = 1, . . . , M. Each MRA is lifted on Ωi by using the parametrization κi to define V j (Ωi ) := {ϕ ◦ κi−1 : ϕ ∈ V jZ i ()} and similarly V˜ j (Ωi ). The M on Ω so that boundary conditions are chosen to fit with a suitable bounded resolution of the identity R = {Ri }i=1 PM PM ∗ V j (Ω ) := i=1 Ri V j (Ωi ) and V˜ j (Ω ) := i=1 Ri V˜ j (Ωi ) will define global biorthogonal multiresolution analysis on Ω with properties MRA1–MRA3.
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
383
Fig. 2. GP functions are illustrated for n = 2 and n = 3 and for different choices of h > n − 1.
Fig. 3. Normalized GP functions adapted to the interval (0, 1) [12, Section 4] for n = 5 and h = 10. ˜ The main ingredient is the construction of biorthogonal refinable bases Φ Zj i := {φ j,0 , . . . , φ j,N j } and Φ˜ Zj i :=
{φ˜ j,0 , . . . , φ j,N j } to define the multiresolution spaces with complementary boundary conditions on = (0, 1)d . In [7] such bases have been derived from integer shifts of B-splines Θ N on the real line. Unfortunately the definition of the projectors Ri makes use of certain extension operators that show increasing norms as soon as higher smoothness is required. This can make the resulting bases on Ω ill conditioned. One way to compensate this drawback is to start at the very beginning with refinable functions that are better conditioned than B-splines. 5.2. GP bases generating a variety of MRA A more general class of positive compactly supported refinable functions on the real line has been Pn+1 n,h (n,h) characterized = k=0 ak ϕ(2x − k), where ak h in [13] as thesolutions i of the refinement equations ϕ(x) = n+1 n−1 −h h−n + 4(2 − 1) k−1 , n ≥ 2, h > n − 1, k = 0, . . . , n + 1. Let us call GP-functions the elements of 2 k this class. In particular B-splines are GP-functions for h = n. This generalization is worth doing for the following reason: It has been pointed out in [12] that they are better conditioned as bases Θ (n,h) := {ϕ(x −k) : k ∈ Z} for increasing values of h. This is due to their remarkable properties that they have a smaller and smaller essential support, see Fig. 2, while preserving the degree of smoothness, as soon as h increases. In particular, in [12] it is shown how such bases can be boundary adapted to intervals, see Fig. 3, together with certain biorthogonal duals as derived in [5] although no complementary boundary conditions have been considered yet. We postpone this second adaptation to a successive contribution. As already mentioned above, composite biorthogonal bases as derived on manifolds in [8] might exhibit either limited smoothness or high condition numbers. As proposed here, the use of better local bases certainly improves this drawback. However an alternative way to overcome this difficulty is to use frames instead of biorthogonal bases, see for example [6]. Such frames are constructed by Overlapping Domain Decompositions where the patches are again smooth images of . The fact that we do not need to implement interfaces through patches to preserve global smoothness reduces the ill conditioning of the global system. The use of frames do not affect, e.g., the core of the proofs of Theorem 2.3 and Proposition 4.1.
384
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
Fig. 4. The solid line represents the uniform logarithmic error in the reconstruction of the polynomial f (x) = (x − 12 )2 by using the B-spline space (5,5)
(5,h)
for h ∈ {4.1, 4.3, 4.5}. Vj . The dashed lines represent the uniform logarithmic error in the reconstruction of f (x) = (x − 21 )2 by means of V j We can see that for h = 4.1, we achieve the machine precision error with a 20% less the number of iterations than for the choice h = n = 5, corresponding to B-splines.
Fig. 5. Segmentation of the brain phantom. The curves are computed with RESTORE from the nonuniform sampling indicated by the white pixels. (n,h) and the quasi-interpolation operator Q Ψ ,X . Different curves are generated by changing the GP space V j
5.3. Implementing GP bases for the sampling problem As stated in the previous subsection GP-bases Θ (n,h) adapted to the interval have better condition numbers for increasing values of h > n − 1. Thus we are tempted to affirm that for any f ∈ V (n,h) := span Θ (n,h) for all h > n − 1, f will be reconstructed with increasing rate of convergence from a fixed ∆-dense sampling set, for increasing values of h. In particular, it has been shown [13] that Pn−3 ⊂ V (n,h) for all h > n − 1. Therefore we can test this claim on polynomials. Numerical experiments on polynomials show instead that the asymptotic rate of (5,h) convergence decreases for increasing values of h. However, this also means that, e.g., as in Fig. 4, by using V j , j = 5, and n − 1 < h < n we can anyway obtain better performances in the reconstruction than implementing (5,5) RESTORE by just using the B-spline space V j . 5.4. Changing parameters in quasi-interpolation M and a fixed space V , it is not ensured that there exists f ∈ V such that For a given sampling set {(x` , y` )}`=0 j j j f j (x` ) = y` for ` = 0, . . . , M. Nevertheless, as soon as the set is dense enough RESTORE will converge anyway and (∞) the resulting function f j will depend on several different parameters: the choice of the quasi-interpolation operator Q Ψ ,X and the choice of the space V j . The former can be tuned according to different partitions of unity Ψ , the latter is parametrized by the particular basis Φ j . Therefore we can tune the choice of Q Ψ ,X and V j in order to optimize (∞) M , see for example certain geometrical properties of the resulting f j with respect to the given data set {(x` , y` )}`=0 Fig. 5.
Acknowledgement Massimo Fornasier acknowledges the financial support provided by the European Union’s Human Potential Programme under contract MOIF-CT-2006-039438.
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
385
References [1] A. Aldroubi, H.G. Feichtinger, Exact iterative reconstruction algorithm for multivariate irregularly sampled functions in spline like spaces: The L p theory, Proc. Amer. Math. Soc. 126 (9) (1998) 2677–2686. [2] A. Aldroubi, K. Gr¨ochenig, Nonuniform sampling and reconstruction in shift-invariant spaces, SIAM Rev. 43 (4) (2001) 585–620. [3] J.J. Benedetto, P.J.S.G. Ferreira, Modern Sampling Theory: Mathematics and Applications, Birkh¨auser, 2000. [4] Z. Ciesielski, T. Figiel, Spline bases in classical function spaces on compact C ∞ manifolds. I and II, Studia Math. 76 (1983) 1–58 and 96–136. [5] M. Cotronei, M.L. Lo Cascio, A method for the construction of families of biorthogonal filters with prescribed properties, Adv. Comput. Math. 17 (3) (2002) 199–210. [6] S. Dahlke, M. Fornasier, T. Raasch, Adaptive frame methods for elliptic operator equations, Adv. Comput. Math. 27 (1) (2007) 27–63. [7] W. Dahmen, R. Schneider, Wavelets with complementary boundary conditions — Function spaces on the cube, Results Math. 34 (3–4) (1998) 255–293. [8] W. Dahmen, R. Schneider, Wavelets on manifolds I. Construction and domain decomposition, SIAM J. Math. Anal. 31 (1999) 184–230. [9] H.G. Feichtinger, K. Gr¨ochenig, Iterative reconstruction of multivariate band-limited functions from irregular sampling values, SIAM J. Math. Anal. 1 (1992) 244–261. [10] H.G. Feichtinger, K. Gr¨ochenig, T. Strohmer, Efficient numerical methods in non-uniform sampling theory, Numer. Math. 69 (1995) 423–440. [11] H.G. Feichtinger, T. Strohmer, Recovery missing segments and lines in images, in: Digital Image Recovery and Synthesis, Opt. Eng. 33 (10) (1994) 3283–3289 (special issue). [12] L. Gori, L. Pezza, F. Pitolli, Recent results on wavelet bases on the interval generated by GP refinable functions, Appl. Numer. Math. 51 (2004) 549–563. [13] L. Gori, F. Pitolli, A class of totally positive refinable functions, Rend. Math. 7 (20) (2000) 305–322. [14] K. Gr¨ochenig, H. Schwab, Fast local reconstruction methods for nonuniform sampling in shift invariant spaces, SIAM J. Matrix Anal. App. 24 (42) (2004) 899–913.
Sampling theorems on bounded domains Massimo Fornasier a , Laura Gori b,∗ a Program in Applied and Computational Mathematics, Princeton University, Fine Hall, Washington Road, 08544-1000 Princeton NJ, USA b Dipartimento di Metodi e Modelli Matematici per le Scienze Applicate, Universit`a di Roma “La Sapienza”, Via A. Scarpa 16,
I-00161 Roma, Italy Received 17 January 2007; received in revised form 5 June 2007
Abstract This paper concerns with iterative schemes for the perfect reconstruction from nonuniform sampling of functions belonging to multiresolution spaces on bounded manifolds. Since the iterations converge uniformly, we can produce the corresponding iterative integration schemes that allow to recover the integral of functions belonging to multiresolution spaces from nonuniform sampling. We present also an error analysis and, in particular, we estimate the L 2 -error we produce in recovering smooth functions in H s , but not necessarily in any multiresolution space. The error analysis extends to integrals. Our results hold regardless of the dimension of the domain and for a variety of multiresolution spaces constructed from certain refinable bases formed by the so-called GP-functions. c 2007 Elsevier B.V. All rights reserved.
MSC: 65D05; 65D15; 65D32; 65T60; 94A20 Keywords: Nonuniform sampling; Quasi-interpolation; Multiresolution analysis; Iterative integration formulas
1. Introduction The reconstruction of a signal from sparse and nonuniform sampling data is a well-known problem [3] and mainly addressed to guessing or learning some relatively small missing part of a function by using the information of the relevant known part. Besides statistical learning, as refinements of the Shannon sampling theorem, mathematical methods and numerical algorithms based on deterministic techniques have been developed to compute missing parts of (band-limited) signals from few and sparse known sampling information. We refer in particular to the classical work of Feichtinger et al. [1,9– 11] and Gr¨ochenig et al. [2,14]. These methods and algorithms are essentially based on an iterative quasi-interpolation of the function defined on the Euclidean space by means of suitable series expansions of irregularly shifted (translated) basic functions or, in its discrete version [10], by discrete finite series expansion of complex exponentials. In this paper we want to present an axiomatic approach to the sampling problem for functions defined on compact subsets Ω of Rd . Because of the boundedness nature of the domain the results are nicely suited for applications to ∗ Corresponding author.
E-mail addresses: [email protected] (M. Fornasier), [email protected] (L. Gori). c 2007 Elsevier B.V. All rights reserved. 0377-0427/$ - see front matter doi:10.1016/j.cam.2007.10.037
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
377
signals in concrete situations. The reconstruction from nonuniform sampling of functions defined on intervals has also been considered in [10,14]. In [10] the functions are modeled as complex trigonometric polynomials and in [14] as restrictions of functions defined on R as linear combinations of shifted compactly supported functions. None of these two approaches can be really extended to the reconstruction of functions defined on manifolds, e.g., on the sphere. The approach that we present in this paper can be indeed generalized to compact manifolds with or without boundary and it is computationally as efficient as the methods in [10,14]. In the following, the functions are assumed to belong to level spaces of a multiresolution analysis directly constructed on the domain (or manifold) [4,7,8], and not as a restriction of functions defined on larger spaces. Moreover, we can associate wavelet bases to these multiresolution analyses to characterize function spaces of smooth functions, e.g., Sobolev spaces, and we derive certain error estimates which would not be available otherwise. The paper is organized as follows: Section 2 describes the abstract setting and the general requirements of the multiresolution spaces we use to model functions, and the corresponding sampling theorem for the reconstruction of functions from nonuniform sampling. Section 3 is devoted to the numerical implementation of the iterative reconstruction formulas derived in Section 2. We show also how the iterative formula can be used to implement iterative numerical integration schemes from nonuniform sampling. We derive error estimators for smooth functions belonging to certain Sobolev spaces in Section 4. The construction on compact manifolds of multiresolution analyses with the properties that allow the application of the sampling theorems is addressed in Section 5. We will show that we can construct a rather large variety of multiresolution spaces from certain refinable functions on the real line chosen in the class described in [13]. 2. Elementary axiomatic sampling theorems In the following Ω is a compact subset or a compact manifold in Rd . Let us assume that (MRA1) is given a multiresolution analysis of finite-dimensional spaces V = {V j } j∈N such that for all p ∈ (1, ∞) χΩ ∈ V0 (Ω ) ⊂ V1 (Ω ) ⊂ . . . V j (Ω ) ⊂ L p (Ω ),
[
V j (Ω )
L p (Ω )
= L p (Ω );
j
(MRA2) each space has a suitable basis of nonnegative continuously differentiable functions Φ j := {φ j,0 , . . . , φ j,N j } ⊂ C 1 (Ω , R+ ). Moreover, we assume that diam(supp(φ j,k )) ≤ C2− j , #{k : supp(φ j,k ) ∩ supp(φ j,k 0 ) 6= ∅} ≤ m for some m ∈ N, uniformly with respect to j, and k∇φ j,k k∞ ≤ C2 j (d/2+1) , for all k, and uniformly with respect to j; (MRA3) P j : C(Ω ) → V j (Ω ) is a bounded linear projector, i.e., P j f = f for all f ∈ V j (Ω ), for all j ∈ N. Since V j (Ω ) ⊂ C(Ω ) ⊂ L q (Ω ) for all q ∈ [1, ∞], we assume that P j are bounded on L q (Ω ) for the full range q ∈ [1, ∞]. S T M of sampling nodes in Ω is called ∆-dense if Ω = ( M B (x )) Definition 2.1. A set X = {x` }`=0 Ω , for `=0 δ ` 0 < δ := ∆−1 , where Bδ (x) is the ball centered at x of radius δ. We have the following immediate result. M ⊂ L ∞ (Ω ) with Lemma 2.2. For a given ∆-dense set X of sampling nodes, there is a partition of unity Ψ = {ψ` }`=0 the following properties
(1) 0 ≤ ψ` ≤ 1 for all ` = 0, . . . , M; (2) supp(ψ P M ` ) ⊂ Bδ (x` ); (3) `=0 ψ` ≡ 1. Associated to a ∆-dense set X of sampling nodes we can define the following quasi-interpolation operator given by Q Ψ ,X f :=
M X
f (x` )ψ` ,
for all f ∈ C(Ω ).
`=0
Let us formulate now the main result of this paper.
(2.1)
378
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
Theorem 2.3. If V = {V j } j∈N is a multiresolution analysis in L p (Ω ) for p ∈ (1, ∞) with properties (MRA1)–(MRA3), then for any j ∈ N there exists ∆ j > 0 large enough such that for any ∆ j -dense M
j sampling set X j := X (∆ j ) any function f j ∈ V j (Ω ) can be exactly recovered from its samples { f j (x j,` )}`=1 by the following iterative algorithm
(n+1)
fj
(n)
(n)
= P j Q Ψ j ,X j ( f j − f j ) + f j ,
n ≥ 1,
(0)
fj
= P j Q Ψ j ,X j f j .
(2.2)
In fact we have (n)
f j := lim f j ,
(2.3)
n→∞
where the convergence is uniform on Ω , and then it is valid in L q (Ω ) for all q ∈ [1, ∞]. For f ∈ C(Ω ) but f 6∈ V j (Ω ) (0) for any j, the iterations (2.2) applied with the initial value f j = P j Q Ψ j ,X j f will converge anyway to a function (∞)
fj
(∞)
with the property P j (Q Ψ j ,X j f j
− Q Ψ j ,X j f ) = 0.
The iterative scheme (2.2) is analogous to that appearing in [1]. Here it is adapted to our abstract setting for the reconstruction of functions defined on bounded domains. Before working out the proof of Theorem 2.3 we need to show the following technical lemma. Lemma 2.4. Under the notations and the assumptions of Theorem 2.3, for any f ∈ V j (Ω ), j ∈ N, and any δ > 0 we define the oscillation function by oscδ ( f )(x) :=
sup y∈Bδ (x)∩Ω
| f (x) − f (y)|,
x ∈ Ω.
(2.4)
Then, we have koscδ ( f )k∞ ≤ C j δ2 j (d/2+1) |Ω |k f k∞ .
(2.5)
Here and in the following |Ω | denotes the measure of the set Ω . PN j Proof. Let us assume f ∈ V j (Ω ). Therefore, f = k=0 ak ( f )φ j,k ∈ C(Ω ) ∩ L 2 (Ω ) and, if y ∈ Bδ (x), then | f (x) − f (y)| ≤
Nj X
|ak ( f )||φ j,k (x) − φ j,k (y)| ≤ k∇φ j,k k∞ |x − y|
k=0
≤ Cδ2 j (d/2+1)
Nj X
Nj X
|ak ( f )|
k=0
|ak ( f )| ≤ C(1 + N j )1/2 δ2 j (d/2+1)
Nj X
1/2 |ak ( f )|2
k=0
k=0
≤ C B(1 + N j )1/2 δ2 j (d/2+1) k f k2 ≤ C j δ2 j (d/2+1) |Ω |k f k∞ . PN j The estimate ( k=0 |ak ( f )|2 )1/2 ≤ Bk f k2 for some B > 0 is valid because Φ j is a basis for V j (Ω ). Here we have denoted C j = C B(1 + N j )1/2 . Hence (2.5) holds. Now we have all the necessary ingredients to prove Theorem 2.3. Proof. We claim that for ∆ j > 0 large enough and for any ∆ j -dense sampling set X j := X (∆ j ) we have k(I − P j Q Ψ j ,X j ) f j k∞ ≤ ηk f j k∞ ,
for all f j ∈ V j (Ω ),
(2.6)
for some 0 < η < 1. This would imply that P j Q Ψ j ,X j is a boundedly invertible operator on V j (Ω ) with inverse P n (P j Q Ψ j ,X j )−1 = ∞ n=0 (I − P j Q Ψ j ,X j ) , so that f j = (P j Q Ψ j ,X j )−1 P j Q Ψ j ,X j f j =
∞ X n=0
(I − P j Q Ψ j ,X j )n P j Q Ψ j ,X j f j ,
(2.7)
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
379
with uniform convergence. Of course, it is immediate to see that (2.2) and (2.3) are just a reformulation of (2.7). Let us show (2.6):
k(I − P j Q Ψ j ,X j ) f j k∞ = kP j (I − Q Ψ j ,X j ) f j k∞ ≤ C 00j (I − Q Ψ j ,X j ) f j ∞
Mj
Mj Mj X X X
00 00
= Cj fj ψ` − f j (x j,` )ψ` ≤ C j | f j (·) − f j (x j,` )|ψ`
`=0
`=0 `=0 ∞ ∞
Mj
X
0 ≤ C 00j |oscδ j f j (x j,` )|ψ`
≤ C j koscδ j f j k∞ .
`=0
∞
RESTORE[n, fjs ] → fcj :
INTEGRATE[n, fjs ] → I j :
fc :=
int := ω Tj M
j PQc (fjs ); j PQs (fjs );
fs := i = 0; While i ≤ n do i := i + 1 j fc := fc + PQc (fjs − fs ); j
fs := fs + PQs (fjs − fs ); enddo fcj := fc .
j f s; ψ φ˜ j j PQs (fjs );
fs := i = 0; While i ≤ n do i := i + 1 j int := int + ω Tj M ˜ (fjs − fs ); j
ψφ
fs := fs + PQs (fjs − fs ); enddo I j := int.
By applying Lemma 2.4, we finally have k(I − P j Q Ψ j ,X j ) f j k∞ ≤ C j δ j 2 j (d/2+1) |Ω |k f j k∞ . Therefore, if ∆ j = δ −1 if f ∈ C(Ω ), but not necessarily f ∈ V j (Ω ) for j > 0 is large enough, then we immediately have (2.6). Moreover, P n any j, then P j Q Ψ j ,X j f ∈ V j (Ω ) and P j Q Ψ j ,X j f = P j Q Ψ j ,X j ( ∞ n=0 (I − P j Q Ψ j ,X j ) P j Q Ψ j ,X j f ). Let us denote P (∞) (∞) n − Q Ψ j ,X j f ) = 0. fj := ( ∞ n=0 (I − P j Q Ψ j ,X j ) P j Q Ψ j ,X j f ) ∈ V j (Ω ). Then we have P j (Q Ψ j ,X j f j 3. Numerical implementation First of all it is important to discuss how to construct a possible projector P j . A natural choice is given by PN j Nj 2 ˜ ˜ ˜ Pj f = k=0 h f, φ j,k iφ j,k , where Φ j := {φ j,k }k= j ⊂ L (Ω ) is a biorthogonal dual basis for Φ j . In particular P j would result as a projection from L 2 (Ω ) onto V j (Ω ). We will assume that it is also bounded from C(Ω ) into V j (Ω ) both being endowed with the sup-norm. Therefore the problem of constructing a good projector is shifted to the construction of suitable dual bases. We will discuss this issue in Section 5. Now that we have a recipe how to construct P j , we should discuss how to implement numerically the general PM j scheme (2.2). Recalling from (2.1) that Q Ψ j ,X j f = for all f ∈ C(Ω ), we immediately `=0 f (x j,` )ψ j,` , PN j PM j PN j PM j ˜ have P j Q Ψ ,X f = h f (x j,` )ψ j,` , φ j,k iφ j,k = ( f (x j,` )hψ j,` , φ˜ j,k i)φ j,k . This suggests j
k=0
j
`=0
k=0
`=0
j to define the notations fsj := [ f j (x j,0 ), . . . , f j (x j,M j )]T , M ˜ := (hψ j,` , φ˜ j,k i)k=0,...,N j ;`=0,...,M j , 8sj := ψφ (φ j,k (x j,` ))`=0,...,M j ;k=0,...,N j , and 8cj := (φ j,k (τ i))i=i∈Zd ,τ i∈Ω ;k=0,...,N j for τ > 0 small, and we can write j
PQs (fs ) := 8sj M
j s f , ψ φ˜ j
j
PQc (fs ) := 8cj M
j s f . ψ φ˜ j
Since the iterative formula (2.2) in Theorem 2.3 converges uniformly to f j we have Z Z (n+1) f j (x)dx = lim fj (x)dx. Ω
n→∞ Ω
(3.1)
(3.2)
380
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
Fig. 1. The original surface belonging to V0 ([0, 3]2 ) = span{B0,2 (x + 2 − i)B0,2 (y + 2 − j) : i, j = 0, . . . , 4} to be reconstructed from random samples is shown on the left. On the right the result of the reconstruction is shown after 6000 iterations. Table 1 Application of the iterative integration Iter.
Approx. integral
Error. (×10−6 )
0 50 100 300 2000
1.07874 0.97298 0.953824 0.947611 0.947643
131 095.0 25 335.3 6 178.9 34.0 1.3
On the other hand Z Z Z (n) (n+1) fj (x)dx = P j Q Ψ j ,X j ( f j − f j )(x)dx + Ω
=
Ω
Mj Nj X X k=0 `=0
=
j ω Tj M ˜ (fsj ψφ
(n) (f j )s
(n)
( f (x j,` ) − f j (x j,` ))hψ j,` , φ˜ j,k i (n) − (f j )s ) +
(n) [ f j (x j,0 ), . . . ,
Z
Z Ω
(n)
Ω
f j (x)dx
Z φ j,k (x)dx +
Ω
(n)
f j (x)dx
(n)
Ω
f j (x)dx,
(n) f j (x j,M j )]T ,
R R and ω j = [ Ω φ j,0 (x)dx, . . . , Ω φ j,N j (x)dx]T . where := If f j ∈ V j (Ω ) is the function that we want to reconstruct from its vector of samples fjs = [ f j (x j,0 ), . . . , f 0 (x j,M j )]T , we can use the procedure RESTORE to approximate the vector fcj = [ f j (τ i) : i ∈ Zd , τ i ∈ Ω ]T of its sampling on a τ −1 -dense regular grid. The convergence of the discrete scheme is of course ensured by the uniform and pointwise convergence of the original iterative algorithm (2.2). See Fig. 1 for an example of the reconstruction of a surface. Moreover, we can formulate a corresponding iterative integration procedure INTEGRATE. In Table 1 we report the numerical results of the iterative integration of the surface of Fig. 1. Proposition 3.1. Under the assumptions of Theorem 2.3, the procedure INTEGRATE has the following properties: (i) It is linear, i.e., INTEGRATE[n, λfs + µgs ] = λ · INTEGRATE[n, fs ] + µ · INTEGRATE[n, gs ], for all λ, µ ∈ R. R (ii) for all f j ∈ V j (Ω ) it is Ω f j (x)dx = limn→∞ INTEGRATE[n, fsj ]. Proof. To show (i), it is sufficient to observe that all the operations executed in the procedure INTEGRATE are linear. The proof of (ii) is a straightforward application of (3.2). 4. Error analysis (∞)
In this section we address the estimate of the error k f − f j k2 , under the assumption that f has some regularity. In particular we will assume that f ∈ H s (Ω ), for s ≥ 0, where H s (Ω ) denotes the Sobolev space of order s. Moreover,
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
381
we also estimate the error that we produce by approximating the integral of f by INTEGRATE[∞, f s ]. The key tool is the characterization of H s (Ω ) by suitable (biorthogonal) wavelets: (W1) there exists a wavelet biorthogonal basis {Ψ j,k } j≥0,k∈J j associated to the MRA {V j (Ω )} j≥0 , so that any function f ∈ L 2 (Ω ) can be written as X X f = h f, φ˜ J,k iφ J,k + h f, Ψ˜ j,k iΨ j,k = PJ ( f ) + f J⊥ ; (4.1) j≥J,k∈J j
k=0,...,N J
(W2) f ∈ H s (Ω ), s ≥ 0, if and only if 1/2
c1 k f k H s (Ω ) ≤ kPJ f k2 +
X
22s j |h f, Ψ˜ j,k i|2
≤ c2 k f k H s (Ω ) ,
(4.2)
j≥J,k∈J j
for c1 , c2 > 0P independent of f . This characterization of H s (Ω ) and a straightforward computation ensure in particular that ( j≥J,k∈J j |h f, Ψ˜ j,k i|2 )1/2 ≤ c2 2−s J k f k H s (Ω ) ; dj
dj 2 (W3) kΨ R j,k k∞ ≤ CΨ 2 and #J j ≤ CJ 2 ; (W4) Ω Ψ j,k (x)dx = 0 for all j ≥ 0 and k ∈ J j .
Proposition 4.1. Let f ∈ C(Ω ) ∩ H s (Ω ) for s ≥ 0, but not necessarily f ∈ V j for any j ≥ 0. Under the assumptions P (∞) n of Theorem 2.3 we have that f j := S j f = ( ∞ n=0 (I − P j Q Ψ j ,X j ) P j Q Ψ j ,X j f ) ∈ V j (Ω ) has the property (∞)
k f − fj
k2 ≤ CkS j k L 2 →L 2 2−s j k f k H s (Ω ) .
(4.3) (∞)
Proof. With similar arguments as in the proofs of Lemma 2.4 and Theorem 2.3, we show that the map S j : f → f j is bounded with the L 2 -norm. Moreover this map coincides with the identity on V j . This together with the properties W1 and W2 imply that (∞)
k f − fj
k2 = k f − S j ( f − P j f + P j f )k2 = kS j ( f − P j f ) − ( f − S j (P j f ))k2 = kS j ( f − P j f ) − ( f − P j f )k2 ≤ (1 + kS j k L 2 →L 2 )k f − P j f k2 ≤ C 0 (1 + kS j k L 2 →L 2 )2−s j k f k H s (Ω ) .
Lemma 4.2. Let f ∈ C(Ω ). Under the assumptions of Theorem 2.3 we have that ∞ ∞ INTEGRATE[∞, f s ] ≤ |Ω |kP j k L →L k f k∞ , 1−η
(4.4)
where 0 < η < 1 is as in formula (2.6). Proof. By Theorem 2.3 and by (3.2), we have Z ! ∞ X s n INTEGRATE[∞, f ] = (I − P j Q Ψ j ,X j ) P j Q Ψ j ,X j f (x) dx Ω n=0
∞
X
|Ω |
kP j Q Ψ j ,X j f k∞ ≤ |Ω | (I − P j Q Ψ j ,X j )n P j Q Ψ j ,X j f ≤
n=0 1−η ∞
≤
|Ω |kP j k L ∞ →L ∞ kQ Ψ j ,X j k L ∞ →L ∞ 1−η
We conclude by observing that kQ Ψ j ,X j k L ∞ →L ∞ = 1. Then we have the following error estimate result.
k f k∞ .
382
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
Proposition 4.3. Assume 3d 2 < s. If the iterative integration formula is executed on a ∆ J -dense set X J , for J ≥ 0 and for ∆ J > 0 large enough, then there exists a constant C > 0 such that, for all f ∈ H s (Ω ), we have Z ( 3d 2 −s)(J +1) ∞ ∞ INTEGRATE[∞, f s ] − ≤ C |Ω |kPJ k L →L k f k H s 2 f (x)dx . 3d 1−η Ω 1 − 2 2 −s Proof. First of all observe that s > d2 . By the Sobolev embedding theorem we have f ∈ C(Ω ) and it makes sense to consider its sampling. Proposition 3.1(i) ensures the linearity of INTEGRATE and by (W1) Z X INTEGRATE[∞, f s ] − ˜ j,k i9 s f (x)dx = INTEGRATE[∞, PJ (f)s ] + INTEGRATE ∞, hf, 9 j,k Ω j≥J,k ! Z Z X − PJ ( f )(x)dx + h f, Ψ˜ j,k i Ψ j,k (x)dx . Ω Ω j≥J,k In the last equality we could exchange the integral with the sum because the sum converges absolutely due to d P P P d (W2) and (W3): j≥J,k |h f, Ψ˜ j,k i||Ψ j,k (x)| ≤ CΨ j≥J,k 2 j |h f, Ψ˜ j,k i| ≤ Ck f k H s j≥J,k 2( 2 −s) j ≤ C 0 k f k H s R P ( 3d 2 −s) j < ∞. Moreover, INTEGRATE is exact on V J . By Proposition 3.1(ii) and j≥J 2 Ω Ψ j,k (x)dx = 0 we have Z X s INTEGRATE[∞, f s ] − ˜ j,k i9 ] . f (x)dx = INTEGRATE[∞, hf, 9 j,k Ω j≥J,k By Lemma 4.2, the latter equation, and by using again (W2) and (W3), we obtain
Z
X |Ω |kPJ k L ∞ →L ∞
s INTEGRATE[∞, f ] − f (x)dx ≤ h f, Ψ˜ j,k iΨ j,k
1 − η Ω j≥J,k ≤ C0
|Ω |kPJ k L ∞ →L ∞ k f kH s 1−η
X
2( 2 −s) j = C 0 3d
j≥J
∞ ( 3d −s)(J +1) 2
2 |Ω |kPJ k L ∞ →L ∞ k f kH s . 3d 1−η 1 − 2 2 −s
5. Construction of MRA on domains and manifolds In this section we want briefly to recall the construction of multiresolution analyses on manifolds with properties MRA1–MRA3 for which the axiomatic sampling theorems illustrated in the previous sections can be applied. It is also ensured that corresponding wavelet bases with properties W1–W3 are also available. 5.1. Composite multiresolution analyses In [4,7,8] the construction of spline-type multiresolution analyses on rather general manifolds has been proposed. Such construction is based on the decomposition of the manifold Ω into submanifolds Ωi smoothly parametrized by M Ω , Ω = κ (), i = 1, . . . , M, where Ω ∩ Ω = ∅ for i 6= j, and κ : Rd → Rd 0 , := (0, 1)d , i.e., Ω = ∪i=1 i i i i j i ˜ Zi 0 d ≤ d are smooth regular functions. The idea is to define biorthogonal multiresolution analyses V () and V˜ Z i () j
j
on adjusted with suitable complementary boundary conditions Z i and Z˜ i [7], for each i = 1, . . . , M. Each MRA is lifted on Ωi by using the parametrization κi to define V j (Ωi ) := {ϕ ◦ κi−1 : ϕ ∈ V jZ i ()} and similarly V˜ j (Ωi ). The M on Ω so that boundary conditions are chosen to fit with a suitable bounded resolution of the identity R = {Ri }i=1 PM PM ∗ V j (Ω ) := i=1 Ri V j (Ωi ) and V˜ j (Ω ) := i=1 Ri V˜ j (Ωi ) will define global biorthogonal multiresolution analysis on Ω with properties MRA1–MRA3.
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
383
Fig. 2. GP functions are illustrated for n = 2 and n = 3 and for different choices of h > n − 1.
Fig. 3. Normalized GP functions adapted to the interval (0, 1) [12, Section 4] for n = 5 and h = 10. ˜ The main ingredient is the construction of biorthogonal refinable bases Φ Zj i := {φ j,0 , . . . , φ j,N j } and Φ˜ Zj i :=
{φ˜ j,0 , . . . , φ j,N j } to define the multiresolution spaces with complementary boundary conditions on = (0, 1)d . In [7] such bases have been derived from integer shifts of B-splines Θ N on the real line. Unfortunately the definition of the projectors Ri makes use of certain extension operators that show increasing norms as soon as higher smoothness is required. This can make the resulting bases on Ω ill conditioned. One way to compensate this drawback is to start at the very beginning with refinable functions that are better conditioned than B-splines. 5.2. GP bases generating a variety of MRA A more general class of positive compactly supported refinable functions on the real line has been Pn+1 n,h (n,h) characterized = k=0 ak ϕ(2x − k), where ak h in [13] as thesolutions i of the refinement equations ϕ(x) = n+1 n−1 −h h−n + 4(2 − 1) k−1 , n ≥ 2, h > n − 1, k = 0, . . . , n + 1. Let us call GP-functions the elements of 2 k this class. In particular B-splines are GP-functions for h = n. This generalization is worth doing for the following reason: It has been pointed out in [12] that they are better conditioned as bases Θ (n,h) := {ϕ(x −k) : k ∈ Z} for increasing values of h. This is due to their remarkable properties that they have a smaller and smaller essential support, see Fig. 2, while preserving the degree of smoothness, as soon as h increases. In particular, in [12] it is shown how such bases can be boundary adapted to intervals, see Fig. 3, together with certain biorthogonal duals as derived in [5] although no complementary boundary conditions have been considered yet. We postpone this second adaptation to a successive contribution. As already mentioned above, composite biorthogonal bases as derived on manifolds in [8] might exhibit either limited smoothness or high condition numbers. As proposed here, the use of better local bases certainly improves this drawback. However an alternative way to overcome this difficulty is to use frames instead of biorthogonal bases, see for example [6]. Such frames are constructed by Overlapping Domain Decompositions where the patches are again smooth images of . The fact that we do not need to implement interfaces through patches to preserve global smoothness reduces the ill conditioning of the global system. The use of frames do not affect, e.g., the core of the proofs of Theorem 2.3 and Proposition 4.1.
384
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
Fig. 4. The solid line represents the uniform logarithmic error in the reconstruction of the polynomial f (x) = (x − 12 )2 by using the B-spline space (5,5)
(5,h)
for h ∈ {4.1, 4.3, 4.5}. Vj . The dashed lines represent the uniform logarithmic error in the reconstruction of f (x) = (x − 21 )2 by means of V j We can see that for h = 4.1, we achieve the machine precision error with a 20% less the number of iterations than for the choice h = n = 5, corresponding to B-splines.
Fig. 5. Segmentation of the brain phantom. The curves are computed with RESTORE from the nonuniform sampling indicated by the white pixels. (n,h) and the quasi-interpolation operator Q Ψ ,X . Different curves are generated by changing the GP space V j
5.3. Implementing GP bases for the sampling problem As stated in the previous subsection GP-bases Θ (n,h) adapted to the interval have better condition numbers for increasing values of h > n − 1. Thus we are tempted to affirm that for any f ∈ V (n,h) := span Θ (n,h) for all h > n − 1, f will be reconstructed with increasing rate of convergence from a fixed ∆-dense sampling set, for increasing values of h. In particular, it has been shown [13] that Pn−3 ⊂ V (n,h) for all h > n − 1. Therefore we can test this claim on polynomials. Numerical experiments on polynomials show instead that the asymptotic rate of (5,h) convergence decreases for increasing values of h. However, this also means that, e.g., as in Fig. 4, by using V j , j = 5, and n − 1 < h < n we can anyway obtain better performances in the reconstruction than implementing (5,5) RESTORE by just using the B-spline space V j . 5.4. Changing parameters in quasi-interpolation M and a fixed space V , it is not ensured that there exists f ∈ V such that For a given sampling set {(x` , y` )}`=0 j j j f j (x` ) = y` for ` = 0, . . . , M. Nevertheless, as soon as the set is dense enough RESTORE will converge anyway and (∞) the resulting function f j will depend on several different parameters: the choice of the quasi-interpolation operator Q Ψ ,X and the choice of the space V j . The former can be tuned according to different partitions of unity Ψ , the latter is parametrized by the particular basis Φ j . Therefore we can tune the choice of Q Ψ ,X and V j in order to optimize (∞) M , see for example certain geometrical properties of the resulting f j with respect to the given data set {(x` , y` )}`=0 Fig. 5.
Acknowledgement Massimo Fornasier acknowledges the financial support provided by the European Union’s Human Potential Programme under contract MOIF-CT-2006-039438.
M. Fornasier, L. Gori / Journal of Computational and Applied Mathematics 221 (2008) 376–385
385
References [1] A. Aldroubi, H.G. Feichtinger, Exact iterative reconstruction algorithm for multivariate irregularly sampled functions in spline like spaces: The L p theory, Proc. Amer. Math. Soc. 126 (9) (1998) 2677–2686. [2] A. Aldroubi, K. Gr¨ochenig, Nonuniform sampling and reconstruction in shift-invariant spaces, SIAM Rev. 43 (4) (2001) 585–620. [3] J.J. Benedetto, P.J.S.G. Ferreira, Modern Sampling Theory: Mathematics and Applications, Birkh¨auser, 2000. [4] Z. Ciesielski, T. Figiel, Spline bases in classical function spaces on compact C ∞ manifolds. I and II, Studia Math. 76 (1983) 1–58 and 96–136. [5] M. Cotronei, M.L. Lo Cascio, A method for the construction of families of biorthogonal filters with prescribed properties, Adv. Comput. Math. 17 (3) (2002) 199–210. [6] S. Dahlke, M. Fornasier, T. Raasch, Adaptive frame methods for elliptic operator equations, Adv. Comput. Math. 27 (1) (2007) 27–63. [7] W. Dahmen, R. Schneider, Wavelets with complementary boundary conditions — Function spaces on the cube, Results Math. 34 (3–4) (1998) 255–293. [8] W. Dahmen, R. Schneider, Wavelets on manifolds I. Construction and domain decomposition, SIAM J. Math. Anal. 31 (1999) 184–230. [9] H.G. Feichtinger, K. Gr¨ochenig, Iterative reconstruction of multivariate band-limited functions from irregular sampling values, SIAM J. Math. Anal. 1 (1992) 244–261. [10] H.G. Feichtinger, K. Gr¨ochenig, T. Strohmer, Efficient numerical methods in non-uniform sampling theory, Numer. Math. 69 (1995) 423–440. [11] H.G. Feichtinger, T. Strohmer, Recovery missing segments and lines in images, in: Digital Image Recovery and Synthesis, Opt. Eng. 33 (10) (1994) 3283–3289 (special issue). [12] L. Gori, L. Pezza, F. Pitolli, Recent results on wavelet bases on the interval generated by GP refinable functions, Appl. Numer. Math. 51 (2004) 549–563. [13] L. Gori, F. Pitolli, A class of totally positive refinable functions, Rend. Math. 7 (20) (2000) 305–322. [14] K. Gr¨ochenig, H. Schwab, Fast local reconstruction methods for nonuniform sampling in shift invariant spaces, SIAM J. Matrix Anal. App. 24 (42) (2004) 899–913.