Blind Deconvolution for Jump-Preserving Curve Estimation

Abstract

In many applications, observed signals are contaminated by both random noise and blur. This paper proposes a blind deconvolution procedure for estimating a regression function with possible jumps preserved, by removing both noise and blur when recovering the signals. Our procedure is based on three local linear kernel estimates of the regression function, constructed from observations in a left-side, a right-side, and a two-side neighborhood of a given point, respectively. The estimated function at the given point is then defined by one of the three estimates with the smallest weighted residual sum of squares. To better remove the noise and blur, this estimate can also be updated iteratively. Performance of this procedure is investigated by both simulation and real data examples, from which it can be seen that our procedure performs well in various cases.

1. Introduction

Nonparametric regression analysis provides us statistical tools for estimating regression functions from noisy data [1]. When the underlying regression function has jumps, the estimated functions by conventional nonparametric regression procedures are not statistically consistent at the jump positions. However, the problem to estimate jump regression functions is important because the true regression functions are often discontinuous in applications [2]. This paper focuses on estimation of the jump regression function when the observed data are contaminated by both random noise and blur.

In the literature, there are some existing procedures for estimating regression curves with jumps preserved in cases when only random noise is present in observed data. These procedures include the one-sided kernel estimation methods (e.g., [3–9]), the local least-squares estimation procedures (e.g., [10–12]), the wavelet transformation method [13], the spline smoothing method [14], and the direct estimation methods (e.g., [15, 16]).

In some applications, our observations are both blurred and contaminated by pointwise noise (e.g., signals of groundwater levels in geothermy). It is, therefore, important to remove both noise and blur when estimating the true regression function. In the nonparametric regression literature, we have not seen any discussion about this problem yet. In the context of image processing, which can be regarded as a two-dimensional nonparametric regression problem [2], there is a related research area concerned about deblurring images. Most existing image deblurring procedures assume that the point-spread function (psf), which describes the blurring mechanism, either is known or has a parametric form, and this function is homogeneous in the entire image (e.g., [17–20]). In many applications, however, it is difficult to specify the psf completely or partially by a parametric model. Our major goal in this paper is to provide a method to estimate the true regression function properly in cases when both noise and blur are present and the psf is unspecified.

The remaining part of the paper is organized as follows. In next section, our proposed method is discussed in detail. Some comparative results are presented in Section 3 with several simulated examples. A real data application is presented in Section 4. Some concluding remarks are included in Section 5.

2. Our Proposed Method

Suppose that the regression model concerned is𝑦_𝑖=ℎ⊗𝑓(𝑥_𝑖)+𝜀_𝑖,𝑖=1,2,…,𝑛, (2.1) where 0<𝑥₁<𝑥₂<⋯<𝑥_𝑛<1 are design points, 𝜀_𝑖 are i.i.d. random noise with mean 0 and variance 𝜎², and 𝑓 is an unknown nonparametric regression function that is continuous in [0,1] except on some jump points 0<𝑠₁<𝑠₂<⋯<𝑠_𝑚<1. In (2.1), ℎ is a psf, and ℎ⊗𝑓 denotes the convolution between ℎ and 𝑓, defined byℎ⊗𝑓(𝑥)=_𝑅ℎ(𝑢)𝑓(𝑥−𝑢)𝑑𝑢. (2.2)

If ℎ is not zero at the origin and zero everywhere else, then there is no blurring in the observed curve. In such cases, model (2.1) is the conventional nonparametric regression model. Generally speaking, the problem described by model (2.1) is ill posed if both ℎ and 𝑓 are unknown, because infinite sets of ℎ and 𝑓 would correspond to the same response in such cases. Therefore, proper estimation of 𝑓 based on (2.1) is a challenging problem. As a demonstration, in Figure 1(a), the solid lines denote a true regression function 𝑓 with one jump point at 𝑥=0.5, the dotted curve denotes the blurred regression function ℎ⊗𝑓, and the small pluses denote a set of observations from model (2.1).

(a)

(b)

Figure 1 

(a) Solid lines denote the true regression function, and dotted curve denotes the blurred regression function. (b) Solid lines denote the estimate of the regression function by the proposed method, dotted curve denotes the blurred true regression function, and dashed curve denotes the conventional local linear kernel estimate of the regression function. In both plots, little “$+$”s denote a set of observations generated from model (2.1).

To estimate 𝑓, by the conventional local linear kernel (LLK) smoothing [1], we would consider the problemmin^𝑛_{𝑎,𝑏𝑖=1}{𝑦_𝑖−[𝑎+𝑏(𝑥_𝑖−𝑥)]}²𝑥𝐾(^𝑖−𝑥ℎ_𝑛), (2.3) where 𝐾 is a density kernel function with support [−1,1] and ℎ_𝑛 is a bandwidth parameter. Then the solution of (2.3) to 𝑎, denoted as 𝑎_𝑐(𝑥), is defined as the LLK estimator of 𝑓(𝑥). In Figure 1(b), the dotted curve denotes the blurred regression function, and the dashed curve denotes the LLK estimate of 𝑓. It can be seen that the LLK estimate removes the noise well; but the blur is not removed and it is actually made more serious around the jump point.

Qiu [15] finds that the major reason why the conventional LLK estimate could not preserve the jump at the jump point is that it uses a “continuous” function (i.e., a linear function) locally to approximate the jump regression function. To overcome this limitation, Qiu suggests fitting a piecewise linear function around 𝑥 as follows: min_{𝑎𝑙,𝑏𝑙;𝑎𝑟,𝑏𝑟}^𝑛_𝑖=1{𝑦_𝑖−[𝑎_𝑙+𝑏_𝑙(𝑥_𝑖−𝑥)]−[(𝑎_𝑟−𝑎_𝑙)𝐼(𝑥_𝑖−𝑥)+(𝑏_𝑟−𝑏_𝑙)(𝑥_𝑖−𝑥)𝐼(𝑥_𝑖−𝑥)]}²𝑥𝐾(^𝑖−𝑥ℎ_𝑛), (2.4) where 𝐼(⋅) is an indicator function defined by 𝐼(𝑢)=1 if 𝑢≥0 and 0 otherwise. The solution to 𝑎_𝑙 and 𝑎_𝑟 are denoted as 𝑎_𝑙(𝑥) and 𝑎_𝑟(𝑥). It is easy to see that 𝑎_𝑙(𝑥) and 𝑎_𝑟(𝑥) are actually LLK estimates of 𝑓(𝑥) constructed from observations in the left-sided neighborhood [𝑥−ℎ_𝑛,𝑥) and the right-sided neighborhood [𝑥,𝑥+ℎ_𝑛], respectively, with kernel functions 𝐾_𝑙(𝑥)=𝐾(𝑥)𝐼(−𝑥) and 𝐾_𝑟(𝑥)=𝐾(𝑥)𝐼(𝑥). Then, Qiu suggests the following jump-preserving estimate of 𝑓(𝑥):𝑓₁𝑎(𝑥)=^𝑙(𝑥),ifWRMS_𝑙(𝑥)<WRMS_𝑟𝑎(𝑥),_𝑟(𝑥),ifWRMS_𝑙(𝑥)>WRMS_𝑟(𝑎(𝑥),_𝑙(𝑥)+𝑎_𝑟(𝑥))2ifWRMS_𝑙(𝑥)=WRMS_𝑟(𝑥), (2.5) where WRMS_𝑙 and WRMS_𝑟 are the Weighted Residual Mean Squares(WRMSs)of the left-sided and right-sided estimates, respectively, defined byWRMS_𝑙∑(𝑥)=^𝑛𝑖=1[𝑌_𝑖−𝑎_𝑙−𝑏̂_𝑙(𝑥_𝑖−𝑥)]²𝐾_𝑙((𝑥_𝑖−𝑥)/ℎ_𝑛)∑^𝑛_𝑖=1𝐾_𝑙(𝑥_𝑖−𝑥/ℎ_𝑛)WRMS_𝑟∑(𝑥)=^𝑛𝑖=1[𝑌_𝑖−𝑎_𝑟−𝑏̂_𝑟(𝑥_𝑖−𝑥)]²𝐾_𝑟((𝑥_𝑖−𝑥)/ℎ_𝑛)∑^𝑛_𝑖=1𝐾_𝑟(𝑥_𝑖−𝑥/ℎ_𝑛). (2.6)

Qiu [15] proves that 𝑓₁ is a consistent estimate of 𝑓 when there is no blurring in the observed data.

Since only part of observations is actually used in 𝑓₁, this estimator would be quite noisy in continuity regions of 𝑓. To overcome this problem, similar to the method in [16], we propose a modification as follows:𝑓₂(𝑥)=𝑎(𝑥),ifWRMS(𝑥)≤min[WWRMS_𝑙(𝑥),WWRMS_𝑟𝑎(𝑥)],_𝑙(𝑥),ifWWRMS_𝑙(𝑥)<WWRMS_𝑟𝑎(𝑥),_𝑟(𝑥),ifWWRMS_𝑙(𝑥)>WWRMS_𝑟(𝑎(𝑥),_𝑙(𝑥)+𝑎_𝑟(𝑥))2ifWWRMS_𝑙(𝑥)=WWRMS_𝑟(𝑥). (2.7) By (2.7), when 𝑥 is far away from any jump points, 𝑓(𝑥) would be estimated by the conventional LLK estimate. It would still be estimated by one of the one-sided estimates around the jump points. Therefore, it is expected that 𝑓₂(𝑥) would be better for estimating 𝑓(𝑥) than 𝑓₁(𝑥). The solid curve in Figure 2 denotes 𝑓₂ constructed from the data shown in either plot of Figure 1. The two dotted curves show 𝑎_𝑙 and 𝑎_𝑟, respectively. It can be seen that 𝑓₂ indeed estimates 𝑓 well in this case.

Figure 2 

The solid curve denotes ${\widehat{f}}_2$ constructed from the data shown in either plot of Figure 1. The two dotted curves show ${\widehat{a}}_l$ and ${\widehat{a}}_r$, respectively.

The estimate defined in (2.7) can also be updated iteratively as follows. The estimated values {𝑓₂(𝑥_𝑖),𝑖=1,2,…,𝑛} can be used as observed data, and the estimate 𝑓₂ can be updated by (2.7) with all quantities on the right-hand side of (2.7) computed from {𝑓₂(𝑥_𝑖),𝑖=1,2,…,𝑛}, and this process can continue iteratively. Numerical results in the next section suggest that a good estimate can usually be generated after about 5 iterations in all the cases considered there.

In our procedure, the bandwidth parameter ℎ_𝑛 should be chosen properly. To this end, we use the following cross-validation (CV) procedure:ℎ_𝑛=argmin_ℎ𝑛1𝑛^𝑛_𝑖=1(𝑦_𝑖−𝑓_−𝑖(𝑥))², (2.8)

where 𝑓_−𝑖(𝑥) is the estimate of 𝑓(𝑥) using all of the data points except the 𝑖th point (𝑥_𝑖,𝑦_𝑖).

3. Simulation Study

In this section, some simulated examples are presented concerning the numerical performance of our proposed procedure. In all numerical examples presented in this paper, the Epanechnikov kernel function 𝐾(𝑥)=3/4(1−𝑥²)𝐼(|𝑥|≤1) is used. We consider the following two true regression functions: 𝑓₁(𝑥)=4𝑥,if0≤𝑥≤0.2,4(0.4−𝑥),if0.2<𝑥≤0.4,exp^−2𝑥2sin(2.5𝜋𝑥)−1,if0.4<𝑥≤0.8,exp^−2𝑥2𝑓sin(2.5𝜋𝑥),if0.8<𝑥≤1,₂(𝑥)=2−2(0.26−𝑥)^0.2,if0≤𝑥≤0.26,2−2(𝑥−0.26)^0.6,if0.26<𝑥≤0.78,2−2(𝑥−0.26)^0.6+1,if0.78<𝑥≤1. (3.1)

The function 𝑓₁ has two jumps at 0.4, 0.8, and a roof discontinuity (i.e., jump in the slope) at 0.2. The function 𝑓₂ has a jump at 0.78, and an unbalanced cusp (i.e., a sharp angle) at 0.26. These functions are shown by the solid curves in Figure 3. Our observations are generated from model (2.1) with random noise from the 𝑁(0,𝜎²) distribution and the psf ℎ=(1−𝑥²)(1+cos(1.5𝑥𝜋))𝐼(|𝑥|≤1). One set of observations from each regression function is shown by the little pluses in Figure 3.

(a)

(b)

Figure 3 

(a) Solid curve denotes the true regression function $f_1$, and little pluses denote a set of observations when $n=400$ and $\sigma =0.2$. (b) Solid curve denotes the true regression function $f_2$, and little pluses denote a set of observations when $n=400$ and $\sigma =0.2$.

For the proposed procedure, its Mean-Squared Error (MSE) values with several different (𝑛,𝜎) combinations are presented in Figure 4 as functions of the number of iterations, where the bandwidth ℎ_𝑛 is chosen by the CV procedure in each iteration. All of the MSE values are computed based on 100 replications. From the plots in Figure 4, it can be seen that, for each (𝑛,𝜎) combination, MSE values first decrease and then increase with the iteration number. The optimal number of iteration is around 5 in each case, which is the number that we recommend to use in applications.

(a)

(b)

(c)

(d)

Figure 4 

MSE values of the estimated regression function. (a) $f=f_1$ and $n=200$, (b) $f=f_1$ and $n=400$, (c) $f=f_2$ and $n=200$, (d) $f=f_2$ and $n=400$.

Next, we compare the proposed procedure (denoted as NEW) with the conventional local linear kernel (LLK) smoothing procedure and the jump-preserving curve estimation (JPCE) procedure by Qiu [15]. Figure 5 presents the estimated regression functions by all three methods from the observed data shown in Figure 3. For procedure NEW, 5 iterations are used. From the plots in Figure 5, it can be seen that LLK blurs the jumps, JPCE preserves the jumps well but its estimates are quite noisy in continuity regions, and our proposed procedure NEW preserves the jumps well and also removes noise efficiently.

(a)

(b)

Figure 5 

Solid curves denote the estimates by the proposed procedure NEW, dotted curves denote the estimates by LLK, and dashed curves denote the estimates by JPCE.

Tables 1 and 2 present the MSE values and the corresponding standard errors of the three methods in various cases. We use 5 iterations in the proposed procedure NEW. From Tables 1 and 2, it can be seen that procedure NEW performs the best in all cases.

Table 1 

Comparison of the MSE values of the three methods in various cases when 𝑓=𝑓₁. The numbers in parentheses are the standard errors.

Table 2 

Comparison of the MSE values of the three methods in various cases when 𝑓=𝑓₂. The numbers in parentheses are the standard errors.

4. An Application

In this section, we apply our proposed method to a groundwater level data. Possible jumps in groundwater level arise from changes in subsurface fluid currents, which has become an important predictor of earthquakes. In Figure 6, little pluses denote the groundwater level observed by the Seismograph Network Stations of China Earthquake Center during January and May 2007. The solid curve denotes the estimated regression curve by our proposed procedure in which all parameters are chosen in the same way as in the simulation examples presented in Section 3. As indicated by the plot, the jumps around Mar 23, Aril 10, and Aril 17 are preserved well by our procedure. We checked the earthquake history and it is confirmed by China Earthquake Center that earthquakes with magnitude of more than 4.0 occurred in these periods in the area of observation.

Figure 6 

Little pluses denote the groundwater level observed by the Seismograph Network Stations of China Earthquake Center during January and May 2007, and solid curve denotes the estimated regression curve by our proposed procedure.

5. Concluding Remarks

We have presented a blind deconvolution procedure for jump-preserving curve estimation when both noise and blur are present in the observed data. Numerical results show that it performs well in various cases. However, theoretical properties of the proposed method are not available yet, which requires much future research. We believe that the proposed method can be generalized to two-dimensional cases to solve problems such as image deblurring and restoration.

Acknowledgments

This research was supported in part by an NSF grant. The authors thank the guest editor of this special issue, Professor Ming Li, for help during the paper submission process. They also thank the three anonymous referees for their careful reading of the paper.