Regularization by Denoising: Clarifications and New Interpretations

I Introduction

One of the most popular approaches to image recovery is the “variational” approach, where one poses and solves an optimization problem of the form

Based on this framework, they proposed several recovery algorithms (based on steepest descent, ADMM, and fixed-point methods, respectively) that yield state-of-the-art performance in deblurring and super-resolution tasks.

In light of these (negative) results, there remains the question of how to explain/understand the RED algorithms from when used with non-symmetric denoisers. In response, we propose a framework called score-matching by denoising (SMD), which aims to match the “score” (i.e., the gradient of the log-prior) rather than to design any explicit regularizer. We then show tight connections between SMD, kernel density estimation , and constrained minimum mean-squared error (MMSE) denoising. In addition, we provide new interpretations of the RED-ADMM and RED-FP algorithms proposed in , and we propose novel RED algorithms with faster convergence. Inspired by , we show that the RED algorithms seek to satisfy a consensus equilibrium condition that allows a direct comparison to the plug-and-play ADMM algorithms from

The remainder of the paper is organized as follows. In Section II we provide more background on RED and related algorithms such as plug-and-play ADMM . In Section III, we discuss the impact of Jacobian symmetry on RED and test whether this property holds in practice. In Section IV, we propose the SMD framework. In Section V, we present new interpretations of the RED algorithms from and new algorithms based on accelerated proximal gradient methods. In Section VI, we perform an equilibrium analysis of the RED algorithms, and, in Section VII, we conclude.

II Background

For use in the sequel, we briefly discuss the Bayesian maximum a posteriori (MAP) estimation framework . The MAP estimate of $\boldsymbol{x}$ from $\boldsymbol{y}$ is defined as

where $p(\boldsymbol{x}|\boldsymbol{y})$ denotes the probability density of $\boldsymbol{x}$ given $\boldsymbol{y}$. Notice that, from Bayes rule $p(\boldsymbol{x}|\boldsymbol{y})=p(\boldsymbol{y}|\boldsymbol{x})p(\boldsymbol{x})/p(\boldsymbol{y})$ and the monotonically increasing nature of $\ln(\cdot)$, we can write

II-B ADMM

A popular approach to solving (2) is through ADMM , which we now review. Using variable splitting, (2) becomes

Using the augmented Lagrangian, problem (7) can be reformulated as

using Lagrange multipliers (or “dual” variables) $\boldsymbol{p}$ and a design parameter $\beta>0$. Using $\boldsymbol{u}\triangleq\boldsymbol{p}/\beta$, (8) can be simplified to

II-C Plug-and-Play ADMM

Image denoising has been studied for decades (see, e.g., the overviews ), with the result that high performance methods are now readily available. Today’s state-of-the-art denoisers include those based on image-dependent filtering algorithms (e.g., BM3D ) or deep neural networks (e.g., TNRD , DnCNN ). Most of these denoisers are not variational in nature, i.e., they are not based on any explicit regularizer $\lambda\rho(\boldsymbol{x})$.

Leveraging the denoising interpretation of ADMM, Venkatakrishnan, Bouman, and Wolhberg proposed to replace line 3 of Algorithm 1 with a call to a sophisticated image denoiser, such as BM3D, and dubbed their approach Plug-and-Play (PnP) ADMM. Numerical experiments show that PnP-ADMM works very well in most cases. However, when the denoiser used in PnP-ADMM comes with no explicit regularization $\rho(\boldsymbol{x})$, it is not clear what objective PnP-ADMM is minimizing, making PnP-ADMM convergence more difficult to characterize. Similar PnP algorithms have been proposed using primal-dual methods and FISTA in place of ADMM.

Approximate message passing (AMP) algorithms also perform denoising at each iteration. In fact, when $\boldsymbol{A}$ is large and i.i.d. Gaussian, AMP constructs an internal variable statistically equivalent to $\boldsymbol{r}$ in (10) . While the earliest instances of AMP assumed separable denoising (i.e., $[\boldsymbol{f}(\boldsymbol{x})]_{n}=f(x_{n})~{}\forall n$ for some $f$) later instances, like , considered non-separable denoising. The paper by Metzler, Maleki, and Baraniuk proposed to plug an image-specific denoising algorithm, like BM3D, into AMP. Vector AMP, which extends AMP to the broader class of “right rotationally invariant” random matrices, was proposed in , and VAMP with image-specific denoising was proposed in . Rigorous analyses of AMP and VAMP under non-separable denoisers were performed in and , respectively.

II-D Regularization by Denoising (RED)

to use within the variational framework (2). The advantage of using an explicit regularizer is that a wide variety of optimization algorithms can be used to solve (2) and their convergence can be tractably analyzed.

In , numerical evidence is presented to show that image denoisers $\boldsymbol{f}(\cdot)$ are locally homogeneous (LH), i.e.,

If $\nabla\rho_{\textsf{red}}(\boldsymbol{x})=\boldsymbol{x}-\boldsymbol{f}(\boldsymbol{x})$, then any minimizer $\hat{\boldsymbol{x}}$ of the variational objective under quadratic loss,

must yield $\nabla C_{\textsf{red}}(\hat{\boldsymbol{x}})=\boldsymbol{0}$, i.e., must obey

Based on this line of reasoning, Romano et al. proposed several iterative algorithms that find an $\hat{\boldsymbol{x}}$ satisfying the fixed-point condition (15), which we will refer to henceforth as “RED algorithms.”

III Clarifications on RED

In this section, we first show that the gradient expression (13) holds if and only if the denoiser $\boldsymbol{f}(\cdot)$ is LH and has Jacobian symmetry (JS). We then establish that many popular denoisers lack JS, such as the median filter (MF) , non-local means (NLM) , BM3D , TNRD , and DnCNN . For such denoisers, the RED algorithms cannot be explained by $\rho_{\textsf{red}}(\cdot)$ in (11). We also show a more general result: When a denoiser lacks JS, there exists no regularizer $\rho(\cdot)$ whose gradient expression matches (13). Thus, the problem is not the specific form of $\rho_{\textsf{red}}(\cdot)$ in (11) but rather the broader pursuit of explicit regularization.

We first state some definitions and assumptions. In the sequel, we denote the $i$th component of $\boldsymbol{f}(\boldsymbol{x})$ by $f_{i}(\boldsymbol{x})$, the gradient of $f_{i}(\cdot)$ at $\boldsymbol{x}$ by

and the Jacobian of $\boldsymbol{f}(\cdot)$ at $\boldsymbol{x}$ by

When $\boldsymbol{J}$ exists, it can be shown [25, p.216] that $\boldsymbol{J}=J\boldsymbol{f}(\boldsymbol{x})$.

III-B The RED Gradient

We first recall a result that was established in .

Suppose that denoiser $\boldsymbol{f}(\cdot)$ is locally homogeneous. Then $[J\boldsymbol{f}(\boldsymbol{x})]\boldsymbol{x}=\boldsymbol{f}(\boldsymbol{x})$.

Our proof is based on differentiability and avoids the need to define a directional derivative. From (18), we have

where (20) follows from local homogeneity (12). Equation (21) implies that $[J\boldsymbol{f}(\boldsymbol{x})]\boldsymbol{x}=\boldsymbol{f}(\boldsymbol{x})~{}\forall\boldsymbol{x}\in\mathcal{X}$. ∎

We now state one of the main results of this section.

For $\rho_{\textsf{red}}(\cdot)$ defined in (11),

For any $\boldsymbol{x}\in\mathcal{X}$ and $n=1,\dots,N$,

using the definition of $J\boldsymbol{f}(\boldsymbol{x})$ from (17). Collecting $\{\frac{\partial\rho_{\textsf{red}}(\boldsymbol{x})}{\partial x_{n}}\}_{n=1}^{N}$ into the gradient vector (13) yields (22). ∎

Note that the gradient expression (22) differs from (13).

Suppose that the denoiser $\boldsymbol{f}(\cdot)$ is locally homogeneous. Then the RED gradient expression (13) holds if and only if $J\boldsymbol{f}(\boldsymbol{x})=[J\boldsymbol{f}(\boldsymbol{x})]^{\top}$.

If $J\boldsymbol{f}(\boldsymbol{x})=[J\boldsymbol{f}(\boldsymbol{x})]^{\top}$, then the last term in (22) becomes $-\frac{1}{2}[J\boldsymbol{f}(\boldsymbol{x})]\boldsymbol{x}$, which equals $-\frac{1}{2}\boldsymbol{f}(x)$ by Lemma 1, in which case (22) agrees with (13). But if $J\boldsymbol{f}(\boldsymbol{x})\neq[J\boldsymbol{f}(\boldsymbol{x})]^{\top}$, then (22) differs from (13). ∎

III-C Impossibility of Explicit Regularization

For denoisers $\boldsymbol{f}(\cdot)$ that lack Jacobian symmetry (JS), Lemma 3 establishes that the gradient expression (13) does not hold. Yet (13) leads to the fixed-point condition (15) on which all RED algorithms in are based. The fact that these algorithms work well in practice suggests that “$\nabla\rho(\boldsymbol{x})=\boldsymbol{x}-\boldsymbol{f}(\boldsymbol{x})$” is a desirable property for a regularizer $\rho(\boldsymbol{x})$ to have. But the regularization $\rho_{\textsf{red}}(\boldsymbol{x})$ in (11) does not lead to this property when $\boldsymbol{f}(\cdot)$ lacks JS. Thus an important question is:

Does there exist some other regularization $\rho(\cdot)$ for which $\nabla\rho(\boldsymbol{x})=\boldsymbol{x}-\boldsymbol{f}(\boldsymbol{x})$ when $\boldsymbol{f}(\cdot)$ is non-JS?

The following theorem provides the answer.

Suppose that denoiser $\boldsymbol{f}(\cdot)$ has a non-symmetric Jacobian. Then there exists no regularization $\rho(\cdot)$ for which $\nabla\rho(\boldsymbol{x})=\boldsymbol{x}-\boldsymbol{f}(\boldsymbol{x})$.

To apply this result to our problem, we define

Thus, if $J\boldsymbol{f}(\boldsymbol{x})$ is non-symmetric, then $J[\boldsymbol{x}-\boldsymbol{f}(\boldsymbol{x})]=\boldsymbol{I}-J\boldsymbol{f}(\boldsymbol{x})$ is non-symmetric, which means that there exists no $\rho$ for which (29) holds. ∎

Thus, the problem is not the specific form of $\rho_{\textsf{red}}(\cdot)$ in (11) but rather the broader pursuit of explicit regularization. We note that the notion of conservative vector fields was discussed in [27, App. A] in the context of PnP algorithms, whereas here we discuss it in the context of RED.

III-D Analysis of Jacobian Symmetry

The previous sections motivate an important question: Do commonly-used image denoisers have sufficient JS?

For some denoisers, JS can be studied analytically. For example, consider the “transform domain thresholding” (TDT) denoisers of the form

where $\boldsymbol{g}(\cdot)$ performs componentwise (e.g., soft or hard) thresholding and $\boldsymbol{W}$ is some transform, as occurs in the context of wavelet shrinkage , with or without cycle-spinning . Using $g_{n}^{\prime}(\cdot)$ to denote the derivative of $g_{n}(\cdot)$, we have

and so the Jacobian of $\boldsymbol{f}(\cdot)$ is perfectly symmetric.

Another class of denoisers with perfectly symmetric Jacobians are those that produce MAP or MMSE optimal $\hat{\boldsymbol{x}}$ under some assumed prior $\widehat{p_{\text{\sf x}}}$. In the MAP case, $\hat{\boldsymbol{x}}$ minimizes (over $\boldsymbol{x}$) the cost $c(\boldsymbol{x};\boldsymbol{r})=\frac{1}{2\nu}\|\boldsymbol{x}-\boldsymbol{r}\|^{2}-\ln\widehat{p_{\text{\sf x}}}(\boldsymbol{x})$ for noisy input $\boldsymbol{r}$. If we define $\phi(\boldsymbol{r})\triangleq\min_{\boldsymbol{x}}c(\boldsymbol{x};\boldsymbol{r})$, known as the Moreau-Yosida envelope of $-\ln\widehat{p_{\text{\sf x}}}$, then $\hat{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{r})=\boldsymbol{r}-\nu\nabla\phi(\boldsymbol{r})$, as discussed in (See also for insightful discussions in the context of image denoising.) The elements in the Jacobian are therefore $[J\boldsymbol{f}(\boldsymbol{r})]_{n,q}=\frac{\partial f_{n}(\boldsymbol{r})}{\partial r_{q}}=\delta_{n-q}-\nu\frac{\partial^{2}\phi(\boldsymbol{r})}{\partial r_{q}\partial r_{n}}$, and so the Jacobian matrix is symmetric. In the MMSE case, we have that $\boldsymbol{f}(\boldsymbol{r})=\boldsymbol{r}-\nabla\rho_{\text{\sf TR}}(\boldsymbol{r})$ for $\rho_{\text{\sf TR}}(\cdot)$ defined in (52) (see Lemma 4), and so $[J\boldsymbol{f}(\boldsymbol{r})]_{n,q}=\delta_{n-q}-\frac{\partial^{2}\rho_{\text{\sf TR}}(\boldsymbol{r})}{\partial r_{q}\partial r_{n}}$, again implying that the Jacobian is symmetric. But it is difficult to say anything about the Jacobian symmetry of approximate MAP or MMSE denoisers.

Now let us consider the more general class of denoisers

sometimes called “pseudo-linear” . For simplicity, we assume that $\boldsymbol{W}(\cdot)$ is differentiable on $\mathcal{X}$. In this case, using the chain rule, we have

and so the following are sufficient conditions for Jacobian symmetry.

$\boldsymbol{W}(\boldsymbol{x})$ is symmetric $\forall\boldsymbol{x}\in\mathcal{X}$,

$\sum_{i=1}^{N}\frac{\partial w_{ni}(\boldsymbol{x})}{\partial x_{q}}x_{i}=\sum_{i=1}^{N}\frac{\partial w_{qi}(\boldsymbol{x})}{\partial x_{n}}x_{i}~{}\forall\boldsymbol{x}\in\mathcal{X}$.

When $\boldsymbol{W}$ is $\boldsymbol{x}$-invariant (i.e., $\boldsymbol{f}(\cdot)$ is linear) and symmetric, both of these conditions are satisfied. This latter case was exploited for RED in . The case of non-linear $\boldsymbol{W}(\cdot)$ is more complicated. Although $\boldsymbol{W}(\cdot)$ can be symmetrized (see ), it is not clear whether the second condition above will be satisfied.

III-E Jacobian Symmetry Experiments

For denoisers that do not admit a tractable analysis, we can still evaluate the Jacobian of $\boldsymbol{f}(\cdot)$ at $\boldsymbol{x}$ numerically via

where $\boldsymbol{e}_{n}$ denotes the $n$th column of $\boldsymbol{I}_{N}$ and $\epsilon>0$ is small ($\epsilon=1\times 10^{-3}$ in our experiments). For the purpose of quantifying JS, we define the normalized error metric

which should be nearly zero for a symmetric Jacobian.

Table I showsMatlab code for the experiments is available at http://www2.ece.ohio-state.edu/~schniter/RED/index.html. the average value of $e^{J}_{\boldsymbol{f}}(\boldsymbol{x})$ for $17$ different image patchesWe used the center $16\times 16$ patches of the standard Barbara, Bike, Boats, Butterfly, Cameraman, Flower, Girl, Hat, House, Leaves, Lena, Parrots, Parthenon, Peppers, Plants, Raccoon, and Starfish test images. of size $16\times 16$, using denoisers that assumed a noise variance of $25^{2}$. The denoisers tested were the TDT from (30) with the 2D Haar wavelet transform and soft-thresholding, the median filter (MF) with a $3\times 3$ window, non-local means (NLM) , BM3D , TNRD , and DnCNN . Table I shows that the Jacobians of all but the TDT denoiser are far from symmetric.

Jacobian symmetry is of secondary interest; what we really care about is the accuracy of the RED gradient expressions (13) and (22). To assess gradient accuracy, we numerically evaluated the gradient of $\rho_{\textsf{red}}(\cdot)$ at $\boldsymbol{x}$ using

and compared the result to the analytical expressions (13) and (22). Table II reports the normalized gradient error

for the same $\epsilon$, images, and denoisers used in Table I. The results in Table II show that, for all tested denoisers, the numerical gradient $\widehat{\nabla\rho_{\textsf{red}}}(\cdot)$ closely matches the analytical expression for $\nabla\rho_{\textsf{red}}(\cdot)$ from (22), but not that from (13). The mismatch between $\widehat{\nabla\rho_{\textsf{red}}}(\cdot)$ and $\nabla\rho_{\textsf{red}}(\cdot)$ from (13) is partly due to insufficient JS and partly due to insufficient LH, as we establish below.

III-F Local Homogeneity Experiments

Recall that the TDT denoiser has a symmetric Jacobian, both theoretically and empirically. Yet Table II reports a disagreement between the $\nabla\rho_{\textsf{red}}(\cdot)$ expressions (13) and (22) for TDT. We now show that this disagreement is due to insufficient local homogeneity (LH).

To do this, we introduce yet another RED gradient expression,

which results from combining (22) with Lemma 1. Here, $\stackrel{{\scriptstyle\text{\tiny\sf LH}}}{{=}}$ indicates that (38) holds under LH. In contrast, the gradient expression (13) holds under both LH and Jacobian symmetry, while the gradient expression (22) holds in general (i.e., even in the absence of LH and/or Jacobian symmetry). We also introduce two normalized error metrics for LH,

which should both be nearly zero for LH $\boldsymbol{f}(\cdot)$. Note that $e^{\text{\sf LH,1}}_{\boldsymbol{f}}$ quantifies LH according to definition (12) and closely matches the numerical analysis of LH in . Meanwhile, $e^{\text{\sf LH,2}}_{\boldsymbol{f}}$ quantifies LH according to Lemma 1 and to how LH is actually used in the gradient expressions (13) and (38).

The middle row of Table II reports the average gradient error of the gradient expression (38), and Table III reports average LH error for the metrics $e^{\text{\sf LH,1}}_{\boldsymbol{f}}$ and $e^{\text{\sf LH,2}}_{\boldsymbol{f}}$. There we see that the average $e^{\text{\sf LH,1}}_{\boldsymbol{f}}$ error is small for all denoisers, consistent with the experiments in . But the average $e^{\text{\sf LH,2}}_{\boldsymbol{f}}$ error is several orders of magnitude larger (for all but the MF denoiser). We also note that the value of $e^{\text{\sf LH,2}}_{\boldsymbol{f}}$ for BM3D is several orders of magnitude higher than for the other denoisers. This result is consistent with Fig. 2, which shows that the cost function associated with BM3D is much less smooth than that of the other denoisers. As discussed below, these seemingly small imperfections in LH have a significant effect on the RED gradient expressions (13) and (38).

Starting with the TDT denoiser, Table II shows that the gradient error on (38) is large, which can only be caused by insufficient LH. The insufficient LH is confirmed in Table III, which shows that the value of $e^{\text{\sf LH,2}}_{\boldsymbol{f}}(\boldsymbol{x})$ for TDT is non-negligible, especially in comparison to the value for MF.

Continuing with the MF denoiser, Table I indicates that its Jacobian is far from symmetric, while Table III indicates that it is LH. The gradient results in Table II are consistent with these behaviors: the $\nabla\rho_{\textsf{red}}(\boldsymbol{x})$ expression (38) is accurate on account of LH being satisfied, but the $\nabla\rho_{\textsf{red}}(\boldsymbol{x})$ expression (13) is inaccurate on account of a lack of JS.

The results for the remaining denoisers NLM, BM3D, TNRD, and BM3D show a common trend: they have non-trivial levels of both JS error (see Table I) and LH error (see Table III). As a result, the gradient expressions (13) and (38) are both inaccurate (see Table II).

In conclusion, the experiments in this section show that the RED gradient expressions (13) and (38) are very sensitive to small imperfections in LH. Although the experiments in suggested that many popular image denoisers are approximately LH, our experiments suggest that their levels of LH are insufficient to maintain the accuracy of the RED gradient expressions (13) and (38).

III-G Hessian and Convexity

From (26), the $(n,j)$th element of the Hessian of $\rho_{\textsf{red}}(\boldsymbol{x})$ equals

where $\delta_{k}=1$ if $k=0$ and otherwise $\delta_{k}=0$. Thus, the Hessian of $\rho_{\textsf{red}}(\cdot)$ at $\boldsymbol{x}$ equals

This expression can be contrasted with the Hessian expression from [1, (60)], which reads

III-H Example RED-SD Trajectory

For a trajectory $\{\boldsymbol{x}_{k}\}_{k=1}^{K}$ produced by the steepest-descent (SD) RED algorithm from , Fig. 1 plots, versus iteration $k$, the RED Cost $C_{\textsf{red}}(\boldsymbol{x}_{k})$ from (14) and the error on the fixed-point condition (15), i.e., $\|\boldsymbol{g}(\boldsymbol{x}_{k})\|^{2}$ with

For this experiment, we used the $3\times 3$ median-filter for $\boldsymbol{f}(\cdot)$, the Starfish image, and noisy measurements $\boldsymbol{y}=\boldsymbol{x}+\mathcal{N}(\boldsymbol{0},\sigma^{2}\boldsymbol{I})$ with $\sigma^{2}=20$ (i.e., $\boldsymbol{A}=\boldsymbol{I}$ in (14)).

Figure 1 shows that, although the RED-SD algorithm asymptotically satisfies the fixed-point condition (15), the RED cost function $C_{\textsf{red}}(\boldsymbol{x}_{k})$ does not decrease with $k$, as would be expected if the RED algorithms truly minimized the RED cost $C_{\textsf{red}}(\cdot)$. This behavior implies that any optimization algorithm that monitors the objective value $C_{\textsf{red}}(\boldsymbol{x}_{k})$ for, say, backtracking line-search (e.g., the FASTA algorithm ), is difficult to apply in the context of RED.

III-I Visualization of RED Cost and RED-Algorithm Gradient

We now show visualizations of the RED cost $C_{\textsf{red}}(\boldsymbol{x})$ from (14) and the RED algorithm’s gradient field $\boldsymbol{g}(\boldsymbol{x})$ from (46), for various image denoisers. For this experiment, we used the Starfish image, noisy measurements $\boldsymbol{y}=\boldsymbol{x}+\mathcal{N}(\boldsymbol{0},\sigma^{2}\boldsymbol{I})$ with $\sigma^{2}=100$ (i.e., $\boldsymbol{A}=\boldsymbol{I}$ in (14) and (46)), and $\lambda$ optimized over a grid (of $20$ values logarithmically spaced between $0.0001$ and $1$) for each denoiser, so that the PSNR of the RED fixed-point $\hat{\boldsymbol{x}}$ is maximized.

Figure 2 plots the RED cost $C_{\textsf{red}}(\boldsymbol{x})$ and the RED algorithm’s gradient field $\boldsymbol{g}(\boldsymbol{x})$ for the TDT, MF, NLM, BM3D, TNRD, and DnCNN denoisers. To visualize these quantities in two dimensions, we plotted values of $\boldsymbol{x}$ centered at the RED fixed-point $\hat{\boldsymbol{x}}$ and varying along two randomly chosen directions. The figure shows that the minimizer of $C_{\textsf{red}}(\boldsymbol{x})$ does not coincide with the fixed-point $\hat{\boldsymbol{x}}$, and that the RED cost $C_{\textsf{red}}(\cdot)$ is not always smooth or convex.

IV Score-Matching by Denoising

As discussed in Section II-D, the RED algorithms proposed in are explicitly based on gradient rule

This rule appears to be useful, since these algorithms work very well in practice. But Section III established that $\rho_{\textsf{red}}(\cdot)$ from (11) does not usually satisfy (47). We are thus motived to seek an alternative explanation for the RED algorithms. In this section, we explain them through a framework that we call score-matching by denoising (SMD).

As a precursor to the SMD framework, we first propose a technique based on what we will call Tweedie regularization.

Recall the measurement model (10) used to define the “denoising” problem, repeated in (48) for convenience:

To avoid confusion, we will refer to $\boldsymbol{r}$ as “pseudo-measurements” and $\boldsymbol{y}$ as “measurements.” From (48), the likelihood of $\boldsymbol{x}^{0}$ is $p(\boldsymbol{r}|\boldsymbol{x}^{0};\nu)=\mathcal{N}(\boldsymbol{r};\boldsymbol{x}^{0},\nu\boldsymbol{I})$.

Now, suppose that we model the true image $\boldsymbol{x}^{0}$ as a realization of a random vector $\boldsymbol{x}$ with prior pdf $\widehat{p_{\text{\sf x}}}$. We write “$\widehat{p_{\text{\sf x}}}$” to emphasize that the model distribution may differ from the true distribution $p_{\text{\sf x}}$ (i.e., the distribution from which the image $\boldsymbol{x}$ is actually drawn). Under this prior model, the MMSE denoiser of $\boldsymbol{x}$ from $\boldsymbol{r}$ is

and the likelihood of observing $\boldsymbol{r}$ is

We will now define the Tweedie regularizer (TR) as

As we now show, $\rho_{\text{\sf TR}}(\cdot)$ has the desired property (47).

For $\rho_{\text{\sf TR}}(\boldsymbol{r};\nu)$ defined in (52),

where $\hat{\boldsymbol{f}}_{\text{\sf mmse},\nu}(\cdot)$ is the MMSE denoiser from (49).

Equation (53) is a direct consequence of a classical result known as Tweedie’s formula . A short proof, from first principles, is now given for completeness.

where (56) used $\frac{\partial}{\partial r_{n}}\mathcal{N}(\boldsymbol{r};\boldsymbol{x},\nu\boldsymbol{I})=\mathcal{N}(\boldsymbol{r};\boldsymbol{x},\nu\boldsymbol{I})(x_{n}-r_{n})/\nu$. Stacking (59) for $n=1,\dots,N$ in a vector yields (53). ∎

Thus, if the TR regularizer $\rho_{\text{\sf TR}}(\cdot;\nu)$ is used in the optimization problem (14), then the solution $\hat{\boldsymbol{x}}$ must satisfy the fixed-point condition (15) associated with the RED algorithms from , albeit with an MMSE-type denoiser. This restriction will be removed using the SMD framework in Section IV-C.

It is interesting to note that the gradient property (53) holds even for non-homogeneous $\hat{\boldsymbol{f}}_{\text{\sf mmse},\nu}(\cdot)$. This generality is important in applications under which $\hat{\boldsymbol{f}}_{\text{\sf mmse},\nu}(\cdot)$ is known to lack LH. For example, with a binary image $\boldsymbol{x}\in\{0,1\}^{N}$ modeled by $\widehat{p_{\text{\sf x}}}(\boldsymbol{x})=\prod_{n=1}^{N}0.5(\delta(x_{n})+\delta(x_{n}-1))$, the MMSE denoiser takes the form $[\hat{\boldsymbol{f}}_{\text{\sf mmse},\nu}(\boldsymbol{x})]_{n}=0.5+0.5\tanh(x_{n}/\nu)$, which is not LH.

IV-B Tweedie Regularization as Kernel Density Estimation

We now show that TR arises naturally in the data-driven, non-parametric context through kernel-density estimation (KDE) .

Recall that, in most imaging applications, the true prior $p_{\text{\sf x}}$ is unknown, as is the true MMSE denoiser $\boldsymbol{f}_{\text{\sf mmse},\nu}(\cdot)$. There are several ways to proceed. One way is to design “by hand” an approximate prior $\widehat{p_{\text{\sf x}}}$ that leads to a computationally efficient denoiser $\hat{\boldsymbol{f}}_{\text{\sf mmse},\nu}(\cdot)$. But, because this denoiser is not MMSE for $\boldsymbol{x}\sim p_{\text{\sf x}}$, the performance of the resulting estimates $\hat{\boldsymbol{x}}$ will suffer relative to $\boldsymbol{f}_{\text{\sf mmse},\nu}$.

Another way to proceed is to approximate the prior using a large corpus of training data $\{\boldsymbol{x}_{t}\}_{t=1}^{T}$. To this end, an approximate prior could be formed using the empirical estimate

but a more accurate match to the true prior $p_{\text{\sf x}}$ can be obtained using

with appropriately chosen $\nu>0$, a technique known as kernel density estimation (KDE) or Parzen windowing . Note that if $\widetilde{p_{\text{\sf x}}}$ is used as a surrogate for $p_{\text{\sf x}}$, then the MAP optimization problem becomes

with $\rho_{\text{\sf TR}}(\cdot;\nu)$ from (50)-(52) constructed using $\widehat{p_{\text{\sf x}}}$ from (60). In summary, TR arises naturally in the data-driven approach to image recovery when KDE is used to smooth the empirical prior.

IV-C Score-Matching by Denoising

A limitation of the above TR framework is that it results in denoisers $\hat{\boldsymbol{f}}_{\text{\sf mmse},\nu}$ with symmetric Jacobians. (Recall the discussion of MMSE denoisers in Section III-D.) To justify the use of RED algorithms with non-symmetric Jacobians, we introduce the score-matching by denoising (SMD) framework in this section.

Let us continue with the KDE-based MAP estimation problem (62). Note that $\hat{\boldsymbol{x}}$ from (62) zeros the gradient of the MAP optimization objective and thus obeys the fixed-point equation

In principle, $\hat{\boldsymbol{x}}$ in (64) could be found using gradient descent or similar techniques. However, computation of the gradient

is too expensive for the values of $T$ typically needed to generate a good image prior $\widetilde{p_{\text{\sf x}}}$.

A tractable alternative is suggested by the fact that

where $\hat{\boldsymbol{f}}_{\text{\sf mmse},\nu}(\boldsymbol{r})$ is the MMSE estimator of $\boldsymbol{x}\sim\widehat{p_{\text{\sf x}}}$ from $\boldsymbol{r}=\boldsymbol{x}+\mathcal{N}(\boldsymbol{0},\nu\boldsymbol{I})$. In particular, if we can construct a good approximation to $\hat{\boldsymbol{f}}_{\text{\sf mmse},\nu}(\cdot)$ using a denoiser $\boldsymbol{f}_{\boldsymbol{\theta}}(\cdot)$ in a computationally efficient function class $\mathcal{F}\triangleq\{\boldsymbol{f}_{\boldsymbol{\theta}}:\boldsymbol{\theta}\in\boldsymbol{\Theta}\}$, then we can efficiently approximate the MAP problem (62).

This approach can be formalized using the framework of score matching , which aims to approximate the “score” (i.e., the gradient of the log-prior) rather than the prior itself. For example, suppose that we want to want to approximate the score $\nabla\ln\widetilde{p_{\text{\sf x}}}(\cdot;\nu)$. For this, Hyvärinen suggested to first find the best mean-square fit among a set of computationally efficient functions $\boldsymbol{\psi}(\cdot;\boldsymbol{\theta})$, i.e., find

and then to approximate the score $\nabla\ln\widetilde{p_{\text{\sf x}}}(\cdot;\nu)$ by $\boldsymbol{\psi}(\cdot;\hat{\boldsymbol{\theta}})$. Later, in the context of denoising autoencoders, Vincent showed that if one chooses

for some function $\boldsymbol{f}_{\boldsymbol{\theta}}(\cdot)\in\mathcal{F}$, then $\hat{\boldsymbol{\theta}}$ from (68) can be equivalently written as

In this case, $\boldsymbol{f}_{\hat{\boldsymbol{\theta}}}(\cdot)$ is the MSE-optimal denoiser, averaged over $\widehat{p_{\text{\sf x}}}$ and constrained to the function class $\mathcal{F}$.

Note that the denoiser approximation error can be directly connected to the score-matching error as follows. For any denoiser $\boldsymbol{f}_{\boldsymbol{\theta}}(\cdot)$ and any input $\boldsymbol{x}$,

where (71) follows from (66) and (72) follows from (69). Thus, matching the score is directly related to matching the MMSE denoiser.

Plugging the score approximation (69) into the fixed-point condition (64), we get

which matches the fixed-point condition (15) of the RED algorithms from . Here we emphasize that $\mathcal{F}$ may be constructed in such a way that $\boldsymbol{f}_{\boldsymbol{\theta}}(\cdot)$ has a non-symmetric Jacobian, which is the case for many state-of-the-art denoisers. Also, $\boldsymbol{\theta}$ does not need to be optimized for (73) to hold. Finally, $\widehat{p_{\text{\sf x}}}$ need not be the empirical prior (60); it can be any chosen prior . Thus, the score-matching-by-denoising (SMD) framework offers an explanation of the RED algorithms from that holds for generic denoisers $\boldsymbol{f}_{\boldsymbol{\theta}}(\cdot)$, whether or not they have symmetric Jacobians, are locally homogeneous, or MMSE. Furthermore, it suggests a rationale for choosing the regularization weight $\lambda$ and, in the context of KDE, the denoiser variance $\nu$.

IV-D Relation to Existing Work

Tweedie’s formula (53) has connections to Stein’s Unbiased Risk Estimation (SURE) , as discussed in, e.g., [41, Thm. 2] and [42, Eq. (2.4)]. SURE has been used for image denoising in, e.g., . Tweedie’s formula was also used in to interpret autoencoding-based image priors. In our work, Tweedie’s forumula is used to provide an interpretation for the RED algorithms through the construction of the explicit regularizer (52) and the approximation of the resulting fixed-point equation (64) via score matching.

Recently, Alain and Bengio studied the contractive auto-encoders, a type of autoencoder that minimizes squared reconstruction error plus a penalty that tries to make the autoencoder as simple as possible. While previous works such as conjectured that such auto-encoders minimize an energy function, Alain and Bengio showed that they actually minimize the norm of a score (i.e., match a score to zero). Furthermore, they showed that, when the coder and decoder do not share the same weights, it is not possible to define a valid energy function because the Jacobian of the reconstruction function is not symmetric. The results in parallel those in this paper, except that they focus on auto-encoders while we focus on variational image recovery. Another small difference is that uses the small-$\nu$ approximation

whereas we use the exact (Tweedie’s) relationship (53), i.e.,

where is $\widetilde{p_{\text{\sf x}}}$ the “Gaussian blurred” version of $\widehat{p_{\text{\sf x}}}$ from (51).

V Fast RED Algorithms

In , Romano et al. proposed several ways to solve the fixed-point equation (15). Throughout our paper, we have been referring to these methods as “RED algorithms.” In this section, we provide new interpretations of the RED-ADMM and RED-FP algorithms from and we propose new RED algorithms based on accelerated proximal gradient methods.

The ADMM approach was summarized in Algorithm 1 for an arbitrary regularizer $\rho(\cdot)$. To apply ADMM to RED, line 3 of Algorithm 1, known as the “proximal update,” must be specialized to the case where $\rho(\cdot)$ obeys (13) for some denoiser $\boldsymbol{f}(\cdot)$. To do this, Romano et al. proposed the following. Because $\rho(\cdot)$ is differentiable, the proximal solution $\boldsymbol{v}_{k}$ must obey the fixed-point relationship

An approximation to $\boldsymbol{v}_{k}$ can thus be obtained by iterating

over $i=1,\dots,I$ with sufficiently large $I$, initialized at $\boldsymbol{z}_{0}=\boldsymbol{v}_{k-1}$. This procedure is detailed in lines 3-6 of Algorithm 2. The overall algorithm is known as RED-ADMM.

V-B Inexact RED-ADMM

Algorithm 2 gives a faithful implementation of ADMM when the number of inner iterations, $I$, is large. But using many inner iterations may be impractical when the denoiser is computationally expensive, as in the case of BM3D or TNRD. Furthermore, the use of many inner iterations may not be necessary.

For example, Fig. 3 plots PSNR trajectories versus runtime for TNRD-based RED-ADMM with $I=1,2,3,4$ inner iterations. For this experiment, we used the deblurring task described in Section V-G, but similar behaviors can be observed in other applications of RED. Figure 3 suggests that $I=1$ inner iterations gives the fastest convergence. Note that also used $I=1$ when implementing RED-ADMM.

With $I=1$ inner iterations, RED-ADMM simplifies down to the 3-step iteration summarized in Algorithm 3. Since Algorithm 3 looks quite different than standard ADMM (recall Algorithm 1), one might wonder whether there exists another interpretation of Algorithm 3. Noting that line 3 can be rewritten as

we see that the $I=1$ version of inexact RED-ADMM replaces the proximal step with a gradient-descent step under stepsize $1/(\lambda+\beta)$. Thus the algorithm is reminiscent of the proximal gradient (PG) algorithm . We will discuss PG further in the sequel.

V-C Majorization-Minimization and Proximal-Gradient RED

We now propose a proximal-gradient approach inspired by majorization minimization (MM) . As proposed in , we use a quadratic upper-bound,

on the regularizer $\rho(\boldsymbol{x})$, in place of $\rho(\boldsymbol{x})$ itself, at the $k$th algorithm iteration. Note that if $\rho(\cdot)$ is convex and $\nabla\rho(\cdot)$ is $L_{\rho}$-Lipschitz, then $\overline{\rho}(\boldsymbol{x};\boldsymbol{x}_{k})$ “majorizes” $\rho(\boldsymbol{x})$ at $\boldsymbol{x}_{k}$ when $L\geq L_{\rho}$, i.e.,

The majorized objective can then be minimized using the proximal gradient (PG) algorithm (also known as forward-backward splitting) as follows. From (82), note that the majorized objective can be written as

where (86) follows from assuming (47), which is the basis for all RED algorithms. The RED-PG algorithm then alternately updates $\boldsymbol{v}_{k}$ as per the gradient step in (86) and updates $\boldsymbol{x}_{k+1}$ according to the proximal step

as summarized in Algorithm 4. Convergence is guaranteed if $L\geq L_{\rho}$; see for details.

We now show that RED-PG with $L=1$ is identical to the “fixed point” (FP) RED algorithm proposed in . First, notice from Algorithm 4 that $\boldsymbol{v}_{k}=\boldsymbol{f}(\boldsymbol{x}_{k})$ when $L=1$, in which case

The RED-PG and inexact RED-ADMM-$I\!=\!1$ algorithms show interesting similarities: both alternate a proximal update on the loss with a gradient update on the regularization, where the latter term manifests as a convex combination between the denoiser output and another term. The difference is that RED-ADMM-$I\!=\!1$ includes an extra state variable, $\boldsymbol{u}_{k}$. The experiments in Section V-G suggest that this extra state variable is not necessarily advantageous.

V-D Dynamic RED-PG

Recalling from (86) that $1/L$ acts as a stepsize in the PG gradient step, it may be possible to speed up PG by decreasing $L$, although making $L$ too small can prevent convergence. If $\rho(\cdot)$ was known, then a line search could be used, at each iteration $k$, to find the smallest value of $L$ that guarantees the majorization of $\rho(\boldsymbol{x})$ by $\overline{\rho}(\boldsymbol{x};\boldsymbol{x}_{k})$ . However, with a non-LH or non-JS denoiser, it is not possible to evaluate $\rho(\cdot)$, preventing such a line search.

We thus propose to vary $L_{k}$ (i.e., the value of $L$ at iteration $k$) according to a fixed schedule. In particular, we propose to select $L_{0}$ and $L_{\infty}$, and smoothly interpolate between them at intermediate iterations $k$. One interpolation scheme that works well in practice is summarized in line 3 of Algorithm 5. We refer to this approach as “dynamic PG” (DPG). The numerical experiments in Section V-G suggest that, with appropriate selection of $L_{0}$ and $L_{\infty}$, RED-DPG can be significantly faster than RED-FP.

V-E Accelerated RED-PG

Another well-known approach to speeding up PG is to apply momentum to the $\boldsymbol{v}_{k}$ term in Algorithm 4 , often known as “acceleration.” An accelerated PG (APG) approach to RED is detailed in Algorithm 6. There, the momentum in line 5 takes the same form as in FISTA . The numerical experiments in Section V-G suggest that RED-APG is the fastest among the RED algorithms discussed above.

By leveraging the principle of vector extrapolation (VE) , a different approach to accelerating RED algorithms was recently proposed in . Algorithmically, the approach in is much more complicated than the PG-DPG and PG-APG methods proposed above. In fact, we have been unable to arrive at an implementation of that reproduces the results in that paper, and the authors have not been willing to share their implementation with us. Thus, we cannot comment further on the difference in performance between our PG-DPG and PG-APG schemes and the one in .

V-F Convergence of RED-PG

Recalling Theorem 1, the RED algorithms do not minimize an explicit cost function but rather seek fixed points of (15). Therefore, it is important to know whether they actually converge to any one fixed point. Below, we use the theory of non-expansive and $\alpha$-averaged operators to establish the convergence of RED-PG to a fixed point under certain conditions.

First, an operator $\boldsymbol{B}(\cdot)$ is said to be non-expansive if its Lipschitz constant is at most $1$ . Next, for $\alpha\in(0,1)$, an operator $\boldsymbol{P}(\cdot)$ is said to be $\alpha$-averaged if

for some non-expansive $\boldsymbol{B}(\cdot)$. Furthermore, if $\boldsymbol{P}_{1}$ and $\boldsymbol{P}_{2}$ are $\alpha_{1}$ and $\alpha_{2}$-averaged, respectively, then [55, Prop. 4.32] establishes that the composition $\boldsymbol{P}_{2}\circ\boldsymbol{P}_{1}$ is $\alpha$-averaged with

Recalling RED-PG from Algorithm 4, let us define an operator called $\boldsymbol{T}(\cdot)$ that summarizes one algorithm iteration:

With Lemma 5, we can prove the convergence of RED-PG.

From (94), we have that Algorithm 4 is equivalent to

We note that similar Mann-based techniques were used in to prove the convergence of PnP-based algorithms. Also, we conjecture that similar techniques may be used to prove the convergence of other RED algorithms, but we leave the details to future work. Experiments in Section V-G numerically study the convergence behavior of several RED algorithms with different image denoisers $\boldsymbol{f}(\cdot)$.

V-G Algorithm Comparison: Image Deblurring

We now compare the performance of the RED algorithms discussed above (i.e., inexact ADMM, FP, DPG, APG, and PG) on the image deblurring problem considered in [1, Sec. 6.1]. For these experiments, the measurements $\boldsymbol{y}$ were constructed using a $9\times 9$ uniform blur kernel for $\boldsymbol{A}$ and using AWGN with variance $\sigma^{2}=2$. As stated earlier, the image $\boldsymbol{x}$ is normalized to have pixel intensities in the range $$.

versus iteration $k$ for the starfish test image. In the figure, the proposed RED-DPG and RED-APG algorithms appear significantly faster than the RED-FP and RED-ADMM-$I\!=\!1$ algorithms proposed in . For example, RED-APG reaches PSNR $=30$ in $15$ iterations whereas RED-FP and inexact RED-ADMM-$I=1$ take about $50$ iterations.

verus iteration $k$. All but the RED-APG and RED-ADMM algorithms appear to converge to the solution set of the fixed-point equation (15). The RED-APG and RED-ADMM algorithms appear to approximately satisfy the fixed-point equation (15), but not exactly satisfy (15), since the fixed-point error does not decay to zero.

Figure 6 shows the update distance $\frac{1}{N}\|\boldsymbol{x}_{k}-\boldsymbol{x}_{k-1}\|^{2}$ vs. iteration $k$ for the algorithms under test. For most algorithms, the update distance appears to be converging to zero, but for RED-APG and RED-ADMM it does not. This suggests that the RED-APG and RED-ADMM algorithms are converging to a limit cycle rather than a unique limit point.

Next, we replace the TNRD denoiser with the TDT denoiser from (30) and repeat the previous experiments. For the TDT denoiser, we used a Haar-wavelet based orthogonal discrete wavelet transform (DWT) $\boldsymbol{W}$, with the maximum number of decomposition levels, and a soft-thresholding function $\boldsymbol{g}(\cdot)$ with threshold value $0.001$. Unlike the TNRD denoiser, this TDT denoiser is the proximal operator associated with a convex cost function, and so we know that it is $\frac{1}{2}$-averaged and non-expansive.

Figure 7 shows PSNR versus iteration with TDT denoising. Interestingly, the final PSNR values appear to be nearly identical among all algorithms under test, but more than $1$ dB worse than the values around iteration $20$. Figure 8 shows the fixed-point error vs. iteration for this experiment. There, the errors of most algorithms converge to a value near $10^{-7}$, but then remain at that value. Noting that RED-PG satisfies the conditions of Theorem 2 (i.e., convex loss, non-expansive denoiser, $L>1$), it should converge to a fixed-point of (15). Therefore, we attribute the fixed-point error saturation in Fig. 8 to issues with numerical precision. Figure 9 shows the normalized distance versus iteration with TDT denoising. There, the distance decreases to zero for all algorithms under test.

We emphasize that the proposed RED-DPG, RED-APG, and RED-PG algorithms seek to solve exactly the same fixed-point equation (15) sought by the RED-SD, RED-ADMM, and RED-FP algorithms proposed in . The excellent quality of the RED fixed-points was firmly established in , both qualitatively and quantitatively, in comparison to existing state-of-the-art methods like PnP-ADMM . For further details on these comparisons, including examples of images recovered by the RED algorithms, we refer the interested reader to .

VI Equilibrium View of RED Algorithms

Like the RED algorithms, PnP-ADMM repeatedly calls a denoiser $\boldsymbol{f}(\cdot)$ in order to solve an inverse problem. In , Buzzard, Sreehari, and Bouman show that PnP-ADMM finds a “consensus equilibrium” solution rather than a minimum of any explicit cost function. By consensus equilibrium, we mean a solution $(\hat{\boldsymbol{x}},\hat{\boldsymbol{u}})$ to

We now show that the RED algorithms also find consensus equilibrium solutions, but with $G\neq G_{\textsf{pnp}}$. First, recall ADMM Algorithm 1 with explicit regularization $\rho(\cdot)$. By taking iteration $k\rightarrow\infty$, it becomes clear that the ADMM solutions must satisfy the equilibrium condition (97) with

where we note that $F_{\textsf{admm}}=F_{\textsf{pnp}}$.

The RED-ADMM algorithm can be considered as a special case of ADMM Algorithm 1 under which $\rho(\cdot)$ is differentiable with $\nabla\rho(\boldsymbol{x})=\boldsymbol{x}-\boldsymbol{f}(\boldsymbol{x})$, for a given denoiser $\boldsymbol{f}(\cdot)$. We can thus find $G_{\textsf{red-admm}}(\cdot)$, i.e., the RED-ADMM version of $G(\cdot)$ satisfying the equilibrium condition (97b), by solving the right side of (101) under $\nabla\rho(\boldsymbol{x})=\boldsymbol{x}-\boldsymbol{f}(\boldsymbol{x})$. Similarly, we see that the RED-ADMM version of $F(\cdot)$ is identical to the ADMM version of $F(\cdot)$ from (100). Now, the $\hat{\boldsymbol{x}}=G_{\textsf{red-admm}}(\boldsymbol{v})$ that solves the right side of (101) under differentiable $\rho(\cdot)$ with $\nabla\rho(\boldsymbol{x})=\boldsymbol{x}-\boldsymbol{f}(\boldsymbol{x})$ must obey

which we note is a special case of (15). Continuing, we find that

where $\boldsymbol{I}$ represents the identity operator and $(\cdot)^{-1}$ represents the functional inverse. In summary, RED-ADMM with denoiser $\boldsymbol{f}(\cdot)$ solves the consensus equilibrium problem (97) with $F=F_{\textsf{admm}}$ from (100) and $G=G_{\textsf{red-admm}}$ from (107).

Next we establish an equilibrium result for RED-PG. Defining $\boldsymbol{u}_{k}=\boldsymbol{v}_{k}-\boldsymbol{x}_{k}$ and taking $k\rightarrow\infty$ in Algorithm 4, it can be seen that the fixed points of RED-PG obey (97a) for

Furthermore, from line 3 of Algorithm 4, it can be seen that the RED-PG fixed points also obey

which matches (97b) when $G=G_{\textsf{red-pg}}$ for

Note that $G_{\textsf{red-pg}}=G_{\textsf{red-admm}}$ when $L=\beta/\lambda$.

VI-B Interpreting the RED Equilibria

The equilibrium conditions provide additional interpretations of the RED algorithms. To see how, first recall that the RED equilibrium $(\hat{\boldsymbol{x}},\hat{\boldsymbol{u}})$ satisfies

or an analogous pair of equations involving $F_{\textsf{red-admm}}$ and $G_{\textsf{red-admm}}$. Thus, from (108), (109), and (114a), we have that

Note that (118) is reminiscent of, although in general not equivalent to,

which was discussed as an “alternative” formulation of RED in [1, Sec. 5.2].

Insights into the relationship between RED and PnP-ADMM can be obtained by focusing on the simple case of

where the overall goal of variational image recovery would be the denoising of $\boldsymbol{y}$. For PnP-ADMM, (90) and (98) imply

and so the equilibrium condition (97a) implies

Meanwhile, (99) and the equilibrium condition (97b) imply

In the case that $\lambda=1/\sigma^{2}$, we have the intuitive result that

which corresponds to direct denoising of $\boldsymbol{y}$. For RED, $\hat{\boldsymbol{u}}_{\textsf{red}}$ is algorithm dependent, but $\hat{\boldsymbol{x}}_{\textsf{red}}$ is always the solution to (15), where now $\boldsymbol{A}=\boldsymbol{I}$ due to (120). That is,

Taking $\lambda=1/\sigma^{2}$ for direct comparison to (126), we find

Thus, whereas PnP-ADMM reports the denoiser output $\boldsymbol{f}(\boldsymbol{y})$, RED reports the $\hat{\boldsymbol{x}}$ for which the denoiser residual $\boldsymbol{f}(\hat{\boldsymbol{x}})-\hat{\boldsymbol{x}}$ negates the measurement residual $\boldsymbol{y}-\hat{\boldsymbol{x}}$. This $\hat{\boldsymbol{x}}$ can be expressed concisely as

VII Conclusion

The RED paper proposed a powerful new way to exploit plug-in denoisers when solving imaging inverse-problems. In fact, experiments in suggest that the RED algorithms are state-of-the-art. Although claimed that the RED algorithms minimize an optimization objective containing an explicit regularizer of the form $\rho_{\textsf{red}}(\boldsymbol{x})\triangleq\frac{1}{2}\boldsymbol{x}^{\top}(\boldsymbol{x}-\boldsymbol{f}(\boldsymbol{x}))$ when the denoiser is LH, we showed that the denoiser must also be Jacobian symmetric for this explanation to hold. We then provided extensive numerical evidence that practical denoisers like the median filter, non-local means, BM3D, TNRD, or DnCNN lack sufficient Jacobian symmetry. Furthermore, we established that, with non-JS denoisers, the RED algorithms cannot be explained by explicit regularization of any form.

None of our negative results dispute the fact that the RED algorithms work very well in practice. But they do motivate the need for a better understanding of RED. In response, we showed that the RED algorithms can be explained by a novel framework called score-matching by denoising (SMD), which aims to match the “score” (i.e., the gradient of the log-prior) rather than design any explicit regularizer. We then established tight connections between SMD, kernel density estimation, and constrained MMSE denoising.

On the algorithmic front, we provided new interpretations of the RED-ADMM and RED-FP algorithms proposed in , and we proposed novel RED algorithms with much faster convergence. Finally, we performed a consensus-equilibrium analysis of the RED algorithms that lead to additional interpretations of RED and its relation to PnP-ADMM.

Acknowledgments

The authors thank Peyman Milanfar, Miki Elad, Greg Buzzard, and Charlie Bouman for insightful discussions.

References