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
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 .
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
Hence, we choose our thresholds of the form
Finally, choose and such that
For J defined in (21), one has
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
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.
Then, the following statement is true.
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 .
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
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).
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.
Proof of Lemma 1. From (1), one has that
This is because on , one has
On one hand, one has
On the other hand,
and can be partitioned to and , where