From Denoising to Compressed Sensing

I Introduction

Oftentimes a signal $x_{o}$ is sparse (or approximately sparse) in some transform domain, i.e., $x_{o}=\mathbf{\Psi}u$ with sparse $u$, where $\mathbf{\Psi}$ represents the inverse transform matrix. In this case we lump the measurement and transformation into a single measurement matrix $\mathbf{A}=\mathbf{\Phi}\mathbf{\Psi}$. When a sparsifying basis is not used $\mathbf{A}=\mathbf{\Phi}$. Future references to the measurement matrix refer to $\mathbf{A}$.

which is known formally as basis pursuit denoising (BPDN). It was first shown in that if $x_{o}$ is sufficiently sparse and $\mathbf{A}$ satisfies certain properties, then (1) can accurately recover $x_{o}$.

The initial work in CS solved (1) using convex programming methods. However, when dealing with large signals, such as images, these convex programs are extremely computationally demanding. Therefore, lower cost iterative algorithms were developed; including matching pursuit, orthogonal matching pursuit , iterative hard-thresholding , compressive sampling matching pursuit, approximate message passing, and iterative soft-thresholding , to name just a few. See for a complete set of references.

Iterative thresholding (IT) algorithms generally take the form

where $\eta_{\tau}(y)$ is a shrinkage/thresholding non-linearity, $x^{t}$ is the estimate of $x_{o}$ at iteration $t$, and $z^{t}$ denotes the estimate of the residual $y-Ax_{o}$ at iteration $t$. When $\eta_{\tau}(y)=(|y|-\tau)_{+}{\rm sign}(y)$ the algorithm is known as iterative soft-thresholding (IST).

AMP extends iterative soft-thresholding by adding an extra term to the residual known as the Onsager correction term:

Here, $\delta=m/n$ is a measure of the under-determinacy of the problem, $\langle\cdot\rangle$ denotes the average of a vector, and $\frac{1}{\delta}\langle\eta^{\prime}_{\tau}(\mathbf{A}^{*}z^{t-1}+x^{t-1})\rangle$, where $\eta^{\prime}_{\tau}$ represents the derivative of $\eta_{\tau}$, is the Onsager correction term. The role of this term is illustrated in Figure 1. This figure compares the QQplotA QQplot is a visual inspection tool for checking the Gaussianity of the data. In a QQplot, deviation from a straight line is an evidence of non-Gaussianity. of $x^{t}+\mathbf{A}^{*}z^{t}-x_{o}$ for IST and AMP. We call this quantity the effective noise of the algorithm at iteration $t$. As is clear from the figure, the QQplot of the effective noise in AMP is a straight line. This means that the noise is approximately Gaussian. This important feature enables the accurate analysis of the algorithm , the optimal tuning of the parameters , and leads to the linear convergence of $x^{t}$ to the final solution . We will employ this important feature of AMP in our work as well.

I-B Main contributions

A sparsity model is accurate for many signals and has been the focus of the majority of CS research. Unfortunately, sparsity-based methods are less appropriate for many imaging applications. The reason for this failure is that natural images do not have an exactly sparse representation in any known basis (DCT, wavelet, curvelet, etc.). Figure 2 shows the wavelet coefficients of the classic signal processing image Barbara. The majority of the coefficients are non-zero and many are far from zero. As a result, algorithms that seek only wavelet-sparsity fail to recover the signal.

In response to this failure, researchers have considered more elaborate structures for CS recovery. These include minimal total variation , block sparsity , wavelet tree sparsity , hidden Markov mixture models , non-local self-similarity , and simple representations in adaptive bases . Many of these approaches have led to significant improvements in imaging tasks.

In this paper, we take a complementary approach to enhancing the performance of CS recovery of non-sparse signals . Rather than focusing on developing new signal models, we demonstrate how the existing rich literature on signal denoising can be leveraged for enhanced CS recovery.In this paper, denoising refers to any algorithm that receives $x_{o}+\sigma z$, where $\sigma z\sim N(0,\sigma^{2}I)$ denotes the noise, as its input and returns an estimate of $x_{o}$ as its output. Refer to Sections III-B and VII-A for more information on denoisers. The idea is simple: Signal denoising algorithms (whether based on an explicit or implicit model) have been developed and optimized for decades. Hence, any CS recovery scheme that employs such denoising algorithms should be able to capture complicated structures that have heretofore not been captured by existing CS recovery schemes.

D-AMP has several advantages over existing CS recovery algorithms: (i) It can be easily applied to many different signal classes. (ii) It outperforms existing algorithms and is extremely robust to measurement noise (our simulation results are summarized in Section VII). (iii) It comes with an analysis framework that not only characterizes its fundamental limits, but also suggests how we can best use the framework in practice.

D-AMP employs a denoiser in the following iteration:

Here, $x^{t}$ is the estimate of $x_{o}$ at iteration $t$ and $z^{t}$ is an estimate of the residual. As we will show later, $x^{t}+\mathbf{A}^{*}z^{t}$ can be written as $x_{o}+v^{t}$, where $v^{t}$ can be considered as i.i.d. Gaussian noise.This conjecture has been validated empirically elsewhere for simpler denoisers. For known $\sigma_{t}$ the conjecture has been proven for scalar denoisers in . By combining the proof of with the proof technique developed in we can prove the above conjecture for the scalar denoisers. Since scalar denoisers are not of our main concern in this paper, we do not include a proof here. We will present empirical evidence that it holds for many of the state-of-the-art image denoising algorithms. $\hat{\sigma}^{t}$ is an estimate of the standard deviation of that noise. ${\rm div}D_{\hat{\sigma}^{t-1}}$ denotes the divergence of the denoiser.In the context of this work the divergence ${\rm div}D(x)$ is simply the sum of the partial derivatives with respect to each element of $x$, i.e., ${\rm div}D(x)=\sum\limits_{i=1}^{n}\frac{\partial D(x)}{\partial x_{i}}$, where $x_{i}$ is the $i^{\rm th}$ element of $x$. The term $z^{t-1}{\rm div}D_{\hat{\sigma}^{t-1}}(x^{t-1}+\mathbf{A}^{*}z^{t-1})/m$ is the Onsager correction term. We will show later that this term has a major impact on the performance of the algorithm. The explicit calculation of this term is not always straightforward: many popular denoisers do not have explicit formulations. However, we will show that it may be approximately calculated without requiring the explicit form of the denoiser.

D-AMP applies an existing denoising algorithm to vectors that are generated from compressive measurements. The intuition is that at every iteration D-AMP obtains a better estimate of $x_{o}$ and that this sequence of estimates eventually converges to $x_{o}$.

To predict the performance of D-AMP, we will employ a novel state evolution framework to theoretically track the standard deviation of the noise, $\hat{\sigma}_{t}$ at each iteration of D-AMP. Our framework extends and validates the state evolution framework proposed in . Through extensive simulations we show that in high-dimensional settings (for the subset of denoisers that we consider in this paper) our state evolution predicts the mean square error (MSE) of D-AMP accurately. Based on the state evolution we characterize the performance of D-AMP and connect the number of measurements D-AMP requires to the performance of the denoiser. We also employ the state evolution to address practical concerns such as the tuning of the parameters of denoisers and the sensitivity of the algorithm to measurement noise. Furthermore, we use the state evolution to explore the optimality of D-AMP. We postpone a detailed discussion to Section III.

Figure 3 compares the performance of the original AMP (which employs sparsity in the wavelet domain) with that of D-AMP based on the non-local means denoising algorithms , called NLM-AMP here. Since NLM is a better denoiser than wavelet thresholding for piecewise constant functions, NLM-AMP dramatically outperforms the original AMP. The details of our simulations are given in Section VII-D.

I-C Related work

Note that in the development of D-AMP we are not concerned about whether the algorithm is approximating a posterior distribution for a certain prior or not. Nor are we concerned about whether or not the denoisers used within D-AMP’s iterations are tied to any prior. Instead, we rely on only one important feature of AMP—that $x^{t}+\mathbf{A}^{*}z^{t}-x_{o}$ behaves similar to i.i.d. Gaussian noise. Our analysis is based on this assumption. We validate this assumption with extensive simulations that are presented in Section VII-C.

Donoho et al. also extended the AMP framework based upon the fact that $x^{t}+\mathbf{A}^{*}z^{t}-x_{o}$ behaves similar to i.i.d. Gaussian noise. In their framework the denoiser D can be any scale-invariant function. There are several major differences between our work and theirs: (i) We do not impose scale-invariance on the denoiser, because this assumption does not hold for many practical denoisers. (ii) We present a far broader validation of our method and state evolution: The empirical validation presented is concerned with very specific simple denoisers and has remained at the level of maximin phase transition . Likewise, the state evolution they employed in their empirical validation is based on the Bayesian framework described above and the validations are restricted to simple distributions. In this paper we consider a deterministic version of the state evolution and for the first time present evidence that such a state evolution can in fact predict the performance of D-AMP. Note that the evidence we present goes far beyond the match in the maximin phase transition. This development is important because the maximin framework employed in is not useful in most practical applications that deal with naturally occurring signals. (iii) We present a signal-dependent parameter tuning strategy for AMP and show that our deterministic state evolution can cope with those situations as well. (iv) We show how practical denoisers whose explicit functional form is not given can be employed in AMP. (v) We investigate the optimality of D-AMP as a means to employ different denoisers in the AMP algorithm.

While writing this paper, we became aware of another relevant paper about extensions to the AMP algorithm. In this work, the authors employ AMP with scalar denoisers that are better adapted to the statistics of natural images. By doing so, they have obtained a major improvement over existing algorithms. In this paper, we consider a much broader class of denoisers. We not only show how the AMP algorithm can be adapted to such denoisers; we also explore the theoretical properties of our recovery algorithms.

I-C2 Model-based CS imaging

Many researchers have noticed the weakness of sparsity-based methods for imaging applications and have therefore explored the use of more complicated signal models. These models can be enforced explicitly, by constraining the solution space, or implicitly, by using penalty functionals to encourage solutions of a certain form.

Initially these model-based methods were restricted to simple concepts like minimal total variation and block sparsity , but they have since been extended to structures such as wavelet trees and mixture models . Furthermore, some researchers have employed more complicated signal models through non-local regularization and the use of adaptive over-complete dictionaries . A non-local regularization method, NLR-CS , represents the current state-of-the-art in CS recovery. Through the use of denoisers, rather than explicit models or penalty functionals, our algorithm outperforms these methods on standard test images.

An additional reconstruction algorithm does not fit into any of the above categories but in many ways relates closely to our own. Egiazarian et al. developed a denoising-based CS recovery algorithm that uses the same research group’s BM3D denoising algorithm to impose a non-parametric model on the reconstructed signal. This method solves the CS problem when the measurement matrix is a subsampled DFT matrix. The method iteratively adds noise to the missing part of the spectra and then applies BM3D to the result. In it was shown that the BM3D-based algorithm performed considerably worse than NLR-CS. Therefore it is not tested here.

Finally we should emphasize another major difference between our work and other approaches designed for imaging applications. D-AMP comes with an accurate analysis that explains the behavior of the algorithm, its optimality properties, and its limitations. Such an accurate analysis does not exist for other methods.

I-D Structure of the paper

The remainder of this paper is structured as follows: Section II introduces our D-AMP algorithm and some of its main features. Section III is devoted to the theoretical analysis of D-AMP and its optimality properties. Section IV establishes a connection between our state evolution and existing state evolutions. Section V explains two different approaches to calculating the Onsager correction term. Section VI explains how to smooth poorly behaved denoisers so that they can be used within our framework. Section VII summarizes our main simulation results: it provides evidence on the validity of our state evolution framework; it provides a detailed guideline on setting and tuning of different parameters of the algorithms; it compares the performance of our D-AMP algorithm with the state-of-the-art algorithms in compressive imaging.

II Denoising-based approximate message passing

Consider a family of denoising algorithms $D_{\sigma}$ for a class of signals $\mathcal{C}$. Our goal is to employ these denoisers to obtain a good estimate of $x_{o}\in\mathcal{C}$ from $y=\mathbf{A}x_{o}+w$, where $w\sim N(0,\sigma_{w}^{2}I)$. We start with the following approach that is inspired by the iterative hard-thresholding algorithm and its extensions for block-based compressive imaging . To better understand this approach, consider the noiseless setting in which $y=\mathbf{A}x_{o}$, and assume that the denoiser is a projection onto $C$. The affine subspace defined by $y=\mathbf{A}x$ and the set $C$ are illustrated in Figure 4. We assume that the point $x_{o}$ is the unique point in the intersection of $y=\mathbf{A}x$ and $C$.

We know that the solution lies in the affine subspace $\{x|y=\mathbf{A}x\}$. Therefore, starting from $x^{0}=0$, we move in the direction that is orthogonal to the subspace, i.e., $\mathbf{A}^{*}y$. $\mathbf{A}^{*}y$ is closer to the subspace however it is not necessarily close to $C$. Hence, we employ denoising (or projection in the figure) to obtain an estimate that satisfies the structure of our signal class $C$. After these two steps we obtain $D(\mathbf{A}^{*}y)$. As is also clear in the figure, by repeating these two steps, i.e., moving in the direction of the gradient and then projecting onto $\mathcal{C}$, our estimate may eventually converge to the correct solution $x_{o}$. This leads us to the following iterative algorithm:Note that if $D$ is a projection operator onto $C$ and $C$ is a convex set, then this algorithm is known as projected gradient descent and is known to converge to the correct answer $x_{o}$.

For ease of notation, we have introduced the vector of estimated residual as $z^{t}$. We call this algorithm denoising-based iterative thresholding (D-IT). Note that if we replace $D$ (that was assumed to be projection onto set $C$ in Figure 4) with a denoising algorithm we implicitly assume that $x^{t}+\mathbf{A}^{*}z^{t}$ can be modeled as $x_{o}+v^{t}$, where $v^{t}\sim N(0,(\sigma^{t})^{2}I)$ and is independent of $x_{o}$. Hence, by applying a denoiser we obtain a signal that is closer to $x_{o}$. Unfortunately, as is shown in Figure 5(a) (we will show stronger evidence in Section III-C), this assumption is not true for D-IT. This is the same phenomenon that we observed in Section I-A for iterative soft-thresholding.

Our proposed solution to avoid the non-Gaussianity of the noise in the case of the iterative thresholding algorithms was to employ message passing/approximate message passing. Following the same path, we propose the following message passing algorithm:

Here, $x_{\cdot\rightarrow a}^{t}=[x^{t}_{1\rightarrow a},x^{t}_{2\rightarrow a},\ldots,x^{t}_{n\rightarrow a}]^{T}$ provides an estimate of $x_{o}$. $\hat{\sigma}^{t}$ denotes the standard deviation of the vector

Our empirical findings, summarized in Section VII-C, show that $v_{\cdot\rightarrow a}^{t}$ closely resembles i.i.d. Gaussian noise in high-dimensional settings (both $m$ and $n$ are large). This result has been rigorously proved for a class of scalar denoisers and can also be proved for a class of block-wise denoisers .There are some subtle differences between our claim regarding the Gaussianity of $v_{\cdot\rightarrow a}^{t}$ and the claims presented in other works. Our claim is made in a deterministic setting, while in existing works the Gaussianity claim is made in regards to stochastic settings. This point will be clarified in Section IV.

Despite their advantage in avoiding the non-Gaussianity of the effective noise vector $v$, message passing algorithms have $m$ (number of measurements) different estimates of $x_{o}$; each $x^{t}_{\cdot\rightarrow a}$ is an estimate of $x_{o}$. Similarly, they have $n$ different estimates of the residual $y-\mathbf{A}x_{o}$. The update of all these messages is computationally demanding. Fortunately, if the problem is high dimensional, we can approximate a message passing algorithm’s iterations and obtain the denoising-based approximate message passing algorithm (D-AMP):

The only difference between D-AMP and D-IT is again in the Onsager correction term $z^{t-1}{\rm div}D_{\hat{\sigma}^{t-1}}(x^{t-1}+\mathbf{A}^{*}z^{t-1})/m$. The derivation of D-AMP from the Denoising-based Message Passing (D-MP) algorithm is similar to the derivation of AMP from Message Passing (MP) which can be found in Chapter 5 of . Similar to D-IT, D-AMP relies on the assumption that the effective noise $v^{t}=x^{t}+\mathbf{A}^{*}z^{t}-x_{o}$ resembles i.i.d. Gaussian noise (independent of the signal $x_{o}$) at every iteration. Our empirical findings confirm this assumption: Figure 5(b) displays the effective noise at iteration 5 of D-AMP with BM3D denoising, which will be briefly explained in Section VII-A (we call this algorithm BM3D-AMP). Notice the clearly Gaussian distribution. Based on this observation, and stronger evidence that we will provide in Section VII-C, we conjecture that $v^{t}$ indeed behaves as additive white Gaussian noise for high dimensional problems. The proof of this property is left for future work. In this paper we not only provide strong empirical evidence to support our conjecture, but also explore its theoretical implications.

III Theoretical analysis of D-AMP

The main objective of this section is to characterize the theoretical properties of the D-AMP framework. In this section (and also in our simulations) we start the algorithm with $x^{0}=0$ and $z^{0}=y$. Our analysis is under the high dimensional setting: $n,m$ are very large while $m/n=\delta<1$ is a fixed number. $\delta$ is called the under-determinacy of the system of equations.

III-B Denoiser properties

To analyze D-AMP, we require the denoiser family to be (near) proper, monotone, and Lipschitz continuous (proper and monotone are defined below). Because most denoisers easily satisfy these first two properties, and can be modified to satisfy the third (see Section VI), the requirements do not overly restrict our analysis.

$D_{\sigma}$ is called a proper family of denoisers of level $\kappa$ ($\kappa\in(0,1)$) for the class of signals $C$ if

for every $\sigma>0$. Note that the expectation is with respect to $\epsilon\sim N(0,I)$.

To clarify the above definition, we consider the following examples:

for every $x_{o}\in C$ and every $\sigma^{2}$. Hence, this family of denoisers is proper of level $k/n$.

First note that since the projection onto a subspace is a linear operator and since $P_{C}(x_{o})=x_{o}$ we have

Next we consider a slightly more complicated example that has been popular in signal processing for the last twenty-five years. Let $\Gamma_{k}$ denote the set of $k$-sparse vectors.

Let $\eta(y;\tau\sigma)=(|y|-\tau\sigma)_{+}{\rm sign}(y)$ denote the family of soft-thresholding denoisers. Then

Similar results can be found in other papers including . But since the proof is short and the result is slightly different from similar existing results, we mention the proof here.

For notational simplicity we assume that the first $k$ coordinates of $x_{o}$ are non-zero and the rest are equal to zero.

Note that the optimal threshold $\tau$ to use within soft-thresholding depends on the sparsity $k/n$ of the signal being denoised. One can optimize the parameter $\tau$ for every value of $k/n$ and obtain an optimized family of denoisers. Figure 6 displays the level $\kappa$ of the optimized soft-thresholding in terms of $k/n$. Note that for sparse signals ($k/n$ small) soft-thresholding is an effective denoiser and thus $\kappa$ is small.

The previous denoisers both utilized prior knowledge about the structure of the signal (its dimensionality and its sparsity) in order to denoise $x_{o}$. When nothing is known about $x_{o}$ a proper denoiser might be too much to ask for. For instance, consider the maximum likelihood estimator.

If $D_{\sigma}(x_{o}+\sigma\epsilon)$ is the maximum likelihood estimate of $x_{o}$ from $x_{o}+\sigma\epsilon$, then

In this example the class of signals we have considered is generic and hence the denoiser cannot employ any specific structure in $x_{o}$.

There are occasions when we want to deal with denoisers that are not proper because of an error/bias term that is independent of the noise level. To deal with scenarios such as these, we introduce the definition near proper.

for every $\sigma>0$. Note that the expectation is with respect to $\epsilon\sim N(0,I)$.

As in Definition 1, the constants $\kappa$ and $B$ determine the quality of the denoiser family. Better denoisers have smaller constants.

for every $x_{o}\in\mathcal{C}_{p}$ and every $\sigma^{2}$. Hence, this family of denoisers is near proper with $\kappa=\frac{k}{n}$ and $B=\frac{(k+1)^{1-2/p}}{n(2/p-1)}$.

Let $\Lambda$ denote the set of indices of the $k$-largest coefficients of $x_{o}$. For a vector $x$, define $x_{\Lambda}$ in the following way: $x_{\Lambda,i}=x_{i}$ if $i\in\Lambda$ and otherwise $x_{\Lambda,i}=0$. Note that $x_{o,\Lambda}$ is the best $k$-term approximation of $x_{o}$. We have

Following the same logic as used in Example 1 we see that

Let $x_{o,(j)}$ denote the $j^{\rm th}$ largest element in absolute value of $x_{o}$. It is clear that $|x_{o,(1)}|\geq|x_{o,(2)}|$, …, $|x_{o,(j-1)}|\geq|x_{o,(j)}|$. Combining this fact with (14) we obtain $j|x_{o,(j)}|^{p}\leq 1$, which in turn implies $|x_{o,(j)}|\leq j^{-1/p}$. Returning to (4), we see that

Substituting (13) and (4) into (4) gives the desired result. ∎

In subsequent sections we assume our signal belongs to a class $C$ for which we have a proper or near proper family of denoisers $D_{\sigma}$. The class and denoiser can be very general. For instance, we may assume $C$ to be the class of natural images and $D_{\sigma}$ to denote the BM3D algorithmWe will review this algorithm briefly in Section VII-A at different noise levels .

We call a denoiser monotone if for every $x_{o}$ its risk function

is a non-decreasing function of $\sigma^{2}$.

We make a few remarks regarding monotone denoisers.

Monotonicity is a natural property to expect from denoisers. Many standard denoisers such as soft-thresholding and group soft-thresholding are monotone if we optimize over the threshold parameter. See Lemma 4.4 in for more information.

If a family of denoisers $D_{\sigma}$ is not monotone, then it is straightforward to construct a new denoiser that outperforms $D_{\sigma}$. Here is a simple proof. Suppose that for $\sigma_{1}<\sigma_{2}$ we have

Then construct a new denoiser for noise level $\sigma_{1}$ in the following way:

In the rest of the paper we consider only monotone denoisers.

III-C State evolution

A key ingredient in our analysis of D-AMP is the state evolution; a series of equations that predict the intermediate MSE of AMP algorithms at each iteration. Here we introduce a new “deterministic” state-evolution to predict the performance of D-AMP. Starting from $\theta^{0}=\frac{\|x_{o}\|_{2}^{2}}{n}$ the state evolution generates a sequence of numbers through the following iterations:

where $(\sigma^{t})^{2}=\frac{\theta^{t}}{\delta}(x_{o},\delta,\sigma_{w}^{2})+\sigma_{w}^{2}$ and the expectation is with respect to $\epsilon\sim N(0,I)$. Note that our notation $\theta^{t+1}(x_{o},\delta,\sigma_{w}^{2})$ is set to emphasize that $\theta^{t}$ may depend on the signal $x_{o}$, the under-determinacy $\delta$, and the measurement noise. Consider the iterations of D-AMP and let $x^{t}$ denote its estimate at iteration $t$. Our empirical findings show that the MSE of D-AMP is predicted accurately by the state evolution. We formally state our finding.

If the D-AMP algorithm starts from $x^{0}=0$, then for large values of $m$ and $n$, state evolution predicts the mean square error of D-AMP, i.e.,

Based on extensive simulations, we believe that this finding is true if the following properties are satisfied: (i) The elements of the matrix $\mathbf{A}$ are i.i.d. Gaussian (or subGaussian) with mean zero and standard deviation $1/m$. (ii) The noise $w$ is also i.i.d. Gaussian. (iii) The denoiser $D$ is Lipschitz continuous.A denoiser is said to be $L$-Lipschitz continuous if for every $x_{1},x_{2}\in C$ we have $\|D(x_{1})-D(x_{2})\|_{2}^{2}\leq L\|x_{1}-x_{2}\|_{2}^{2}.$ Many advanced image denoisers have no closed form expression, thus it is very hard to verify whether or not they are Lipschitz continuous. That said, every advanced denoisers we tested was found to closely follow our state evolution equations (Finding 1), suggesting they are in fact Lipschitz. In Section VI we show examples in which Lipschitz continuity is violated and propose a simple approach for dealing with discontinuous denoisers. In all our simulations the elements of $A$ are i.i.d. Gaussian. The same is true for the elements of $w$.

Figure 7 compares the state evolution predictions of D-AMP (based on the BM3D denoising algorithm ) with the empirical performance of D-AMP and D-IT. As is clear from this figure, the state evolution is accurate for D-AMP but not for D-IT. We have checked the validity of the above finding for the following denoising algorithms: (i) BM3D, (ii) BLS-GSM, (iii) Non-local means, (iv) AMP with soft-wavelet-thresholding. We report some of our simulations on this phenomenon in Section VII-C. We have posted our code onlinehttp://dsp.rice.edu/software/DAMP-toolbox to enable other researchers to verify our findings in more general settings and explore the validity of this conjecture on a wider range of denoisers.

In the following sections we assume that the state evolution is accurate for D-AMP and derive some of the main features of D-AMP based on this assumption.

III-D Analysis of D-AMP in the absence of measurement noise

In this section we consider the noiseless setting $\sigma_{w}^{2}=0$ and characterize the number of measurements D-AMP requires (under the validity of the state evolution framework) to recover the signal $x_{o}$ exactly. We consider monotone denoisers, as defined in section III-B. Consider the state evolution equation under the noiseless setting $\sigma_{w}^{2}=0$:

where $(\sigma^{t})^{2}=\frac{\theta^{t}(x_{o},\delta,0)}{\delta}$. Starting with $\theta^{0}(x_{o},\delta,0)=\frac{\|x_{o}\|_{2}^{2}}{n}$, depending on the value of $\delta$ there are two conceivable scenarios for the state evolution equation:

$\theta^{t}(x_{o},\delta,0)\rightarrow 0$ as $t\rightarrow\infty$.

$\theta^{t}(x_{o},\delta,0)\nrightarrow 0$ as $t\rightarrow\infty$.

$\theta^{t}(x_{o},\delta,0)\rightarrow 0$ implies the success of D-AMP algorithm, while $\theta^{t}(x_{o},\delta,0)\nrightarrow 0$ implies its failure in recovering $x_{o}$. The main goal of this section is to study the success and failure regions.

For monotone denoisers, if for $\delta_{0}$, $\theta^{t}(x_{o},\delta_{0},0)\rightarrow 0$, then for any $\delta>\delta_{0}$, $\theta^{t}(x_{o},{\delta},0)\rightarrow 0$ as well.

Define $(\sigma^{t})^{2}=\frac{\theta^{t}(x_{o},\delta,\sigma_{w}^{2})}{\delta}$. Clearly, since $\theta^{t}(x_{o},\delta_{0},\sigma_{w}^{2})\rightarrow 0$ so does $\sigma^{t}$. Our first claim is that for every $\sigma^{2}<\frac{\|x_{o}\|_{2}^{2}}{n\delta_{0}}=(\sigma^{0})^{2}$ (this is where D-AMP is initialized) we have

By using the monotonicity of the denoiser we have for every $\sigma\geq\sigma_{*}$

This (through simple induction) implies that for every $t$,

Furthermore according to the definition of $\sigma_{*}^{2}$ and the fact that $\sigma^{t}>\sigma_{*}$, we have

Therefore, $\sigma^{t+1}$ is a decreasing sequence with lower bound $\sigma_{*}$. Hence, $\sigma^{t}$ converges to $\sigma^{\infty}\geq\sigma_{*}$. The last step is to show that $\sigma^{\infty}=\sigma_{*}$. If this is not the case, then $\sigma^{\infty}>\sigma_{*}$. But according the definition of $\sigma_{*}$ and the supposition that $\sigma^{\infty}>\sigma_{*}$, we have

which is a contradiction to $\sigma^{\infty}$ being a fixed point. Hence $\sigma^{\infty}=\sigma_{*}$. Since $\sigma^{\infty}=0$, we conclude that $\sigma_{*}=0$ and we have

Since, $\delta>\delta_{0}$ we can conclude that

Note that for very small values of $\delta$, it is straightforward to see that $\theta^{t}(x_{o},\delta,0)\nrightarrow 0$ as $t\rightarrow\infty$. If we combine this result with Lemma 1 we conclude the following simple result: For small values of $\delta$ D-AMP fails in recovering $x_{o}$. As $\delta$ increases, after a certain value of $\delta$ D-AMP will successfully recover $x_{o}$ from its undersampled measurements. Define

$\delta^{*}(x_{o})$ denotes the minimum number of measurements required for the successful recovery of $x_{o}$. Our goal is to characterize $\delta^{*}(x_{o})$ in terms of the performance (we will clarify what we mean by performance) of the denoising algorithm. However, since the number of measurements $\delta^{*}(x_{o})$ depends on the signal $x_{o}$, a more natural question in the design of a system is the following: How many measurements does D-AMP require to recover every signal $x_{o}\in C$? The following result addresses this question.

Suppose that for signal class $C$ the denoiser $D_{\sigma}$ is proper at level $\kappa$. Then

The proof of this proposition is a simple application of the state evolution equation. Similar to the proof of Lemma 1 define

Also for notational simplicity we use the notation $\sigma^{t}$ instead of $\sigma^{t}(x_{o},\delta,0)$ in the equation below. According to state evolution we have

Hence, if $\delta>\kappa$, then $(\sigma^{t}(x_{o},\delta,0))^{2}\rightarrow 0$ as $t\rightarrow\infty$. ∎

We can apply Proposition 1 to the examples of Section III-B and derive some well-known results, such as the phase transition of AMP with the soft-threshold denoiser .

If our denoiser is only nearly proper, perfect recovery may not be possible. However, we can use the same technique to bound the recovery error of D-AMP.

Let $D_{\sigma}$ denote a near proper family of denoisers with levels $\kappa$ and $B$, as defined in Definition 2. Then, if $\delta>\kappa$, the error of D-AMP is upper bounded by

The proof of this result is much like the one used for proper denoisers. Again define $\sigma^{t}(x_{o},\delta,\sigma_{w}^{2})=\frac{\theta^{t}(x_{o},\delta,\sigma_{w}^{2})}{\delta}$. Using the state evolution and the definition of near proper we have

For $\delta>\kappa$, the limit of this sequence is as follows

Note that the proof techniques employed above was first developed in and was later employed to establish the phase transition of AMP extensions . There are some minor differences between our derivation and the derivations presented in the other papers since we have not adopted the minimax setting.

III-E Noise sensitivity of D-AMP

In Section III-D we considered the performance of D-AMP in the noiseless setting where $\sigma_{w}^{2}=0$. This section will be devoted to the analysis of D-AMP in the presence of the measurement noise. Here we assume that the denoiser is near proper at levels $\kappa$ and $B$, i.e.,

The following result shows that D-AMP is robust to the measurement noise. Let $\theta^{\infty}(x_{o},\sigma_{w}^{2},\delta)$ denote the fixed point of the state evolution equation. Since there is measurement noise, $\theta^{\infty}(x_{o},\sigma_{w}^{2},\delta)\neq 0$, i.e., D-AMP will not recover $x_{o}$ exactly. We define the noise sensitivity of D-AMP as

The following proposition provides an upper bound for the noise sensitivity as a function of the number of measurements and the variance of the measurement noise.

Let $D_{\sigma}$ denote a near proper family of denoisers at levels $\kappa$ and $B$. Then, for $\delta>\kappa$, the noise sensitivity of D-AMP satisfies

Note that $\theta^{\infty}(x_{o},\delta,\sigma_{w}^{2})$ is a fixed point of the state evolution equation and hence it satisfies

where $\sigma^{\infty}(x_{o},\delta,\sigma_{w}^{2})=\sqrt{\theta^{\infty}(x_{o},\delta,\sigma_{w}^{2})/\delta+\sigma_{w}^{2}}$. Therefore,

A simple calculation completes the proof. ∎

Substituting in $B=0$ into the above result gives the noise sensitivity for proper denoisers.

There are several interesting features of this proposition that we would like to emphasize.

The bound we presented in Proposition 2 is a worst case analysis. The bound may be achieved for certain signals in $C$ and certain noise variances. However, for most signals in $C$ and most noise variances D-AMP will perform better than what is predicted by the bound. Figure 8 shows the performance of BM3D-AMP in terms of the standard deviation of the noise.

The technique we employed above was first developed in . The result we derived in Proposition 2 can be considered as a generalization of the result of to much broader class of denoisers.

As an aside, upper and lower bounds were recently derived for the minimax noise sensitivity of any recovery algorithm when the measurement matrix is i.i.d. Gaussian and the compressively sampled signal is sparse . Note that while our results can be applied to sparse signals, they have been derived under much more general setting. In this section we discussed upper bounds on the noise sensitivity. See Section III-G for some preliminary results on the lower bound.

III-F Tuning the parameters of D-AMP

Practical denoisers typically have a few free parameters and the denoisers’ performance relies on the effective tuning of these parameters. One of the simplest examples of a denoiser with parameters is soft-thresholding (introduced in Example 2), for which the threshold can be regarded as a parameter. There exists extensive literature on tuning the free parameters of denoisers . Diverse and powerful algorithms such as SURE (Stein’s Unbiased Risk Estimation) have been proposed for this purpose .

D-AMP can employ any of these tuning schemes. However, once we use a denoising algorithm in the D-AMP framework the problem of tuning the free parameters of the denoiser seems to become dramatically more difficult: to produce good performance from D-AMP the parameters must be tuned jointly across different iterations. To state this challenge we overload our notation of a denoiser to $D_{\sigma,\tau}$, where $\tau$ denotes the denoiser’s parameters. According to this notation the state evolution is given by

where $(\sigma^{t})^{2}=\frac{o^{t}(\tau^{0},\tau^{1},\ldots,\tau^{t-1})}{\delta}+\sigma_{w}^{2}$. Note that we have changed our notation for the state evolution variables to emphasize the dependence of $o^{t+1}$ on the choice of the parameters we pick at at the previous iterations. The first question that we ask is the following: What does the optimality of $\tau^{0},\tau^{1},\ldots,\tau^{t}$ mean? Suppose that the sequence of parameters $\tau^{t}$ is bounded.

A sequence of parameters $\tau_{*}^{1},\ldots,\tau_{*}^{t}$ is called optimal at iteration $t+1$ if

Note that $\tau_{*}^{0},\ldots,\tau_{*}^{t}$ is optimal in the sense that they produce the smallest mean square error D-AMP can achieve after $t$ iterations. This definition was first given in for the AMP algorithm based on soft-thresholding.

It seems from our formulation that we should solve a joint optimization on $\tau^{0},\ldots,\tau^{t}$ to obtain the optimal values of these parameters. However, it turns out that in D-AMP the optimal parameters can be found much more easily. Consider the following greedy algorithm for setting the parameters:

Tune $\tau^{0}$ such that $o^{1}(\tau^{0})$ is minimized. Call the optimal value $\tau_{*}^{0}$.

If $\tau^{0},\ldots,\tau^{t-1}$ are set to $\tau^{0}_{*},\ldots,\tau^{t-1}_{*}$, then set $\tau^{t}$ such that it minimizes $o^{t+1}(\tau_{*}^{0},\ldots,\tau^{t-1}_{*},\tau^{t})$.

Note that the above strategy is a greedy parameter selection. The following result proves that in the context of D-AMP this greedy strategy is optimal:

Our proof is based on an induction. According to the first step of our procedure we know that

Now suppose that the claim of the theorem is true for every $t\leq T$. We would like to prove that the result also holds for $t=T+1$, i.e.,

Suppose that it is not true and for $\tau_{o}^{0},\ldots,\tau_{o}^{T}$ we have

where $(\sigma_{*}^{T})^{2}=\frac{o^{T}(\tau_{*}^{0},\ldots,\tau_{*}^{T-1})}{\delta}+\sigma_{w}^{2}$. If we define $(\sigma_{o}^{T})^{2}=\frac{o^{T}(\tau_{o}^{0},\ldots,\tau_{o}^{T-1})}{\delta}+\sigma_{W}^{2}$, then according to the induction assumption $\sigma_{*}^{T}\leq\sigma_{o}^{T}$. Therefore, according to the monotonicity of the denoiser

This is in contradiction with (21). Hence,

To summarize the above discussion, greedy parameter tuning is optimal for D-AMP, thus the tuning of D-AMP is as simple (or as difficult) as the tuning of the denoising algorithm that is employed in D-AMP. Many researchers in the area of signal denoising have optimized the parameters of state-of-the-art denoisers, such as BM3D. Lemma 3 implies that optimally tuned denoisers will induce the best possible performance from D-AMP. Therefore, the tuning of D-AMP has already been thoroughly addressed in the denoising literature.

III-G Optimality of D-AMP

D-AMP is a framework by which to employ denoisers to solve linear inverse problems. But is D-AMP optimal? In other words, given a family of denoisers, $D_{\sigma}$, for a set $C$, can we come up with an algorithm for recovering $x_{o}$ from $y=\mathbf{A}x_{o}+w$ that outperforms D-AMP? Note that this problem is ill-posed in the following sense: the denoising algorithm might not capture all the structures that are present in the signal class $C$. Hence, a recovery algorithm employs extra structures not used by the denoiser (and thus not used by D-AMP) clearly might outperform D-AMP. In the following sections we use two different approaches to analyze the optimality of D-AMP.

III-G2 Uniform optimality

Let $\mathcal{{E}}_{\kappa}$ denote the set of all classes of signals $C$ for which there exists a family of denoisers $D_{\sigma}^{C}$ that satisfies

We know from Proposition 1 that for any $C\in\mathcal{E}_{k}$, D-AMP recovers all the signals in $C$ from $\delta>\kappa$ measurements.

We now ask our uniform optimality question: Does there exist any other signal recovery algorithm that can recover all the signals in all these classes with fewer measurements than D-AMP? If the answer is affirmative, then D-AMP is sub-optimal in the uniform sense, meaning there exists an approach that outperforms D-AMP uniformly over all classes in $\mathcal{E}_{\kappa}$. The following proposition shows that any recovery algorithm requires at least $m=\kappa n$ measurements for accurate recovery, i.e., D-AMP is optimal in this sense.

If $m^{*}$ denotes the minimum number of measurements required (by any recovery algorithm) for a set $C\in\mathcal{E}_{\kappa}$, then

According to this simple result, D-AMP is optimal for at least certain classes of signals and certain denoisers. Hence, it cannot be uniformly improved.

III-G3 Single class optimality

The uniform optimality framework we introduced above considers a set of signal classes and measures the performance of an algorithm on every class in this set. However, in many applications such as imaging we are interested in the performance of D-AMP on a specific class of signals, such as images. Unfortunately, the uniform optimality framework does not provide any conclusion in such cases. Therefore, in this section we introduce another framework for evaluating the optimality of D-AMP that we call single class optimality.

A family of denoisers $D^{*}_{\sigma}$ is called minimax optimal for D-AMP at noise level $\sigma_{w}^{2}$, if it achieves

Note that according to our definition, the optimal denoiser may depend on both $\sigma_{w}^{2}$ and $\delta$ and it is not necessarily unique. We call the version of D-AMP that employs $D^{*}_{\sigma}$, ${\rm D}^{*}$-AMP.

Armed with this definition, we formally ask the single class optimality question: Can we provide a new algorithm that can recover signals in class $C$ with fewer measurements than ${\rm D}^{*}$-AMP? If negative, it means that if we employ the optimal denoiser for D-AMP algorithm no other algorithm can outperform D-AMP. Unfortunately, we will show that there are signal classes for which D-AMP is not optimal in this sense. Our proof requires the following standard definition from statistics text books :

The minimax risk of a set of signals $C$ at the noise level $\sigma^{2}$ is defined as

where the expected value is with respect to $\epsilon\sim N(0,I)$. If $D^{M}_{\sigma}$ achieves $R_{MM}(C,\sigma^{2})$, then it will be called the family of minimax denoisers for the set $C$ under the square loss.

The family of minimax denoisers for $C$ is a family of optimal denoisers for D-AMP. Furthermore, in order to recover every $x_{o}\in C$, ${\rm D}^{*}$-AMP requires at least $n\kappa_{MM}$ measurements:

Since the proof of this result is slightly more involved, we postpone it to Appendix -A. ∎

Based on this result, we can simplify the single class optimality question: Does there exist any recovery algorithm that can recover every $x_{o}\in C$ from fewer observations than $n\kappa_{MM}$? Unfortunately, the answer is affirmative.

Consider the following extreme example. Let $B^{n}_{k}$ denote the class of signals that consist of $k$ ones and $n-k$ zeros. Define $\rho=k/n$ and let $\phi(z)$ denote the density function of a standard normal random variable.

For very high dimensional problems, there are recovery algorithms that can recover signals in $B_{k}$ accurately from $1$ measurement. On the other hand, $D^{*}$-AMP requires at least $n(\kappa_{MM}-o(1))$ measurement to recover signals from this class, where

The proof of this result is slightly more involved and hence is postponed to Appendix -B. According to this proposition, since $\kappa_{MM}$ is non-zero, the number of measurements $D^{*}$-AMP requires is proportional to the ambient dimension $n$, while the actual number of measurements that is required for recovery is equal to $1$. Hence, in such cases ${\rm D}^{*}$-AMP is sub-optimal.

However, it is also important to note that while D-AMP is sub-optimal for this class, according to Proposition 3 D-AMP is optimal for other classes. Characterizing the classes of signals for which D-AMP is optimal is left as an open direction for future research. Despite this sub-optimality result, we will show in Section VII that D-AMP provides impressive results for the class of natural images and outperforms state-of-the-art recovery algorithms.

III-H Additional miscellaneous properties of D-AMP

This intuitive result is a key feature of D-AMP. We formalize it below.

Let a family of denoisers $D_{\sigma}^{1}$ be a better denoiser than a family $D^{2}_{\sigma}$ for signal $x_{o}$ in the following sense:

Also, let $\theta^{\infty}_{D_{i}}(x_{o},\delta,\sigma_{w}^{2})$ denote the fixed point of state evolution for denoiser $D_{i}$. Then,

The proof of this result is straightforward. Since, the state evolution of $D^{1}$ is uniformly lower than $D^{2}$, its fixed point is lower as well. ∎

III-H2 D-AMP as a regularization technique

Since in many cases $J(x)$ is non-convex and non-differentiable, iterative heuristic methods have been proposed for solving the above optimization problem.Many of these methods solve (24) accurately when $J$ is convex. D-AMP provides another heuristic approach for solving (24). It has two main advantages over the other heuristics: (i) D-AMP can be analyzed by the state evolution theoretically. Hence, we can theoretically predict the number of measurements required and the noise sensitivity of D-AMP. (ii) The performance of most heuristic methods depend on their free parameters. As discussed in Section III-F there are efficient approaches for tuning the parameters of D-AMP optimally. Below we briefly review the application of D-AMP for solving (24).

Assume that there exists a computationally efficient scheme for solving the optimization problem $\chi_{J}(u;\lambda)=\arg\min\frac{1}{2}\|u-x\|_{2}^{2}+\lambda J(x)$. $\chi_{J}(u,\lambda)$ is called the proximal operator for the function $J$. The D-AMP algorithm for solving (24) is given by

Considering $\chi_{J}$ as a denoiser, this algorithm has exactly the same interpretation as our generic D-AMP algorithm. Furthermore, if the explicit calculation of the Onsager correction term is challenging we can employ the Monte Carlo technique that will be discussed in Section V-B.

IV Connection with other state evolutions

In Section III we introduced a new, “deterministic” state evolution (SE) and used it to analyze D-AMP. Here we review this SE and compare it with AMP’s Bayesian SE, which was first introduced in .

The deterministic SE assumes that $x_{o}$ is an arbitrary but fixed vector in $C$. Starting from $\theta^{0}=\frac{\|x_{o}\|_{2}^{2}}{n}$, the deterministic SE generates a sequence of numbers through the following iterations:

where $(\sigma^{t})^{2}=\frac{1}{\delta}\theta^{t}(x_{o},\delta,\sigma_{w}^{2})+\sigma_{w}^{2}$ and $\epsilon\sim N(0,I)$.

IV-B Bayesian state evolution

The Bayesian SE assumes that $x_{o}$ is a vector drawn from a probability density function (pdf) $p_{x}$, where the support of $p_{x}$ is a subset of $C$. Starting from $\bar{\theta}^{0}=\frac{\|x_{o}\|_{2}^{2}}{n}$, the Bayesian SE generates a sequence of numbers through the following iterations:

where $(\bar{\sigma}^{t})^{2}=\frac{1}{\delta}\bar{\theta}^{t}(p_{x},\delta,\sigma_{w}^{2})+\sigma_{w}^{2}$. We have used the notation $\bar{\theta}$ to distinguish the Bayesian SE from its deterministic counterpart. In Definition 6 we presented a definition of the optimal denoiser under the deterministic framework. One can do the same for the Bayesian framework.

While the deterministic and Bayesian SEs are different, we can establish a connection between them by employing standard results in theoretical statistics regarding the connection between the minimax risk and the Bayesian risk. Next section briefly discusses this connection.

IV-C Connection between the two state evolutions

In this section we would like to establish a connection between the fixed points of the Bayes-optimal denoisers and the minimax-optimal denoisers for D-AMP. Let $\bar{\theta}^{\infty}(p_{x},\delta,\sigma_{w}^{2})$ denote the fixed point of the Bayesian SE (26) associated with the family of Bayes-optimal denoisers from Definition 7. Also, let $\theta^{\infty}(x_{o},\delta,\sigma_{w}^{2})$ denote the fixed point of the deterministic SE (16) for the family of minimax denoisers from Definition 5.

Let $\mathcal{P}$ denote the set of all distributions whose support is a subset of $C$. Then,

For an arbitrary family of denoisers $D_{\sigma}$ we have

If we take the minimum with respect to $D_{\sigma}$ on both sides of (27), we obtain the following inequality

This inequality implies that $\bar{\theta}^{\infty}(p_{x},\delta,\sigma_{w}^{2})$ is below the fixed point of the deterministic SE using $D^{MM}$ at $x_{mm}$. Therefore, because $\sup_{x_{o}\in C}\theta^{\infty}(x_{o},\delta,\sigma_{w}^{2})$ will be equal to or above the fixed point of $D^{MM}$ at $x_{mm}$, it will satisfy $\sup_{p_{x}\in\mathcal{P}}\bar{\theta}^{\infty}(p_{x},\delta,\sigma_{w}^{2})\leq\sup_{x_{o}\in C}\theta^{\infty}(x_{o},\delta,\sigma_{w}^{2})$. ∎

Under some general conditions it is possible to prove that

then (29) holds as well. Since we work with square loss in the SE, swapping the infimum and supremum is permitted under quite general conditions on $\mathcal{P}$. For more information, see Appendix A of . If (29) holds, then we can follow similar steps as in the proof of Theorem 2 to prove that under the same set of conditions we can have

In words, the supremum of the fixed point of the Bayesian SE with the Bayes-optimal denoiser is equivalent to the supremum of the fixed point of the deterministic SE with the minimax denoiser.

IV-D Why bother?

Considering that the deterministic and Bayesian SEs look so similar, and under certain conditions have the same suprememums, it is natural to ask why we developed the deterministic SE at all. That is, what is gained by using SE (16) rather than (26)?

The deterministic SE is useful because it enables us to deal with signals with poorly understood distributions. Take, for instance, natural images. To use the Bayesian SE on imaging problems, we would first need to characterize all images according to some generalized, almost assuredly inaccurate, pdf. In contrast, the deterministic SE deals with specific signals, not distributions. Thus, even without knowledge of the underlying distribution, so long as we can come up with representative test signals, we can use the deterministic SE. Because the SE shows up in the parameter tuning, noise sensitivity, and performance guarantees of AMP algorithms, being able to deal with arbitrary signals is invaluable.

V Calculation of the Onsager correction term

So far, we have emphasized that the key to the success of approximate message passing algorithms is the Onsager correction term, $z^{t-1}{\rm div}D_{\hat{\sigma}^{t-1}}(x^{t-1}+\mathbf{A}^{*}z^{t-1})/m$, but we have not yet addressed how one can calculate it for an arbitrary denoiser. In this section we provide some guidelines on the calculation of this term. If the input-output relation of the denoiser is known explicitly, then calculating the divergence, ${\rm div}D(x)$, and thus the Onsager correction term, is usually straightforward.In the context of this work the divergence ${\rm div}D(x)$ is simply the sum of the partial derivatives with respect to each element of $x$, i.e., ${\rm div}D(x)=\sum\limits_{i=1}^{n}\frac{\partial D(x_{i})}{\partial x_{i}}$. We will review some popular denoisers and calculate the corresponding Onsager correction terms in the next section. However, most state-of-the-art denoisers are complicated algorithms for which the input-output relation is not explicitly known. In Section V-B we show that, even without an explicit input-output relationship, we can calculate a good approximation for the Onsager correction term.

Three of the most popular signal classes in the literature are sparse, group-sparse, and low-rank signals (when the signal has a matrix form). The most popular denoisers for these signals are soft-thresholding, block soft-thresholding, and singular value soft-thresholding, respectively. The goal of this section is to derive the Onsager correction term for each of these denoisers. Most of the results mentioned in this section have been derived elsewhere. We summarize these results to help the reader understand the steps involved in explicitly computing the Onsager correction term.

V-A2 Block soft-thresholding

V-A3 Singular value thresholding

in which $\lambda$ is a regularization parameter that can be optimized for the best performance. Again this denoiser can be employed in the D-AMP framework to recover low-rank matrices from their underdetermined set of linear equations. To calculate the Onsager correction term we should compute ${\rm div}\mathbf{SVT}_{\lambda}(\mathbf{\Phi})$. According to the divergence of singular value thresholding is given by

V-B Monte Carlo method

While simple denoisers often yield a closed form for their divergence, high-performance denoisers are often data dependent; making it very difficult to characterize their input-output relation explicitly. Here we explain how a good approximation of the divergence can be obtained in such cases. This method relies on a Monte Carlo technique first developed in . The authors of that work showed that given a denoiser $D_{\sigma,\tau}(x)$, using an i.i.d. random vector $b\sim N(0,I)$, we can estimate the divergence with

The only challenge in using this formula is calculating the expected value. This can be done efficiently using Monte Carlo simulation. We generate $M$ i.i.d., $N(0,I)$ vectors $b^{1},b^{2},\ldots,b^{M}$. For each vector $b^{i}$ we obtain an estimate of the divergence $\widehat{\rm div}^{i}$. We then obtain a good estimate of the divergence by averaging

VI Smoothing a denoiser

The denoiser used within D-AMP can take on almost any form. However, the state evolution predictions are not necessarily accurate if the denoiser is not Lipschitz continuous. This requirement seems to disallow some popular denoisers with discontinuities, such as hard-thresholding. Figure 9 compares the state evolution predictions for the hard thresholding denoiser alongside the actual performance of D-AMP based on hard thresholding; the state evolution predictions fail entirely. One simple idea to resolve this issue is to “smooth” the denoisers. The smoothed version should behave nearly the same as the original denoiser but, because it has no discontinuities, should satisfy the state evolution equations. The concept of smoothing simple denoisers and this process’s effects on the performance of simple denoisers has been analyzed in . Here we explain how smoothing can be applied in practice.

where $d\zeta=d\zeta_{1}d\zeta_{2}\ldots d\zeta_{n}$.

Suppose that ${\eta}(x_{1},x_{2},\ldots,x_{n})$ satisfies the following condition:

Note that this condition implies that $\eta_{i}$ is not growing very fast as $\zeta_{1},\zeta_{2},\ldots,\zeta_{n}\rightarrow\infty$.

where the last equality is due to the fact that the $j^{\rm th}$ element of $b^{j}$ is equal to one. Also, note that for $|\gamma|<1$ we have

where to obtain the last equality we used the fact that all the element of $b_{j}$ except the $j^{\rm th}$ one are zero. Define

It is straightforward to use (31) and check that

It is straightforward to conclude that this derivative is bounded. Proving the continuity of the derivative employs the same lines of reasoning and hence we skip it. ∎

where $h^{i}\sim N(0,r^{2}\mathbf{I})$ for all $i$.

Figure 11 compares the input-output relationship of the hard-thresholding denoiser before and after it has been smoothed using this method. Notice the smoothing process completely removes the discontinuities but otherwise leaves the function intact.

The above discussion does not provide any suggestion on how we should pick the smoothing parameter $r$. In fact, rigorous study of the effect of $r$ in AMP requires the evaluation of the difference $|\frac{\|x^{t}-x_{o}\|_{2}^{2}}{N}-\theta^{t}(x_{o},\delta,\sigma_{w}^{2})|$, in terms of the dimension. We leave it as an open problem for future research. Nevertheless, from a practical perspective $r$ can be considered as just another denoiser parameter. The problem of optimizing denoiser parameters has been extensively studied in the field of image processing .

Figures 9 and 10 demonstrate the benefits of smoothing a denoiser using this approach. Unlike D-AMP using the original hard-thresholding denoiser, D-AMP using the smoothed denoiser closely follows the state evolution. This change is significant because it allows us to take advantage of the theory and tuning strategies developed in Section III. More importantly, Figure 10 illustrates how D-AMP based on smoothed-hard-thresholding dramatically outperforms its discontinuous counterpart.

Before proceeding, we would like to emphasize that the above process is not needed for any of the advanced denoisers that we explored in this paper. We found that advanced denoisers satisfy the state evolution and perform exceptionally in D-AMP without any smoothing. We believe this finding implies they are sufficiently smooth to begin with.

VII Simulation results for imaging applications

To demonstrate the efficacy of the D-AMP framework, we evaluate its performance on imaging applications.

As we have discussed so far, D-AMP employs a denoising algorithm for signal recovery problems. In this section, we briefly review some well-known image denoising algorithms that we would like to use in D-AMP. We later demonstrate that any of these denoisers, as well as many others, can be used within our D-AMP algorithm in order to reconstruct various compressively sampled signals. As we discussed in Section III-H, theory says that if denoising algorithm $\cal{M}$ outperforms denoising algorithm $\cal{N}$, then D-AMP based on $\cal{M}$ will outperform D-AMP based on $\cal{N}$. We will see this behavior in our simulations as well.

Below we represent a noisy image with the vector $f$; $f=x+\sigma z$ where $x$ is the noise-free version of the image, $\sigma$ is the standard deviation of the noise, and the elements of $z$ follow an i.i.d. Gaussian distribution.

Gaussian kernel regression: One of the simplest and oldest denoisers is Gaussian kernel regression, which is implemented via a Gaussian filter. As the name suggests, it simply applies a filter whose coefficients follow a Gaussian distribution to the noisy image. It takes the form:

where $G$ and $\star$ denote the Gaussian kernel and the convolution operator, respectively. The Gaussian filter operates under the model that a signal is smooth. That is, neighboring pixels should have similar values. Note that this assumption is violated on image edges and hence this denoiser tends to over-smooth them. Compared to other approaches Gaussian kernel regression has a very low implementation cost. However, it does not remove noise as well as other denoisers.

Bilateral filter: Similar to kernel regression, the bilateral filter sets each pixel value according to a weighted average of neighboring pixels. However, whereas the Gaussian filter computes weights based on how close to one another two pixels are, the bilateral filter computes weights based on the similarity of the pixel values (in addition to their spatial proximity). The estimate produced by the bilateral filter can be written as

where $f(j)$ is the value of the $j^{th}$ pixel, $\Omega_{i}$ is a search window around pixel $i$, and $h$ is a smoothing parameter set according to the amount of noise in the signal. Note that the bilateral filter tries to avoid averaging together light and dark pixels on opposite sides of an edge. The bilateral filter has generally proven much more effective than the Gaussian filter. However, it fails entirely when a very large amount of noise is present and the denoiser cannot determine which pixels should be alike.

Non-local means (NLM): Non-local means extends the bilateral filter concept of averaging pixels with similar values to pixels with similar neighborhoods. NLM’s original implementation takes the same form as the bilateral filter (38) but with the following weights:

where $N(i)$ represents a patch of pixels neighboring pixel $i$ and $h$ is a smoothing parameter set according to the variance of the noise. Because the true value of a pixel is better reflected by the noisy value of its neighborhood than by just its noisy pixel value, NLM better recognizes which pixels should be alike and thus outperforms the bilateral filter. However, because two pixels on opposite sides of an edge usually have very similar neighborhoods, NLM still produces artifacts around edges.

Wavelet thresholding: Wavelet thresholding denoises natural images by assuming they are sparse in the wavelet domain. It transforms signals into a wavelet basis, thresholds the coefficients, and then inverses the transform. Hence if $\Psi^{\prime}$ and $\Psi$ denote the wavelet transform and its inverse, respectively, then the denoised image is given by

BLS-GSM: Bayes least squares Gaussian scale mixtures extends simple wavelet thresholding by using an overcomplete wavelet basis and computing denoised coefficient values not with a thresholding function, but via a Bayesian least squares estimate. This estimate is computed by considering a neighborhood around every coefficient and then modeling the distribution of the coefficients within that neighborhood as the product of a Gaussian random vector and a random scalar, each with a carefully defined prior. The algorithm uses these priors to compute the expected value of the noiseless coefficient value. Because the distributions of the wavelet coefficients of natural images are highly dependent on one another, a Bayesian least squares estimate can remove noise while retaining far more structure than coefficient thresholding alone. Accordingly, BLS-GSM significantly outperforms wavelet thresholding. Its performance relative to NLM depends on the statistics of the image being denoised.

BM3D: Block matching 3D collaborative filtering can be considered a combination of NLM and wavelet thresholding. The algorithm begins by comparing patches around the pixels in an image and then grouping similar patches into stacks. It then performs 2D and 1D transforms on the group. These transforms are a 2D DCT and a 1D Haar transform or a 2D bi-orthogonal spline wavelet (Bior) and a 1D Haar transform. Which pair is used depends on the amount of noise in the image. Next the algorithm shrinks the coefficients of these groups and performs an inverse transform to estimate each pixel. It performs this process twice; once by hard-thresholding the coefficients and a second time using Wiener filtering based on the spectra of the initial estimate. In practice BM3D significantly outperforms NLM and wavelet thresholding techniques. It does a great job at removing noise and produces fewer artifacts than competing methods. Additionally, the authors of BM3D have provided well optimized code that makes this complicated algorithm quite efficient.

BM3D-SAPCA: BM3D with shape adaptive principal component analysis combines two extensions to the original BM3D algorithm; block matching using shape adaptive patches and thresholding/filtering in a PCA derived basis. Using adaptive patches helps ensure that the algorithm groups only similar patches. The use of an adaptive basis means that features not well captured by the DCT/Bior and Haar bases of BM3D will be retained. The performance of BM3D-SAPCA tends to be incrementally better than BM3D. Unfortunately, this small increase in performance comes at a huge increase in computational cost.

Table I provides a comparison among the above denoising algorithms. The parameters of the Gaussian filter, the bilateral filter, non-local means, and wavelet thresholding were all experimentally tuned so as to maximize PSNR.PSNR stands for peak signal-to-noise ratio and is defined as $10\log_{10}(\frac{255^{2}}{{\rm mean}((\hat{x}-x_{o})^{2})})$ when the pixel range is 0 to 255. It is a measure of how closely a signal estimate $\hat{x}$ is to the true signal $x_{o}$. In this paper we use PSNR to measure both the denoising algorithms’ and CS recovery algorithms’ rescaled MSE. The parameters for the other 3 algorithms were set automatically using their respective packages. The BM3D, BM3D-SAPCA, and BLS-GSM packages are available online.http://www.cs.tut.fi/~foi/GCF-BM3D/ http://www.io.csic.es/PagsPers/JPortilla/software

VII-B Implementation details of D-AMP and D-IT

Our goal is to plug each of the denoising algorithms that we reviewed in Section VII-A into our D-AMP algorithm. In the rest of the paper we use the following terminology: If denoising method $\cal{M}$ is employed in D-AMP, then we call the reconstruction algorithm $\cal{M}$-AMP. For instance, if we use NLM, the resulting algorithm will be called NLM-AMP and if we use BM3D, the resulting algorithm will be called BM3D-AMP.

We use the same terminology for D-IT: If we use the BM3D denoiser then we call the resulting algorithm BM3D-IT.

VII-B2 Denoising parameters

One of the main challenges in comparing different recovery algorithms is the tuning of each algorithm’s free parameters. As discussed in Section III-F, the parameters of D-AMP can be tuned efficiently with a greedy strategy. In other words, at every iteration we optimize the parameters to obtain the minimum MSE at that iteration. Toward this goal, we can employ different strategies that have been proposed in the literature for setting the parameters of denoising algorithms .

A variety of techniques exist to estimate the standard deviation of the noise in an image; however, we tackled this problem by using a convenient feature of AMP algorithms: $||z^{t}||_{2}^{2}/m\approx(\sigma^{t})^{2}$ . Additionally, the packages provided with many of the state-of-the-art denoising algorithms , work with just two inputs; the noisy signal and an estimate of the standard deviation of the Gaussian noise. The packages then tune all other parameters internally so as to minimize the MSE. Thus, for the BM3D, BM3D-SAPCA, and BLS-GSM variants of D-AMP we use $(\hat{\sigma}^{t})^{2}=||z^{t}||_{2}^{2}/m$ along with the packages and skip the parameter tuning problem.

For denoisers without self-tuning packages, such as NLM, the tuning problem is challenging because at early iterations the effective noise has a large standard deviation but at later iterations the effective noise has a small standard deviation. This means the best parameters for early iterations are very different than the best parameters for later iterations. To get around this problem we use look-up-tables to set the parameters according to $\hat{\sigma}^{t}$. We naively generated these tables by first constructing artificial denoising problems with varying amounts of additive white Gaussian noise and then sweeping through the tuning parameters at each noise level. Figure 12 presents how we chose the parameter $h$ used in NLM. At each noise level we simply chose the parameter values that maximized the PSNR of the denoising problem. For example, for NLM our look-up-table set $h$ to .9 for $\hat{\sigma}^{t}$ between 15 and 30. The same parameters were applied to all images; we did not optimize our code for individual images.

Recall that the state evolution comparison (Figure 7) showed that the MSE of BM3D-IT rose as the number of iterations increased. We attribute this to non-Gaussian effective noise and correct for this behavior by over-smoothing BM3D-IT at each iteration. The over-smoothing was set by using parameters optimized for $2\hat{\sigma}$ rather than $\hat{\sigma}$. The scalar 2 was chosen as it provided the best MSE among the scalar values we tested.

VII-B3 Stopping criterion

AMP is typically designed to stop after some number of iterations or when $\frac{\|x^{t}-x^{t-1}\|_{2}}{\|x^{t}\|_{2}}$ is less than a threshold. Figure 13 demonstrates the PSNR evolution of BM3D-AMP (as a function of iterations) for different sampling rates of the $128\times 128$ Barbara test image. As the figure suggests, after about 10 iterations the PSNR has generally approached its maximum, but the variance of the estimates remains very high. After 30 iterations the variance is quite low. Therefore to reduce variation in our results, we decided to run BM3D-AMP for $30$ iterations. The other D-AMP algorithms, as well as D-IT, IST, and AMP, exhibited similar behavior and were also run for 30 iterations.

VII-B4 Onsager correction

In all our implementations of D-AMP (except for the original AMP for which we used the closed form solution) we have used the Monte Carlo method for calculating the Onsager correction term, as reviewed in Section V-B. While the algorithm seems to be insensitive to the exact value of $\epsilon$ and works for a wide range of values of $\epsilon$, we used $\epsilon=\frac{\|x\|_{\infty}}{1000}$. We found this value was small enough for the approximation to be effective while not so small as to result in rounding errors. In the case of the original AMP, we have used the calculations we described in Section V-A.

VII-C State evolution of D-AMP

Because the effective noise within D-AMP iterations is Gaussian, as further illustrated in Figure 14, state evolution serves as an effective predictor of D-AMP’s performance. As the first step in our simulations, we would like to provide evidence of this prediction accuracy. To do so we compare the predicted and observed performance of D-AMP with NLM, wavelet thresholding, BLS-GSM, and BM3D.

Recall that the state evolution of D-AMP is defined by

where $(\sigma^{t})^{2}=\frac{\theta^{t}}{\delta}(x_{o},\delta,\sigma_{w}^{2})+\sigma_{w}^{2}$. To compute this value in practice, at every iteration we added white Gaussian noise with standard deviation $\sigma^{t}$ to $x_{o}$, denoised the signal with denoiser $D_{\sigma^{t}}$ (using the true, rather than estimated, $\sigma^{t}$), and then computed the MSE.

Figure 15 displays the state evolutions alongside the true intermediate MSEs of the four aforementioned denoising-based algorithms when applied to a $\delta=0.4$ sampled $128\times 128$ House test image with no measurement noise. The average true MSEs at iteration 29 of AMP, NLM-AMP, BLS-GSM-AMP, and BM3D-AMP are all within 1.2% of the MSEs predicted by their respective state evolutions. We have posted our code onlinehttp://dsp.rice.edu/software/DAMP-toolbox to enable other researchers to verify and explore the validity of our state evolution predictions for these and other D-AMP algorithms.

VII-D One-dimensional synthetic test

As a simple demonstration of the improvements that can be achieved by employing better denoising algorithms in AMP, we compare the performance of the original AMP (that employs sparsity in the wavelet domain) with the performance of NLM-AMP on a piecewise constant signal. Within the test AMP used a Haar basis for wavelet thresholding and used the max-min optimal threshold as determined by . The Haar basis was chosen because it well captures signal discontinuities. NLM-AMP used a length 11 patch, $|N(i)|=11$, a length 21 search window, $|\Omega_{i}|$=21, and a smoothing parameter of 1.5, $h=1.5$. These settings were chosen because they allow NLM to effectively denoise piecewise constant signals at a variety of noise levels. The results of our simulation are shown in Figure 3. As is clear from the figure, NLM-AMP significantly outperforms the original AMP. Even though the signal is relatively sparse in the wavelet domain, NLM captures its structure far more effectively. Hence NLM-AMP outperforms the standard AMP that employs sparsity in the wavelet domain.

VII-E Imaging tests

In this section we compare the performance of D-AMP, using a variety of denoisers, with other CS reconstruction algorithms. In particular, we compare the performance of our D-AMP algorithm with turbo-AMP http://www2.ece.ohio-state.edu/~schniter/turboAMPimaging/, which is a hidden Markov tree model-based AMP algorithm, and ALSB http://idm.pku.edu.cn/staff/zhangjian/ALSB/ and NLR-CS http://see.xidian.edu.cn/faculty/wsdong/NLR_Exps.htm, which both utilize non-local group-sparsity. NLR-CS represents the current state-of-the-art in CS image reconstruction algorithms. We compare these 3 algorithms to D-AMP based on the NLM, BLS-GSM, BM3D, and BM3D-SAPCA denoisers. The performance of D-AMP using the Gaussian filter and the bilateral filter was not competitive and has been omitted from the results. We include comparisons with AMP, using a wavelet basis. We also include comparisons with the BM3D-IT algorithm to illustrate the importance of the Onsager correction term in the performance of D-AMP. Other D-IT algorithms demonstrated considerably worse performance and are therefore omitted from the results.

VII-E2 Test Settings

ALSB uses rows drawn from a $32^{2}\times 32^{2}$ orthonormalized Gaussian measurement matrix to perform block-based compressed sensing, as described in . All other tests used an $m\times n$ measurement matrix that was generated by first using Matlab’s randn(m,n) command and then normalizing the columns. All simulations were conducted on a 3.16 GHz Xeon quad-core processor with 32GB of memory.

For the AMP algorithm we used Daubechies 4 wavelets as the sparsifying basis and set its threshold optimally according to . The parameters of D-AMP and D-IT were set following the methods described in section VII-B2. All D-IT and D-AMP algorithms were run for 30 iterations. AMP was run for 30 iterations as well. Turbo-AMP was run for 10 iterations. We experimented with running turbo-AMP for 30 iterations but found that this yielded no improvement in performance while nearly tripling the computation time. Because the DCT-sparsity-based iterative soft-thresholding method used to generate an initial estimate in NLR-CS’s provided source code failed for Gaussian measurement matrices, we generated the initial estimates used by NLR-CS by running BM3D-AMP for 8 iterations for noiseless tests and 4 iterations for noisy tests. Only 4 iterations of BM3D-AMP were used during noisy tests because if run for 8 iterations the initial estimates from BM3D-AMP were often better than the final estimates from NLR-CS. Turbo-AMP, ALSB, and NLR-CS were otherwise tested under their default settings.

VII-E3 Image database

The data was generated using six standard image processing images drawn from Javier Portilla’s dataset:http://www.io.csic.es/PagsPers/JPortilla/software Lena, Barbara, Boat, Fingerprint, House, and Peppers. The images each have a pixel range of roughly $0-255$. Each of these images, except the examples presented in Figures 16 and 17, were rescaled to $128\times 128$ for testing. Restricting the tests to $128\times 128$ enabled the entire measurement matrix $\mathbf{A}$ to be stored in memory. We also created a version of D-AMP that does not store $\mathbf{A}$ but instead generates sections of $\mathbf{A}$ as required. This version can handle images of arbitrarily large size but is extremely slow.

VII-E4 Noiseless image recovery

While matching the denoiser to the signal produces impressive results in one-dimensional settings (as summarized in Section VII-D), the results in 2D are even more pronounced. We begin this section with a visual comparison of three algorithms: Figure 16 illustrates the image recovery performance of AMP, NLR-CS, and our BM3D-SAPCA-AMP algorithm. BM3D-SAPCA-AMP outperformed NLR-CS slightly; 29.96 dB vs 29.31 dB. Both of these algorithms dramatically outperformed the wavelet sparsity-based AMP algorithm; 20.07 dB.

We also present a more complete comparison of D-AMP with other algorithms in Table II.Model-CoSaMP and other model-based techniques have not been included in our simulation results. First and foremost these methods were too slow for us to gather data before finishing the report. Additionally, we found they were not competitive: In the original Model-CoSaMP paper the authors reported a RMSE of 11.1 (PSNR of 27.22 dB) from a reconstruction of a $128\times 128$ pepper test image using 5000 Gaussian measurements. By comparison, BM3D-AMP returns a RMSE of 5.1 (PSNR of 33.98 dB) on the same test. As is clear from this table, BM3D-AMP or BM3D-SAPCA-AMP outperform all the other algorithms in a large majority of the tests. In the next section we demonstrate that the denoising based-AMP algorithms perform far better than competing methods when in the presence of measurement noise.

VII-E5 Imaging in the presence of measurement noise

In realistic settings compressive samples are subject to measurement noise. Noisy sampling can be modeled by $y=\mathbf{A}x+w$ where $w$ represents additive white Gaussian noise (AWGN). In Figure 17 we provide a visual comparison between the reconstructions of BM3D-SAPCA-AMP (26.86 dB) and NLR-CS (25.30 dB) in the presence of measurement noise. In Table III we compare the performance of the BM3D variant of D-AMP to NLR-CS and ALSB when varying amounts of measurement noise are present. As one might expect from a denoising-based algorithm, D-AMP was found to be exceptionally robust to noise. It outperformed the other methods in almost all tests and in some tests by as much as 7.4 dB.

VII-E6 Computational complexity

Table IV demonstrates that, depending on the denoiser in use, D-AMP can be quite efficient: The BM3D variant of D-AMP is dramatically faster than NLR-CS and ALSB. The table also illustrates how using different denoisers within D-AMP presents not only a means of capturing different signal models, but also a way to balance performance and run times.

VIII Conclusions

Through extensive testing we have demonstrated that the approximate message passing (AMP) compressed sensing recovery algorithm can be extended to use arbitrary denoisers to great effect. Variations of this denoising-based AMP algorithm (D-AMP) deliver state-of-the-art compressively sampled image recovery performance while maintaining a low computational footprint. Our theoretical results and simulations show that the performance of D-AMP can be predicted accurately by state evolution. We have also proven that the problem of tuning the parameters of D-AMP is no more difficult than the tuning of the denoiser that is used in the algorithm. Finally, we have shown that D-AMP is extremely robust to measurement noise. D-AMP represents a plug and play method to recover compressively sampled signals of arbitrary class; simply choose a denoiser well matched to the signal model and plug it in the AMP framework. Since designing denoising algorithms that employ complicated structures is usually much easier than designing recovery algorithms, D-AMP can benefit many different application areas.

A significant amount of work remains to be done. First and foremost, all of the theory we developed for D-AMP relies upon the assumption that residual signals follow Gaussian distributions. In this paper we supported this assumption with state evolution and QQplot experiments. Theoretical validation of this assumption is left for future research. Likewise, all theory and results have been for i.i.d. Gaussian (or subGaussian) measurement matrices. Extension to other measurement matrices such as Fourier samples is another open direction that is left for future research.

where $\sigma^{2}=\frac{\theta^{\infty}_{D^{M}}(x_{o},\delta,\sigma_{w}^{2})}{\delta}+\sigma_{w}^{2}$. Define

Again for notational simplicity assume that the supremum is achieved at $x_{o}^{m*}$. The following lemma will be useful in our proof. It also has a nice interpretation that we describe after proving it.

If $\theta^{\infty}_{D}(x_{o},\delta,\sigma_{W}^{2})$ denotes the fixed point of the state evolution with denoiser $D$ at signal $x_{o}$, then

This result has an interesting interpretation. The least favorable signal for D-AMP, i.e., the signal that leads to the highest fixed point, is one of the least favorable signals for the denoiser $D_{\sigma}$. While we proved this result for a specific denoiser $D_{\sigma}^{M}$, the proof can be easily extended to any denoiser $D_{\sigma}$.

We may now return to the proof of Proposition 4. Similar to Lemma 5 define

Hence the fixed point of $D_{\sigma}^{M}$ is less than or equal to the fixed point of $D^{*}_{\sigma}$. Hence the proof is complete.

-B Proof of Proposition 5

Let $B_{k}$ denote the class of $k$-sparse signals with zero-one elements. Suppose that we have observed $y=Ax_{o}$ ($x_{o}\in B_{k}$) and the goal is to recover $x_{o}$ from $y$. Consider the following recovery algorithm that is a special form of compressible signal pursuit proposed in :

When the algorithm incorrectly estimates $x_{o}$, $\|\hat{x}_{o}-x_{o}\|_{2}^{2}$ is at most $2k$. Since the error is bounded our result is established. ∎

This is essentially the proof of the second part of the theorem. We now prove the first part of the theorem.

We first prove that the samples we draw from $\underline{\pi}^{*}$ belong to $B_{k}$ with high probability.

We employ the result of step one to derive a lower bound for the minimax risk.

Step (i) is a simple application of Hoeffding inequality. Let $\underline{x}_{o}$ be a sample from this distribution. By using Hoeffding inequality we obtain:

In other words, with very high probability the samples that are generated from $\underline{\pi}^{*}$ belong to $B_{k}$. Set $\gamma=\frac{\sqrt{2}}{n^{1/4}}$, define the event $\mathcal{A}$ as $\|x_{o}\|_{0}\leq k$ and let $\pi^{**}$ denote the distribution of $x_{o}$ conditioned on event $\mathcal{A}$. Note that the support of $\pi^{**}$ is a subset of $B_{k}$. Now we can discuss step (ii), i.e., deriving a lower bound for minimax risk. Since the support of $\pi^{**}$ is a subset of $B_{k}$, for every denoiser $D_{\sigma}$ we have

By taking the infimum over $D_{\sigma}$ from both sides, since the optimal denoiser on the left is the Bayes denoiser, we obtain

Define $\phi(z_{i})=\frac{1}{\sqrt{2\pi}}{\rm e}^{-z_{i}^{2}/2}$ and $\underline{\phi}(z)=\phi(z_{1})\times\phi(z_{2})\times\ldots\times\phi(z_{n})$.

Finally, by the dominated convergence theorem we prove that

References