The Little Engine that Could: Regularization by Denoising (RED)

Introduction

We open this paper with a bold and possibly controversial statement: To a large extent, removal of zero-mean white additive Gaussian noise from an image is a solved problem in image processing.

Before justifying this statement, let us describe the basic building block that will be the star of this paper: the image denoising engine. From the humble linear Gaussian filter to the recently developed state-of-the-art methods using convolutional neural networks, there is no shortage of denoising approaches. In fact, these algorithms are so widely varied in their definition and underlying structure that a concise description will need to be made carefully. Our story begins with an image x corrupted by zero-mean white additive Gaussian noise,

In our notation, we consider an image as a vector of length $n$ (after lexicographic ordering). In the above description, the noise vector is normally distributed, $\textbf{e}\sim{\cal N}\left(\underline{0},\sigma^{2}\textbf{I}\right)$. In the most general terms, the image denoising engine is a function $f:^{n}\longrightarrow^{n}$ that maps an image y to another image of the same size ${\widehat{\textbf{x}}}=f(\textbf{y})$, with the hope to get as close as possible to the original image x. Ideally, such functions operate on the input image y to remove the deleterious effect of the noise while maintaining edges and textures beneath.

The claim made above about the denoising problem being solved is based on the availability of algorithms proposed in the past decade that can treat this task extremely effectively and stably, getting very impressive results, which also tend to be quite concentrated (see for example the work reported in ). Indeed, these documented results have led researchers to the educated guess that these methods are getting very close to the optimally possible denoising performance . This aligns well with the unspoken observation in our community in recent years that investing more work to improve image denoising algorithms seems to lead to diminishing returns.

While the above may suggest that work on denoising algorithms is turning to a dead-end avenue, a new opportunity emerges from this trend: Seeking ways to leverage the vast progress made on the image denoising front in order to treat other tasks in image processing, bringing their solutions to new heights. One natural path towards addressing this goal is to take an existing and well-performing denoising algorithm, and generalize it to handle a new problem. This has been the logical path that has led to contributions such as , and many others. These papers, and others like them, offer an exhaustive manual adaptation of existing denoising algorithms, carefully re-tailored to handle specific alternative problems. This line of work, while often successful, is quite limited, as it offers no flexibility and no general scheme for diverting image denoising engines to treat new image processing tasks.

Could one offer a more systematic way to exploit the abundance of high-performing image-denoising algorithms to treat a much broader family of problems? The recent work by Venkatakrishnan, Bouman and Wohlberg provides a positive and tantalizing answer to this question, in the form of the Plug-and-Play Prior ($P^{3}$) method . This technique builds on the use of an implicit prior for regularizing general inverse problems. When handling the obtained optimization task via the ADMM optimization scheme , the overall problem decomposes into a sequence of image denoising tasks, coupled with simpler $L_{2}$-regularized inverse problems that are much easier to handle.

While the $P^{3}$ scheme may sound like the perfect answer to our prayers, reality is somewhat more complicated. First, this method is not always accompanied by a clear definition of the objective function, since the regularization being effectively used is only implicit, implied by the denoising algorithm. Indeed, it is not clear at all that there is an underlying objective function behind the $P^{3}$ scheme, if arbitrary denoising engines are used . Second, parameter tuning of the ADMM scheme is a delicate matter, and especially so under a non-provable convergence regime, as is the case when using sophisticated denoising algorithms. Third, being intimately coupled with the ADMM, the $P^{3}$ scheme does not offer easy and flexible ways of replacing the iterative procedure. Because of these reasons, the $P^{3}$ scheme is not a turn-key tool, nor is it free from emotional-involvement. Nevertheless, the $P^{3}$ method has drawn much attention (e.g., ), and rightfully so, as it offers a clear path towards harnessing a given image denoising engine for treating more general inverse problems, just as described above.

Is there a more general alternative to the $P^{3}$ method that could be simpler and more stable? This paper puts forward such a framework, offering a systematic use of such denoising engines for regularization of inverse problems. We term the proposed method “Regularization by Denoising” (RED), relying on a general structured smoothness penalty term harnessed to regularize any desired inverse problem. More specifically, the regularization term we propose in this work is of the following

in which the denoising engine itself is applied on the candidate image x, and the penalty induced is proportional to the inner-product between this image and its denoising residual. This defined smoothness regularization is effectively using an image-adaptive Laplacian, which in turn draws its definition from the arbitrary image denoising engine of choice, $f(\cdot)$. Surprisingly, under mild assumptions on $f(\cdot)$, it is shown that the gradient of the regularization term is manageable, given as the denoising residual, $\textbf{x}-f(\textbf{x})$. Therefore, armed with this regularization expression, we show that any inverse problem can be handled while calling the denoising engine iteratively.

RED, the newly proposed framework, is much more flexible in the choice of the optimization method to use, not being tightly coupled to one specific technique, as in the case of the $P^{3}$ scheme (relying on ADMM). Another key difference w.r.t. the $P^{3}$ method is that our adaptive Laplacian-based regularization functional is explicit, making the overall Bayesian objective function clearer and better defined. RED is capable of incorporating any image denoising algorithm, and can treat general inverse problems very effectively, while resulting in an overall algorithm with very simple structure.

An important advantage of RED over the $P^{3}$ scheme is the flexibility with which one can choose the denoising engine $f(\cdot)$ to plug in the regularization term. While most of the discussion in this paper keeps focusing on White Gaussian Noise (WGN) removal, RED can actually deploy almost any denoising engine. Indeed, we define a set of two mild conditions that $f(\cdot)$ should satisfy, and show that many known denoising methods obey these properties. As an example, in our experiments we show how the median filter can become an effective regularizer. Last but not least, we show that the defined regularization term is a convex function, implying that in most cases, in which the log-likelihood term is convex too, the proposed algorithms are guaranteed to converge to a global optimum solution. We demonstrate this scheme, showing state-of-the-art results in image deblurring and single image super-resolution.

This paper is organized as follows: In the next section we present the background material for this work, discussing the general form of inverse problems as optimization tasks, and presenting the Plug-and-Play Prior scheme. Section 3 focuses on the image denoising engine, defining it and its properties clearly, so as to enable its use in the proposed Laplacian paradigm. Section 4 serves the main part of this work – introducing RED: a new way to use an image denoising engine to handle general structured inverse problems. In Section 5 we analyze the proposed scheme, discussing convexity, an alternative formulation, and a qualitative comparison to the $P^{3}$ scheme. Results on the image deblurring and single-image super-resolution problems are brought in Section 6, demonstrating the strength of the proposed scheme. We conclude the paper in Section 7 with a summary of the open questions that we identify for future work.

Preliminaries

In this section we provide background material that serves as the foundation to this work. We start by presenting the breed of optimization tasks we will work on throughout the paper for handling the inverse problems of interest. We then introduce the $P^{3}$ method and discuss its merits and weaknesses.

Bayesian estimation of an unknown image x given its measured version y uses the posterior conditional probability, $P(\textbf{x}|\textbf{y})$, in order to infer x. The most popular estimator in this regime is the Maximum aposteriori Probability (MAP), which chooses the mode (x for which the maximum probability is obtained) of the posterior. Using Bayes’ rule, this implies that the estimation task is turned into an optimization problem of the form

In the above derivations we exploited the fact that $P(\textbf{y})$ is not a function of x and thus can be omitted. We also used the fact that the $-\log$ function is monotonic decreasing, turning the maximization into a minimization problem.

The term $-\log\{P(\textbf{y}|\textbf{x})\}$ is known as the log-likelihood term, and it encapsulates the probabilistic relationship between the desired image x and the measurements y, under the assumption that the desired image is known. We shall rewrite this term as

The second term in Equation (2.1), $-\log\{P(\textbf{x})\}$, refers to the prior, bringing in the influence of the statistical nature of the unknown. This term is also referred to as the regularization, as it helps in better conditioning the overall optimization task in cases where the likelihood alone cannot lead to a unique or stable solution. We shall rewrite this term as

where $\lambda$ is a scalar that encapsulates the confidence in this term.

What is $\rho(\textbf{x})$ and how is it chosen? This is the holy grail of image processing, with a progressive advancement over the years of better modeling the image statistics and leveraging this for handling various tasks in image processing. Indeed, one could claim that almost everything done in our field surrounds this quest for choosing a proper prior, from the early smoothness prior $\rho(\textbf{x})=\lambda\textbf{x}^{T}\textbf{L}\textbf{x}$ using the classic Laplacian , through total variation and wavelet sparsity , all the way to recent proposals based on patch-based GMM and sparse-representation modeling . Interestingly, the work we report here builds on the surprising comeback of the Laplacian regularization in a much more sophisticated form, as reported in .

Armed with a clear definition of the relation between the measurements and the unknown, and with a trusted prior, the MAP estimation boils down to the optimization problem of the form

This defines a wide family of inverse problems that we aim to address in this work, which includes tasks such as denoising, deblurring, super-resolution, demosaicing, tomographic reconstruction, optical-flow estimation, segmentation, and many other problems. The randomness in these problems is typically due to noise contamination of the measurements, and this could be Gaussian, Laplacian, Gamma-distributed, Poisson, and other noise models.

For completeness of this exposition, we briefly review the $P^{3}$ approach. Aiming to solve the problem posed in Equation (6), the ADMM technique suggests to handle this by variable splitting, leading to the equivalent problem

The constraint is turned into a penalty term, relying on the augmented Lagrangian method (in its scaled dual form ), leading to

where u serves as the Lagrange multiplier vector for the set of constraints. ADMM addresses the resulting problem by updating x, v, and u sequentially in a block-coordinate-descent fashion, leading to the following series of sub-problems:

Update of x: When considering v (and u) as fixed, the term $\rho(\textbf{v})$ is omitted, and our task becomes

which is a far simpler inverse problem, where the regularization is an $L_{2}$ proximity one, which is easy to solve in most cases.

Update of v: In this stage we freeze x (and u), and thus the log-likelihood term drops, leading to

This stage is nothing but a denoising of the image $\textbf{x}+\textbf{u}$, assumed to be contaminated by a white additive Gaussian noise of power $\sigma^{2}=1/\beta$. This is easily verified by returning to Equation (6) and plugging the log-likelihood term $\|\textbf{v}-\textbf{x}-\textbf{u}\|_{2}^{2}/2\sigma^{2}$ referring to this case. Indeed, this is the prime observation in , as they suggest to replace the direct solution of (10) by activating an image denoising engine of choice. This way, we do not need to define explicitly the regularization $\rho(\cdot)$ to be used, as it is implied by the engine chosen.

Update of u: We complete the algorithm description by considering the update of the Lagrange multiplier vector u, which is done by ${\widehat{\textbf{u}}}=\textbf{u}+\textbf{x}-\textbf{v}$.

Although the above algorithm has a clear mathematical formulation and only two parameters, denoted by $\beta$ and $\lambda$, it turns out that tuning these is not a trivial task. The source of complexity emerges from the fact that the input noise-level to the denoiser is equal to $\sqrt{\lambda/\beta}$. The confidence in the prior is determined by $\lambda$, and the penalty on the distance between x and v is affected by $\beta$. Empirically, setting a fixed value of $\beta$ does not seize the potential of this algorithm; following previous work (e.g. ), a common practical strategy to achieve a high-quality estimation is to increase the value of $\beta$ as a function of the iterations: Starting from a relatively small value, i.e allowing an aggressive regularization, then proceeding to a more conservative one that limits the smoothing effect, up-to a point where $\beta$ should be large enough to ensure convergence and to avoid an undesired over-smoothed outcome. As one can imagine, it is cumbersome to choose the rate in which $\beta$ should be increased, especially because the corrupted image $\textbf{x}+\textbf{u}$ is a function of the Lagrange multiplier, which varies through the iterations as well.

In terms of convergence, the $P^{3}$ scheme has been shown to be well-behaved under some conditions on the denoising algorithm used. While the work reported in requires the denoiser to be a symmetric smoothing and non-expansive filter, the later work in relaxes this condition to much simpler boundedness of the denoising effect. However, both these prove at best a convergence to a steady-state outcome, which is very far from the desired claim of getting to the global minimizer of the overall objective function. The work reported in offers clear conditions for a global convergence of $P^{3}$, requiring the denoiser to be non-expansive, and emerging as the minimizer of a convex functional. A recently released paper extends the above by using a specific GMM-based denoiser, showing that these two conditions are met, thus guaranteeing global convergence of their ADMM scheme .

Indeed, in that respect, a delicate matter with the $P^{3}$ approach is the fact that given a choice of a denoising engine, it does not necessarily refer to a specific choice of a prior $\rho(\cdot)$, as not every such engine could have a MAP-oriented interpretation. This implies a fundamental difficulty in the $P^{3}$ scheme, as in this case we will be activating a denoising algorithm while departing from the original setting we have defined, and having no underlying cost function to serve. Indeed, the work reported in addresses this very matter in a narrower setting, by studying the identity of the effective prior obtained from a chosen denoising engine. The author chooses to limit the answer to symmetric smoothing filters, showing that even in this special case, the outcome is far from being trivial. As we are about to see in the next section, this shortcoming can be overcome by adopting a different regularization strategy.

The Image Denoising Engine

Image denoising is a special case of the inverse problem posed in Equation (6), referring to the case $\textbf{y}=\textbf{x}+\textbf{e}$, where e is white Gaussian noise contamination of variance $\sigma^{2}$. In this case, the MAP problem becomes

The image denoising engine, which is the focal point of this work, is any candidate solver to the above problem, under a specific choice of a prior. In fact, in this work we choose to widen the definition of the image denoising engine to be any function $f:^{n}\longrightarrow^{n}$ that maps an image y to another image $f(\textbf{y})$ of the same size, and which aims to treat the denoising problem by the operation ${\widehat{\textbf{x}}}=f(\textbf{y})$, be it MAP-based, MMSE-based, or any other approach.

Below, we accompany the definition of a denoiser with few basic conditions on the function $f$. Just before doing so, we make the following broad observation: Among the various degradations that inverse problems come to remedy, removal of noise is fundamentally different. Consider the set of all reasonable “natural” images living on a manifold ${\cal M}$. If we blur any given image or down-scale it, it is still likely to live in ${\cal M}$. However, if the image is contaminated by an additive noise, it pops out of the manifold along the normal to ${\cal M}$ with high probability. Denoising is therefore a fundamental mechanism for an orthogonal “projection” of an image back onto ${\cal M}$In the context of optimization, a smaller class of the general denoising algorithms we define are characterized as “proximal operators” . These operators are in fact direct generalizations of orthogonal projections.. This may explain why denoising is such a central operation, which has been so heavily studied. In the context of this work, in any given step of our iterations, this projection would allow us to project the temporary result back onto ${\cal M}$, so as to increase chances of getting a good-quality restored version of our image.

We pose the following two necessary conditions on $f(\textbf{x})$ that will facilitate our later derivations. Both these conditions rely on the differentiabilityA discussion on this requirement and possible ways to relax it appear in appendix D. of the denoiser $f(\textbf{x}$).

Condition 1: (Local) Homogeneity. A denoiser applied to a positively scaled image $f(c\textbf{x})$ should result in a scaled version of the original image. More specifically, for any scalar $c\geq 0$ we must have $f(c\textbf{x})=cf(\textbf{x})$. In this work we shall relax this condition and demand its satisfaction for $|c-1|\leq\epsilon$ for a very small $\epsilon$.

A direct implication of the above property refers to the behavior of the directional derivative of the denoiser $f(\textbf{x})$ along the direction x. This derivative can be evaluated as

for a very small $\epsilon$. Invoking the homogeneity condition this leads to

Thus, the filter $f(\textbf{x})$ can be written asThis result is sometimes known as Euler’s homogeneous function theorem .

Condition 2: Strong Passivity. The Jacobian $\nabla_{\textbf{x}}f(\textbf{x})$ of the denoising algorithm is stable, satisfying the condition

where $\eta(\textbf{A})$ is the spectral radius of the matrix A. We interpret this condition as the restriction of the denoiser not to magnify the norm of an input image since

Here we have relied on the relation $f(\textbf{x})=\nabla_{\textbf{x}}f(\textbf{x})\textbf{x}$ that has been established in Equation (14). Note that we have chosen this specific condition over the natural weaker alternative, $\|f(\textbf{x})\|\leq\|\textbf{x}\|$, since the strong condition implies the weak one.

A reformulation of the denoising engine that will be found useful throughout this work is the one suggested in , where we assume that the algorithm is built of two phases – a first in which highly non-linear decisions are made, and a second in which these decisions are used to adapt a linear filter to the raw noisy image in order to perform the actual noise removal. Algorithms such as the NLM, kernel-regression, K-SVD, and many others admit this structure, and for them we can write

The matrix W is an $n\times n$ matrix, representing the (pseudo-) linear filter that multiplies the $n\times 1$ noisy image vector y. This matrix is image dependent, as it draws its content from the pixels in y. Nevertheless, this pseudo-linear format provides a convenient description of the denoising engine for our later derivations. We should emphasize that while this notation is true for only some of the denoising algorithms, the proposed framework we outline in this paper is general and applies to any denoising filter that satisfies the two conditions posed above. Indeed, the careful reader will observe that this pseudo-linear form is closely related to the directional derivative relation shown above: $f(\textbf{y})=\nabla_{\textbf{y}}f(\textbf{y})\>\textbf{y}$. In this form, the right hand side is now reminding us of the pseudo-linear form where the matrix $\textbf{W}(\textbf{y})$ is replaced by the Jacobian matrix $\nabla_{\textbf{y}}f(\textbf{y})$.

An interesting consequence of the homogeneity property is the following stability of the pseudo-linear operator W. Starting with a first-order Taylor expansion of $f(\textbf{x}+\textbf{h})$, and invoking the directional derivative relation $f(\textbf{y})=\nabla_{\textbf{y}}f(\textbf{y})\>\textbf{y}$, we get

This result implies that while $\textbf{W}(\textbf{y})=\nabla_{\textbf{y}}f(\textbf{y})$ may indeed depend on y, its sensitivity to its perturbation is negligible, rendering it as an essentially constant linear operator on the perturbed image $\textbf{y}+\textbf{h}$.

2 Denoisers Obeying the Above Conditions

We cannot conclude this section without answering the key question: Which are the denoising engines to which we are constantly referring? While these could include any of the thousands of denoising algorithms published over the years, we obviously focus on the best performing ones, such as the Non-Local Means (NLM) and its advanced variants , the K-SVD denoising method that relies on sparse representation modeling of image patches and its non-local extension , the kernel-regression method that exploits local orientation , the well-known BM3D that combines sparsity and self-similarity of patches , the EPLL scheme that suggests patch-modeling based on the GMM model , CSR and NCSR, which cluster the patches and sparsifies them jointly , the group-Wiener filtering applied on patches , the multi-layer Perceptron method trained to clean an image or the more recent CNN-based alternative called Trainable Nonlinear Reaction-Diffusion (TNRD) algorithm , more recent work that proposed low-rank modeling of patches and the use of the weighted nuclear-norm , non-local sparsity with GSM model , and the list goes on and on. Each and every one of these options (and many others) is a candidate engine that could be fit into our scheme.

A fair and necessary question is whether the denoisers we work with obey the two conditions we have posed above (homogeneity and passivity), and whether the preliminary requirement of differentiability is met. We choose to defer the discussion on the differentiability to Appendix D, due to its relevance to several spread parts of this paper and focus here on the homogeneity and passivity.

Starting with the homogeneity property, we give an experimental evidence, accompanied by a theoretical analysis, to substantiate the fulfillment of this property by a series of well-known denoisers. Figure 1 shows $f((1+\epsilon)\textbf{x})$ versus $(1+\epsilon)f(\textbf{x})$ as a scatter-plot, tested for K-SVD, BM3D, NLM, EPLL, and the TNRD. In all these experiments, the image x is set to be Peppers, $\epsilon=0.01$ and $\sigma=5$ (level of noise assumed within $f$). As can be seen, a tendency to an equality $f((1+\epsilon)\textbf{x})\approx(1+\epsilon)f(\textbf{x})$ is obtained, suggesting that all these are indeed satisfying the homogeneity property. The deviation from exact equality in each of these tests has been evaluated as the standard deviation of the difference $f((1+\epsilon)\textbf{x})-(1+\epsilon)f(\textbf{x})$, leading to $2.95e-4,~{}3.38e-4,~{}1.38e-4,~{}1.46e-4,~{}9.51e-5$, respectively. A further discussion on the homogeneity property from a theoretical perspective is given in Appendix C.

1italic-ϵxf((1+\epsilon)\textbf{x}) versus $(1+\epsilon)f(\textbf{x})$ as a scatter-plot for K-SVD, BM3D, NLM, EPLL, and the TNRD. Equality implies satisfaction of the homogeneity, and the numbers in the brackets provide the STD of the difference. Note that these results were observed on various test images, but shown here for the image Peppers. Turning to the passivity condition, a conceptual difficulty is the need to explicitly obtain the Jacobian of the denoiser in question. Assuming that we overcame this problem somehow and got $\nabla_{\textbf{x}}f(\textbf{x})$, its spectral radius would be evaluated using the Power-Method that applies iterations of the form

The spectral radius itself is easily obtained as

In order to bypass the need to explicitly obtain $\nabla_{\textbf{x}}f(\textbf{x})\textbf{h}_{k}$, we rely on the first order Taylor expansion again,

implying that $\nabla_{\textbf{x}}f(\textbf{x})\cdot\textbf{h}\approx f(\textbf{x}+\textbf{h})-f(\textbf{x})$, which holds true if $\|\textbf{h}\|_{2}$ is small enough. Thus, our alternative Power-Method activates one denoising step per iteration,

The vector $\textbf{h}_{k}$ is normalized in each iteration, and thus $\|\textbf{h}_{k}\|_{2}=1$. This vector is an image, and thus the gray values in it must be very small ($\textbf{h}_{k}(j)\ll 1$), so as to lead to a sum of squares to be equal to 1. This agrees with the need for the perturbation $\textbf{x}+\textbf{h}_{k}$ to be small.

This algorithm has been applied to K-SVD, BM3D, NLM, EPLL, and the TNRD (x set to be the image Cameraman, $\sigma=5$, number of iterations set to give an accuracy of $1e-5$), resulting all with values smaller or equal to $1$, verifying the passivity of these filters.

Regularization by Denoising (RED)

The new and alternative framework we propose relies on a form of an image-adaptive Laplacian which builds a powerful (empirical) prior that can be used to regularize a variety of inverse problems. As a place to start and motivate this definition, let’s go back to the description of the denoiser given in Equation (17), namelyNote that we conveniently assume that the prior is applied to the clean image x, a matter that will be clarified as we dive into our explanations. $\textbf{W}(\textbf{x})\textbf{x}$. We may think of this pseudo-linear filter as one where a set of coefficients (depending on x) are first computed in the matrix W, and then applied to the image x. From this we can construct the Laplacian form,

This definition by itself is not novel, as it is similar to ideas brought up in a series of recent contributions . This expression relies on using an image-adaptive Laplacian – one that draws its definition from the image itself.

Observing the obtained expression, we note that it can be interpreted as the unnormalized cross-correlation between the image x and its corresponding residual $\textbf{x}-\textbf{W}(\textbf{x})\textbf{x}$. As a prior expression should give low values for likely images, in our case this would be achieved in one of two ways (or their combination):

A small value is obtained for $\rho_{L}(\textbf{x})$ if the residual is very small, implying that the image x serves as a near fixed-point of the denoising engine, $\textbf{x}\approx\textbf{W}(\textbf{x})\textbf{x}$.

A small value is obtained for $\rho_{L}(\textbf{x})$ if the cross-correlation of the residual to the image itself is small, a feature that implies that the residual behaves like white noise, or alternatively, if it does not contain elements from the image itself. Interestingly, this concept has been harnessed successfully by some denoising algorithms such as the Dantzig-Selector and by image denoising boosting techniques . Indeed, enforcing orthogonality between the signal and its treated residual is the underlying force behind the Normal equations in statistical estimation (e.g. Least Squares and Kalman Filtering).

Given the above prior, we return to the general inverse-problem posed in Equation (6), and define our new objective,

The prior expression, while exhibiting a possibly complicated dependency on the unknown x, is well-defined and clear. Nevertheless, an attempt to apply any gradient-based algorithm for solving the above minimization task encounters an immediate problem, due to the need to differentiate $\textbf{W}(\textbf{x})$ with respect to x. We overcome this problem by observing that $\textbf{W}(\textbf{x})\textbf{x}$ is in fact the activation of the image denoising engine on x, i.e., $f(\textbf{x})=\textbf{W}(\textbf{x})\textbf{x}$. This observation inspires the following more general definition of the Laplacian regularizer, which is the prime message of this paper:

This is the Regularization by Denoising (RED) paradigm that this work advocates. In this expression, the residual is defined more generally for any filter $f(\textbf{x})$ even if it can not be written in the familiar (pseudo-)linear form. Note that all the preceding intuition about the meaning of this prior remains intact; namely, the value is low if the cross-correlation between the image and its denoising residual is small, or if the residual itself is small due to x being a fixed point of $f$.

Surprisingly, while this expression is more general, it leads to a better-managed optimization problem due to the careful properties we have outlined in Section 3 on our denoising engines $f$. The overall energy functional to minimize is

and the gradient of this expression is readily available by

Based on our prior assumption regarding the availability of a directional derivative for the denoising engine, the term $\nabla_{\textbf{x}}f(\textbf{x})\textbf{x}$ can be replacedA better approximation can be applied in which we replace $\nabla_{\textbf{x}}f(\textbf{x})\textbf{x}$ by the difference $(f((1+\epsilon)\textbf{x})-f(\textbf{x}))/\epsilon$, but this calls for two activations of the denoising engine per gradient evaluation. by $\nabla_{\textbf{x}}f(\textbf{x})\>\textbf{x}=f(\textbf{x})$, based on Equation (14), implying that the gradient expression is further simplified to be

requiring only one activation of the denoising engine for the gradient evaluation. Interestingly, if we bring back now the pseudo-linear interpretation of the denoising engine, the gradient would be the residual, just as posed above, implying that

Observe that this is a non-trivial derivation of the gradient of the original penalty function posed in Equation (24).

2 Deploying the Denoising Engine for Solving Inverse Problems

Gradient Descent Methods: Given the gradient of the energy function $E(\textbf{x})$, the Steepest-Descent (SD) is the simplest option that can be considered, and it amounts to the update formula

Figure 2 describes this algorithm in more details.

A line-search can be proposed in order to set $\mu$ dynamically per iteration, but this is necessarily more involved. For example, in the case of the Armijo rule, it requires a computation of the above gradient $\textbf{g}_{k}$ and then assessing the energy $E({\widehat{\textbf{x}}}_{k}-\mu\textbf{g}_{k})$ for different values of $\mu$ in a retracting fashion, each of which calling for a computation of the denoising engine once.

One could envision using the Conjugate-Gradient (CG) to speed this method, or better yet, applying the Sequential Subspace Optimization (SESOP) algorithm . SESOP holds the current gradient and the last several update directions as the columns of a matrix $\textbf{V}_{k}$ (referring to the $k^{th}$ iteration), and seeks the best linear combination of these columns as an update direction to the current solution, namely $\textbf{x}_{k+1}=\textbf{x}_{k}+\textbf{V}_{k}\boldsymbol{\alpha}_{k}$. When restricted to have only one column, this reduces to a simple SD with line-search. When using two columns, it has the flavor (and strength) of CG, and when using more columns, this method can lead to much faster convergence in non-quadratic problems. The key points of SESOP are (i) The matrix V is updated easily from one iteration to another by discarding the last direction, bringing in the last one, and adding the new gradient; and (ii) The unknown weights vector $\boldsymbol{\alpha}_{k}$ is low-dimensional, and thus updating it can be done using a Newton method. Naturally, one should evaluate the first and second derivatives of the penalty function w.r.t. $\boldsymbol{\alpha}_{k}$, and these will leverage the relations established above. We shall not dive deeper into this option because it will not be included in our experiments.

One possible shortcoming of the gradient approach (in all its manifestations) is the fact that per activation of the denoising engine, the likelihood is updated rather mildly as a simple step toward the current log-likelihood gradient. This may imply that the overall algorithm will require many iterations to converge. The next two methods propose a way to overcome this limitation, by treating the log-likelihood term more “aggressively”.

ADMM: Addressing the optimization task given in Equation (26), we can imitate the path taken by the $P^{3}$ scheme, and apply variable splitting and ADMM. The steps would follow the description given in Section 2 almost exactly, with one major difference – the prior in our case is explicit, and therefore, the stage termed “update of v” would become

Rather than applying an arbitrary denoising engine to compute ${\widehat{\textbf{v}}}$ as a replacement to the actual minimization, we should target this minimization directly by some iterative scheme. For example, setting the gradient of the above expression to zero leads to the equation

which can be solved iteratively using the fixed-point strategy, by

This means that our approach in this case is computationally more expensive, as it will require several activations of the denoising engine. However, a common approach to speed up the convergence (in terms of runtime) of the ADMM is called “early termination” , suggesting to approximate the solution of the v-update stage. We found this approach useful for our setting, especially because the application of a denoiser is computationally expensive. To this end, we may choose to apply only one iteration of the iterative process described in Equation (33), which amounts to one operation of a denoising algorithm. Figure 3 describes this specific algorithm in more details. If one changes all Part 2 (in Figure 3) with the computation ${\widehat{\textbf{v}}}_{k}=f_{1/\sqrt{\beta}}(\textbf{z}^{*})$, we obtain the $P^{3}$ scheme for the same choice of the denoising engine. While this difference is quite delicate, we should remind the reader that (i) this bridge between the two approaches is valid only when we deploy ADMM on our scheme, and (ii) as opposed to the $P^{3}$ method, our method is guaranteed to converge to the global optimum of the overall penalty function, as will be described hereafter.

We should point out that when using the ADMM, the update of x applies an aggressive inversion of the log-likelihood term, which is followed by the above optimization task. Thus, the shortcoming mentioned above regarding the lack of balance between the treatments given to the likelihood and the prior is mitigated.

Fixed-Point Strategy: An appealing alternative to the above exists, obtained via the fixed-point strategy. As our aim is to find x that nulls the gradient, this could be posed as an implicit equation to be solved directly,

Using the fixed-point strategy, this could be handled by the iterative formula

As an example, in the case of linear degradation model and Gaussian white additive noise, this equation would be

This formula suggests one activation of the denoising per iteration, followed by what seems to be a plain Wiener filtering computationNote that Equation (33) refers to the same problem posed here under the choice $\textbf{H}=\textbf{I}$ and $\beta=1/\sigma^{2}$.. The matrix inversion itself could be done in the Fourier domain for block-circulant H, or iteratively using for example, the Richardson algorithm: Defining

our goal is to solve the linear system $\textbf{A}\textbf{x}=\textbf{b}$. This is achieved by a variant of the SD methodAll this refer to a specific iteration $k$ within which we apply inner iterations to solve the linear system, and thus the use of the different index $j$., $\textbf{x}_{j+1}=\textbf{x}_{j}-\mu(\textbf{A}\textbf{x}_{j}-\textbf{b})=\textbf{x}_{j}-\mu\textbf{e}_{j}$, where we have defined $\textbf{e}_{j}=\textbf{A}\textbf{x}_{j}-\textbf{b}$. By setting the step size to be $\mu_{j}=\textbf{e}_{j}^{T}\textbf{A}\textbf{e}_{j}/\textbf{e}_{j}^{T}\textbf{A}^{T}\textbf{A}\textbf{e}_{j}$, we greedily optimize the potential of each iteration.

Convergence of the above algorithm is guaranteed since

This approach, similarly to the ADMM, has the desired balance mentioned above between the likelihood and the regularization terms, matching the efforts dedicated to both. A pseudo-code describing this algorithm appears in Figure 4.

A basic question that has not been discussed so far is how to set the parameters of $f(\textbf{x})$ in defining the regularization term. More specifically, assuming that the denoising engine depends on one parameter – the noise standard-deviation $\sigma_{f}$ – the question is which value to use. While one could envision using varying values as the iterations progress and the outcome improves, the approach we take in this work is to set this parameter to be a small and fixed value. Our intuition for this choice is the desire to have a clear and fixed regularization term, which in turn implies a clear cost function to work with. Furthermore, the prior we propose should encapsulate in it our desire to get to a final image that is a stable point of such a weak denoising engine, $\textbf{x}\approx f(\textbf{x})$. Clearly, more work is required to better understand the influence of this parameter and its automatic setting.

Analysis

Is our proposed regularization function $\rho_{L}(\textbf{x})$ convex? At first glance, this may seem like too much to expect. Nevertheless, it appears that for reasonably performing denoising engines obeying the conditions posed in Section 3, this is exactly the case. For the function $\rho_{L}(\textbf{x})=\textbf{x}^{T}(\textbf{x}-f(\textbf{x}))$ to be convex, we should demand that the second derivative is a positive semi-definite matrix . We have already seen that the first derivative is simply $\textbf{x}-f(\textbf{x})$, which leads to the conclusion that the second derivative is given by $\textbf{I}-\nabla_{\textbf{x}}f(\textbf{x})$.

As already mentioned earlier, in the context of some algorithms such as the NLM and the K-SVD, this is associated with the Laplacian $\textbf{I}-\textbf{W}(\textbf{x})$, and it is positive semi-definite if W has all its eigenvalues in the rangeWe could in fact allow negative eigenvalues for W, but this is unnatural in the context of denoising. $$. This is indeed the case for the NLM filter , the K-SVD-denoising algorithm , and many other denoising engines.

In the wider context of general image denoising engines, convexity is assured if the Jacobian $\nabla_{\textbf{x}}f(\textbf{x})$ of the denoising algorithm is stable, as indeed required in Condition 2 in Section 3, $\eta(\nabla_{\textbf{x}}f(\textbf{x}))\leq 1$. In this case we have that $\rho_{L}(\cdot)$ is convex, and this implies that if the log-likelihood expression is convex as well, the proposed scheme is guaranteed to converge to the global optimum of our cost function in Equation (6). In this respect the proposed algorithm is superior to the $P^{3}$ scheme in its most general form, which at best is known to get to a stable-point . Furthermore, this result may seem similar to the one posed in , as our two conditions for global convergence are homogeneity and passivity of $f$, while in these papers the requirements are passivity of $f$, along with a convex energy functional whom $f$ minimizes. The second condition – having a convex origin to derive $f(\textbf{x})$ – is in fact more restrictive than demanding homogeneity, as it is unclear which of the known denoisers meet this requirement.

2 An Alternative Prior

In Section 4 we motivated the choice of the proposed prior by the desire to characterize the unknown image x as one that is not affected by the denoising algorithm, namely, $\textbf{x}\approx f(\textbf{x})$. Rather than taking the route proposed, we could have suggested a prior of the form

This prior term makes sense intuitively, being based on the same desire to see the denoising residual being small. Indeed, this choice is somewhat related to the option we chose since

suggesting a symmetrization of our own expression.

In order to understand the deeper meaning of this alternative, we resort again to the pseudo-linear denoisers, for which this prior is nothing but

This means that rather than regularizing with the Laplacian, we do so with its square. While this is a worthy possibility which has been considered in the literature under the term “fourth order regularization” , it is known to be more delicate. We leave this and other possibilities of formulating the regularization with the use of $f(\textbf{x})$ for future work.

3 When is Plug-and-Play-Prior = RED ?

In Section 4 we described the use of ADMM as one of the possible avenues for handling our proposed regularization. When handling the inverse problem posed in Equation (26) with ADMM, we have shown that the only difference between this and the $P^{3}$ scheme resides in the update stage for v. Here we aim to answer the following question: Assuming that the numerical algorithm used is indeed the ADMM, under what conditions would the two methods ($P^{3}$ and ours) become equivalent? The answer to this question resides in the optimization task for updating v, which is a denoising task. Thus, purifying this question, our goal is to find conditions on $f(\cdot)$ and $\lambda$ such that the two treatments of this update stage coincide. Starting from our approach, we would seek the solution of

or, putting it in terms of nulling the gradient of this energy, require

The ${\widehat{\textbf{x}}}$ that is the solution of this equation is our updated image. On the other hand, the $P^{3}$ scheme would propose to simply computeA delicate matter not considered here is that $P^{3}$ may apply $\frac{1}{c}f(c\textbf{y})$ in order to tune to a specific noise level. We assume $c=1$ for simplicity. ${\widehat{\textbf{x}}}=f(\textbf{y})$ as a replacement to this minimization task. Therefore, for the two methods to coincide, we should demand that the gradient expression posed above is solved for the choice of the $P^{3}$ scheme, namely,

This means that the denoising residual should remain the same (up to a constant) for the first activation of the denoising engine $\textbf{y}-f(\textbf{y})$, and the second one applied on the filtered image $f(\textbf{y})$.

In order to get a better intuition towards this result, let’s return to the pseudo-linear case, $f(\textbf{y})=\textbf{W}\textbf{y}$ with the assumption that W is a fixed and diagonalizable matrix. Plugged into the above condition, this gives

As the above equation should hold true for any image y, we require

Without loss of generality, we can assume that W is diagonal, after multiplying the above equation from the left and right by the diagonalizing matrix. With this simplification in mind, we now consider the eigenvalues of W, and the above equation implies that exact equivalence between our scheme and the $P^{3}$ one is obtained only if our denoising engine has eigenvalues that are purely $1$’s, or $\beta/\lambda$. Clearly, this is a very limiting case, which suggests that for all other cases, the two methods are likely to differ.

Interestingly, the above analysis is somewhat related to the one given in . Both and our treatment assume that the actual applied denoising engine is $f(\textbf{y})=\textbf{W}\textbf{y}$ within the ADMM scheme. While we ask for the conditions on W to fit our regularization term $\textbf{x}^{T}(\textbf{x}-\textbf{W}\textbf{x})$, the author of seeks the actual form of the prior to match this step, reaching the conclusion that the prior should be $\textbf{x}^{T}(\textbf{I}-\textbf{W})\textbf{W}^{\dagger}\textbf{x}$. Bearing in mind that the conditions we get for the equivalence between the two methods are too restricting and rarely met, the result in shows the actual gap between the two methods: While we regularize with the expression $\textbf{x}^{T}(\textbf{I}-\textbf{W})\textbf{x}$, an equivalence takes place only if the $P^{3}$ modifies this to involve $\textbf{W}^{\dagger}$, getting a far more complicated and less natural term.

Just before we conclude this section, we turn briefly to discuss the computational complexity of the proposed algorithm and its relation to the complexity of the $P^{3}$ scheme. Put very simply, RED and $P^{3}$ are roughly of the same computational cost. This is the case when RED is deployed via ADMM and assuming only one iteration in the update of v, as shown above. Similarly, when using the fixed-point option, RED has the same cost as $P^{3}$ per iteration.

To conclude, we must state that this paper is about a more general framework rather than a comparison to the $P^{3}$. Indeed, one could consider this work as an attempt to provide more solid mathematical foundations for methods like the $P^{3}$. In addition, when comparing $P^{3}$ and RED, one can identify several major differences that are far more central than the complexity issue, such as (1) a lack of a clear objective function that $P^{3}$ serves, while our scheme has a very well-defined penalty; (2) the inability to claim much in terms of convergence of the $P^{3}$, while our penalty is shown to be convex; (3) the complications of tuning the $P^{3}$ algorithm, which is very different from the experience we show with RED.

Results

In this section we compare the performance of the proposed framework to the $P^{3}$ approach, along with various other leading algorithms that are designed to tackle the image deblurring and super-resolution problems. To this end, we plug two substantially different denoising algorithms into the proposed scheme. The first is the (simple) median filter, which surprisingly turns out to act as a reasonable regularizer to our ill-posed inverse problems. This option is brought as a core demonstration of the idea that an arbitrary denoiser can be deployed in RED without difficulties. The second denoising engine we use is the state-of-the-art Trainable Nonlinear Reaction Diffusion (TNRD) method. This algorithm trains a nonlinear reaction-diffusion model in a supervised manner. As such, in order to treat different restoration problems, one should re-train the underlying model for every specific task – something we aim to avoid. In the experiments below we build upon the published pre-trained model by the authors of TNRD, tailored to denoise images that are contaminated by white Gaussian noise with a fixedIn order to handle an arbitrary noise-level, $\sigma_{f}$, we rely on the following relation $f_{\sigma_{f}}(\textbf{y})~{}=~{}\frac{1}{c}f_{5}(c\cdot\textbf{y})$, where $c=5/\sigma_{f}$. noise-level, which is equal to $5$. Leveraging this, we show how state-of-the-art deblurring and super-resolution results can be achieved simply by integrating the TNRD denoiser in RED. In all the experiments that follow, the parameters were manually set in order to enable each method to get its best possible results over the subset of images tested.

In order to have a fair comparison to previous work, we follow the synthetic non-blind deblurring experiments conducted in the state-of-the-art work that introduced the Non-locally Centralized Sparse Representation (NCSR) algorithm , which combines the self-similarity assumption with the sparsity-inspired model . More specifically, we degrade the test images, supplied by the authors of NCSR, by convolving them with two commonly used point spread functions (PSF); the first is a $9\times 9$ uniform blur, and the second is a 2D Gaussian function with a standard deviation of $1.6$. In both cases, an additive Gaussian noise with $\sigma=\sqrt{2}$ is then added to the blurred images. Similarly to NCSR, restoring an RGB image is done by converting it to the YCbCr color-space, applying the deblurring algorithm on the luminance channel only, and then converting the result back to the RGB domain.

Table 1 provides the restoration performance of the three RED schemes – the steepest-descent (SD), the ADMM, and the fixed-point (FP) methods – along with the results of theWe note that $P^{3}$ using TNRD has never appeared in an earlier publication, and it is brought here in order to let $P^{3}$ perform as best as it possibly can. $P^{3}$, the state-of-the-art NCSR and IDD-BM3D , and two additional baseline deblurring methods . For brevity, only the steepest-descent scheme is presented when considering the basic median filter as a denoiser. The performance is evaluated using the Peak Signal to Noise Ratio (PSNR) measure, higher is better, computed on the luminance channel of the ground-truth and the estimated image. The parameters of the proposed approach, as well as the ones of the $P^{3}$, are tuned to achieve the best performance on this dataset; in the case of the TNRD denoiser, these are depicted in Table 2 and 3, respectively. In the setting of the median filter, which extracts the median value of a $3\times 3$ window, we choose to run the suggested steepest-descent scheme for $N=400$ iterations with $\lambda=0.12$ for the uniform PSF, and $N=200$ with $\lambda=0.225$ for the Gaussian PSF.

Several remarks are to be made with regard to the obtained results. When the image is degraded by a Gaussian blur kernel, integrating the median filter in the proposed framework leads to a surprising restoration performance that is similar to the total variation deblurring . Furthermore, by choosing the state-of-the-art TNRD to be our denoising engine we achieve results that are competitive with the impressive NCSR and IDD-BM3D methods, which are specifically designed to tackle the deblurring task. Notice that the three versions of the proposed framework obtain a similar PSNR score. However, while the ADMM and the fixed-point variants are of similar complexity, the steepest-descent requires many more steps to converge and thereby more applications of the denoiser. As our last observation, based on the obtained results we conclude that the proposed approach is equivalent in quality to the alternative $P^{3}$ framework. However, tuning the parameters of the proposed algorithm is significantly simpler than the ones of the $P^{3}$; while in the $P^{3}$ the parameters should be modified throughout the iterations, in our approach these are always fixed (refer to Section 6.3 for a broader discussion).

An illustration of the convergence of the proposed approach using the three different numerical schemes (steepest-descent, ADMM, and fixed-point) and the two denoisers (median filter, and TNRD) is given in Figure 5. As can be observed, the three algorithms indeed converge, but at different rates; the steepest-descent is the slowest, while the ADMM is the fastest. The fixed-point strategy is slightly slower than the ADMM alternative, but requires only one application of the denoiser per iteration while the ADMM demands $m_{2}=3$ such operations. However, when comparing the fixed-point to the ADMM with the setting $m_{2}=1$ (now both have the same computational cost) we observe that the fixed-point is faster.

The above discussion is also supported visually in Figure 6 and 7, comparing the proposed method to the $P^{3}$ and the NCSR both for uniform and Gaussian PSF. As can be seen, by plugging the median filter into RED we improve the quality of the blurry image, yet the gap in performance between this simplistic approach and the state-of-the-art is substantial. Once relying on the TNRD, we obtain an efficient deblurring machine that has comparable results to $P^{3}$, and both are competitive or even slightly better than the NCSR.

2 Image Super-Resolution

Similarly to the previous subsection, we imitate the super-resolution experiments done in . To this end, a low-resolution image is generated by blurring the ground-truth high-resolution one with a $7\times 7$ Gaussian blur kernel with standard deviation $1.6$, followed by down-sampling by a factor of $3$ in each axis. Next, white Gaussian noise of standard deviation $5$ is added to the low-resolution images. The upscaling of an RGB image is done by transforming it first to the YCbCr color-space, super-resolving the luminance channel using the proposed approach (or the baseline methods), while the chroma channels are upscaled by bicubic interpolation. Lastly, the outcome is converted back to the RGB color-space.

In terms of PSNR, Table 4 presents the restoration performance of the three variants of the proposed approach in addition to the ones of the $P^{3}$, the NCSR and the ASDS-Reg algorithms. Similarly to the deblurring case, the PSNR is computed on the luminance channel only. In the case of the TNRD denoiser, the parameters that are used in our approach and in the $P^{3}$ are listed in Table 5 and 6, respectively. In the simpler case of the median filter (defined by a $3\times 3$ window), the number of iterations of the proposed steepest-descent algorithm is set to $N=50$ with a parameter $\lambda=0.0325$.

Interestingly, when setting the median filter to be our denoising engine we get a 2.19dB improvement (on average) over the bicubic interpolation. Alternatively, when choosing a stronger denoiser such as the TNRD, we achieve state-of-the-art results. Notably, the $P^{3}$ and the three variants of the proposed approach lead to similar restoration performance, consistent with the observation that was made in the context of the deblurring problem. These support once again the claim that our framework is a tangible alternative to the $P^{3}$. Figure 8 and 9 visually compare the proposed method to the $P^{3}$ and also to the state-of-the-art NCSR. As shown, the three algorithms offer an impressive restoration with sharp and clear edges, complying with the quantitative results which are given in Table 4.

3 Robustness to the Choice of Parameters

In this subsection we test the robustness of RED to the choice of its parameters, and contrast it to the $P^{3}$. To this end, we choose the single image super-resolution problem as a case-study, described in the previous subsection. Figure 10 (a) plots the average PSNR obtained by the different approaches as a function of the outer iterations. One can observe that RED (in all its forms) converges to a similar PSNR value. Also, an increase of $m_{2}$ (the number denoising steps within each iteration) leads to an improved rate of convergence of the ADMM. On the other hand, the curve describing the $P^{3}$ shows an unstable behavior and tendency to decrease in the PSNR after the first 200 iterations. Note that no tool has been suggested in the literature so far to automatically stop the $P^{3}$ for extracting the best performing outcome.

This unstable nature of the $P^{3}$ appears again as we modify the values of $\alpha$ and $\beta_{0}$. Figure 10 (b) shows the behavior of $P^{3}$ for several settings of these two parameters, clearly exhibiting an erratic behavior. One could observe that for specific choices of these two parameters, convergence is obtained, as manifested by the flattened curves. However, this is a fake convergence, caused by a large enough value of $\beta_{k}$. We stress that, in principle, a change in $\beta_{0}$ is expected to modify the convergence rate of the ADMM, but the steady state outcome should remain the same. However, when observing the curves in Figure 10 (b), it is clear that this is not the case in the $P^{3}$.

Back to RED, we repeat a similar experiment and test the sensitivity (or better yet the robustness) of the ADMM to the choice of $\beta$. Figure 10 (c) shows that different values of $\beta$ indeed affect the convergence rate (as expected), but the PSNR of the final outcome is always the same. This figure also indicates that the more accurate the solution of Part II in Figure 3 (obtained by increasing the value of $m_{2}$), the better the convergence rate.

The sensitivity of RED to the choice of $\sigma_{f}$ – the input noise-level to the denoiser – is depicted in Figure 10 (d). Notice that the choice of $\sigma_{f}$ affects directly the proposed regularizer, given by $\lambda\rho(\textbf{x},\sigma_{f})=\frac{\lambda}{2}\>\textbf{x}^{T}\left[\textbf{x}-f_{\sigma_{f}}(\textbf{x})\right]$. As such, a change in $\sigma_{f}$ is expected to modify the objective and thereby the resulting PSNR, as shown empirically in Figure 10 (d) for the FP methodOne expects that we could tune this parameter using existing techniques such as SURE, but we leave this topic for future research.. Clearly, a similar behavior is expected to occur when modifying the weight of the regularizer, $\lambda$, in which we choose to omit from this experimental part for brevity.

Conclusions

The idea put forward in this paper is strange – using a denoising engine within the regularization term in order to handle general inverse problems. A surprising outcome of this proposal is the fact that differentiation of the regularization terms remains tractable, still using the very same denoiser engine and not its derivative. This led us to the proposed scheme, termed Regularization by Denoising (RED). We have shown and discussed various appealing properties of this approach, such as convexity and its relation to advanced Laplacian smoothing. We have contrasted this scheme with the plug-and-play-prior method , and we have provided experiments that demonstrate the validity of the regularization strategy and the resulting competitive performance.

Is there a wider message in this work? Could it be that one could use more general $f$ functions and still form the proposed regularization? We have seen that all that it takes is the availability of the directional derivative of this function in order to follow all through. Strangely enough, we have shown how even a median filter could fit into this scheme, despite the fact that it does not have any relation to Gaussian denoising. More work is required to investigate the ability to further generalize RED to other and perhaps more daring regularization functionals.

A key question we have left open at this stage is the setting of the parameter $\sigma_{f}$. We chose this to be a fixed value, but we did not address the question of its influence on the overall performance, or whether a varying value strategy could lead to a benefit. More work is required here as well.

Appendices

Appendix A Can We Mimic any Prior?

So far we have shown that denoisers $f(\textbf{x})$ can be used to define the powerful family of priors $\textbf{x}^{T}(\textbf{x}-f(\textbf{x}))$ that can regularize a wide variety of inverse problem. Now, let’s consider the reverse question: For a given $\rho(\textbf{x})$, what is the effective $f(\textbf{x})$ behind it (if any)? More specifically, given $\rho(\textbf{x})$, we wish to find $f(\textbf{x})$ such that

Recall that one of the key conditions we assumed for the denoiser $f(\textbf{x})$ is that it is homogenous of degree $1$,

This immediately notifies us that the $\rho(\textbf{x})$ we wish to mimic in (50) must be $2$-homogenous because

Of course, not all priors satisfy this condition, but the class of priors that do is wide. For instance, while the total-variation function TV($\textbf{x})=\|\nabla x\|_{1}$ is only $1$-homogeneous, its square $\|\nabla x\|_{1}^{2}$ keeps us in business. More generally, since any norm is by definition $1$-homogenous, all regularizers of the type $\rho(\textbf{x})=\|\textbf{A}x\|_{q}^{2}$ for $q=1,2,\cdots,\infty$ are $2$-homogeneousThe sparsity-inspired $L_{0}$ is an exception that can not be treated this way.:

Squared norms are not the only functions at our disposal. Beyond norms, any $k$-homogenous function (such as homogeneous polynomials of degree $k$) can be made $2$-homogeneous by raising to power $2/k$. As well, order-statistics such as maximum, minimum and median are $1$-homonegous, and can be squared to the same end.

With the above discussion, we move forward with the assumption that we have a $2$-homogeous prior $\rho(\textbf{x})$ in our hands. Let’s recall that the task is to find $f(\textbf{x})$ so that

We proceed by differentiating both sides:

where in the last step we have invoked the directional derivative expression developed in (14). The solution, $f(\textbf{x})$, is a denoiser explicitly defined in terms of the prior,

This is quite a natural result in retrospect – it resembles a steepest descent stepIt is interesting to note the resemblance of this expression to a similar expression that holds for proximal operators (see ). In this context, the class of denoisers we have described using conditions 1 and 2 is a more general form of proximal mappings.. To see this, consider the denoising problem regularized by $\rho(\textbf{x})$:

One step of steepest descent with step-size of $1$ and initial condition $\textbf{x}_{0}=\textbf{y}$, would read

Appendix B Kernelizing Priors

We showed in the previous section that a prior with the right properties implies a denoiser beneath; and that this denoiser is directly given by the gradient of the prior. Next, let’s address a related, but more specific question: Given a prior $\rho(\textbf{x})$, does it imply a denoising filter of the pseudo-linear form $f(\textbf{x})=\textbf{W}(\textbf{x})\textbf{x}$? These types of filters are of course of great interest because they reveal the kind of weighted averaging done by the denoiser, and they are easy to implement. Given x, we can compute the weights $\textbf{W}(\textbf{x})$ in one step, and apply them as $\textbf{W}(\textbf{x})\textbf{x}$ to the image in another step. Some of the most popular filters to date, such as bilateral, NLM, and K-SVD, are of this convenient form.

Before we go about finding the hidden filter matrix $\textbf{W}(\textbf{x})$, let’s illustrate a useful property of the pseudo-linear filters. Take $f(\textbf{x})=\textbf{W}(\textbf{x})\textbf{x}$, and again invoke the expression $f(\textbf{x})=\nabla f(\textbf{x})\textbf{x}$ developed in (14). Substitution gives

for all x. From this we conclude that the gradient of pseudo-linear filters is in fact the weight matrix

Now we can go after the weight matrix by taking the analysis from the previous section one step further and computing the second derivative of (50). Starting with the expression (55) that arose from the first derivative, we differentiate again,

Replacing $\nabla f(\textbf{x})=\textbf{W}(\textbf{x})$, we obtain the pleasing result that the weight matrix implied by the prior is the identity matrix minus the Hessian of the prior,

or posed another way, $\textbf{L}(\textbf{x})=\textbf{I}-\textbf{W}(\textbf{x})=\textbf{H}(\rho(\textbf{x}))$ (i.e. the Laplacian filter is directly given by the Hessian of the prior, which is not surprising, bearing in mind that we seek a relation of the form $\rho(\textbf{x})=\textbf{x}^{T}\textbf{L}(\textbf{x})\textbf{x}$). What we have done here is to “kernelize” the regularizer and find an explicit expression for the weights of the implied filter. Several observations about this result are in order:

Convexity: If $\rho(\textbf{x})$ is convex, then its Hessian is symmetric positive semi-definite (PSD), and therefore $\textbf{L}(\textbf{x})$ is PSD. Furthermore, if $\eta\left(\textbf{L}(\textbf{x})\right)\leq 1$, i.e., we get that (i) $\textbf{W}(\textbf{x})$ has spectral radius equal or smaller than $1$, implying that strong passivity condition (Condition 2) is guaranteed; and (ii) $\textbf{W}(\textbf{x})=\textbf{I}-\textbf{L}(\textbf{x})$ is also PSD – a desirable property .

Homogeneity: Since $\rho(\textbf{x})$ is 2-homogenous, its gradient is $1$-homogeneous. We can show this by differentiating:

Similarly the Hessian $\textbf{H}(\rho(\textbf{x}))$ is invariantOr -homogenous to scaling of the image x, and so is the implied filter matrix $\textbf{W}(\textbf{x})$. Consequently, the applied filter $\textbf{W}(\textbf{x})\textbf{x}$ is $1$-homogenous, which again is consistent with our earlier conditions.

Row Stochasticity: Let’s recall the expression we derived above for the filter in terms of the Hessian of the prior,

This relation does not hold in general. However, consider defining the prior in terms of the gradient of the image instead of the image itself. Namely, this involves a change of variables in the prior $\rho(\textbf{x})$ from x to Dx, where D is the gradient (e.g. difference) operator. For instance, instead of $\rho(\textbf{x})=\|\textbf{x}\|_{1}^{2}$, consider $\rho_{\textbf{D}}(\textbf{x})=\|\textbf{D}\textbf{x}\|_{1}^{2}$. The Hessian of the prior under this linear transformation is given by the chain rule as

Appendix C More on Homogeneity

The concept of homogeneity of a filter played an important role in the development of the key results of the paper. So it is worth saying a bit more about it and to question whether this condition is satisfied for some popular and familiar filters.

A general construction of a denoising filter could be based on a symmetric positive semi-definite kernel $\textbf{K}_{i,j}(\textbf{x})=\textbf{K}(x_{i},x_{j})\geq 0$ from which the filter matrix $\textbf{W}(\textbf{x})$ is constructed by normalization. More specifically,

Whether such a denoiser is homogenous very much depends on the choice of the kernel K. For instance, if the kernel is homogeneous of any degree $p$, then the resulting filter matrix is invariant to scaling through cancellation,

Examples of homogeneous kernels include (homogeneous) polynomials, and several others which are not in common use in image processing. Most commonly used kernels are the exponentials (Gaussian to be exact), which are used in the bilateral or non-local means (NLM) cases. The Gaussian function is not homogeneous, but as we will show below, the resulting pseudo-linear filter is nearly so. We will show that for $c=1+\epsilon$ with very small $\epsilon$ we have $1$-homogeneity for the NLM-type filters, namely:

The $i,j$-th element of the filter weight matrix for the NLM filterThe bilateral filter, which includes a spatial distance weight can be treated similarly, since the spatial weights are invariant to scaling of the values of the image in any case. is

where $\textbf{R}_{i}\textbf{x}$ is a patch centred at pixel position $i$, extracted from the image; the normalization constant $d_{j}(\sigma)$ is given by summing across the rows of the kernel matrix.

Do these weights change much when the image x is replaced by a scaled version $c\textbf{x}$? First, note that if $\sigma$ is nearly zero, then all weights are essentially equal to $1/n$ and therefore they are automatically invariant to scaling of the image. Next, let’s consider the other extreme where the value of $\sigma$ is away from zero. Now, note that the scaling in x can be absorbed in the parameter $\sigma$ as follows:

Second, the effect of this (multiplicative) scaling can be approximated by an additive perturbation,

We now compute an approximation to $\textbf{W}_{i,j}(\textbf{x},\sigma+\delta)$ using a Taylor series:

The derivative of the weight values will be calculated in terms of the functions $e_{ij}(\sigma)$ and $d_{j}(\sigma)$ as follows:

Replacing $\delta=-\epsilon\sigma$, we obtain

We can simplify further, but this is not necessary since $\phi(\sigma)$ does not depend on $\epsilon$. For its part, $\phi(\sigma)$ behaves like $n/\sigma^{2}$. To see this note that the first term in the definition of $\phi$ is at worst $n/\sigma^{2}$ since $\|\textbf{R}_{i}\textbf{x}-\textbf{R}_{j}\textbf{x}\|^{2}$ is bounded by $n$, given that the pixel values are in the range $$.

Similarly, $d$ is on the order of $n\cdot\exp(\|\textbf{R}_{i}\textbf{x}-\textbf{R}_{j}\textbf{x}\|^{2}/2\sigma^{2})$, and its derivative $d^{\prime}$ is on the order $n\cdot\|\textbf{R}_{i}\textbf{x}-\textbf{R}_{j}\textbf{x}\|^{2}\cdot d/\sigma^{3}$. Consequently, the second term $\sigma\cdot d^{\prime}/d$ behaves like $\|\textbf{R}_{i}\textbf{x}-\textbf{R}_{j}\textbf{x}\|^{2}/\sigma^{2}$ is also on the order $n/\sigma^{2}$. Therefore, choosing $\epsilon=1/n$, for sufficiently large $\sigma$, the term $\epsilon\phi(\sigma)$ becomes negligible.

What we have shown is that the filter weights change very little as a result of the scaling $(1+\epsilon)\textbf{x}$ as long as $\epsilon$ is very small. Therefore the NLM (and bilateral) filters are (almost exactly) $1$-homogeneous as we had hoped.

C.2 Tikhonov Regularizer and Wiener Filtering

We now turn to show that Tikhonov regularization obeys the homogeneity condition. In this case, the denoised image is the solution of

where B, for example, can be a discrete approximation of a derivative operator. The closed-from expression of the above minimization is given by

When we feed this linear denoiser with the scaled image $(1+\epsilon)\textbf{y}$, it gives

which is trivially the same as $(1+\epsilon)f(\textbf{y})$.

A more challenging case is obtained when the denoiser adapts to the input such that when we apply $f((1+\epsilon)\textbf{y})$, the denoiser modifies $\lambda$ to be $\lambda(1+\epsilon)^{2}$. In what follows we return to Equation (74), but this time with the modified $\lambda$, and study the behavior of the filter for infinitesimal change in $\epsilon$, i.e.,

By relying on the first-order Taylor expansion of the above expression and taking the limit $\epsilon\rightarrow 0$ we get the following

where $\textbf{W}=\textbf{W}(0)$. Therefore, we obtain that the adaptive filter $\textbf{W}(\epsilon)$ changes linearly with $\epsilon$. Moreover, $\|\sigma^{2}\textbf{W}\textbf{B}^{T}\textbf{B}\|<1$, and thus, if $\lambda\ll 1$ the second term becomes negligible (similar to what we have shown for the NLM filter). Therefore, we conclude that the image-adaptive Wiener filtering satisfies the homogeneity condition under mild conditions.

C.3 Patch-Based Denoising

We now turn to discuss state-of-the-art patch-based denoising algorithms. These methods clean the input noisy image by (i) breaking it into small overlapping patches, (ii) applying a local denoising step, and (iii) reconstructing the output by averaging the denoised overlapping patches to form the final image. The wisdom of these algorithms relies on the choice of the local non-linear prior. Moreover, in most cases, the denoising process can be divided into two parts; the first contains all the non-linear decisions made by the prior, while the second is nothing but a linear filter that cleans the noisy patches followed by a patch-averaging step. Clearly, the latter fulfills the homogeneity condition due to its linearity. In what follows we argue that the non-linear part of various popular denoisers is stable to an infinitesimal change in scale, thus leading to an overall satisfaction of the homogeneity property.

We start by the EPLL which can be considered as an iterative GMM denoising algorithm. The non-linear part of the GMM prior is the choice of a (pre-trained) Gaussian model that best fits the input noisy patch, and its linear part is the subsequent Wiener filtering step. Consider a Gaussian model for an $n$-dimensional patch $\textbf{x}_{i}\sim N(\textbf{0},\mathbf{\Sigma}_{k})$, where $\mathbf{\Sigma}_{k}\in\mathbf{R}^{n\times n}$ is the $k$-th Gaussian, taken from the mixture. Following the derivations in , the MAP estimate is formulated by

which is nothing but the Wiener filter that cleans the $i$-th noisy patch $\textbf{y}_{i}$. The best model for the $i$-th patch, $k_{i}^{*}$, is the one that maximizing the MAP over all the possible models, given by

By plugging Equation (C.3) into the above we get

When we feed the denoiser with $(1+\epsilon)\textbf{y}_{i}$, the above can be written as

By relying on the relation $1/(1+\epsilon)^{2}\approx 1-2\epsilon$ we further simplify the above and obtain

Now we turn to compare Equation (C.3) to the one derived above, and get the following condition on $\epsilon$ that guarantees that $k_{i}^{*,\epsilon}=k_{i}^{*}$:

Since $\Psi(\textbf{y},\sigma,\mathbf{\Sigma}_{k})>\Psi(\textbf{y},\sigma,\mathbf{\Sigma}_{k^{*}})$ one can always choose $\epsilon\rightarrow 0$ which will keep this inequality intact. In an extremely rare case, when $\Psi(\textbf{y},\sigma,\mathbf{\Sigma}_{k})=\Psi(\textbf{y},\sigma,\mathbf{\Sigma}_{k^{*}})$, we can modify the GMM denoiser and propose a simple rule for choosing the model that have a smaller $\log|\mathbf{\Sigma}_{k}|$ term, ensuring that $k_{i}^{*,\epsilon}=k_{i}^{*}$. To conclude, the non-linear part of the GMM denoising algorithm is stable to small scaling of the input.

Sparsity-Inspired Denoisers: K-SVD

Given a dictionary, the non-linear mechanism of the K-SVD is the Orthogonal Matching Pursuit (OMP) algorithm, which estimates the sparse representation of an input noisy patch. This is a greedy method, aiming to approximate the solution of

where $n$ is the size of the patch, and $\textbf{D}\in\mathbf{R}^{n\times m}$ is a (possibly redundant, $m>n$, and non-orthogonal) dictionary. At each step, denoted by $k$, the OMP picks a new atom (a column from D) that minimizes the residual. Formally, the rule for choosing the first atom $d_{j_{1}}$ can be written as

while in the $t$-th step of the OMP, it is the one that maximizing

where $\textbf{r}_{i}^{t}=\textbf{y}_{i}-\textbf{D}_{S_{i}^{t}}\alpha_{i}^{t}$ is the residual. We denoted by $S_{i}^{t}$ the set of chosen atoms, obtained in the previous steps, and by $\textbf{D}_{S_{i}^{t}}\in\mathbf{R}^{n\times|{S_{i}^{t}}|}$ the corresponding dictionary – a matrix having the chosen atoms as its columns. This expression can be further simplified by relying on the fact that the representation is the outcome of a least-squares solution, given by

Substituting the above in Equation (85) results in

This process is repeated until reaching to the error constraint.

Trivially, scaling the input patch would not modify the result of Equations (84) and (87). However, scaling the input may modify the stopping rule, and thereby the number of atoms that the OMP picks. This will happen only when $\|\textbf{y}_{i}-\textbf{D}\alpha_{i}\|_{2}^{2}=n\sigma^{2}$, which is an extremely rare case. Notice that given the final set of chosen atoms, $S_{i}^{*}$, the cleaned patch is given by $\hat{\textbf{x}}_{i}=\textbf{D}_{S_{i}^{*}}\left(\textbf{D}_{S_{i}^{*}}^{T}\textbf{D}_{S_{i}^{*}}\right)^{-1}\textbf{D}_{S_{i}^{*}}^{T}\textbf{y}_{i}$, which clearly satisfies the homogeneity condition. To conclude, we showed that the non-linear part of the OMP is stable to an infinitesimal change in scale, and since the cleaned patch is simply obtained by a linear projection onto the chosen atoms we have that with high probability the OMP (and thus the K-SVD) fulfills our hope for homogeneity.

Appendix D The Differentiability Requirement

In the opening of Section 3.1, we required the denoising engine $f(\textbf{x})$ to be differentiable. Why? There are several benefits for having this behavior:

The directional derivative property, $\nabla_{\textbf{x}}f(\textbf{x})\textbf{x}=f(\textbf{x})$, which emerges from the homogeneity condition, becomes possible.

The passivity condition, which refers to the spectral radius of $\nabla_{\textbf{x}}f(\textbf{x})$, stands on solid grounds.

The convergence of the fixed-point algorithm discussed in Section 4.2 is guaranteed.

The convexity of the regularization term $\rho_{L}(\textbf{x})$, discussed in Section 5.1, and emerging from Condition 2, becomes possible.

All these are good reasons for demanding differentiability of the denoiser $f(\textbf{x})$, and yet, the question is whether it is too limiting or whether this requirement could be circumvented.

Observe that Point 1 above could rely on the availability of a far weaker requirement of the availability of all directional derivatives of the form $\nabla_{\textbf{x}}f(\textbf{x})\textbf{x}$. Indeed, the passivity mentioned in Point 2 could also be posed in terms of a directional derivative, replacing Equation (15) by the somewhat weaker requirement

under the assumption that $\nabla_{\textbf{x}}f(\textbf{x})$ is symmetric.

The question we leave as open at this stage is whether Points 3 and 4 above (convergence of the FP algorithm and the convexity of our regularization) could rely on the existence of all directional derivatives. At worst, if indeed the denoising engine is not differentiable but has all its directional derivatives, Points 3 and 4 are lost, and the behavior of the proposed algorithm is not clear from a theoretical standpoint.

A different, and more practical question is whether differentiability assumption on $f(\textbf{x})$ is feasible in existing algorithms. The NLM and the Bilateral filter are clearly differentiable. TNRD is differentiable as well since its non-linearities are smooth (linear combination of Gaussian RBF’s). EPLL, BM3D, are K-SVD more challenging since their non-linear parts include sharp-decisions. Each of these methods could be $\epsilon$-modified to have a fuzzy decision, thus rendering all of them differentiable with hardly any change in their actual behavior.

References