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

TL;DR

This paper develops an adaptive wavelet-based method for anisotropic two-dimensional deconvolution with long-memory fractional Gaussian noise, achieving near-optimal convergence rates and extending to higher dimensions without curse of dimensionality.

Contribution

It introduces a wavelet-vaguelette approach for anisotropic fractional Gaussian noise deconvolution, providing asymptotic optimality and extending results to higher dimensions.

Findings

The estimator attains asymptotically quasi-optimal convergence rates.

Convergence rates depend on Besov parameters, ill-posedness, and noise parameters.

Simulation studies support theoretical results.

Abstract

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 \leq p < \infty$, 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.

Figures

Tables

DIP=0.5		(${\alpha }_1$,${\alpha }_2$)
signal	$\sigma$	(0.6,0.8)	(0.4,0.8)	(0.2,0.8)	(0.4,0.6)	(0.2,0.6)	(0.2,0.4)
Lidar (10dB)	0.20	0.0760 (3)	0.0851 (3)	0.0956 (3)	0.1038 (3)	0.1293 (3)	0.2251 (3)
Lidar (20dB)	0.06	0.0423 (4)	0.0453 (4)	0.0486 (4)	0.0536 (4)	0.0612 (4)	0.0791 (4)
Lidar (30dB)	0.02	0.0120 (5)	0.0132 (5)	0.0149 (5)	0.0180 (5)	0.0210 (5)	0.0331 (5)
Doppler (10dB)	0.15	0.2601 (3)	0.2640 (3)	0.2749 (3)	0.2836 (3)	0.3090 (3)	0.3922 (3)
Doppler (20dB)	0.05	0.1247 (4)	0.1264 (4)	0.1292 (4)	0.1331 (4)	0.1376 (4)	0.1578 (4)
Doppler (30dB)	0.015	0.0513 (5)	0.0510 (5)	0.0522 (5)	0.0527 (5)	0.0541 (5)	0.0591 (5)

Equations

$$dY(t,x)=q(t,x)dtdx+{\epsilon }^{\overline{\alpha }}d{B}^{\alpha }(t,x),q(t,x)={\int }_0^{1}f(s,x)g(t-s,x)ds.$$

$$B^{\alpha }(t,x)=\frac{1}{C_H}{\int }_0^t{\int }_0^{x}Q(s,y)dW(s,y),$$

$$Q(s,y)={\left[(t-s{)}_{+}\right]}^{H_1-1/2}{\left[(x-y{)}_{+}\right]}^{H_2-1/2},$$

$$\text{mbox}Cov(B^{\alpha }(t_1,s_1),B^{\alpha }(t_2,s_2))=c\left[|t_1{|}^{2H_1}+|t_2{|}^{2H_1}-|t_1-t_2{|}^{2H_1}\right]\left[|s_1{|}^{2H_2}+|s_2{|}^{2H_2}-|s_1-s_2{|}^{2H_2}\right],$$

$$Y(t_i,x_l)=q(t_i,x_l)+\sigma {\xi }_{il},(i,l)=1,2,\cdots ,N,$$

$$\gamma (h_1,0)\asymp h_1^{-{\alpha }_1},\gamma (0,h_2)\asymp h_2^{-{\alpha }_2}and\gamma (h_1,h_2)\asymp h_1^{-{\alpha }_1}{h}_2^{-{\alpha }_2}.$$

$$X_N\stackrel{D}{\to }{B}^{\alpha },$$

$$X_N(t,x)=\frac{1}{N^{2-{\alpha }_1/2-{\alpha }_2/2}}\sum _{i=1}^{\lceil Nt\rceil }\sum _{l=1}^{\lceil Nx\rceil }{\xi }_{il},(t,x)\in U^2,$$

$${\Psi }_{\omega }(t,x)={\psi }_{j_1,k_1}(t){\psi }_{j_2,k_2}(x),(t,x)\in U,$$

$$\Omega =\left\{\omega =(j_1,k_1;j_2,k_2):j_i=m_0-1,\cdots ,\infty ;k_i=0,1,\cdots ,2^{j_i}-1,i=1,2\right\}.$$

$$f(t,x)=\sum _{\omega \in \Omega }{\beta }_{\omega }{\Psi }_{\omega }(t,x),$$

$$W_{j_i}=\left\{m_i:{\tilde{\psi }}_{j_i,k_i}(m_i)\ne 0\right\}\subseteq 2\pi /3\left[-2^{j_i+2},-2^{j_i}\right]\cup \left[2^{j_i},2^{j_i+2}\right].$$

$$\tilde{Y}(m_1,m_2)=\overline{\tilde{f}(m_1,m_2)}\tilde{g}(m_1,m_2)+{\epsilon }^{\overline{\alpha }}{\tilde{B}}^{\alpha }(m_1,m_2).$$

$${\beta }_{\omega }=\sum _{m_1\in W_{j_1}}\sum _{m_2\in W_{j_2}}\overline{{\tilde{\psi }}_{j_1,k_1}(m_1)}{\tilde{\psi }}_{j_2,k_2}(m_2)\tilde{f}(m_1,m_2).$$

$${\tilde{\beta }}_{\omega }=\sum _{m_1\in W_{j_1}}\sum _{m_2\in W_{j_2}}\overline{{\tilde{\psi }}_{j_1,k_1}(m_1)}{\tilde{\psi }}_{j_2,k_2}(m_2)\frac{\tilde{Y}(m_1,m_2)}{\tilde{g}(m_1,m_2)}.$$

$${\hat{f}}_{\epsilon }(t,x)=\sum _{\omega \in \Omega (J_1,J_2)}{\tilde{\beta }}_{\omega }\mathbb{I}\left(|{\tilde{\beta }}_{\omega }|>{\lambda }_{j;\epsilon }^{\alpha }\right){\Psi }_{\omega }(t,x),$$

$$\Omega (J_1,J_2)=\left\{\omega =(j_1,k_1;j_2,k_2):j_i=m_0-1,\cdots ,J_i;k_i=0,1,\cdots ,2^{j_i}-1,i=1,2\right\},$$

$$C_1|m_1{|}^{-2\nu }\le |\tilde{g}(m_1,m_2){|}^2\le C_2|m_1{|}^{-2\nu }.$$

$$\mathbb{E}{\left|{\tilde{\beta }}_{\omega }-{\beta }_{\omega }\right|}^{2p}\le K_1{\epsilon }^{2p\overline{\alpha }}{2}^{j_{1}p(2\nu +{\alpha }_1-1)+j_{2}p({\alpha }_2-1)}.$$

$$\mathbb{E}{\left|{\tilde{\beta }}_{\omega }-{\beta }_{\omega }\right|}^p\le K_2{\epsilon }^{p\overline{\alpha }}{2}^{j_1\frac{p}{2}(2\nu +{\alpha }_1-1)+j_2\frac{p}{2}({\alpha }_2-1)}.$$

$${\lambda }_{j;\epsilon }^{\alpha }=\gamma {\epsilon }^{\overline{\alpha }}\sqrt{|\ln (\epsilon )|}{2}^{\frac{j_1}{2}(2\nu +{\alpha }_1-1)}{2}^{\frac{j_2}{2}({\alpha }_2-1)},$$

$$2^{J_2}={\left[\frac{{\epsilon }^{2\overline{\alpha }}}{A^2}\right]}^{-\frac{1}{{\alpha }_2}},2^{J_1}={\left[\frac{{\epsilon }^{2\overline{\alpha }}}{A^2}\right]}^{-\frac{1}{2\nu +{\alpha }_1}}.$$

$$s_i^{\ast }$$

$$s_i^{\prime }$$

$$s_i^{\prime \prime }$$

$$B_{\pi ,q}^{s_1,s_2}(A)=\left\{f\in L^2(U):{\left(\sum _{j_1,j_2}{2}^{(j_1{s}_1^{\ast }+j_2{s}_2^{\ast })q}{\left(\sum _{k_1,k_2}|{\beta }_{j_1,k_1,j_2,k_2}{|}^{\pi }\right)}^{q/\pi }\right)}^{1/q}\le A\right\}$$

$$R(V)=\underset{\tilde{f}}{\inf }\underset{f\in V}{\sup }\mathbb{E}\Vert {\tilde{f}}_{\epsilon }-f{\Vert }_p^p,$$

$$R(B^{s_{1},s_{2}}_{\pi,q}(A))\geq CA^{p}\left\{\begin{array}[]{ll}\left[\frac{\varepsilon^{2\overline{\alpha}}}{A^{2}}\right]^{\frac{ps_{2}}{2s_{2}+\alpha_{2}}},&\mbox{if}\ \ s_{1}>\frac{s_{2}}{\alpha_{2}}(2\nu+\alpha_{1}),\ \&\ \frac{s_{2}}{\alpha_{2}}\geq\frac{p}{2}(\frac{1}{p^{\prime}}-\frac{1}{p})\\ \left[\frac{\varepsilon^{2\overline{\alpha}}}{A^{2}}\right]^{\frac{ps_{1}}{2s_{1}+2\nu+\alpha_{1}}},&\mbox{if}\ \ \ \frac{p}{2}(2\nu+\alpha_{1})(\frac{1}{\pi}-\frac{1}{p})\leq s_{1}\leq\frac{s_{2}}{\alpha_{2}}(2\nu+\alpha_{1}),\\ \left[\frac{\varepsilon^{2\overline{\alpha}}}{A^{2}}\right]^{\frac{p(s_{1}+\frac{1}{p}-\frac{1}{\pi})}{2s^{*}_{1}+2\nu+\alpha_{1}-1}},&\mbox{if}\ s_{1}<(2\nu+\alpha_{1})\frac{p}{2}(\frac{1}{\pi}-\frac{1}{p}),\ \&\ \frac{s_{2}}{\alpha_{2}}\geq\frac{p}{2}(\frac{1}{p^{\prime}}-\frac{1}{p}).\end{array}\right.$$

$$\sum _{k_1=0}^{2^{j_1}-1}\sum _{k_2=0}^{2^{j_2}-1}\mid {\beta }_{j_1,k_1,j_2,k_2}{\mid }^p\le A^p{2}^{-p[(j_1{s}_1+j_2{s}_2)+(\frac{1}{2}-\frac{1}{p^{\prime }})(j_1+j_2)]},\forall j_1,j_2\ge 0.$$

$${\Theta }_{\omega ,\eta }=\left\{\Theta :|{\tilde{\beta }}_{\omega }-{\beta }_{\omega }|>\eta {\lambda }_{j;\epsilon }^{\alpha }\right\}.$$

$$\mathrm{Pr}\left({\Theta }_{\omega ,\eta }\right)=O\left(\frac{{\left[{\epsilon }^{2\overline{\alpha }}\right]}^{\tau }}{|\ln (\epsilon ){|}^{\frac{1}{2}}}\right),$$

Full text

**Keywords and phrases: Anisotropic functional deconvolution, Besov space, anisotropic fractional Brownian sheet, minimax convergence rate

**

**AMS (2000) Subject Classification: 62G05, 62G20, 62G08 **

1 Introduction.

Consider the problem of estimating a periodic two-dimensional function, $f(\cdot,\cdot)\in L^{2}(U)$, based on noisy convolutions that are described by the equations

[TABLE]

Here, $U=[0,1]^{2}$, $g(\cdot,\cdot)$ is the convolution kernel and it is assumed to be known, $B^{{\alpha}}(t,x)$ is an anisotropic two-dimensional fractional Brownian sheet (fBs), $\overline{\alpha}={(\alpha_{1}+\alpha_{2})}/{2}$, and $\alpha_{i}=2-2H_{i}\in(0,1]$, $i=1,2$, are the parameters of the long-memory in the direction of $t$ and $x$, respectively. A two-dimensional fractional Brownian sheet is defined by the formula

[TABLE]

where $W=\left\{W(t,x),(t,x)\in U\right\}$ is a two-dimensional standard Brownian sheet, $C_{H}$ is some explicit constant and

[TABLE]

with $Y_{+}=\max(0,Y)$. In addition, the two-parameter fBs is characterized by a covariance function of the form

[TABLE]

for some $c\in R$, $(t_{i},s_{i})\in U$, $i=1,2$, where $H_{i}=1-\frac{\alpha_{i}}{2}$ are the Hurst parameters in the direction of $t$ and $s$, respectively (see Kamont (1996)).

The discrete, and the more realistic, version of model (1) is given by

[TABLE]

where $t_{i}=\frac{i}{N}$, $x_{l}=\frac{l}{N}$, $\sigma$ is a positive variance constant, and $\xi_{il}$ are zero-mean second-order stationary Gaussian random variables satisfying the auto-covariance structure

[TABLE]

**Assumption A.1. ** The error structure $\{\xi_{il}:i,l\in Z\}$ in model (5) is a zero-mean second-order stationary process satisfying

[TABLE]

where

[TABLE]

and $B^{{\alpha}}$ 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 $\varepsilon^{{\alpha_{1}+\alpha_{2}}}\asymp n^{\alpha_{1}/2+\alpha_{2}/2}$, models (1) and (5) are asymptotically equivalent. Therefore, for sufficiently large sample size $n=N^{2}$, 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 $f(t,x)\equiv f(t)$, Pensky and sapatinas (2009, 2010, 2011) pioneered into the formulation and further development of the problem.

Functional deconvolution problem of type (1) with $\alpha_{1}=\alpha_{2}=1$ corresponds to the white noise case and it was investigated under the $L^{2}$-risk in Benhaddou et al. (2013), and under $L^{p}$-risk, $1\leq p<\infty$, 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 $\varepsilon$ 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 $x=x_{l}$. 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 $f(\cdot)$ 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 $L^{p}$-risk, with $1\leq p<\infty$, when $f(t,x)$ belongs to an anisotropic Besov ball and the blurring function $g(t,x)$ 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 $\varepsilon$ 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 $\alpha_{1}$ and $\alpha_{2}$ 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 $r$-dimensional function, $r>2$, 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 $L^{p}$-risk of estimators of $f$ 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 $r$-dimensional case. Finally, Section 6 contains the proofs of the theoretical results.

2 Estimation Algorithm.

In what follows, denote the complex conjugate of $a$ by $\bar{a}$. Let $\tilde{Y}(m_{1},m_{2})$, $\tilde{B^{\alpha}}(m_{1},m_{2})$, $\tilde{g}(m_{1},m_{2})$ and $\tilde{f}(m_{1},m_{2})$ be the two-dimensional Fourier coefficients of functions $dY(t,x)$, $dB^{\alpha}(t,x)$, $g(t,x)$, and $f(t,x)$, respectively.

Consider a bandlimited wavelet basis $\psi_{j,k}(y)$ (e.g., Meyer-type) on $[0,1]$, and form the product wavelet basis

[TABLE]

where $\omega\in\Omega$, and

[TABLE]

Let $m_{0}$ be the lowest resolution level for the Meyer basis and denote the scaling function for the wavelet by $\psi_{m_{0}-1,k_{i}}$, $i=1,2$. Since the functions (9) form an orthonormal basis of the $L^{2}(U)$-space, the function $f(t,x)$ can be expanded over these bases with coefficients $\beta_{\omega}$ into a wavelet series

[TABLE]

where $\beta_{\omega}=\int^{1}_{0}\int^{1}_{0}f(t,x)\Psi_{\omega}(t,x)dtdx$. Denote the two-dimensional Fourier transform of $\Psi_{\omega}(t,x)$ by $\tilde{\Psi}_{\omega}(m_{1},m_{2})=\tilde{\psi}_{j_{1},k_{1}}(m_{1})\tilde{\psi}_{j_{2},k_{2}}(m_{2})$. It is well-known (see, e.g, Johnstone et al (2004), section 3.1) that under the Fourier domain and for any $j_{i}\geq m_{0}$, $i=1,2$, one has

[TABLE]

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

[TABLE]

Now, applying Plancherel formula in both directions, the wavelet coefficients $\beta_{\omega}$ in (11) can be rewritten as

[TABLE]

Therefore, an unbiased estimator for $\beta_{\omega}$ is given by

[TABLE]

Finally, consider the hard thresholding estimator for $f(t,x)$

[TABLE]

where

[TABLE]

and the quantities $J_{1}$, $J_{2}$ and $\lambda^{\alpha}_{j;\varepsilon}$ are to be specified.

**Assumption A.2. ** In the Fourier domain, the convolution kernel $g(t,x)$, for some positive constants $\nu$, and $C_{1}$ and $C_{2}$, independent of $m_{1}$ and $m_{2}$, is such that

[TABLE]

Let us now evaluate the variance of (15).

Lemma 1

Let $\tilde{\beta}_{\omega}$ be defined in (15). Then, under the condition (18), for $1\leq p<\infty$, one has

[TABLE]

In addition,

[TABLE]

where $K_{1}$ and $K_{2}$ are some positive constants independent of $\varepsilon$.

According to Lemma 1, choose the thresholds $\lambda^{\alpha}_{j;\varepsilon}$ of the form

[TABLE]

where $\gamma$ is some positive constant independent of $\varepsilon$. Furthermore, the finest resolution levels $J_{1}$ and $J_{2}$ are such that

[TABLE]

3 Convergence rates and asymptotic optimality.

Denote

[TABLE]

where $\pi^{\prime}=\min\{2,\pi\}$ and $p^{\prime}=\min\{p,\pi\}$.

Assumption A.3. The function $f(t,x)$ belongs to an anisotropic two-dimensional Besov space. In particular, its wavelet coefficients $\beta_{j_{1},k_{1},j_{2},k_{2}}$ satisfy

[TABLE]

To construct minimax lower bounds for the $L^{p}$-risk, we define the $L^{p}$-risk over the set $V$ as

[TABLE]

where $\|g\|_{p}$ is the $L^{p}$-norm of a function $g$ and the infimum is taken over all possible estimators $\tilde{f}$ of $f$.

Theorem 1

Let $\min\{s_{1},s_{2}\}\geq\max\{\frac{1}{\pi},\frac{1}{2}\}$ with $1\leq\pi,q\leq\infty$, and $A>0$. Then, under conditions (18) and (26), for $1\leq p<\infty$, as $\varepsilon\rightarrow 0$,

[TABLE]

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 $f(t,x)$ is dense in both dimensions (the dense-dense case) and the case when $f(t,x)$ is dense in $x$ but sparse in $t$. Lemma $A.1$ 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. $\Box$

The derivation of upper bounds of the $L^{p}$-risk relies on the following two lemmas.

Lemma 2

Under condition (26), one has

[TABLE]

Lemma 3

Let $\tilde{\beta}_{\omega}$ and $\lambda^{\alpha}_{j;\varepsilon}$ be defined by (15) and (21), respectively. Define, for some positive constant $\eta$, the set

[TABLE]

Then, under condition (18) and as $\varepsilon\rightarrow 0$, one has

[TABLE]

where $\tau=\frac{\eta^{2}\gamma^{2}C^{2}_{1}}{8\overline{\alpha}}\left(\frac{3}{8\pi}\right)^{2\nu+\alpha_{1}-1+\alpha_{2}-1}$, $C_{1}$ and $\gamma$ appear in (18) and (21), respectively, and $\overline{\alpha}=\frac{\alpha_{1}+\alpha_{2}}{2}$.

Theorem 2

Let $\widehat{f}(.,.)$ be the wavelet estimator in (16), with $J_{1}$ and $J_{2}$ given by (22). Let $\min\{s_{1},s_{2}\}\geq\max\{\frac{1}{\pi},\frac{1}{2}\}$ with $1\leq\pi,q\leq\infty$, and let conditions (18) and (26) hold. If $\gamma$ in (21) is large enough, then, for $1\leq p\leq\infty$ and as $\varepsilon\rightarrow 0$, one has

[TABLE]

where $\xi_{1}$ and $\xi_{2}$ are defined as

[TABLE]

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

Remark 1

Notice that the finest resolution levels $J_{1}$ and $J_{2}$ in (22) and the thresholds $\lambda^{\alpha}_{j;\varepsilon}$ 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 $L^{p}$-risk, the estimator in (16) is asymptotically quasi-optimal within a logarithmic factor of $\varepsilon$, over a wide range of anisotropic Besov balls $B^{s_{1},s_{2}}_{\pi,q}(A)$.

Remark 3

The convergence rates depend on a delicate balance between the parameters of the Besov ball, smoothness of the convolution kernel $\nu$ and the parameters of the anisotropic fractional Brownian sheet, $\alpha_{1}$ and $\alpha_{2}$. These rates deteriorate as the values of $\alpha_{1}$ and $\alpha_{2}$ 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 $\alpha_{1}=\alpha_{2}=1$, the rates of convergence match exactly those in Benhaddou (2017), and, with $p=2$, those in Benhaddou et al. (2013) in their bivariate white noise case.

Remark 6

If we hold the second dimension $x$ 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) ${\mathbb{E}}\|\hat{f}-f\|^{2}$ of the functional deconvolution estimator.

We generated the data using equation (5) with convolution kernel $g(t,x)=0.5\exp(-|t|(1+(x-0.5)^{2}))$, and $f(t,x)=f(t)f(x)$. In particular, we chose $f(t)$ to be a LIDAR or Doppler signal over the grid $t_{i}=\frac{i}{N}$ with $i=1,2,...,N$, and we chose $f(x)=\exp(-|x-0.5|x^{3})$ over the grid $x_{i}=\frac{i}{N}$ with $i=1,2,...,N$, and $N=2^{10}=1024$. 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 $f(t,x)$ 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 $g(t,x)$, the degree of ill-posedness (DIP) of the convolution is DIP=0.5. 2. 2.

We simulated the LM error in the direction of both $t$ and $x$ 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$(0,d,0)$ (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 $\alpha=(\alpha_{1},\alpha_{2})$, we used different combinations of the levels $0.8$, $0.6$, $0.4$ and $0.2$. To create a dependence structure that is anisotropic, we only used combinations such that $\alpha_{1}\neq\alpha_{2}$. Note that the fractional differencing parameter $d$ is obtained from $\alpha$ by the relation $d=\frac{1-\alpha}{2}$ . 3. 3.

The choice of $\sigma$ in (5) was determined by the blurred signal-to-noise ratio (SNR), where $\text{SNR}=10\log_{10}\left(\frac{\|q(t,x)\|^{2}}{\sigma^{2}}\right).$ We considered three choices, SNR=10dB (high noise), 20dB (medium noise) and 30dB (low noise). 4. 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 $t$-direction by FWaveD, with no thresholding by setting thr to be a zero vector. Then, we applied the second thresholded Meyer wavelet transform in $x$-direction by FWaveD, however, the level dependent threshold were modified according to (21). For the tuning parameter $\gamma$, we tried different values ranging from $\gamma=\sqrt{2}$ to $\gamma=\sqrt{6}$, but the performance of the estimator reached its best at the latter value. Therefore, our default choice of $\gamma$ was $\sqrt{6}$. The finest resolution level $J_{1}$ was estimated using the default method in the waved package, together with the choice $J_{2}=5$.

Table 1 reports the averages of the errors over 2000 simulation runs using a sample size of $n=N^{2}=(2^{10})^{2}$. Figure (1) illustrates the function $f(t.x)$ and its estimator $\widehat{f}(t,x)$ and shows a relatively good precision that deteriorates as the anisotropic pair of LM parameters $(\alpha_{1},\alpha_{2})$ 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 $x$ than it is in the direction of $t$. Table 1 complements Figure (1) and confirms our theoretical results in the sense that as the level of LM increases ($\alpha$ gets smaller), the mean integrated squared error increases (the performance deteriorates).

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

Consider the $r$-dimensional case of model (1)

[TABLE]

Here, $U=[0,1]^{r}$, ${\bf x}=(x_{1},x_{2},\cdots,x_{r-1})$, $g(t,{\bf x})$ is the convolution kernel, $B^{{\alpha}}(t,{\bf x})$ is an anisotropic $r$-dimensional fBs, $\overline{\alpha}={\sum^{r}_{k=1}\alpha_{k}}/{r}$, and $\alpha_{k}=2-2H_{k}\in(0,1]$, $k=1,2,\cdots,r$, is the parameter of the long-memory in the direction of the $k^{th}$ variable. An $r$-parameter fBs is defined as a zero-mean Gaussian process $\{B^{{\alpha}}({\bf t}):{\bf t}\in[0,1]^{r}\}$ whose covariance function is of the form

[TABLE]

for some $c\in\mathbb{R}$, $(t_{k1},t_{k2})\in[0,1]^{2}$, $k=1,2,\cdots,r$, where $H_{k}=1-\alpha_{k}/2$ is the Hurst parameter in the direction of the $k^{th}$ dimension.

Consider a bandlimited wavelet basis $\psi_{j,k}(y)$ (e.g., Meyer-type) on $[0,1]$, and form the product wavelet basis

[TABLE]

where $\omega_{r}\in\Omega_{r}$, and

[TABLE]

Then, the function $f(t,{\bf x})$ can be expanded over these bases with coefficients $\beta_{\omega_{r}}$ into a wavelet series

[TABLE]

where $\beta_{\omega_{r}}=\int_{U}f(t,{\bf x})\Psi_{\omega_{r}}(t,{\bf x})dtd{\bf x}$. Apply the $r$-dimensional Fourier transform to equation (35) to obtain

[TABLE]

where ${\bf m_{r}}=(m_{1},m_{2},\cdots,m_{r})$. Then, applying Plancherel formula in all directions, the wavelet coefficients $\beta_{\omega_{r}}$ in (39) can be rewritten as

[TABLE]

Therefore, an unbiased estimator for $\beta_{\omega_{r}}$ is given by

[TABLE]

Finally, consider the hard thresholding estimator for $f(t,{\bf x})$

[TABLE]

where

[TABLE]

with ${\bf j}=(j_{1},j_{2},\cdots,j_{r})^{T}$, and the quantities $J_{1}$, $J_{2}$, $\cdots$, $J_{r}$ and $\lambda^{\alpha}_{{\bf j};\varepsilon}$ are to be determined.

**Assumption A.4. ** In the Fourier domain, the convolution kernel $g(t,{\bf x})$, for some positive constants $\nu$, and $C_{r1}$ and $C_{r2}$, independent of ${\bf m_{r}}=(m_{1},m_{2},\cdots,m_{r})$, is such that

[TABLE]

Lemma 4

Let $\tilde{\beta}_{\omega_{r}}$ be defined in (42). Then, under the condition (45), for $1\leq p<\infty$, one has

[TABLE]

In addition,

[TABLE]

where $K_{3}$ and $K_{4}$ are some positive constants independent of $\varepsilon$.

Consequently, choose the thresholds $\lambda^{\alpha}_{{\bf j};\varepsilon}$ of the from

[TABLE]

where $\rho$ is some positive constant that is independent of $\varepsilon$, $\overline{\alpha}=\frac{1}{r}{\sum^{r}_{k=1}\alpha_{k}}$ and $J_{1}$ and $J_{i}$, $i=2,3,\cdots,r$, are such that

[TABLE]

Lemma 5

Let $\tilde{\beta}_{\omega_{r}}$ and $\lambda^{\alpha}_{{\bf j};\varepsilon}$ be defined by (42) and (48), respectively. Define, for some positive constant $\eta$, the set

[TABLE]

Then, under condition (45), as $\varepsilon\rightarrow 0$, one has

[TABLE]

where $\tau_{r}=\frac{\eta^{2}\rho^{2}C^{2}_{r1}}{8\overline{\alpha}}\left(\frac{3}{8\pi}\right)^{2\nu+\sum^{r}_{i=1}(\alpha_{i}-1)}$, and $C_{r1}$ and $\rho$ appear in (45) and (48), respectively.

**Assumption A.5. ** The $r$-dimensional function $f(t,{\bf x})$ belongs to an anisotropic Besov space. In particular, its wavelet coefficients $\beta_{j_{1},k_{1};j_{2},k_{2};\cdots;j_{r},k_{r}}$ satisfy

[TABLE]

where ${\bf j}=(j_{1},j_{2},\cdots,j_{r})^{T}$, ${\bf s^{*}}=(s^{*}_{1},s^{*}_{2},\cdots,s^{*}_{r})^{T}$ and $s^{*}_{i}$ are defined in (23).

Define the quantity

[TABLE]

and denote the dimension of $f$ that corresponds to (53) by $i_{o}$.

Theorem 3

Let $\min\{s_{1},s_{2}\}\geq\max\{\frac{1}{\pi},\frac{1}{2}\}$ with $1\leq\pi,q\leq\infty$, and $A>0$. Then, under conditions (45) and (52), for $1\leq p<\infty$, as $\varepsilon\rightarrow 0$,

[TABLE]

Theorem 4

Let $\widehat{f}(.,.)$ be the wavelet estimator in (43), with $J_{1}$ and $J_{i}$ with $i=2,\cdots,r$ given by (49). Let $\min\{s_{1},s_{2},\cdots,s_{r}\}\geq\max\{\frac{1}{\pi},\frac{1}{2}\}$ with $1\leq\pi,q\leq\infty$, and let conditions (45) and (52) hold. If $\rho$ in (48) is large enough, then, for $1\leq p\leq\infty$, as $\varepsilon\rightarrow 0$,

[TABLE]

where $\xi_{1}$ and $\xi_{2}$ are defined as

[TABLE]

Proof of Theorem 4. The proof is very similar to that of Theorem 5 in Benhaddou et al. (2013) so we skip it. $\Box$

Remark 7

Notice that the convergence rates in Theorems 3 and 4 depend on a delicate balance between $s_{1}$, $\pi$, $\nu$, $\alpha_{1}$, $\min\left\{\frac{s_{i}}{\alpha_{i}},i=2,3,\cdots,r\right\}$ and the average of the parameters of the anisotropic fractional Brownian sheet, $\alpha_{i}$, $i=1,2,\cdots,r$. Such rates are comparable to those in Benhaddou et al. (2013) in their white noise case when $\alpha_{1}=\alpha_{2}=\cdots=\alpha_{r}=1$. In addition, the rates are independent of the dimension $r$ as expected.

6 Proofs.

In the proofs of Theorems 1 and 2 we will use the properties described in Petsa and Sapatinas (2009) associated with Meyer wavelet basis. Namely, the properties of concentration, unconditionality and Temlyakov.

6.1 Proof of the lower bounds.

In order to prove Theorem 1, we use Lemma A.1 of [Bunea, Tsybakov & Wegkamp (2007)].

Lemma 6

*Let $\Theta$ be a set of functions of cardinality $\mbox{card}(\Theta)\geq 2$ such that

(i) $\|f-g\|_{p}^{p}\geq 4\delta^{p},\ for\ f,g\in\Theta,\ f\neq g,$

(ii) the Kullback divergences $K(P_{f},P_{g})$ between the measures $P_{f}$ and $P_{g}$ satisfy the inequality $K(P_{f},P_{g})\leq\log(\mbox{card}(\Theta))/16,\ for\ f,\ g\in\Theta$.

Then, for some absolute positive constant $C_{1}$, one has*

[TABLE]

where $\inf_{f_{n}}$ denotes the infimum over all estimators.

Proof of Theorem 1. Let $P_{f}$ be the probability law of the process $\{q(t,x)+\varepsilon^{\overline{\alpha}}dB^{\alpha}(t,x),(t,x)\in U\}$ when $f$ is true. Let us first transform this process into a process that involves the standard Brownian sheet, taking into account relation (2). This can be accomplished using the following functions

[TABLE]

where $C_{H}$ is some explicit constant, $\Gamma(\cdot)$ is the Gamma function. Consider $0<u<t<1$, then (58) reduces to

[TABLE]

Now define the operator $K_{\Gamma}$ such that

[TABLE]

Taking into account relation (2), it is well known that operator (60) converts an anisotropic two-dimensional fractional Brownian sheet with Hurst parameters $(H_{1},H_{2})=(1-{\alpha_{1}}/{2},1-\alpha_{2}/2)$ into a standard Brownian sheet. Let us then apply operator (60) to equation (1) to obtain,

[TABLE]

where $B(t,x)$ are standard Brownian sheet. Finally, differentiating both sides of (61) with respect to $t$ and $x$, yields

[TABLE]

where $K^{\prime}_{\Gamma}q(t,x)$ is

[TABLE]

Notice that in the Fourier domain, by the convolution theorem and condition (18), (63) has the form

[TABLE]

The dense-dense case.** Let $\eta$ be the matrix with components $\eta_{k_{1}k_{2}}=\{0,1\}$, $k_{i}=0,1,\cdots,2^{j_{i}}-1$, $i=1,2$, and denote the set of all possible values $\eta$ by $H$. Define the functions**

[TABLE]

Note that the matrix $\eta$ has $N=2^{j_{1}+j_{2}}$ components and therefore, $\mbox{card}(H)=2^{N}$. To guarantee that functions (65) satisfy (26), choose $\gamma_{j_{1}j_{2}}=A2^{-j_{1}(s_{1}+1/2)}2^{-j_{2}(s_{2}+1/2)}$. Now, take another function of the form (65), $\tilde{f}_{j_{1}j_{2}}$, with $\tilde{\eta}_{k_{1}k_{2}}\in H$ instead of $\eta_{k_{1}k_{2}}$ and compute the $L^{p}$-norm of the difference between $f_{j_{1}j_{2}}$ and $\tilde{f}_{j_{1}j_{2}}$, with the application of the Varshamov-Gilbert Lemma ([Tsybakov (2008)], p 104), to obtain

[TABLE]

In addition, define $\widehat{P}_{f}$, the probability law of the process $\left\{K^{\prime}_{\Gamma}q(t,x)+\varepsilon^{\overline{\alpha}}dW(t,x),(t,x)\in U\right\}$ when $f$ is true, where $W(t,x)$ is a two-parameter standard Wiener sheet. Note that the probability measures $\widehat{P}_{f}$ and ${P}_{f}$ are stochastically equivalent (see, e.g., Huang et al. (2006), Theorem 3.1). Hence, by the multi-parameter Girsanov formula (see, e.g., Dozzi (1989), p.89), the Kullback divergence takes the form

[TABLE]

Since $\left|\tilde{\eta}_{k_{1}k_{2}}-\eta_{k_{1}k_{2}}\right|\leq 1$, plugging $f_{j_{1}j_{2}}$ and $\tilde{f}_{j_{1}j_{2}}$ into (67), applying Plancherel’s formula, $|\psi_{j,k,m}|\leq 2^{-j/2}$ and (64), it can be shown that

[TABLE]

Part $(ii)$ of Lemma 6 gives the constraint

[TABLE]

Now, define

[TABLE]

Therefore, we need to find combination $\{j_{1},j_{2}\}$ which is the solution to the following optimization problem

[TABLE]

It is easy to check that the solution is $\{j_{1},j_{2}\}=\{\frac{\tau_{\varepsilon}}{2s_{1}+2\nu+\alpha_{1}},0\}$, if $\alpha_{2}s_{1}<s_{2}(2\nu+\alpha_{1})$, and $\{j_{1},j_{2}\}=\{0,\frac{\tau_{\varepsilon}}{2s_{2}+\alpha_{2}}\}$, if $\alpha_{2}s_{1}\geq s_{2}(2\nu+\alpha_{1})$. Hence, the lower bounds are

[TABLE]

The sparse-dense case.** Using the same test functions $\tilde{f}_{j_{1}j_{2}}$ and $f_{j_{1}j_{2}}$, as in Benhaddou et al. (2013), with a finitely supported basis $\eta_{j_{2}k_{2}}$ replaced by a band-limited basis $\psi_{j_{2}k_{2}}$, and following the same procedure as in the dense-dense case, it can be shown that the lower bounds are**

[TABLE]

To complete the proof, notice that the highest of the lower bounds corresponds to

[TABLE]

$\Box$

**Proof of Theorem 3.

The dense-dense case. Let $\eta$ be the matrix with components $\eta_{{\bf k}}=\{0,1\}$, ${\bf k}=(k_{1},k_{2},\cdots,k_{r})$, $k_{i}=0,1,\cdots,2^{j_{i}}-1$, $i=1,2,\cdots,r$, and denote the set of all possible values $\eta$ by $H$. Define the functions**

[TABLE]

Note that here $\mbox{card}(H)=2^{N}$, with $N=2^{j_{1}+j_{2}+\cdots+j_{r}}$. Hence, following the same steps as in the proof of Theorem 1, it can be shown that Note that the lower bounds are

[TABLE]

**and $i=2,3,\cdots,r$.

The sparse-dense case. We apply similar approach to test functions**

[TABLE]

Here $N=2^{j_{2}+\cdots+j_{r}}$ and $\mbox{card}(H)=2^{N}$ and it can be shown that the lower bounds are

[TABLE]

and $i=2,3,\cdots,r$. Hence, the lower bounds of the $L^{p}$-risk correspond to

[TABLE]

and the pair $(s_{2o},\alpha_{2o})$ is defined in (53). $\Box$

6.2 Proof of the upper bounds.

Proof of Lemma 1. Define the quantities

[TABLE]

Consider the Riesz poly-potential operator

[TABLE]

where $\gamma_{2}(\alpha)=4\Gamma\left(\frac{1}{2}-\frac{\alpha_{1}}{2}\right)\Gamma\left(\frac{1}{2}-\frac{\alpha_{2}}{2}\right)cos((1-\alpha_{1})\frac{\pi}{4})cos((1-\alpha_{2})\frac{\pi}{4})$. Then, the anisotropic two-dimensional fractional Brownian sheet $Z^{\alpha}(t,x)=dW^{\alpha}(t,x)$ allows the wavelet-vaguelette representation

[TABLE]

where $\zeta_{j_{1}k_{1};j_{2}k_{2}}$ are white noise processes, and $\psi_{jk}$ is a Meyer-type wavelet basis. Then applying the two dimensional Fourier transform to (82) yields

[TABLE]

where ${\bf{m}}=(m_{1},m_{2})$. Let us evaluate the covariance of (83). Indeed,

[TABLE]

Consequently,

[TABLE]

Now, evaluating the magnitude of (85), and using Hölder’s Inequality along with the fact that $|\tilde{\psi}_{jk}(m)|\leq 2^{-j/2}$, yields

[TABLE]

Hence,

[TABLE]

Finally, let us evaluate the variance of $\tilde{\beta}_{\omega}$. Indeed, using (80) and (86), the variance has the form

[TABLE]

**Results (19) and (20) follow from properties of Gaussian random variables. **$\Box$

Proof of Lemma 2. First, note that under Assumption A.3., one has

[TABLE]

If $p\geq\pi$, one has

[TABLE]

If $p<\pi$, then by Hölder’s Inequality, one has

[TABLE]

**Combining the above results completes the proof. **$\Box$

Proof of Lemma 3. Since the quantities (80) are zero mean Gaussian random variables having variance of form (87), by the Gaussian tail probability inequality, one has

[TABLE]

**where $Z$ is a standard normal, $C$ is some explicit positive constant and $\tau$ appears in Lemma 2. **$\Box$

Proof of Theorem 2. The proof follows a procedure which is similar to that in Benhaddou et al. (2013) with the adaptation to the $L^{p}$-norm case. Denote

[TABLE]

and observe that with the choice of $J_{1}$ and $J_{2}$ given by (22), the $L^{p}$-risk allows the decomposition

[TABLE]

where

[TABLE]

For $R_{1}$ and $R_{4}$, using (19) in $R_{1}$ and (26) in $R_{4}$, as $\varepsilon\rightarrow 0$, yields

[TABLE]

Using Minkowski’s Inequality for $p\geq 2$ and Jensen’s Inequality for $p<2$, and Temlyakov property, $R_{2}$ and $R_{3}$ can be partitioned as $R_{2}\leq R_{21}+R_{22}$ and $R_{3}\leq R_{31}+R_{32}$, where

[TABLE]

Combining (93) and (95), similar calculations as in Benhaddou et al. (2013), yield

[TABLE]

Now, combining (94) and (96), and using (19) and (21), gives

[TABLE]

Finally, $\Delta$ can be decomposed into the following components

[TABLE]

where $\Xi=\left\{j_{1},j_{2}:2^{\frac{p}{2}(j_{1}(2\nu+\alpha_{1})+\alpha_{2}j_{2})}\leq\left[\chi_{\varepsilon}\right]^{d-\frac{p}{2}}\right\}$.

Case 1: $s_{1}\geq\frac{s_{2}}{\alpha_{2}}(2\nu+\alpha_{1})$ and $\frac{s_{2}}{\alpha_{2}}>\frac{p}{2}(\frac{1}{p^{\prime}}-\frac{1}{p})$.** In this case, $d=\frac{ps_{2}}{2s_{2}+\alpha_{2}}$, and**

[TABLE]

Case 2: $(2\nu+\alpha_{1})\frac{p}{2}(\frac{1}{\pi}-\frac{1}{p})<s_{1}<\frac{s_{2}}{\alpha_{2}}(2\nu+\alpha_{1})$****. In this case, $d=\frac{ps_{1}}{2s_{1}+2\nu+\alpha_{1}}$, and

[TABLE]

Case 3: $s_{1}\leq(2\nu+\alpha_{1})\frac{p}{2}(\frac{1}{\pi}-\frac{1}{p})$ and $\frac{s_{2}}{\alpha_{2}}>\frac{p}{2}(\frac{1}{p^{\prime}}-\frac{1}{p})$****. In this case, $d=\frac{p(s_{1}+\frac{1}{p}-\frac{1}{\pi})}{2s^{*}_{1}+2\nu+\alpha_{1}-1}$, and

[TABLE]

**Combining the results from (92) to (104) completes the proof. $\Box$ **

Proof of Lemma 4. The proof will be similar to the proof of Lemma 1. Define the quantities

[TABLE]

Consider the Riesz poly-potential operator

[TABLE]

where $\gamma_{r}(\alpha)=2^{r}\Pi^{r}_{i=1}\left[\Gamma\left(\frac{1}{2}-\frac{\alpha_{i}}{2}\right)cos((1-\alpha_{i})\frac{\pi}{4})\right]$. Then, the anisotropic $r$-dimensional fractional Brownian sheet $Z^{\alpha}(t,x_{1},\ldots,x_{r-1})=dW^{\alpha}(t,x_{1},\ldots,x_{r-1})$ allows the wavelet-vaguelette representation

[TABLE]

Applying the two dimensional Fourier transform to (107), yields

[TABLE]

The covariance of (108) is given by

[TABLE]

Evaluating the magnitude of (109), yields

[TABLE]

Using (105) and (110), similar to (87), the variance has the form

[TABLE]

**Finally, results (46) and (47) follow from properties of Gaussian random variables. **$\Box$

Proof of Lemma 5. The proof is very similar to that of Lemma 3, so we skip it. $\Box$