A) PE risk for pointwise shrinkage

Since the estimate of s _i does not depend on x _j, for j ≠ i, in the pointwise shrinkage scenario, we drop the index i for brevity of notation. The PE risk is defined as

(1)\begin{equation} \mathcal{R} = \mathbb{P}\left( \left\vert \,\widehat{s}-s \right\vert>\epsilon \right), \tag{1}\end{equation}

where ε > 0 is a predefined tolerance parameter. The risk $\mathcal {R}$ quantifies the estimation error using the probability measure and implicitly takes into account the noise distribution beyond the first- and second-order statistics. On the contrary, the MSE relies only on the first- and second-order statistics of noise for shrinkage estimators. Substituting $\widehat {s}=ax=a(s+w)$, the risk $\mathcal {R}$ evaluates to

\begin{align*} \mathcal{R}\left(s;a\right) &= \mathbb{P}\left( \left\vert a(s+w) -s \right\vert > \epsilon \right) \\ &= 1-F \left( \frac{\epsilon-(a-1)s}{a} \right) + F\left( - \frac{\epsilon+(a-1)s}{a} \right), \end{align*}

where F( · ) is the cumulative distribution function (c.d.f.) of the noise. Since $\mathcal {R}$ depends on s, which is the parameter to be estimated, it is impractical to optimize it directly over a. To circumvent the problem, we minimize an estimate of $\mathcal {R}$, which is obtained by replacing s with x (which is also the ML estimate of s). Such an estimate $\widehat { \mathcal {R}}=\mathcal {R}(x;a)$ takes the form

(2)\begin{equation} \widehat{\mathcal{R}}=1-F \left(\frac{\epsilon-(a-1)x}{a} \right) + F\left(-\frac{\epsilon+(a-1)x}{a}\right), \tag{2}\end{equation}

and correspondingly, the optimal shrinkage parameter is obtained as $a_{{\rm opt}} = {\rm arg}\,\underset {0 \leq a \leq 1}{\min }\, \widehat {\mathcal {R}}$. A grid search is performed to optimize $\widehat {\mathcal {R}}$ over a ∈ [0, 1], and the clean signal is obtained as $\widehat {s}=a_{{\rm opt}}x$. We next derive explicit formulae for the risk function for Gaussian, Laplacian, and Student's-t noise distributions.

(i) Gaussian distribution: In this case, the noisy observation x also follows a Gaussian distribution, and therefore, $\widehat {s}-s$ is distributed as $\mathcal {N}((a-1)s,a^2\sigma ^2 )$. The PE risk estimate is given as

(3)\begin{equation} \widehat{\mathcal{R}} = Q \left(\frac{\epsilon-(a-1)x}{a\sigma} \right) + Q \left( \frac{\epsilon+(a-1)x}{a\sigma} \right), \tag{3}\end{equation}

where $Q(u)=({1}/{\sqrt {2 \pi }})\int _{u}^{\infty }e^{-{t^2}/{2}}\,\mathrm {d}t$.

(ii) Student's-t distribution: Consider the case where the noise follow a Student's-t distribution with parameter λ > 2 and the probability density function (p.d.f.) of noise is given by

\[ f(w) = \frac{\Gamma \left( \frac{\lambda+1}{2} \right) }{\sqrt {\lambda \pi } \Gamma \left( \frac{\lambda}{2} \right) } \left( 1+\frac{w^2}{\lambda} \right)^{-{(\lambda+1)}/{2}}.\]

The variance of w is σ² = λ/(λ − 2). The expression for $\widehat {\mathcal {R}}$ is the one given in (2) with

\[ F(w) = \frac{1}{2}+ w\Gamma \left( \frac{\lambda+1}{2} \right) \frac{G_1 \left( \frac{1}{2},\frac{\lambda+1}{2};\frac{3}{2};- \frac{w^2}{\lambda} \right)} {\sqrt {\lambda \pi } \Gamma\left( \frac{\lambda}{2} \right) },\]

where G ₁ is the hypergeometric function defined as

\[ G_1 \left( a,b;c;z \right) = \sum_{k=0}^{\infty}\frac{(a)_k (b)_k}{(c)_k}\frac{z^k}{k!},\]

and (q)_k denotes the Pochhammer symbol:

\[ (q)_k \stackrel{\Delta}{=}\begin{cases} 1 & {\rm for}\ k=0, \\ q(q+1)(q+2)\cdots(q+k-1), & \mbox{for}\ k>0. \end{cases}\]

(iii) Laplacian distribution: Considering the noise to be i.i.d. Laplacian with zero-mean and parameter b (variance σ² = 2b ²), with the p.d.f. f(w) = (1/2b)exp ( − |w|/b), the PE risk can be obtained by using the following expression for F(w) in (2):

(4)\begin{equation} F(w) = \frac{1}{2} + \frac{1}{2} {\rm sgn}(w) \left( 1-\exp \left(-\frac{\vert w\vert }{b} \right) \right). \tag{4}\end{equation}

B) Closeness of $\widehat {\mathcal {R}}$ to $\mathcal {R}$

To measure the closeness of $\widehat {\mathcal {R}}$ to $\mathcal {R}$, consider the example of estimating a scalar s = 4 from a noisy observation x. The PE risk estimate $\widehat {\mathcal {R}}$ is obtained by setting s = x. In Figs 1(a),1(b), and 1(c), we show the variation of the actual risk $\mathcal {R}$ and its estimate $\widehat {\mathcal {R}}$ with a, averaged over 100 independent trials, for Gaussian, Student's-t, and Laplacian noise distributions, respectively. The noise has zero mean, and the variance is taken as σ² = 1 for Gaussian and Laplacian models, whereas for Student's-t model, the variance is σ² = 2. The value of ε is set equal to σ while computing the PE risk. We observe that $\widehat {\mathcal {R}}$ is a good approximation to $\mathcal {R}$, particularly in the vicinity of the minima. The deviation of the shrinkage parameter a _opt(x), obtained by minimizing $\widehat {\mathcal {R}}$, with respect to its true value a _opt(s) resulting from the minimization of $\mathcal {R}$, is shown in Fig. 1(d) for the three noise models under consideration. The central red lines in Fig. 1(d) indicate the medians, whereas the black lines on the top and bottom denote the 25 and the 75 percentile points, respectively. We observe that a _opt(x) is concentrated around a _opt(s), especially for Gaussian and Laplacian noise, barring a few outliers.

Fig. 1. (Color online) The PE risk averaged over 100 realizations for: (a) Gaussian, (b) Student's-t, and (c) Laplacian noise, versus the shrinkage parameter a; and (d) the percentiles of error in minima.

C) Perturbation probability of the location of minimum

The location of the minimum of the PE risk determines the shrinkage parameter. Therefore, one must ensure that it does not deviate too much from its actual value, with high probability, when s is replaced by x in the original risk $\mathcal {R}$. Let $a_{{\rm opt}}(s)=\arg \underset {0\leq a \leq 1}{\min } \mathcal {R}(s;a)$ denote the argument that minimizes the true risk $\mathcal {R}$. Consider the probability of deviation given by

(5)\begin{equation} P_e^{{\rm MPE}}=\mathbb{P}\left( \left\vert a_{{\rm opt}}\left(s\right) - a_{{\rm opt}}\left(x\right) \right\vert \geq \delta \right), \tag{5}\end{equation}

for some δ > 0. Using a first-order Taylor-series approximation of a _opt(x) about s, and substituting x = s + w, we obtain a _opt(x) ≈ a _opt(s) + w a′_opt(s), where ′ denotes the derivative. We would like to point out that the Taylor-series approximation is not rigorously justified here but it is reasonable to guide our analysis. The deviation probability $P_e^{{\rm MPE}}$ in (5) simplifies to $P_e^{{\rm MPE}}=\mathbb {P}( \left \vert w \right \vert \geq {\delta }/{\left \vert a'_{{\rm opt}}(s) \right \vert } ).$ For additive Gaussian noise w with zero mean and variance σ², placing the Chernoff bound [Reference Mitzenmacher and Upfal36,Reference Chang, Cosmon and Milstein37] on $P_e^{{\rm MPE}}$ leads to

\[ P_e^{{\rm MPE}} \leq 2\exp\left( -\frac{\delta^2}{2\sigma^2\left\vert a'_{{\rm opt}} \left(s\right) \right\vert ^2} \right).\]

To ensure that $P_e^{{\rm MPE}}$ is less than α, for a given α ∈ (0, 1), it suffices to have

(6)\begin{equation} \left\vert a'_{{\rm opt}}\left(s\right) \right\vert ^2 \leq \frac{\delta^2}{2\sigma^2 \log \left(\frac{2}{\alpha}\right)}, \tag{6}\end{equation}

which translates to a lower-bound on the input SNR. Since there is no closed-form expression available for $a'_{\rm opt}(s)$ in the context of the PE risk, we empirically obtain the range of input SNR values s ²/σ², for which (6) is satisfied.

Analogously, to satisfy an upper bound on the deviation probability $P_e^{{\rm SURE}}$ of the minimum in the case of SURE, for a given deviation δ > 0, one must ensure that

(7)\begin{equation} \frac{s^6}{8\sigma^6} \left( \delta - \frac{\sigma^4}{\left( s^2+\sigma^2 \right)s^2}\right)^2 \geq \log \left(\frac{2}{\alpha}\right). \tag{7}\end{equation}

The proof of (7) is given in Appendix A.

The minimum input SNR required to ensure P _e ≤ α for both SURE- and MPE-based shrinkage estimators is shown in Fig. 2, for different values of α and δ. The PE risk estimate is obtained by replacing s with x and setting ε = σ. We observe that reducing the amount of deviation δ for a given probability α, or vice versa, leads to a higher input SNR requirement for both SURE and MPE. We also observe from Fig. 2 that, for given δ and α, SURE requires a higher input SNR than MPE to keep the δ-deviation probability under α. Also, for a given input SNR, the δ-deviation probability of the estimated shrinkage parameter a _opt(x) from the optimum a _opt(s) is smaller for MPE than SURE, thereby indicating that the MPE-based shrinkage is comparatively more reliable and robust than the SURE-based one at low input SNRs.

Fig. 2. (Color online) Input SNR requirement for SURE (black) and MPE (blue) to ensure that the probability of δ-perturbation of the minima is less than or equal to α.

D) Unknown noise distributions

In practical applications, the distribution of noise may not be known in a parametric form and may also be multimodal. At best, one would have access to realizations of the noise, from which the distribution has to be estimated. In such cases, approximation of the noise p.d.f. using a GMM is a viable alternative [Reference Plataniotis and Hatzinakos30,Reference Sorenson and Alspach31], wherein one can estimate the parameters of the GMM using the expectation-maximization algorithm [Reference Redner and Walker38]. Gaussian mixture modeling is attractive as it comes with the guarantee that, asymptotically, as the number of Gaussians increases, the GMM approximation to a p.d.f. that has a finite number of discontinuities converges uniformly except at the points of the discontinuity [Reference Sorenson and Alspach31]. The GMM approximation can be used even for non-Gaussian, unimodal distributions. For the GMM-based noise p.d.f.

(8)\begin{equation} f(w)=\sum_{m=1}^{M}\frac{\alpha_m}{\sigma_m \sqrt{2\pi}}\exp \left(-\frac{\left( w-\theta_m \right)^2}{2\sigma_m^2} \right), \tag{8}\end{equation}

the PE risk turns out to be

(9)\begin{align} \widehat{\mathcal{R}} &= \sum_{m=1}^{M}\alpha_m \left[Q \left(\frac{\epsilon-(a-1)s-\theta_m}{a\sigma_m} \right)\right. \nonumber\\ &\quad + \left.Q \left(\frac{\epsilon+(a-1)s+\theta_m}{a\sigma_m} \right)\right], \tag{9}\end{align}

using (3). For illustration, consider the estimation of a scalar s = 4 in the transform domain from its noisy observation x. The additive noise is Laplacian distributed with zero mean and variance σ² = 1. The noise distribution is modeled using a GMM with M = 4 components and the corresponding PE risk estimate is obtained using (9) by setting s = x. In Fig. 3(a), we show a Laplacian p.d.f. and its GMM approximation. Fig. 3(b) shows the GMM approximation to a multimodal distribution. Figure 4(a) shows the PE risk based on the original Laplacian distribution as well as the GMM approximation, as a function of the shrinkage parameter a. The close match between the two indicates that the GMM is a viable alternative when the noise distribution is unknown or follows a complicated model. In Fig. 4(b), we plot the GMM-based PE risk and its estimate averaged over 100 independent trials. We observe that the locations of the minima of the actual risk and its estimate are in good agreement, thereby justifying the minimization of $\widehat {\mathcal {R}}$. The PE risk and its estimate are shown in Fig. 4(c) for the multimodal p.d.f. of Fig. 3(b).

Fig. 3. (Color online) Original noise distribution and a GMM approximation: (a) Laplacian p.d.f. and its approximation using a four-component GMM; and (b) A multimodal p.d.f. and its three-component GMM approximation.

Fig. 4. (Color online) The PE risk estimate versus the shrinkage parameter a: (a) PE risk for Laplacian noise, considering the Laplacian p.d.f. and its GMM approximation; (b) GMM-based PE risk estimate for Laplacian noise; and (c) PE risk estimate for multimodal noise; for ε = σ. The risk estimates are averaged over 100 Monte Carlo realizations.

E) PE risk for subband shrinkage

Let a _J be the shrinkage factor applied to the set of coefficients {x _i, i ∈ J} in subband J. The estimate $\,\widehat {\boldsymbol {s}}_J$ of the clean signal is obtained by $\,\widehat {\boldsymbol {s}}_J= a_J \boldsymbol {x}_J$, where x_J ∈ ℝ^|J| and a _J ∈ [0, 1]. For notational brevity, we drop the subscript J, as we did for pointwise shrinkage, and express the estimator as $\widehat {\boldsymbol {s}} = a \boldsymbol {x}$, where boldface letters indicate vectors. Analogous to pointwise shrinkage, the PE risk for subband shrinkage is defined as $\mathcal {R} = \mathbb {P}( \left \Vert \widehat {\boldsymbol {s}}-\boldsymbol {s} \right \Vert _2>\epsilon ),$ which, for $\,\widehat {\boldsymbol {s}} = a \boldsymbol {x}$, becomes $\mathcal {R} = \mathbb {P}( \left \Vert a \boldsymbol {w} + (a-1)\boldsymbol {s} \right \Vert _2>\epsilon )$. For $\boldsymbol {w} \sim \mathcal {N}(0,\sigma ^2 I)$,

(10)\begin{equation} \mathcal{R} = 1-F(\theta\vert k,\lambda), \tag{10}\end{equation}

where k = |J|, $\lambda = \sum _{j=1}^{k}({(1-a)^2 s_j^2}/{a^2\sigma ^2})$, θ = (ε/aσ)², and F(θ | k, λ) is the c.d.f. of the non-central χ² distribution, given by

\[ F(\theta\,\vert \,k,\lambda)=\sum_{m=0}^{\infty} \frac{\lambda^m e^{-{\lambda}/{2}}}{2^mm!} \frac{\gamma\left(\frac{\theta}{2},\frac{k+2m}{2}\right)}{\Gamma\left(\frac{k+2m}{2}\right)},\]

where $\Gamma (a)\,{=}\,\int _{0}^{\infty } \exp (-t) t^{a-1}\,\mathrm {d} t$ and $\gamma (x,a)\,{=}\, \int _{0}^{x} \exp\break (-t) t^{a-1}\,\mathrm {d} t$.

Similar to pointwise shrinkage, we propose to obtain an estimate $\widehat {\mathcal {R}}$ of $\mathcal {R}$ for subband shrinkage estimators by replacing s _j with x _j. The optimum subband shrinkage factor is obtained by minimizing $\widehat {\mathcal {R}}$.

Fig. 5 shows the subband PE risk and its estimate versus a, where the underlying clean signal s ∈ ℝ^|J| is corrupted by Gaussian noise and the subband size is chosen to be |J| = k = 8. The clean signal s is generated by drawing samples from $\mathcal {N}(2\times \bold 1_k, I_k)$, where $\bold 1_k$ and I _k denote a k-length vector of all ones and a k × k identity matrix, respectively. The observation x is obtained by adding zero-mean i.i.d. Gaussian noise to s, with an input SNR of 5 dB, where the input SNR is defined as ${\rm SNR}_{{\rm in}}= 10 \log _{10}(({1}/{k\sigma ^2}) \sum _{n=1}^{k}s_n^2)\,{\rm dB}$. The PE risk estimate is obtained by replacing s with x in (10), which does not drastically alter the minimum (cf. Fig. 5).

Fig. 5. (Color online) The PE risk and its estimate averaged over 100 Monte Carlo trials for the subband shrinkage estimator versus a; where $\epsilon = \sqrt {k}\sigma$, with k=8. The additive noise is Gaussian with SNR_in = 5 dB. In each trial, s is generated by drawing samples from $\mathcal {N}(2\times \bold 1_k, I_k)$.

A) Performance of pointwise-shrinkage estimator

1) Harmonic signal denoising

Consider the signal

(11)\begin{equation} s_n = \cos\left( \frac{5\pi n}{2048} \right)+ 2 \sin\left( \frac{10\pi n}{2048} \right),\quad 0\leq n \leq 2047, \tag{11}\end{equation}

in additive white Gaussian noise, with zero mean and variance σ². Since the denoising is carried out in the DCT [Reference Rao and Yip39] domain, the Gaussian noise statistics remain unaltered. For the purpose of illustration, we assume that σ² is known. In practice, σ² may not be known a priori and could be replaced by the robust median estimate [Reference Mallat40] or the trimmed estimate [Reference Pastor and Socheleau41]. The clean signal is estimated using inverse DCT after applying the optimum shrinkage. The denoising performance of the MPE and SURE-based approaches is compared in Table 1. In case of the Wiener filter, the power spectrum of the clean signal is estimated using the standard spectral subtraction technique [Reference Boll42,Reference Loizou43]. We observe that MPE-based shrinkage with ε = 3.5σ is superior to SURE and Wiener filter (WF) by 8–12 dB. The comparison also shows that the performance of the MPE depends critically on the choice of ε.

Table 1. Comparison of MPE, SURE-based shrinkage estimator and Wiener filter (WF) for different input SNRs. The output SNR values are averaged over 100 noise realizations.

2) Piece-Regular signal denoising

We consider several noisy realizations of the Piece-Regular signal, taken from the Wavelab toolbox [33], under Gaussian, Student's-t, multimodal (Fig. 3(b)), and Laplacian noise contaminations. The noise variance is assumed to be known. Notably, the Gaussian and Student's-t distributions of noise are preserved by an orthonormal transform [Reference Kibria and Joarder44], unlike the Laplacian distribution. Therefore, the PE estimates for Laplacian noise and multimodal noise are computed based on a four-component and a three-component GMM approximation, respectively, in the DCT domain. The denoised signal corresponding to Laplacian noise is shown in Fig. 6 to serve as an illustration. The MPE estimates are better than SURE estimates. The SNR plots in Fig. 7 indicate that the MPE outperforms SURE for the noise types under consideration and that the gains are particularly high in the input SNR range of −5 to 20 dB and tend to reduce beyond 20 dB.

Fig. 6. (Color online) Denoising performance of the MPE- and SURE-based pointwise shrinkage estimators for the Piece-Regular signal corrupted by Laplacian noise. The PE risk is calculated by using a GMM approximation and setting ε = 3σ.

Fig. 7. (Color online) Output SNR versus input SNR corresponding to the MPE- and SURE-based pointwise shrinkages, under various noise distributions. The output SNR values are calculated by averaging over 100 independent noise realizations. (a) Gaussian noise, (b) Laplacian noise, (c) Student's-t noise, and (d) Multimodal noise.

3) Effect of ε on the denoising performance of MPE

Obtaining a closed-form expression for the ε that maximizes the output SNR is not straightforward. We determine the optimum ε empirically by measuring the SNR gain as a function of ε (cf. Fig. 8) for i.i.d. Gaussian noise. We observe that the output SNR exhibits a peak approximately at β = ε/σ = 3.5 for the harmonic signal in (11) and at β = 3 for the Piece-Regular signal. As a rule of thumb, we recommend ε = 3σ for pointwise shrinkage estimators. Some insights into the role of ε and its choice will be presented in Section 4-4.3.

Fig. 8. (Color online) Average output SNR of the pointwise MPE shrinkage as a function of β = ε/σ, for different values of the input SNR. The output SNR curves peak when β ≈ 3.

B) Performance of subband MPE shrinkage

To validate the performance of the MPE-based subband shrinkage estimator (cf. Section 2-2.5), we consider denoising of the Piece-Regular signal in additive Gaussian noise. The clean signal and its noisy measurement are shown in Fig. 9(a). Denoising is carried out by grouping k adjacent DCT coefficients to form a subband. The denoised signals obtained using SURE and MPE are shown in Fig. 9(b) and 9(c), respectively. The subband size k is chosen to be 16 and the parameter ε is set equal to $1.75\sqrt {k}\sigma$, a value that was determined experimentally and found to be nearly optimal. We observe that the MPE gives 1 dB improvement in SNR than the SURE approach.

Fig. 9. (Color online) Comparison of denoising performance of the subband shrinkage estimators using MPE and SURE, for the Piece-Regular signal corrupted by additive Gaussian noise. The subband size is taken as k = 16 and the value of ε is $1.75\sqrt {k}\sigma$.

The variation of the output SNR is also studied as a function of k (cf. Fig. 10). We experimented with ε = 3σ, $\epsilon =1.75\sqrt {k}\sigma$, and $\epsilon =1.25\sqrt {k}\sigma$ corresponding to subband sizes k = 1, k ∈ [2, 16], and k>16, respectively. For both SURE and MPE, as k increases, the output SNR also increases and eventually saturates for k ≥ 40. For input SNR below 15 dB, MPE gives a comparatively higher SNR than SURE, and the margin diminishes with increase in input SNR or the subband size k. The degradation in performance of SURE for low SNRs is due to the large error in estimating the MSE at such SNRs. The SURE-based estimate of MSE becomes increasingly reliable as k increases, thereby leading to superior performance.

Fig. 10. (Color online) Output SNR versus subband size k, averaged over 100 noise realizations, for different input SNRs. The output SNR of MPE is consistently superior to that obtained using SURE, especially when k ≤ 40 and the input SNR is below 15 dB.

A) Expected ℓ₁ risk using GMM approximation

For the GMM p.d.f. in (8), the expected ℓ₁ distortion evaluates to (cf. Appendix B for the derivation)

(15)\begin{equation} \mathcal{R}_{\ell_1}=\sum_{m=1}^{M} a \alpha_m \sigma_{m} \left(\sqrt{\frac{2}{\pi}}e^{-{\mu_{m}^2}/{2}} -2 \mu_{m} Q \left( \mu_{m} \right) + \mu_{m} \right), \tag{15}\end{equation}

where μ_m = −((a − 1)s + θ_m)/aσ_m. The expected ℓ₁ risk and its estimate for a multimodal (cf. Fig. 3(b)) and Laplacian noise p.d.f.s are shown in Figs 11(b) and 11(c), respectively. We observe that, in both cases, the locations of the minima of the true risk and its estimate are in good agreement.

B) Optimum shrinkage versus posterior SNR

We next study the behavior of a _opt for different input SNRs to compare the denoising capabilities of the MPE and the expected ℓ₁-distortion-based shrinkage estimators. The optimal pointwise shrinkage parameter a _opt for Gaussian noise statistics, obtained by minimizing SURE, PE risk estimate, and the estimated ℓ₁ risk, for different values of the a posteriori SNR x ²/σ² is plotted in Fig. 12(a). To illustrate the effect of ε, the variation of a _opt versus a posteriori SNR for MPE corresponding to Gaussian noise is shown in Fig. 12(b), for different ε. We observe that the shrinkage profiles are characteristic of a reasonable denoising algorithm, as Fig. 12(a) and 12(b) exhibit that the shrinkage parameters increase as the a posteriori SNR increases. Whereas in the case of MPE, the choice of ε is crucial, the expected ℓ₁ distortion does not require tuning such a parameter.

Fig. 12. (Color online) Shrinkage parameter profiles as a function of the a posteriori SNR, corresponding to different risk functions: (a) MPE, SURE, expected ℓ₁ distortion; and (b) MPE for different values of ε. The shrinkage factor a _opt is plotted on a log scale mainly to highlight the fine differences among various attenuation profiles.

C) Shrinkage profiles of estimators based on MPE/ℓ₁ vis-à-vis SURE

To gain insight into the denoising behavior of the MPE and the ℓ₁-minimization-based shrinkage estimators, it is instructive to study the resulting thresholding functions shown in Fig. 13. Figures 13(a) and 13(b) correspond to Gaussian noise. The comparison of the shrinkage functions based on MPE, ℓ₁ risk, and SURE (cf. Fig. 13(a)) reveals that MPE results in a hard-thresholding-type shrinkage, whereas the ℓ₁ risk and SURE produce shrinkage profiles similar to a soft-thresholding function (the threshold operator is called “soft” as it is a continuous function of the input data [Reference Johnstone27]). The shrinkage functions resulting from MPE with different ε are shown in Fig. 13(b). We observe that increasing ε is tantamount to applying a larger attenuation on the noisy signal. In Fig. 13(c), we demonstrate the variation of the MPE-based shrinkage with respect to two different noise distributions, namely, Gaussian and Laplacian, having the same first- and second-order moments. We observe that the MPE-based estimator corresponding to Laplacian noise attenuates larger amplitudes more than its Gaussian counterpart. Such an attenuation profile helps suppress Laplacian noise more effectively, considering the heavier tail of the Laplacian distribution. A distribution-specific shrinkage profile also distinguishes MPE from the conventional MSE-based estimators. In contrast, the MSE-based shrinkage estimator, given by [Reference Sadasivan, Mukherjee and Seelamantula28]: $\hat {s}_{{\rm MSE}}=\max \{0,1-{\sigma ^2}/{x^2}\}x,$ relies only on moments up to second order.

Fig. 13. (Color online) Thresholding functions corresponding to different risk functions: (a) MPE, SURE, expected ℓ₁ distortion; (b) MPE for various values of ε; and (c) MPE, MSE. In (a) and (b), the noise considered is Gaussian. In (c), both Laplacian and Gaussian noise types have been considered with σ = 1 and ε = 3σ.

An inspection of the shrinkage functions corresponding to different ε values in Fig. 13(b) reveals a particularly interesting property: Transform coefficients with magnitudes smaller than ε are set to zero by the MPE-based shrinkage function. Thus, the PE risk essentially results in a distribution-specific threshold operator parameterized by ε. In this context, we recall an observation made by Johnstone et al. in [Reference Johnstone and Silverman45]: A threshold value of 3σ to 4σ or higher works well for signals that are sparse in the transform domain, whereas a threshold value of 2σ or lower is more appropriate for a dense signal. The optimal ε values obtained empirically for the harmonic and the Piece-Regular signals (cf. Section 3-3.1.3.1.3, Fig. 8) also support this statement. Since the harmonic signal is sparser than the Piece-Regular signal in the DCT domain, the optimal ε tends to be higher for the harmonic signal as compared with that for the Piece-Regular signal. Thus, the shrinkage functions shown in Fig. 13(b) provide us with valuable insights on how to fix ε. If the underlying signal is sparse in the transform considered, one should fix ε in the range [3σ, 4σ], whereas for dense signals, one should set ε to about 2σ.

A) Iterative minimization of the expected ℓ₁-risk

When M = 1, the ML-ℓ₁ estimator is obtained by minimizing $\mathcal {R}_{\ell _1}(x,a)$, where x is the noisy version of s. We refer to this estimate as the non-iterative ℓ₁-based shrinkage estimator. Following Algorithm 1, one could iteratively refine the estimate, starting from x. We compare the non-iterative ℓ₁-based estimator with its iterative counterpart (ITER-ℓ₁), and present the results in Figs 15,16, and 17, corresponding to Gaussian, multimodal (Fig. 3(b)), and a GMM approximation to the Laplacian noise, respectively. The output SNR obtained using the oracle-ℓ₁ estimator, calculated by minimizing $\mathcal {R}_{\ell _1}(s,a)$, is also shown for benchmarking the performance.

Fig. 15. (Color online) Performance of ℓ₁-risk-minimization-based pointwise shrinkage estimator (Piece-Regular signal in Gaussian noise): (a) Output SNR versus iterations (input SNR 5 dB); and (b) Output SNR versus input SNR, averaged over 100 independent noise realizations. The number of iterations in Algorithm 1 is fixed at N _iter = 20.

Fig. 16. (Color online) Performance of pointwise shrinkage estimator based on ℓ₁ risk minimization: (a) Output SNR versus iterations, corresponding to the Piece-Regular signal contaminated by noise (SNR 5 dB) whose p.d.f. is given in Fig. 3(b); and (b) Output SNR versus input SNR (averaged over 100 independent noise realizations). The number of iterations in Algorithm 1 is set to N _iter = 20.

Fig. 17. (Color online) Performance of pointwise shrinkage estimator obtained by ℓ₁ risk minimization: (a) Output SNR versus iterations for the Piece-Regular signal in Laplacian noise with an input SNR of 5 dB; and (b) Output SNR versus input SNR for N _iter = 20. In (a) and (b), the Laplacian distribution is modeled using a four-component GMM to calculate the ℓ₁-risk estimate. The results displayed in (b) are obtained after averaging over 100 realizations.

We make the following observations from Figs 15,16, and 17: (i) the output SNR increases with iterations, albeit marginally after about 10 iterations; (ii) the iterative method consistently dominates the non-iterative one, with an overall SNR improvement of about 2–3 dB, for input SNR in the range −5 dB to 20 dB; and (iii) the SNR gain of the iterative technique also reduces for higher input SNR, similar to other denoising algorithms.