 # Anisotropic functional deconvolution with long-memory noise: the case of a multi-parameter fractional Wiener sheet

We look into the minimax results for the anisotropic two-dimensional functional deconvolution model with the two-parameter fractional Gaussian noise. We derive the lower bounds for the L^p-risk, 1 ≤ p < ∞, and taking advantage of the Riesz poly-potential, we apply a wavelet-vaguelette expansion to de-correlate the anisotropic fractional Gaussian noise. We construct an adaptive wavelet hard-thresholding estimator that attains asymptotically quasi-optimal convergence rates in a wide range of Besov balls. Such convergence rates depend on a delicate balance between the parameters of the Besov balls, the degree of ill-posedness of the convolution operator and the parameters of the fractional Gaussian noise. A limited simulations study confirms theoretical claims of the paper. The proposed approach is extended to the general r-dimensional case, with r> 2, and the corresponding convergence rates do not suffer from the curse of dimensionality.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction.

Consider the problem of estimating a periodic two-dimensional function, , based on noisy convolutions that are described by the equations

 dY(t,x)=q(t,x)dtdx+ε¯¯¯αdBα(t,x), q(t,x)=∫10f(s,x)g(t−s,x)ds. (1)

Here, , is the convolution kernel and it is assumed to be known, is an anisotropic two-dimensional fractional Brownian sheet (fBs), , and , , are the parameters of the long-memory in the direction of and , respectively. A two-dimensional fractional Brownian sheet is defined by the formula

 Bα(t,x)=1CH∫t0∫x0Q(s,y)dW(s,y), (2)

where is a two-dimensional standard Brownian sheet, is some explicit constant and

 Q(s,y)=[(t−s)+]H1−1/2[(x−y)+]H2−1/2, (3)

with . In addition, the two-parameter fBs is characterized by a covariance function of the form

 Cov (Bα(t1,s1),Bα(t2,s2))=c[|t1|2H1+|t2|2H1−|t1−t2|2H1][|s1|2H2+|s2|2H2−|s1−s2|2H2], (4)

for some , , , where are the Hurst parameters in the direction of and , respectively (see Kamont (1996)).
The discrete, and the more realistic, version of model (1) is given by

 Y(ti,xl)=q(ti,xl)+σξil, (i,l)=1,2,⋯,N, (5)

where , ,

is a positive variance constant, and

are zero-mean second-order stationary Gaussian random variables satisfying the auto-covariance structure

 γ(h1,0)≍h−α11,  γ(0,h2)≍h−α22 and  γ(h1,h2)≍h−α11h−α22. (6)

Assumption A.1. The error structure in model (5) is a zero-mean second-order stationary process satisfying

 XND−→Bα, (7)

where

 XN(t,x)=1N2−α1/2−α2/2⌈Nt⌉∑i=1⌈Nx⌉∑l=1ξil, (t,x)∈U2, (8)

and is a two-parameter fractional Brownian sheet.
Assumption A.1 is valid under auto-covariance structure (6) and it appeared in Adu and Richardson (2018). It guarantees that, with the calibration , models (1) and (5) are asymptotically equivalent. Therefore, for sufficiently large sample size , model (1) can be used to approximate model (5).

Deconvolution model has been the subject of a great deal of papers since late 1980s, but the most significant contribution was that of Donoho (1995) who was the first to devise a wavelet solution to the problem. Other attempts include, Abramovich and Silverman (1998), Walter and Shen (1999), Johnstone et al. (2004), Donoho and Raimondo (2004), among others. In the case of functional deconvolution model with , Pensky and sapatinas (2009, 2010, 2011) pioneered into the formulation and further development of the problem.

Functional deconvolution problem of type (1) with

corresponds to the white noise case and it was investigated under the

-risk in Benhaddou et al. (2013), and under -risk, , in Benhaddou (2017), where they constructed an adaptive hard-thresholding wavelet estimator, and showed that it is asymptotically quasi-optimal within a logarithmic factor of over a wide range of Besov balls of mixed regularity. This model is motivated by experiments in which one needs to recover a two-dimensional function using observations of its convolutions along profiles . This situation occurs, for example, in seismic inversions (see Robinson (1999)). In these articles, it is assumed that the error terms are white noise processes or i.i.d noise. However, empirical evidence has shown that, even at large lags, the correlation structure in the errors takes the power-like form (6). This phenomenon is referred to as long-memory (LM) or long-range dependence (LRD). The presence of long memory in oceanic seismic data was pointed out in Wood et al. (2014) for instance, and it may be due for example to sea floor temperature variations.

Long-memory has been investigated quite considerably in the standard deconvolution model. One can list; Wang (1996, 1997), Wishart (2013), Benhaddou et al. (2014) and Kulik et al. (2015). In a few other relevant contexts, LM was also investigated in density deconvolution in Comte et al. (2008), and in the Laplace deconvolution in Benhaddou (2018) where the unknown response function is non-periodic and defined on the entire positive real half-line.

The objective of the paper is to extend the work of Benhaddou et al. (2013) to the case when the noise is an anisotropic multi-parameter fractional Brownian sheet. Following a standard procedure, we derive minimax lower bounds for the -risk, with , when belongs to an anisotropic Besov ball and the blurring function is regular smooth. Taking advantage of the Riesz poly-potential operator, and inspired by the work of Wang (1996, 1997), we apply the wavelet-vaguelette expansion (WVD) to de-correlate the anisotropic fBs in two directions. In addition, we take advantage of the flexibility of the Meyer wavelet basis in the Fourier domain to construct a wavelet hard-thresholding estimator that is adaptive and asymptotically quasi-optimal within a logarithmic factor of over a large array of Besov balls. Furthermore, the proposed estimator attains convergence rates that depend on a delicate balance between the parameters of the Besov ball, the degree of ill-posedness of the convolution, and the parameters and associated with the fBs. Similar to the white noise case studied in Benhaddou et al. (2013), the proposed approach is easily extended to recovering an -dimensional function, , and the corresponding convergence rates turn out to be dimension-free. Finally, a simulation study is carried out to further confirm the results of our theoretical findings.

The rest of the paper is organized as follows. Section 2 introduces some notation as well as the estimation algorithm. Section 3 describes the derivation of the lower bounds for the -risk of estimators of as well as the upper bounds and establishes the asymptotic optimality of the estimator. Section 4 presents a limited simulations study to complement the theoretical results from Section 3. Section 5 extends the results in Sections 2 and 3 to the general -dimensional case. Finally, Section 6 contains the proofs of the theoretical results.

## 2 Estimation Algorithm.

In what follows, denote the complex conjugate of by . Let , , and be the two-dimensional Fourier coefficients of functions , , , and , respectively.
Consider a bandlimited wavelet basis (e.g., Meyer-type) on , and form the product wavelet basis

 Ψω(t,x)=ψj1,k1(t)ψj2,k2(x), (t,x)∈U, (9)

where , and

 Ω={ω=(j1,k1;j2,k2):ji=m0−1,⋯,∞;ki=0,1,⋯,2ji−1,i=1,2}. (10)

Let be the lowest resolution level for the Meyer basis and denote the scaling function for the wavelet by , . Since the functions (9) form an orthonormal basis of the -space, the function can be expanded over these bases with coefficients into a wavelet series

 f(t,x)=∑ω∈ΩβωΨω(t,x), (11)

where

. Denote the two-dimensional Fourier transform of

by . It is well-known (see, e.g, Johnstone et al (2004), section 3.1) that under the Fourier domain and for any , , one has

 Wji={mi:~ψji,ki(mi)≠0}⊆2π/3[−2ji+2,−2ji]∪[2ji,2ji+2]. (12)

Apply the two-dimensional Fourier transform to equation (1) to obtain

 ~Y(m1,m2)=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯~f(m1,m2)~g(m1,m2)+ε¯¯¯α~Bα(m1,m2). (13)

Now, applying Plancherel formula in both directions, the wavelet coefficients in (11) can be rewritten as

 βω=∑m1∈Wj1∑m2∈Wj2¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯~ψj1,k1(m1)~ψj2,k2(m2)~f(m1,m2). (14)

Therefore, an unbiased estimator for

is given by

 ~βω=∑m1∈Wj1∑m2∈Wj2¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯~ψj1,k1(m1)~ψj2,k2(m2)~Y(m1,m2)~g(m1,m2). (15)

Finally, consider the hard thresholding estimator for

 ˆfε(t,x)=∑ω∈Ω(J1,J2)~βωI(|~βω|>λαj;ε)Ψω(t,x), (16)

where

 Ω(J1,J2)={ω=(j1,k1;j2,k2):ji=m0−1,⋯,Ji;ki=0,1,⋯,2ji−1,i=1,2}, (17)

and the quantities , and are to be specified.
Assumption A.2. In the Fourier domain, the convolution kernel , for some positive constants , and and , independent of and , is such that

 C1|m1|−2ν≤|~g(m1,m2)|2≤C2|m1|−2ν. (18)

Let us now evaluate the variance of (15).

###### Lemma 1

Let be defined in (15). Then, under the condition (18), for , one has

 E∣∣~βω−βω∣∣2p≤K1ε2p¯¯¯α2j1p(2ν+α1−1)+j2p(α2−1). (19)

 E∣∣~βω−βω∣∣p≤K2εp¯¯¯α2j1p2(2ν+α1−1)+j2p2(α2−1). (20)

where and are some positive constants independent of .

According to Lemma 1, choose the thresholds of the form

 λαj;ε=γε¯¯¯α√|ln(ε)|2j12(2ν+α1−1)2j22(α2−1), (21)

where is some positive constant independent of . Furthermore, the finest resolution levels and are such that

 2J2=[ε2¯¯¯αA2]−1α2,   2J1=[ε2¯¯¯αA2]−12ν+α1. (22)

## 3 Convergence rates and asymptotic optimality.

Denote

 s∗i = si+1/2−1/π, (23) s′i = si+1/2−1/π′, (24) s′′i = si+1/p−1/p′, (25)

where and .
Assumption A.3. The function belongs to an anisotropic two-dimensional Besov space. In particular, its wavelet coefficients satisfy

 Bs1,s2π,q(A)=⎧⎪ ⎪⎨⎪ ⎪⎩f∈L2(U):⎛⎜⎝∑j1,j22(j1s∗1+j2s∗2)q⎛⎝∑k1,k2|βj1,k1,j2,k2|π⎞⎠q/π⎞⎟⎠1/q≤A⎫⎪ ⎪⎬⎪ ⎪⎭ (26)

To construct minimax lower bounds for the -risk, we define the -risk over the set as

 R(V)=inf~fsupf∈VE∥~fε−f∥pp, (27)

where is the -norm of a function and the infimum is taken over all possible estimators of .

###### Theorem 1

Let with , and . Then, under conditions (18) and (26), for , as ,

 R(Bs1,s2π,q(A))≥CAp⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩[ε2¯¯¯αA2]ps22s2+α2,if  s1>s2α2(2ν+α1), & s2α2≥p2(1p′−1p)[ε2¯¯¯αA2]ps12s1+2ν+α1,if   p2(2ν+α1)(1π−1p)≤s1≤s2α2(2ν+α1),[ε2¯¯¯αA2]p(s1+1p−1π)2s∗1+2ν+α1−1,if s1<(2ν+α1)p2(1π−1p), & s2α2≥p2(1p′−1p). (28)

Proof of Theorem 1. In order to prove the theorem, first we establish the equivalence between a process that involves the fractional Wiener sheet (eq. (1)) and another that involves the standard Wiener sheet by applying an appropriate integral operator, taking advantage of relation (2). Then, we consider two cases, the case when is dense in both dimensions (the dense-dense case) and the case when is dense in but sparse in . Lemma of Bunea et al. (2007) is then applied to find such lower bounds using conditions (18) and (26), combined with some useful properties of Meyer wavelet basis. To complete the proof, we choose the highest of the lower bounds.
The derivation of upper bounds of the -risk relies on the following two lemmas.

###### Lemma 2

Under condition (26), one has

 2j1−1∑k1=02j2−1∑k2=0∣βj1,k1,j2,k2∣p≤Ap2−p[(j1s1+j2s2)+(12−1p′)(j1+j2)], ∀j1,j2≥0. (29)
###### Lemma 3

Let and be defined by (15) and (21), respectively. Define, for some positive constant , the set

 Θω,η={Θ:|~βω−βω|>ηλαj;ε}. (30)

Then, under condition (18) and as , one has

 Pr(Θω,η)=O⎛⎜ ⎜⎝[ε2¯¯¯α]τ|ln(ε)|12⎞⎟ ⎟⎠, (31)

where , and appear in (18) and (21), respectively, and .

###### Theorem 2

Let be the wavelet estimator in (16), with and given by (22). Let with , and let conditions (18) and (26) hold. If in (21) is large enough, then, for and as , one has

 R(Bs1,s2π,q(A))≤CAp⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩[ε2¯¯¯αA2]ps22s2+α2[|ln(ε)|]ξ1+ps22s2+α2,if  s1≥s2α2(2ν+α1) & s2α2≥p2(1p′−1p),[ε2¯¯¯αA2]ps12s1+2ν+α1[|ln(ε)|]ps12s1+2ν+α1,if  p2(2ν+α1)(1π−1p)

where and are defined as

 ξ1 = I(s1=s2α2(2ν+α1)), (33) ξ2 = I(s1=(2ν+α1)p2(1π−1p))+I(s2α2=p2(1p′−1p)). (34)

The proof of Theorem 2. The proof is very similar to that of Theorem 2 in Benhaddou. (2017).

###### Remark 1

Notice that the finest resolution levels and in (22) and the thresholds in (21), are independent of the unknown parameters of the Besov ball (26) and therefore estimator (16) is adaptive with respect to those parameters.

###### Remark 2

Theorems 1 and 2 imply that, for the -risk, the estimator in (16) is asymptotically quasi-optimal within a logarithmic factor of , over a wide range of anisotropic Besov balls .

###### Remark 3

The convergence rates depend on a delicate balance between the parameters of the Besov ball, smoothness of the convolution kernel and the parameters of the anisotropic fractional Brownian sheet, and . These rates deteriorate as the values of and get closer and closer to zero.

###### Remark 4

To de-correlate the two-dimensional fractional Gaussian noise a wavelet-vaguelette expansion based on Meyer wavelet bases has been applied to the fractional Brownian sheet. The validity of a wavelet-based series representation for the fractional Brownian motion has been established in Wang (1997), Meyer, Sellan and Taqqu (1999) and Ayache and Taqqu (2002). In particular, Ayache and Taqqu (2002) established the optimality of wavelet-based series representation of the fractional Brownian sheets.

###### Remark 5

For , the rates of convergence match exactly those in Benhaddou (2017), and, with , those in Benhaddou et al. (2013) in their bivariate white noise case.

###### Remark 6

If we hold the second dimension fixed, our rates are comparable to those in Wang (1997) and Wishart (2013) in their standard deconvolution with the one-parameter fractional Gaussian noise case.

## 4 Simulation Study

In order to investigate the effect of the long-memory on the performance of our estimator, we carried out a limited simulation study. We implemented an estimation algorithm for the model (5) using a modification of WaveD method of Raimondo and Stewart (2007). We evaluated mean integrated square error (MISE) of the functional deconvolution estimator.

1. We generated the data using equation (5) with convolution kernel , and . In particular, we chose to be a LIDAR or Doppler signal over the grid with , and we chose over the grid with , and . The LIDAR and Doppler signals were generated from the waved package. All test functions were scaled to have a unit norm. Bear in mind that, though is a product of two univariate functions, the method does not “recognize” this and, therefore, cannot take advantage of this information. Also, notice that with our choice of , the degree of ill-posedness (DIP) of the convolution is DIP=0.5.

2. We simulated the LM error in the direction of both and using the fracdiff package of R available from CRAN. In particular, the fracdiff.sim command was used to simulate two one-dimensional fractionally differenced ARIMA (fARIMA) sequences which behave similar to a fractional Gaussian noise, and then multiply them together to generate a two-dimensional error structure. For the LM parameters , we used different combinations of the levels , , and . To create a dependence structure that is anisotropic, we only used combinations such that . Note that the fractional differencing parameter is obtained from by the relation .

3. The choice of in (5) was determined by the blurred signal-to-noise ratio (SNR), where We considered three choices, SNR=10dB (high noise), 20dB (medium noise) and 30dB (low noise).

4. To compute the estimator, we used the function FWaveD and IWaveD from the waved package. First, we applied the Meyer wavelet transform to the data in -direction by FWaveD, with no thresholding by setting thr

to be a zero vector. Then, we applied the second thresholded Meyer wavelet transform in

-direction by FWaveD, however, the level dependent threshold were modified according to (21). For the tuning parameter , we tried different values ranging from to , but the performance of the estimator reached its best at the latter value. Therefore, our default choice of was . The finest resolution level was estimated using the default method in the waved package, together with the choice .

Table 1 reports the averages of the errors over 2000 simulation runs using a sample size of . Figure (1) illustrates the function and its estimator and shows a relatively good precision that deteriorates as the anisotropic pair of LM parameters decreases from 0.8 to 0.2. Notice here that according to Table 1, the influence of LM is a little more pronounced in the direction of the dimension than it is in the direction of . Table 1 complements Figure (1) and confirms our theoretical results in the sense that as the level of LM increases ( gets smaller), the mean integrated squared error increases (the performance deteriorates). Figure 1: The function f(t,x) (top) and its estimate ˆf(t,x) (bottom) at different combinations of (α1,α2).

## 5 The r-dimensional case: estimation, convergence rates and asymptotic optimality.

Consider the -dimensional case of model (1)

 dY(t,x)=q(t,x)dtdx+ε¯¯¯αdBα(t,x), q(t,x)=∫10f(s,x)g(t−s,x)ds. (35)

Here, , , is the convolution kernel, is an anisotropic -dimensional fBs, , and , , is the parameter of the long-memory in the direction of the variable. An -parameter fBs is defined as a zero-mean Gaussian process whose covariance function is of the form

 Cov (Bα(t1),Bα(t2))=cΠrk=1[|tk1|2Hk+|tk2|2Hk−|tk1−tk2|2Hk], (36)

for some , , , where is the Hurst parameter in the direction of the dimension.
Consider a bandlimited wavelet basis (e.g., Meyer-type) on , and form the product wavelet basis

 Ψωr(t,x)=ψj1,k1(t)Πri=2ψji,ki(xi−1), (t,x)∈U, (37)

where , and

 Ωr={ωr=(j1,k1;j2,k2;⋯;jr,kr):ji=m0−1,⋯,∞;ki=0,1,⋯,2ji−1,i=1,2,⋯,r}. (38)

Then, the function can be expanded over these bases with coefficients into a wavelet series

 f(t,x)=∑ωr∈ΩrβωrΨωr(t,x), (39)

where . Apply the -dimensional Fourier transform to equation (35) to obtain

 ~Y(mr)=~f(mr)~g(mr)+ε¯¯¯α~Bα(mr), (40)

where . Then, applying Plancherel formula in all directions, the wavelet coefficients in (39) can be rewritten as

 βωr=∑m1∈Wj1∑m2∈Wj2⋯∑mr∈Wjr~Ψωr(mr)~f(mr). (41)

Therefore, an unbiased estimator for is given by

 ~βωr=∑m1∈Wj1∑m2∈Wj2⋯∑mr∈Wjr~Ψωr(mr)~Y(mr)~g(mr). (42)

Finally, consider the hard thresholding estimator for

 ˆfε(t,x)=∑ωr∈Ωr(J)~βωrI(|~βωr|>λαj;ε)Ψωr(t,x), (43)

where

 Ωr(J)={ωr=(j1,k1;j2,k2;⋯;jr,kr):ji=m0−1,⋯,Ji;ki=0,1,⋯,2ji−1,i=1,2,⋯,r}, (44)

with , and the quantities , , , and are to be determined.
Assumption A.4. In the Fourier domain, the convolution kernel , for some positive constants , and and , independent of , is such that

 Cr1|m1|−2ν≤|~g(mr)|2≤Cr2|m1|−2ν. (45)
###### Lemma 4

Let be defined in (42). Then, under the condition (45), for , one has

 E∣∣~βωr−βωr∣∣2p≤K3εp2¯¯¯α2j1p(2ν+α1−1)+j2p(α2−1)+⋯+jrp(αr−1). (46)