Minimax adaptive wavelet estimator for the anisotropic functional deconvolution model with unknown kernel

11/22/2018 ∙ by Rida Benhaddou, et al. ∙ Ohio University 0

In the present paper, we consider the estimation of a periodic two-dimensional function f(·,·) based on observations from its noisy convolution, and convolution kernel g(·,·) unknown. We derive the minimax lower bounds for the mean squared error assuming that f belongs to certain Besov space and the kernel function g satisfies some smoothness properties. We construct an adaptive hard-thresholding wavelet estimator that is asymptotically near-optimal within a logarithmic factor in a wide range of Besov balls. The proposed estimation algorithm implements a truncation to estimate the wavelet coefficients, in addition to the conventional hard-thresholds. A limited simulations study confirms theoretical claims of the paper.



There are no comments yet.


page 1

page 2

page 3

page 4

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.

We consider the estimation problem of an unknown function from observations

contaminated by Gaussian white noise in the following convolution model:


where , is the unknown blurring function with observations , and and are two independent two-dimensional Gaussian white noises with covariance function


and is the Dirac delta function. The parameters and are positive and satisfy asymptotically.
The discrete version of model (1) when and are observed at points , , is as follows


where and are two positive constants independent of and , and . The quantities , with

, are zero-mean i.i.d. normal random variables with

. In addition, and are independent of each other.

Deconvolution model has witnessed a considerable number of publications since late 1980s and Donoho (1995) was the first to devise a wavelet solution to the problem. The list also includes Abramovich and Silverman (1998), Pensky and Vidakovic (1999), Walter and Shen (1999), Johnstone, Kerkyacharian, Picard and Raimondo (2004), Donoho and Raimondo (2004), and Pensky and Sapatinas (2009). Functional deconvolution has been investigated in Benhaddou et al. (2013), where they considered model (1) with , which corresponds to the case when the kernel is known. 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)). Recently, Benhaddou, Pensky and Rajapakshage (2019) investigated functional Laplace deconvolution where the function under consideration is not periodic.

In the present setting, the convolution kernel is unknown, but observations are available. This problem is referred to as the blind deconvolution. Inverse problem with unknown operators in its general aspect was studied by Hoffmann and Reiss (2008) where they proposed two nonlinear methods of estimating a one dimensional function. Delattre et al. (2012) considered the blind deconvolution problem and applied the block SVD procedure to construct a wavelet estimator for a one dimensional function that belongs to Sobolev space. Benhaddou (2018a) considered the blind deconvolution model under fractional Gaussian noise, where the function under consideration is univariate and periodic. The common feature between Hoffmann and Reiss (2008) and Benhaddou (2018a) is that they implement preliminary thresholding procedures that eliminate the estimated wavelet coefficients that are judged to be too large.

The objective of the paper is to construct an adaptive hard-thresholding wavelet estimator for model (1). We focus on the regular-smooth convolution. A preliminary stabilizing thresholding procedure is applied to the functional Fourier coefficients of the ”data” to estimate the wavelet coefficients, taking advantage of the flexibility of the Meyer wavelet basis in the Fourier domain. More specifically, we apply the Meyer wavelet transform in the Fourier domain, and for each resolution level , we truncate the estimated wavelet coefficients at values that are close to zero. We show that the proposed approach is asymptotically near-optimal over a wide range of Besov balls under the -risk. In addition, we show that the convergence rates are expressed as the maxima between two terms, taking into account both the noise sources. Similar behavior has been pointed out in Hoffmann and Reiss (2008), Vareschi (2015) and Benhaddou (2018a, 2018b). It should be noted that with , our convergence rates coincide with those in Benhaddou et al. (2013), and with and , our convergence rates match those in Benhaddou (2017). Finally, with , and with and , our rates are comparable to those in Benhaddou (2018a) and Benhaddou (2018b), respectively.

The rest of the paper is organized as follows. In section 2, we describe in details the estimation procedure. In section 3, we study the asymptotic performance of the proposed estimator in terms of its minimax squared loss. In Section 4, we consider a limited simulations study to assess the goodness of our estimator when the sample size is finite. Finally, Section 5 contains the proofs of our theoretical findings.

2 Estimation algorithm

Let denote the inner product in the Hilbert space , i.e., for , where is the complex conjugate of . Denote as a Fourier basis on the interval . Let , , , , be functional Fourier coefficients of functions

, respectively. Applying Fourier transform to equation (

1) we get


Consider a bounded bandwidth periodized wavelet basis ( Meyer-type) and a finitely supported periodized -regular wavelet basis (Daubechies-type). Denote the wavelet functions of these two bases by and respectively, where and correspond to the lowest resolution levels for the two bases. Then, the function has the wavelet series representation given by




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


due to the fact that Meyer wavelets are band-limited, and the cardinality of is . If were known, the problem would reduce to the regular deconvolution model studied in Benhaddou et al. (2013). However, since is unknown and contaminated with gaussian white noise, a preliminary thresholding procedure is to be applied to provide the estimate of . Let us define the quantities


where is a positive constant independent of and . Then, using (6) and (8), and by Plancherel formula, a truncated estimator for is given by


Bear in mind that the above thresholding on the Fourier coefficients enables us to have a stable inversion. Now, define


Then, consider the hard-thresholding estimator for given by


where and will be determined later, and

will be chosen based on the variance of

Assumption A.1. Assume that the functional Fourier coefficients of are uniformly bounded from below and above, that is, there exists positive constants and , all independent of and , such that


The next step is to evaluate the variance of (9). Define for and positive constant , the sets and such that


and denote .

Lemma 1

Let the constant defined in (14) be such that . Then on , one has

Lemma 2 (Upper bound for )

Let be defined by (9), and let constant be such that . Then on , one has


and on ,


Now, to determine the choice of the thresholds , let us define the quantity


Then by (7) and (12), we obtain


Hence, we choose our thresholds of the form


Finally, choose and such that

Lemma 3

For J defined in (21), one has

Remark 1

Notice that our choices of the threshold , and the finest resolution levels and do not depend on any unknown parameter, and therefore estimator (11) is adaptive. It remains to investigate how the estimator performs. Next, we evaluate the lower and upper bound for the -risk.

3 Minimax rates of convergence

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


where is the -norm of a function and the infimum is taken over all possible estimators of .
Assumption A.2. Let . Assume that belongs to a two-dimensional Besov ball, and its wavelet coefficients satisfy

Theorem 1

Let with , let and define , for . Then under conditions (12) and (25), as , simultaneously,




Before we derive the minimax upper bounds for the -risk, we first prove the next lemma which provides the large deviation inequalities for the wavelet coefficients estimators.

Lemma 4

Let and be defined by (9) and (20) respectively. Define for the set


Then, under and the condition (12), and as , simultaneously, one has




and appear in (12), and , appear in (14).

Then, the following statement is true.

Theorem 2

Let be the wavelet estimator defined in (11), with and given by (21) and (22). Let conditions (12) and (25) hold and , with , and choose and in (29) large enough. Then as , simultaneously,


where is defined in (27) and

Remark 2

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

Remark 3

The convergence rates are expressed as a maxima between two terms, taking into account both noise sources. Similar behavior has been pointed out in Hoffmann and Reiss (2008), Vareschi (2015) and Benhaddou (2018a, 2018b). These rates depend on a delicate balance between the parameters of the Besov ball, smoothness of the convolution kernel and the noise parameters; and .

Remark 4

With , our convergence rates coincide with those in Benhaddou et al. (2013), and with and , our convergence rates match those in Benhaddou (2017). In addition, if we fix the variable , and with , our rates are comparable to those in Benhaddou (2018a) in their univariate case, and with and are convergence rates match those in Benhaddou (2018b), in their univariate but multichannel case.

4 Simulation Study

In this section, we carried out a limited simulation study in order to investigate the finite sample performance of our estimator. The first step though is to provide the sample equivalent to equations (4), (8), (9), (13)-(14) and (20)-(22). Indeed, apply the Fourier transform to (3) to obtain


Then, the discrete version of (9) is given by


In addition, equations (3) and (1) are equivalent by setting


The sets and are now of the form


and (18) has the sample counterpart


Finally, the threshold and the finest resolution levels are given by


The simulation was implemented based on the above equations. In particular, it was formatted through MATLAB using the Wavelab toolbox. Similar to Benhaddou, et al. (2013), degree 3 Meyer wavelets family and degree 6 Daubechies wavelets were utilized for the wavelet transform. We generated our data by equation (33) with testing kernel , and with various testing functions and different combinations of values , and , . In particular, we generate from with being a quadratic function scaled to have a unit norm, and being the routinely used testing functions blip, bumps, heavy sine, and doppler. We rescale to have a unit norm. As noted in Benhaddou, et al. (2013), our method did not know that were generated by a product of two functions, thus can not take advantage of this information. For the choices of and , they are dependent on the signal-to-noise ratio. In particular, , where , and is the signal-to-noise ratio of the first equation in (33); and where is the signal-to-noise ratio of the second equation in (33). We used and .
We evaluated the mean integrated square error (MISE)

of the functional deconvolution estimator. Figure 1 reports the averages of those errors over 100 simulation repetitions together with their standard deviations (in the parentheses).

N= 512 1024 512 1024
3,7 7,7 3,7 7,7 3,7 7,7 3,7 7,7
M=128 0.1441 0.1019 0.0588 0.0138 0.0559 0.0511 0.0211 0.0204
(0.00126) (0.000603) (0.00181) (0.000242) (0.000256) (0.000275) (0.000119) (0.0000909)
M=256 0.1435 0.1010 0.0550 0.02461 0.0555 0.0509 0.0204 0.0204
(0.000966) (0.000437) (0.000932) (0.000205) (0.000187) (0.000173) (0.0000748) (0.0000626)
M=128 0.01585 0.00297 0.00518 0.00132 0.0306 0.0168 0.0131 0.0105
(0.000315) (0.0000727) (0.000291) (0.0000252) (0.000409) (0.0000903) (0.000363) (0.0000884)
M=256 0.01561 0.00280 0.00402 0.00273 0.0302 0.01648 0.0116 0.00999
(0.000280) (0.0000453) (0.000161) (0.0000556) (0.000270) (0.0000631) (0.000214) (0.0000661)
M=128 0.1569 0.1103 0.0811 0.0276 0.0323 0.0179 0.0145 0.00916
(0.00139) (0.000412) (0.00190) (0.000331) (0.000411) (0.0000725) (0.000365) (0.0000530)
M=256 0.1529 0.1103 0.1113 0.0353 0.0320 0.0175 0.0128 0.00862
(0.000861) (0.000346) (0.00119) (0.000262) (0.000306) (0.0000521) (0.000237) (0.0000344)
M=128 0.1520 0.1068 0.0686 0.0187 0.0341 0.0194 0.0170 0.0112
(0.00133) (0.000350) (0.00163) (0.000213) (0.000455) (0.000103) (0.000345) (0.0000879)
M=256 0.1516 0.1064 0.0634 0.0169 0.0336 0.0192 0.0153 0.0108
(0.000927) (0.000219) (0.00110) (0.000134) (0.000271) (0.0000844) (0.000242) (0.0000591)
Table 1: MISE averaged over 100 repetitions

The simulation is aimed to study the effect of two components; the sample budget and the noise levels , . The simulation results confirm the theory and are consistent with previous results in the literature. That is, as the sample size increases ( decreases), the performance of estimation increases. At the same time, as the noise level increases ( increases), the performance deteriorates. The simulation results confirm that as the noise level decreases ( increases), the convergence rate improves.

5 Proofs

Proof of Lemma 1. From (1), one has that




This is because on , one has


Thus, using (43) and (44), yields


On one hand, one has


On the other hand,


Combining (14) and (46) with , one has


Now, using (47) and (48), one obtains


Hence, combining (45) and (49) completes the proof.
Proof of Lemma 2. Note that , thus on , the variance of (9) is


and can be partitioned to and , where