text:
\section{Introduction}

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

\noindent Before justifying this statement, let us describe the basic building block that will be the star of this paper: {\em 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 $\textbf{x}$ corrupted by zero-mean white additive Gaussian noise,
\begin{equation}\label{eqn:model}
\textbf{y} = \textbf{x} + \textbf{e}.
\end{equation}
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:[0,255]^n \longrightarrow [0,255]^n$ that maps an image $\textbf{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 $\textbf{x}$. Ideally, such functions operate on the input image $\textbf{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 \cite{NLM3,KSVDden,BM3D,KerReg,Mairal,GMM1,CSR,NCSR,ChattPeyman,NNden,Lebrun1,Lebrun2,DDID,Hossein,WNNM,NL-GSM,Morel}). Indeed, these documented results have led researchers to the educated guess that these methods are getting very close to the optimally possible denoising performance \cite{PeymanDenoisingDead,LevinNadler1,LevinNadler2}. 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 \cite{DenGen1,DenGen2,DenGen3,DenGen4,DenGen5,DenGen6}, 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 \cite{P3,Follows2,Follows1,Follows11}. 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 \cite{ADMM1,ADMM2,ADMM3}, 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 \cite{Follows1}. 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., \cite{Follows1,Follows11,Follows2,Follows3,Follows4,Follows5,Follows6,PPP-Convergence}), 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
\begin{equation}
\rho(\textbf{x}) = \frac{1}{2} \: \textbf{x}^T\left[\textbf{x}-f(\textbf{x})\right],
\end{equation}
in which the denoising engine itself is applied on the candidate image $\textbf{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 \ref{sec:Engine} focuses on the image denoising engine, defining it and its properties clearly, so as to enable its use in the proposed Laplacian paradigm. Section \ref{sec:Algorithm} 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 \ref{sec:Analysis} 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 \ref{sec:Results}, demonstrating the strength of the proposed scheme. We conclude the paper in Section \ref{sec:Conclusions} with a summary of the open questions that we identify for future work.

\section{Preliminaries}\label{sec:Background}

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.

\subsection{Inverse Problems as Optimization Tasks}
Bayesian estimation of an unknown image $\textbf{x}$ given its measured version $\textbf{y}$ uses the posterior conditional probability, $P(\textbf{x} \vert \textbf{y})$, in order to infer $\textbf{x}$. The most popular estimator in this regime is the Maximum aposteriori Probability (MAP), which chooses the mode ($\textbf{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
\begin{eqnarray}\label{eqn:BayesianMAP}
{\widehat \textbf{x}}_{MAP} = \argmax{\textbf{x}} P(\textbf{x} \vert \textbf{y}) & = & \argmax{\textbf{x}} \frac{P(\textbf{y} \vert \textbf{x}) P(\textbf{x})}{P(\textbf{y})}
= \argmax{\textbf{x}} P(\textbf{y} \vert \textbf{x}) P(\textbf{x}) \\ \nonumber
& = & \argmin{\textbf{x}} -\log\{P(\textbf{y} \vert \textbf{x})\} -\log\{ P(\textbf{x})\}.
\end{eqnarray}
In the above derivations we exploited the fact that $P(\textbf{y})$ is not a function of $\textbf{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} \vert \textbf{x})\}$ is known as the log-likelihood term, and it encapsulates the probabilistic relationship between the desired image $\textbf{x}$ and the measurements $\textbf{y}$, under the assumption that the desired image is known. We shall rewrite this term as
\begin{eqnarray}
\ell(\textbf{y}, \textbf{x}) = - \log\{P(\textbf{y} \vert \textbf{x})\}.
\end{eqnarray}
As a classic example for the log-likelihood that will accompany us throughout this paper, the expression $\ell(\textbf{y}, \textbf{x}) = \frac{1}{2\sigma^2}\|\textbf{H} \textbf{x} - \textbf{y}\|_2^2$ refers to the case of $\textbf{y} = \textbf{H} \textbf{x}+ \textbf{e} $, where $\textbf{H}$ is a linear degradation operator and $\textbf{e}$ is white Gaussian noise contamination of variance $\sigma^2$. Naturally, if the noise distribution changes, we depart form the comfortable $L_2$ form.

The second term in Equation (\ref{eqn:BayesianMAP}), $-\log\{ P(\textbf{x})\}$, refers to the {\em 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
\begin{eqnarray}
\lambda \rho(\textbf{x}) = - \log\{P( \textbf{x})\},
\end{eqnarray}
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 \cite{LaplacianReg}, through total variation \cite{TV} and wavelet sparsity \cite{Wavelet}, all the way to recent proposals based on patch-based GMM \cite{GMM1,GMM2} and sparse-representation modeling \cite{SparseModel}. Interestingly, the work we report here builds on the surprising comeback of the Laplacian regularization in a much more sophisticated form, as reported in \cite{LapReg1,LapReg2,LapReg3,LapReg4,Milanfar_SPM,LapReg5,LapReg6,LapReg7,LapReg8,LapReg9}.

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
\begin{eqnarray}\label{eqn:BayesianMAP0}
{\widehat \textbf{x}}_{MAP} = \argmin{\textbf{x}} \ell(\textbf{y} , \textbf{x}) + \lambda \rho(\textbf{x}).
\end{eqnarray}
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.

\subsection{The Plug-and-Play Prior ($P^3$) Approach}

For completeness of this exposition, we briefly review the $P^3$ approach. Aiming to solve the problem posed in Equation (\ref{eqn:BayesianMAP0}), the ADMM technique \cite{ADMM1,ADMM2,ADMM3} suggests to handle this by variable splitting, leading to the equivalent problem
\begin{eqnarray}\label{eqn:VarSplit1}
\left\{{\widehat \textbf{x}}_{MAP},{\widehat \textbf{v}} \right\} = \argmin{\textbf{x}, \textbf{v}} \ell(\textbf{y} , \textbf{x}) + \lambda \rho(\textbf{v})~~~~s.t.~~~\textbf{x}=\textbf{v}.
\end{eqnarray}
The constraint is turned into a penalty term, relying on the augmented Lagrangian method (in its scaled dual form \cite{ADMM1}), leading to
\begin{eqnarray}\label{eqn:VarSplit2}
\left\{{\widehat \textbf{x}}_{MAP},{\widehat \textbf{v}} \right\} = \argmin{\textbf{x}, \textbf{v}} \ell(\textbf{y} , \textbf{x}) + \lambda \rho(\textbf{v}) + \frac{\beta}{2}\|\textbf{x} - \textbf{v} +\textbf{u}\|_2^2,
\end{eqnarray}
where $\textbf{u}$ serves as the Lagrange multiplier vector for the set of constraints. ADMM addresses the resulting problem by updating $\textbf{x}$, $\textbf{v}$, and $\textbf{u}$ sequentially in a block-coordinate-descent fashion, leading to the following series of sub-problems:
\begin{enumerate}
\item \textbf{Update of $\textbf{x}$:} When considering $\textbf{v}$ (and $\textbf{u}$) as fixed, the term $\rho(\textbf{v})$ is omitted, and our task becomes
\begin{eqnarray}\label{eqn:UpdateX}
{\widehat \textbf{x}} = \argmin{\textbf{x}} \ell(\textbf{y} , \textbf{x}) + \frac{\beta}{2}\|\textbf{x} - \textbf{v} +\textbf{u}\|_2^2,
\end{eqnarray}
which is a far simpler inverse problem, where the regularization is an $L_2$ proximity one, which is easy to solve in most cases.

\item \textbf{Update of $\textbf{v}$:} In this stage we freeze $\textbf{x}$ (and $\textbf{u}$), and thus the log-likelihood term drops, leading to
\begin{eqnarray}\label{eqn:UpdateV}
{\widehat \textbf{v}} = \argmin{\textbf{v}}  \lambda \rho(\textbf{v}) + \frac{\beta}{2}\|\textbf{x} - \textbf{v} + \textbf{u}\|_2^2.
\end{eqnarray}
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 (\ref{eqn:BayesianMAP0}) 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 \cite{P3}, as they suggest to replace the direct solution of (\ref{eqn:UpdateV}) 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.

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

\end{enumerate}

\noindent
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 $ \textbf{x} $ and $ \textbf{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. \cite{Follows6,Follows11}), 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 \cite{Follows11} 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 \cite{Follows1} requires the denoiser to be a symmetric smoothing and non-expansive filter, the later work in \cite{Follows11} 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 \cite{Follows2} 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 \cite{PPP-Convergence}.

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 \cite{Follows1} 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.

\section{The {\em Image Denoising Engine}}\label{sec:Engine}

Image denoising is a special case of the inverse problem posed in Equation (\ref{eqn:BayesianMAP0}), referring to the case $\textbf{y} = \textbf{x}+ \textbf{e} $, where $\textbf{e}$ is white Gaussian noise contamination of variance $\sigma^2$. In this case, the MAP problem becomes
\begin{eqnarray}\label{eqn:BayesianMAP1}
{\widehat \textbf{x}}_{Denoise} = \argmin{\textbf{x}} \frac{1}{2\sigma^2} \|\textbf{y} - \textbf{x}\|_2^2 + \lambda \rho(\textbf{x}).
\end{eqnarray}
The {\em 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:[0,255]^n \longrightarrow [0,255]^n$ that maps an image $\textbf{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}$\footnote{In the context of optimization, a smaller class of the general denoising algorithms we define are characterized as ``proximal operators'' \cite{Proximal}. 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.

\subsection{Conditions and Properties of $f(\textbf{x})$}

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

\begin{itemize}

\item \textbf{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\geq0$ we must have $f(c\textbf{x}) = c f(\textbf{x})$. In this work we shall relax this condition and demand its satisfaction for $|c-1|\le \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 $\textbf{x}$. This derivative can be evaluated as
\begin{equation}
\label{eq:gradient}
\nabla_{\textbf{x}} f(\textbf{x})\: \textbf{x} = \frac{f(\textbf{x}+\epsilon \textbf{x})- f(\textbf{x})}{\epsilon}
\end{equation}
for a very small $\epsilon$. Invoking the homogeneity condition this leads to
\begin{eqnarray}\nonumber
\nabla_{\textbf{x}} f(\textbf{x})\: \textbf{x} & = &  \frac{(1+\epsilon) f(\textbf{x})- f(\textbf{x})}{\epsilon}  \\
& = & f(\textbf{x}).
\end{eqnarray}
Thus, the filter $f(\textbf{x})$ can be written as\footnote{This result is sometimes known as Euler's homogeneous function theorem \cite{Euler}.}
\begin{equation}
\label{eq:Euler3}
f(\textbf{x}) = \nabla_{\textbf{x}} f(\textbf{x})\: \textbf{x}.
\end{equation}

\item \textbf{Condition 2: Strong Passivity.} The Jacobian $\nabla_{\textbf{x}} f(\textbf{x})$ of the denoising algorithm is stable, satisfying the condition
\begin{eqnarray}
\eta \left(\nabla_{\textbf{x}} f(\textbf{x}) \right) \le 1,
\end{eqnarray}
where $\eta(\textbf{A})$ is the spectral radius of the matrix $\textbf{A}$. We interpret this condition as the restriction of the denoiser not to magnify the norm of an input image since
\begin{eqnarray}
\|f(\textbf{x})\| = \|\nabla_{\textbf{x}} f(\textbf{x}) \textbf{x}\| \le \eta\left(\nabla_{\textbf{x}} f(\textbf{x}) \right) \cdot \|\textbf{x}\|\le \|\textbf{x}\|.
\end{eqnarray}
Here we have relied on the relation $f(\textbf{x}) = \nabla_{\textbf{x}} f(\textbf{x}) \textbf{x}$ that has been established in Equation (\ref{eq:Euler3}). Note that we have chosen this specific condition over the natural weaker alternative, $\|f(\textbf{x})\| \le \|\textbf{x}\|$, since the strong condition implies the weak one.

\end{itemize}

\noindent A reformulation of the denoising engine that will be found useful throughout this work is the one suggested in \cite{Milanfar_SPM}, 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
\begin{eqnarray}\label{eqn:linearfilter}
\widehat{\textbf{x}}_{\text{Denoise}} = f(\textbf{y}) = \textbf{W} (\textbf{y}) \textbf{y}.
\end{eqnarray}
The matrix $\textbf{W}$ is an $n\times n$ matrix, representing the (pseudo-) linear filter that multiplies the $n\times 1$ noisy image vector $\textbf{y}$. This matrix is image dependent, as it draws its content from the pixels in $\textbf{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 {\em 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})$.

As a side note we mention that yet another natural requirement on $\textbf{W}$ (or $\nabla_{\textbf{y}} f(\textbf{y})$ in a wider perspective) is that it is row-stochastic, implying that (i) this matrix is (entry-wise) non-negative, and that (ii) the vector $\mathbb{1}$ is an eigenvector of $\textbf{W}$. This would imply a constancy behavior -- a denoiser $f(\textbf{y})$ does not alter a constant image. More specifically, defining $\mathbb{1}$ as an $n$-dimensional column vector of all ones, for any scalar $c\geq0$ we have $f(c\mathbb{1}) = c\mathbb{1}$. Algorithms such as the NLM \cite{NLM1} and its many variants all lead to such a row-stochastic Jacobian. We note that this property, while nice to have, is not required for the derivations in this paper.

An interesting consequence of the homogeneity property is the following stability of the pseudo-linear operator $\textbf{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
\begin{eqnarray}\label{eqn:stabilityW}
f(\textbf{y} + \textbf{h}) & \approx &  f(\textbf{y}) + \nabla_{\textbf{y}} f(\textbf{y}) \cdot \textbf{h} \\
& \approx & \nabla f(\textbf{y}) \cdot \textbf{y} + \nabla_{\textbf{y}} f(\textbf{y}) \cdot \textbf{h} \nonumber  \\
& \approx & \nabla_{\textbf{y}} f(\textbf{y}) \cdot (\textbf{y} + \textbf{h}). \nonumber
\end{eqnarray}
This result implies that while $\textbf{W}(\textbf{y}) = \nabla_{\textbf{y}} f(\textbf{y})$ may indeed depend on $\textbf{y}$, its sensitivity to its perturbation is negligible, rendering it as an essentially constant linear operator on the perturbed image $\textbf{y}+\textbf{h}$.

\subsection{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 \cite{NLM1,NLM2,NLM3}, the K-SVD denoising method that relies on sparse representation modeling of image patches \cite{KSVDden} and its non-local extension \cite{Mairal}, the kernel-regression method that exploits local orientation \cite{KerReg}, the well-known BM3D that combines sparsity and self-similarity of patches \cite{BM3D}, the EPLL scheme that suggests patch-modeling based on the GMM model \cite{GMM1}, CSR and NCSR, which cluster the patches and sparsifies them jointly \cite{CSR,NCSR}, the group-Wiener filtering applied on patches \cite{ChattPeyman}, the multi-layer Perceptron method trained to clean an image or the more recent CNN-based alternative called Trainable Nonlinear Reaction-Diffusion (TNRD) algorithm \cite{NNden,chen2015trainable}, more recent work that proposed low-rank modeling of patches and the use of the weighted nuclear-norm \cite{WNNM}, non-local sparsity with GSM model \cite{NL-GSM}, 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 \ref{fig:Homogeneity-Demo} 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 $\textbf{x}$ is set to be \textsf{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.

\begin{figure*}
	\begin{center}
		\subfigure[K-SVD ($2.95e-4$)]{\includegraphics[width=0.3\linewidth]{./scatter_homog_sigma_5_ksvd.pdf}}
		\subfigure[BM3D ($3.38e-4$)]{\includegraphics[width=0.3\linewidth]{./scatter_homog_sigma_5_bm3d.pdf}}
		\subfigure[NLM ($1.38e-4$)]{\includegraphics[width=0.3\linewidth]{./scatter_homog_sigma_5_nlm.pdf}}\\
		\subfigure[EPLL ($1.46e-4$)]{\includegraphics[width=0.3\linewidth]{./scatter_homog_sigma_5_epll.pdf}}
		\subfigure[TNRD ($9.51e-5$)]{\includegraphics[width=0.3\linewidth]{./scatter_homog_sigma_5_tnrd.pdf}}
	\end{center}
	\caption{An empirical evaluation of the homogeneity property. These graphs show $f((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 \textsf{Peppers}.} }
	\label{fig:Homogeneity-Demo}
\end{figure*}

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 \cite{Power-Method} that applies iterations of the form
\begin{eqnarray}
\textbf{h}_{k+1} = \frac{\nabla_{\textbf{x}} f(\textbf{x})\cdot \textbf{h}_{k}}{\| \nabla_{\textbf{x}} f(\textbf{x}) \cdot \textbf{h}_{k}\|_2}.
\end{eqnarray}
The spectral radius itself is easily obtained as
\begin{eqnarray}
\eta \left(\nabla_{\textbf{x}} f(\textbf{x})\right) = \frac{\textbf{h}_{k+1}^T \textbf{h}_{k}}{\textbf{h}_{k}^T\textbf{h}_{k}}.
\end{eqnarray}
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,
\begin{eqnarray}
f(\textbf{x}+\textbf{h}) = f(\textbf{x})+ \nabla_{\textbf{x}} f(\textbf{x}) \cdot \textbf{h},
\end{eqnarray}
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,
\begin{eqnarray}
\textbf{h}_{k+1} \approx \frac{ f(\textbf{x}+\textbf{h}_k) - f(\textbf{x})}{\| f(\textbf{x}+\textbf{h}_k) - f(\textbf{x})\|_2}.
\end{eqnarray}
{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 ($\textbf{x}$ set to be the image \textsf{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.

\section{Regularization by Denoising (RED)}\label{sec:Algorithm}

\subsection{The Image-Adaptive Laplacian}

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 (\ref{eqn:linearfilter}), namely\footnote{Note that we conveniently assume that the prior is applied to the clean image $\textbf{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 $\textbf{x}$) are first computed in the matrix $\textbf{W}$, and then applied to the image $\textbf{x}$. From this we can construct the Laplacian form,
\begin{equation}
\rho_L(\textbf{x}) = \frac{1}{2} \: \textbf{x}^T \textbf{L}(\textbf{x}) \textbf{x} = \frac{1}{2} \: \textbf{x}^T (\textbf{I} - \textbf{W}(\textbf{x})) \textbf{x} = \frac{1}{2} \: \textbf{x}^T \left[\textbf{x}- \textbf{W} (\textbf{x}) \textbf{x} \right].
\end{equation}
This definition by itself is not novel, as it is similar to ideas brought up in a series of recent contributions
\cite{LapReg1,LapReg2,LapReg3,LapReg4,Milanfar_SPM,LapReg5,LapReg6,LapReg7,LapReg8,LapReg9}. 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 {\em unnormalized cross-correlation} between the image $\textbf{x}$ and its corresponding residual $\textbf{x}-\textbf{W}(\textbf{x}) \textbf{x} $. As a {\em prior} expression should give low values for likely images, in our case this would be achieved in one of two ways (or their combination):
\begin{itemize}
\item A small value is obtained for $\rho_L(\textbf{x})$ if the residual is very small, implying that the image $\textbf{x}$ serves as a near fixed-point of the denoising engine, $\textbf{x} \approx \textbf{W}(\textbf{x}) \textbf{x}$.
\item 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 \cite{Dantzig} and by image denoising boosting techniques \cite{Milanfar_SPM,Boosting1,Boosting2}. 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)
\end{itemize}

\noindent Given the above prior, we return to the general inverse-problem posed in Equation (\ref{eqn:BayesianMAP0}), and define our new objective,
\begin{eqnarray}\label{eqn:BayesianMAPLap1}
{\widehat \textbf{x}} = \argmin{\textbf{x}} \ell(\textbf{y} , \textbf{x}) + \frac{\lambda}{2} \: \textbf{x}^T \left[\textbf{x}- \textbf{W} (\textbf{x}) \textbf{x}\right].
\end{eqnarray}
The prior expression, while exhibiting a possibly complicated dependency on the unknown $\textbf{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 $\textbf{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 $\textbf{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:
\begin{equation}
\rho_L(\textbf{x}) = \frac{1}{2} \: \textbf{x}^T\left[\textbf{x}-f(\textbf{x})\right].
\end{equation}
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 $\textbf{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 \ref{sec:Engine} on our denoising engines $f$. The overall energy functional to minimize is
\begin{eqnarray}\label{eqn:BayesianMAPLap2}
E(\textbf{x}) = \ell(\textbf{y} , \textbf{x}) + \frac{\lambda}{2} \: \textbf{x}^T\left(\textbf{x}- f(\textbf{x})\right),
\end{eqnarray}
and the gradient of this expression is readily available by
\begin{eqnarray}\label{eqn:Gradient1}
\nabla_{\textbf{x}} E(\textbf{x})  & = & \nabla_{\textbf{x}} \ell(\textbf{y} , \textbf{x}) + \frac{\lambda}{2} \: \nabla_{\textbf{x}} \left\{ \textbf{x}^T\left(\textbf{x}-f(\textbf{x})\right)\right\} \\ \nonumber
& = & \nabla_{\textbf{x}} \ell(\textbf{y} , \textbf{x}) + \frac{\lambda}{2} \: \nabla_{\textbf{x}} \textbf{x}^T\textbf{x} - \frac{\lambda}{2} \: \nabla_{\textbf{x}} \left[\textbf{x}^T f(\textbf{x}) \right]\\ \nonumber
& = & \nabla_{\textbf{x}} \ell(\textbf{y} , \textbf{x}) + \lambda \textbf{x} - \frac{\lambda}{2}\left[ f(\textbf{x}) + \nabla_{\textbf{x}} f(\textbf{x}) \textbf{x} \right].
\end{eqnarray}
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 replaced\footnote{A 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 (\ref{eq:Euler3}), implying that the gradient expression is further simplified to be
\begin{eqnarray}\label{eqn:Gradient3}
\nabla_{\textbf{x}} E(\textbf{x})  = \nabla_{\textbf{x}} \ell(\textbf{y} , \textbf{x}) + \lambda \textbf{x} - {\lambda} f(\textbf{x}) = \nabla_{\textbf{x}} \ell(\textbf{y} , \textbf{x}) + \lambda (\textbf{x} - f(\textbf{x})),
\end{eqnarray}
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
\begin{eqnarray}\label{eqn:Gradient4}
\nabla_{\textbf{x}} E(\textbf{x})  = \nabla_{\textbf{x}} \ell(\textbf{y} , \textbf{x}) + \lambda (\textbf{x} - \textbf{W}(\textbf{x}) \textbf{x}).
\end{eqnarray}
Observe that this is a non-trivial derivation of the gradient of the original penalty function posed in Equation (\ref{eqn:BayesianMAPLap1}).

\subsection{Deploying the Denoising Engine for Solving Inverse Problems}

In the discussion above we have seen that the gradient of the energy functional to minimize (given in Equation (\ref{eqn:BayesianMAPLap2})) is easily computable, given in Equation (\ref{eqn:Gradient3}). We now turn to show several options for using this in order to solve a general inverse problem. Common to all these methods is the fact that the eventual algorithm is iterative, in which each step is composed of applying the denoising engine (once or more), accompanied by other simpler calculations. In Figures \ref{fig:algoSD}, \ref{fig:algoADMM}, and \ref{fig:algoFP} we present pseudo-code for several such algorithms, all in the context of handling the case in which the likelihood function is given by $\ell(\textbf{y},\textbf{x}) = \|\textbf{H}\textbf{x} - \textbf{y}\|_2^2/2\sigma^2$.

\begin{itemize}

\item \textbf{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
\begin{eqnarray}\label{eqn:SD}
{\widehat \textbf{x}}_{k+1} & = & {\widehat \textbf{x}}_{k} - \mu \nabla_{\textbf{x}} \left. E(\textbf{x})\right|_{{\widehat \textbf{x}}_{k}} \\ \nonumber & = & {\widehat \textbf{x}}_{k} - \mu \left[\nabla_{\textbf{x}} \left. \ell(\textbf{y} , \textbf{x})\right|_{{\widehat \textbf{x}}_{k}} + \lambda ({\widehat \textbf{x}}_{k} -  f({\widehat \textbf{x}}_{k}))\right].
\end{eqnarray}
Figure \ref{fig:algoSD} 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 \cite{SESOP}. 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 $\textbf{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''.

\begin{figure}[hbtp]
\begin{center}
\begin{tabular}{|c|}
\hline
\begin{minipage}[b]{0.8\linewidth}
    \renewcommand{\baselinestretch}{0.5}
    \small \vspace{0.1in}
\scriptsize{
    \begin{description}
    \item {\textbf{f} Objective:} Minimize $E(\textbf{x}) = \frac{1}{2\sigma^2}\|\textbf{H}\textbf{x} -\textbf{y}\|_2^2 + \frac{\lambda}{2} \textbf{x}^T\left(\textbf{x}-f(\textbf{x})\right)$
    \item {\textbf{f} Input:} Supply the following ingredients and parameters:
    \begin{itemize}
            \item Regularization parameter $\lambda$,
            \item Denoising algorithm $f(\cdot)$ and its noise level $\sigma_f$,
            \item Log-Likelihood parameters: $\textbf{H}$ and $\sigma$, and
            \item Number of iterations $N$.
    \end{itemize}
    \item {\textbf{f} Initialization:}
    \item ~~~$\bullet$ Set ${\widehat \textbf{x}}_0 = \textbf{y}$.
    \item ~~~$\bullet$ Set $\mu = \frac{2}{1/\sigma^2 +\lambda}$.
    \item {\textbf{f} Loop:} For $k=1,2 ~\ldots~, N$ do:
    \item ~~~$\bullet$ Apply the denoising engine, ${\tilde \textbf{x}}_k=f_{\sigma_f}({\widehat \textbf{x}}_{k-1})$.
    \item ~~~$\bullet$ Update the solution by ${\widehat \textbf{x}}_k = {\widehat \textbf{x}}_{k-1} - \mu\left[\frac{1}{\sigma^2}\textbf{H}^T (\textbf{H} {\widehat \textbf{x}}_{k-1}- \textbf{y}) +\lambda ({\widehat \textbf{x}}_{k-1}-{\tilde \textbf{x}}_{k})\right]$.
    \item {\textbf{f} End of Loop}
    \item {\textbf{f} Result:} The output of the above algorithm is ${\widehat \textbf{x}}_N$.

    \end{description}
}
\end{minipage}
\\ \hline
\end{tabular}
\\ \vspace{0.1in}
\caption{The proposed scheme (RED) via the steepest descent method.} \label{fig:algoSD}
\end{center}
\end{figure}

\item \textbf{ADMM:} Addressing the optimization task given in Equation  (\ref{eqn:BayesianMAPLap2}), 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 \ref{sec:Background} almost exactly, with one major difference -- the prior in our case is explicit, and therefore, the stage termed ``update of $\textbf{v}$'' would become
\begin{eqnarray}\label{eqn:UpdateVours}
{\widehat \textbf{v}} = \argmin{\textbf{v}}  \frac{\lambda}{2} \textbf{v}^T\left(\textbf{v}-f(\textbf{v})\right) + \frac{\beta}{2}\|\textbf{v} - \textbf{x} - \textbf{u}\|_2^2.
\end{eqnarray}
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
\begin{eqnarray}\label{eqn:UpdateVours1}
\lambda (\textbf{v}-f(\textbf{v})) + \beta(\textbf{v} - \textbf{x} - \textbf{u})= \textbf{0},
\end{eqnarray}
which can be solved iteratively using the fixed-point strategy, by
\begin{eqnarray}\label{eqn:UpdateVours2}
\lambda (\textbf{v}_j-f(\textbf{v}_{j-1})) + \beta(\textbf{v}_{j} - \textbf{x} - \textbf{u})= \textbf{0}
\\ \nonumber \rightarrow~~
\textbf{v}_j = \frac{1}{\beta+\lambda}\left(  \lambda f(\textbf{v}_{j-1} + \beta (\textbf{x} + \textbf{u}) \right).
\end{eqnarray}
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'' \cite{ADMM1}, suggesting to approximate the solution of the $ \textbf{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 \eqref{eqn:UpdateVours2}, which amounts to one operation of a denoising algorithm.} Figure \ref{fig:algoADMM} describes this specific algorithm in more details. If one changes all {\textbf{f} Part 2} (in Figure \ref{fig:algoADMM}) 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 $\textbf{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.

\begin{figure}[hbtp]
\begin{center}
\renewcommand{\baselinestretch}{1}
\begin{tabular}{|c|}
\hline
\begin{minipage}[b]{0.8\linewidth}
    \small \vspace{0.1in}
\scriptsize{
    \begin{description}
    \item {\textbf{f} Objective:} Minimize $E(\textbf{x}) = \frac{1}{2\sigma^2}\|\textbf{H}\textbf{x} -\textbf{y}\|_2^2 + \frac{\lambda}{2} \textbf{x}^T\left(\textbf{x}-f(\textbf{x})\right)$
    \item {\textbf{f} Input:} Supply the following ingredients and parameters:
    \begin{itemize}
            \item Regularization parameter $\lambda$,
            \item Denoising algorithm $f(\cdot)$ and its noise level $\sigma_f$,
            \item Log-Likelihood parameters: $\textbf{H}$ and $\sigma$,
            \item Number of outer and inner iterations, $N$, $m_1$, and $m_2$, and
            \item ADMM coefficient $\beta$.
    \end{itemize}
    \item {\textbf{f} Initialization:} Set ${\widehat \textbf{x}}_0 = \textbf{y}$, ${\widehat \textbf{v}}_0 = \textbf{y}$, and ${\widehat \textbf{u}}_0 = \textbf{0}$.
    \item {\textbf{f} Outer Loop:} For $k=1,2 ~\ldots~, N$ do:

    \item ~~~~~~ \textbf{Part 1:} Solve ${\widehat \textbf{x}}_k = \argmin{\textbf{z}} \frac{1}{2\sigma^2}\|\textbf{H}\textbf{z} - \textbf{y}\|_2^2+ \frac{\beta}{2}\|\textbf{z}-{\widehat \textbf{v}}_{k-1}+{\widehat \textbf{u}}_{k-1} \|_2^2$ by

        \item ~~~~~~~~~$\bullet$ Initialization: ${\textbf{z}}_0 = {\widehat \textbf{x}}_{k-1}$, and define $\textbf{z}^*={\widehat \textbf{v}}_{k-1}-{\widehat \textbf{u}}_{k-1}$.
        \item ~~~~~~~~~$\bullet$ Inner Iteration: For $j=1,2 ~\ldots~, m_1$ do:
        \begin{description}
            \item ~~~~$-$ Compute the gradient $\textbf{e}_j = \frac{1}{\sigma^2}\textbf{H}^T(\textbf{H}\textbf{z}_{j-1}-\textbf{y})+\beta(\textbf{z}_{j-1}-\textbf{z}^*)$.
            \item ~~~~$-$ Compute $\textbf{r}_j= \frac{1}{\sigma^2}\textbf{H}^T\textbf{H}\textbf{e}_j+\beta\textbf{e}_j$.
            \item ~~~~$-$ Compute the step size $\mu= \textbf{e}_j^T \textbf{e}_j/\textbf{e}_j^T \textbf{r}_j$.
            \item ~~~~$-$ Update the solution by $\textbf{z}_j=\textbf{z}_{j-1}+\mu \textbf{e}_j$.
            \item ~~~~$-$ Project the result to the interval $[0,255]$.
        \end{description}
        \item ~~~~~~~~~$\bullet$ End of Inner Loop
        \item ~~~~~~~~~$\bullet$ Set ${\widehat \textbf{x}}_k = \textbf{z}_{m_1}$.

    \item ~~~~~~ \textbf{Part 2:} Solve  ${\widehat \textbf{v}}_k = \argmin{\textbf{z}} \lambda \textbf{z}^T (\textbf{z}-f_{\sigma_f}(\textbf{z}))+ \frac{\beta}{2}\|\textbf{z}-{\widehat \textbf{x}}_{k}-{\widehat \textbf{u}}_{k-1} \|_2^2$ by

        \item ~~~~~~~~~$\bullet$ Initialization: ${\textbf{z}}_0 = {\widehat \textbf{v}}_{k-1}$, and define $\textbf{z}^*={\widehat \textbf{x}}_{k}+{\widehat \textbf{u}}_{k-1}$.
        \item ~~~~~~~~~$\bullet$ Inner Iteration: For $j=1,2 ~\ldots~, m_2$ do:
            \begin{description}
                \item ~~~~$-$ Apply the denoising engine, ${\tilde \textbf{z}}_j=f_{\sigma_f}({\widehat \textbf{z}}_{j-1})$.
                \item ~~~~$-$ Compute the gradient $\textbf{z}_j = \frac{1}{\beta+\lambda} \left(\lambda {\tilde \textbf{z}}_j +\beta \textbf{z}^*\right)$.
            \end{description}
        \item ~~~~~~~~~$\bullet$ End of Inner Loop
        \item ~~~~~~~~~$\bullet$ Set ${\widehat \textbf{v}}_k = \textbf{z}_{m_2}$.

    \item ~~~~~~ \textbf{Part 3:} Update  ${\widehat \textbf{u}}_k = {\widehat \textbf{u}}_{k-1}+ {\widehat \textbf{x}}_{k}-{\widehat \textbf{u}}_{k}$.

    \item {\textbf{f} End of Outer Loop}
    \item {\textbf{f} Result:} The output of the above algorithm is ${\widehat \textbf{x}}_N$.

    \end{description}
}
\end{minipage}
\\ \hline
\end{tabular}
\\ \vspace{0.1in}
\caption{The proposed scheme (RED) via the ADMM method.} \label{fig:algoADMM}
\end{center}
\end{figure}

\item \textbf{Fixed-Point Strategy:} An appealing alternative to the above exists, obtained via the fixed-point strategy. As our aim is to find $\textbf{x}$ that nulls the gradient, this could be posed as an implicit equation to be solved directly,
\begin{eqnarray}
\nabla_{\textbf{x}} \ell(\textbf{y} , \textbf{x}) + \lambda (\textbf{x} -  f(\textbf{x}))=0.
\end{eqnarray}
Using the fixed-point strategy, this could be handled by the iterative formula
\begin{eqnarray}
\nabla_{\textbf{x}} \ell(\textbf{y} , \textbf{x}_{k+1}) + \lambda (\textbf{x}_{k+1} -  f(\textbf{x}_k))=0.
\end{eqnarray}
As an example, in the case of linear degradation model and Gaussian white additive noise, this equation would be
\begin{eqnarray}
\frac{1}{\sigma^2} \textbf{H}^T (\textbf{y} - \textbf{H}\textbf{x}_{k+1}) + \lambda (\textbf{x}_{k+1} -  f(\textbf{x}_k))=0,
\end{eqnarray}
leading to the recursive update relation
\begin{eqnarray}
\textbf{x}_{k+1}= \left[\frac{1}{\sigma^2} \textbf{H}^T \textbf{H} + \lambda \textbf{I}\right]^{-1} \left[ \frac{1}{\sigma^2}\textbf{H}^T \textbf{y} + \lambda f(\textbf{x}_k)\right].
\end{eqnarray}
This formula suggests one activation of the denoising per iteration, followed by what seems to be a plain Wiener filtering computation\footnote{Note that Equation (\ref{eqn:UpdateVours2}) 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 $\textbf{H}$, or iteratively using for example, the Richardson algorithm: Defining
\begin{eqnarray}
\textbf{A} = \frac{1}{\sigma^2} \textbf{H}^T \textbf{H} + \lambda \textbf{I}~~~\mbox{and}~~~
\textbf{b}= \frac{1}{\sigma^2}\textbf{H}^T \textbf{y} + \lambda f(\textbf{x}_k),
\end{eqnarray}
our goal is to solve the linear system $\textbf{A}\textbf{x} = \textbf{b}$. This is achieved by a variant of the SD method\footnote{All 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
\begin{eqnarray}
\left\| \left[\frac{1}{\sigma^2} \textbf{H}^T \textbf{H} + \lambda \textbf{I}\right]^{-1} \lambda \nabla_{\textbf{x}} f(\textbf{x}_k)\right\|<1.
\end{eqnarray}

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 \ref{fig:algoFP}.

\end{itemize}

\noindent 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.

\begin{figure}[hbtp]
\begin{center}
\begin{tabular}{|c|}
\hline
\begin{minipage}[b]{0.8\linewidth}
    \renewcommand{\baselinestretch}{1}
    \small \vspace{0.1in}
\scriptsize{
    \begin{description}
    \item {\textbf{f} Objective:} Minimize $E(\textbf{x}) = \frac{1}{2\sigma^2}\|\textbf{H}\textbf{x} -\textbf{y}\|_2^2 + \frac{\lambda}{2} \textbf{x}^T\left(\textbf{x}-f(\textbf{x})\right)$
    \item {\textbf{f} Input:} Supply the following ingredients and parameters:
    \begin{itemize}
            \item Regularization parameter $\lambda$,
            \item Denoising algorithm $f(\cdot)$ and its noise level $\sigma_f$,
            \item Log-Likelihood parameters: $\textbf{H}$ and $\sigma$, and
            \item Number of outer and inner iterations, $N$ and $m$.
    \end{itemize}
    \item {\textbf{f} Initialization:} ${\widehat \textbf{x}}_0 = \textbf{y}$.
    \item {\textbf{f} Outer Loop:} For $k=1,2 ~\ldots~, N$ do:
    \item ~~$\bullet$ Apply the denoising engine, ${\tilde \textbf{x}}_k=f_{\sigma_f}({\widehat \textbf{x}}_{k-1})$.
    \item ~~$\bullet$ Solve $\textbf{A} \textbf{z} = \textbf{b}$ for $\textbf{A} = \frac{1}{\sigma^2} \textbf{H}^T \textbf{H} + \lambda \textbf{I}~~\mbox{and}~~\textbf{b}= \frac{1}{\sigma^2}\textbf{H}^T \textbf{y} + \lambda {\tilde \textbf{x}}_k$
    \item ~~~~~~- Initialization: ${\textbf{z}}_0 = {\tilde \textbf{x}}_k$.
    \item ~~~~~~- Iterate: For $j=1,2 ~\ldots~, m$ do:
               \begin{itemize}
                  \item Compute the residual $\textbf{r}_j = \textbf{A} \textbf{z}_{j-1} - \textbf{b}$.
                  \item Compute the vector $\textbf{e}_j = \textbf{A} \textbf{r}_j$.
                  \item Compute the optimal step $\mu = \textbf{r}_j^T \textbf{e}_j/\textbf{e}_j^T \textbf{e}_j$.
                  \item Update the solution by $\textbf{z}_j=\textbf{z}_{j-1}+\mu\cdot \textbf{r}_j$.
                  \item Project the result to the interval $[0,255]$.
               \end{itemize}
    \item ~~~~~~- End of Inner Loop
    \item ~~~~~~-  Set ${\widehat \textbf{x}}_k = \textbf{z}_m$.
    \item {\textbf{f} End of Outer Loop}
    \item {\textbf{f} Result:} The output of the above algorithm is ${\widehat \textbf{x}}_N$.

    \end{description}
}
\end{minipage}
\\ \hline
\end{tabular}
\\ \vspace{0.1in}
\caption{The Proposed Scheme (RED) via the fixed-point method.} \label{fig:algoFP}
\end{center}
\end{figure}

\section{Analysis}\label{sec:Analysis}

\subsection{Convexity}

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 \ref{sec:Engine}, 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 \cite{OptBook}. 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 $\textbf{W}$ has all its eigenvalues in the range\footnote{We could in fact allow negative eigenvalues for $\textbf{W}$, but this is unnatural in the context of denoising.} $[0,1]$. This is indeed the case for the NLM filter \cite{NLM1}, the K-SVD-denoising algorithm \cite{LapReg9}, 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 \ref{sec:Engine}, $\eta(\nabla_{\textbf{x}} f(\textbf{x}) ) \le 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 (\ref{eqn:BayesianMAP0}). 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 \cite{P3,Follows11}. Furthermore, this result may seem similar to the one posed in  \cite{Follows2,PPP-Convergence}, 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.

\subsection{An Alternative Prior}

In Section \ref{sec:Algorithm} we motivated the choice of the proposed prior by the desire to characterize the unknown image $\textbf{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
\begin{eqnarray}
\rho_Q(\textbf{x}) = \| \textbf{x} - f(\textbf{x}) \|_2^2.
\end{eqnarray}
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
\begin{eqnarray}
\rho_Q(\textbf{x}) = \| \textbf{x} - f(\textbf{x}) \|_2^2 = \rho_L(\textbf{x}) + f(\textbf{x})^T ( f(\textbf{x})- \textbf{x}),
\end{eqnarray}
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
\begin{eqnarray}
\rho_Q(\textbf{x}) = \| \textbf{x} - \textbf{W}(\textbf{x}) \textbf{x} \|_2^2 = \textbf{x}^T (\textbf{I} - \textbf{W}(\textbf{x}))^T(\textbf{I} - \textbf{W}(\textbf{x})) \textbf{x}.
\end{eqnarray}
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'' \cite{FourthOrder}, 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.

\subsection{When is Plug-and-Play-Prior = RED ?}

In Section \ref{sec:Algorithm} 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 (\ref{eqn:BayesianMAPLap2}) with ADMM, we have shown that the only difference between this and the $P^3$ scheme resides in the update stage for $\textbf{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 $\textbf{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
\begin{eqnarray}
{\hat \textbf{x}} = \argmin{\textbf{x}}  \frac{\beta}{2}\|\textbf{x} - \textbf{y}\|_2^2 + \frac{\lambda}{2} \textbf{x}^T\left(\textbf{x}-f(\textbf{x})\right),
\end{eqnarray}
or, putting it in terms of nulling the gradient of this energy, require
\begin{eqnarray}
 \beta (\textbf{x} - \textbf{y}) + \lambda (\textbf{x}-f(\textbf{x})) =\textbf{0}.
\end{eqnarray}
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 compute\footnote{A 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,
\begin{eqnarray}
 \beta (f(\textbf{y}) - \textbf{y})  + \lambda (f(\textbf{y})-f(f(\textbf{y}))) =\textbf{0},
\end{eqnarray}
or posed slightly different,
\begin{eqnarray}
f(\textbf{y})-f(f(\textbf{y})) = \frac{\beta}{\lambda} ( \textbf{y} - f(\textbf{y})).
\end{eqnarray}
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 $\textbf{W}$ is a fixed and diagonalizable matrix. Plugged into the above condition, this gives
\begin{eqnarray}
\textbf{W} \textbf{y} -\textbf{W}^2\textbf{y} = \frac{\beta}{\lambda} ( \textbf{y} - \textbf{W}\textbf{y} ),
\end{eqnarray}
or posed differently,
\begin{eqnarray}
\left(\frac{\beta}{\lambda}\textbf{I} - \textbf{W}\right) \left(\textbf{I} -\textbf{W}\right)\textbf{y} =\textbf{0}.
\end{eqnarray}
As the above equation should hold true for any image $\textbf{y}$, we require
\begin{eqnarray}
\left(\frac{\beta}{\lambda}\textbf{I} - \textbf{W}\right) \left(\textbf{I} -\textbf{W}\right)=\textbf{0}.
\end{eqnarray}
Without loss of generality, we can assume that $\textbf{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 $\textbf{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 \cite{Follows1}. Both \cite{Follows1} 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 $\textbf{W}$ to fit our regularization term $\textbf{x}^T (\textbf{x} - \textbf{W} \textbf{x})$, the author of \cite{Follows1} 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 \cite{Follows1} 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 $\textbf{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.

\section{Results}\label{sec: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) \emph{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) \cite{chen2015trainable} method. This algorithm trains a nonlinear reaction-diffusion model in a \emph{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 fixed\footnote{In 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.

\subsection{Image Deblurring}

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 \cite{NCSR}, which combines the self-similarity assumption \cite{NLM1} with the sparsity-inspired model \cite{Elad_Book}. 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 \ref{tab:deb_results} 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 the\footnote{We 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 \cite{DenGen4}, and two additional  baseline deblurring methods \cite{beck2009fast,dong2011image}. 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 \ref{tab:params_laplacian_deblurring} and \ref{tab:params_ppp_deblurring}, respectively. In the setting of the median filter, which extracts the median value of a $ 3\times3 $ 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.

\begin{table}[htbp]
	\centering
	\tiny
	\setlength{\tabcolsep}{4.5pt}
	\renewcommand{\arraystretch}{0.7}
	\renewcommand{\tabcolsep}{6.1pt}
	\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|c|}
		\hline
		\textbf{Image} & \textsf{Butterfly} & \textsf{Boats} & \textsf{C. Man} & \textsf{House} & \textsf{Parrot} & \textsf{Lena}  & \textsf{Barbara} & \textsf{Starfish} & \textsf{Peppers} & \textsf{Leaves} & \textbf{Average} \\
		\hline
		\multicolumn{12}{|c|}{\textbf{Deblurring: Uniform kernel, $ \sigma=\sqrt{2} $}} \\
		\hline
		Total Variation \cite{beck2009fast} & 28.37 & 29.04 & 26.82 & 31.99 & 29.11 & 28.33 & 25.75 & 27.75 & 28.43 & 26.49 & 28.21 \\
		IDD-BM3D \cite{DenGen4} & 29.21 & \textbf{31.20} & 28.56 & \textbf{34.44} & 31.06 & 29.70 & \textbf{27.98} & 29.48 & 29.62 & 29.38 & 30.06 \\
		ASDS-Reg \cite{dong2011image} & 28.70  & 30.80  & 28.08 & 34.03 & 31.22 & 29.92 & 27.86 & 29.72 & 29.48 & 28.59 & 29.84 \\
		NCSR  & 29.68 & 31.08 & 28.62 & 34.31 & \textbf{31.95} & 29.96 & 28.10 & 30.28 & 29.66 & 29.98 & 30.36 \\
		$ P^3 $-TNRD & 30.32 & \textbf{31.19} & 28.73 & 33.90 & 31.86 & \textbf{30.13} & 27.21 & 30.27 & \textbf{30.11} & 30.08 & 30.38 \\
		RED: SD-Median Filter & 26.10 & 28.03 & 25.57 & 29.81 & 28.67 & 27.29 & 25.62 & 27.84 & 27.40 & 25.45 & 27.18 \\
		RED: SD-TNRD & 30.20 & \textbf{31.20} & 28.67 & 33.83 & 31.62 & 29.98 & 27.35 & 30.47 & 30.10 & 29.72 & 30.31 \\
		RED: ADMM-TNRD & \textbf{30.40}	& 31.12	& 28.71	& 33.77	& 31.86	& 30.03	& 27.27	& \textbf{30.58}	& \textbf{30.11}	& \textbf{30.12}	& \textbf{30.40} \\
		RED: FP-TNRD & \textbf{30.41} & 31.12 & \textbf{28.76} & 33.76 & 31.83 & 30.02 & 27.27 & \textbf{30.57} & \textbf{30.12} & \textbf{30.13} & \textbf{30.40} \\
		\hline
		\multicolumn{12}{|c|}{\textbf{Deblurring: Gaussian kernel, $ \sigma=\sqrt{2} $}} \\
		\hline
		Total Variation \cite{beck2009fast} & 30.36 & 29.36 & 26.81 & 31.50 & 31.23 & 29.47 & 25.03 & 29.65 & 29.42 & 29.36 & 29.22 \\
		IDD-BM3D \cite{DenGen4} & 30.73 & \textbf{31.68} & 28.17 & \textbf{34.08} & 32.89 & 31.45 & 27.19 & 31.66 & 29.99 & 31.4  & 30.92 \\
		ASDS-Reg \cite{dong2011image} & 29.83 & 30.27 & 27.29 & 31.87 & 32.93 & 30.36 & 27.05 & 31.91 & 28.95 & 30.62 & 30.11 \\
		NCSR \cite{NCSR} & 30.84 & 31.49 & 28.34 & 33.63 & 33.39 & 31.26 & \textbf{27.91} & 32.27 & 30.16 & 31.57 & 31.09 \\
		$ P^3 $-TNRD & \textbf{31.73} & \textbf{31.67} & 28.08 & {33.95} & \textbf{33.43} & \textbf{31.52} & 27.11 & \textbf{32.71} & \textbf{30.94} & \textbf{32.18} & \textbf{31.33} \\
		RED: SD-Median Filter & 29.02 & 30.01 & 26.45 & 31.59 & 31.32 & 30.00 & 25.02 & 30.29 & 28.53 & 28.69 & 29.09 \\
		RED: SD-TNRD & 31.57 & 31.53 & 28.31 & 33.71 & 33.19 & 31.47 & 26.62 & 32.46 & 29.98 & 31.95 & 31.08 \\
		RED: ADMM-TNRD &  31.66	& 31.55	& 28.31	& 33.73	& 33.33	& 31.40	& 26.76	& 32.49	& 30.48	& 31.93	& 31.16 \\
		RED: FP-TNRD & 31.66 & 31.55 & \textbf{28.38} & 33.74 & 33.33 & 31.39 & 26.76 & 32.49 & 30.51 & 31.93 & 31.17 \\
		\hline
	\end{tabular}
	\caption{Deblurring results measured in PSNR [dB] and evaluated on the set of images provided by the authors of NCSR \cite{NCSR}. The $ P^3 $ and RED build upon the TNRD \cite{chen2015trainable} as the denoising engine. We also provide the results obtained by integrating the median filter with the steepest-descent RED scheme. {PSNR scores being less than $ 0.01 $[dB] away from the highest result are highlighted. Note that the $ P^3 $ does not converge when setting the TNRD to be the denoising algorithm. Therefore, we run the $ P^3 $ for a fixed number of iterations, chosen to achieve the best PSNR on average (otherwise the restoration quality would be significantly inferior). Please refer to Figure \ref{fig:red_vs_ppp} and Section \ref{sec:Parameters} for more details regarding the sensitivity of the $ P^3 $ to the choice of parameters. }}
	\label{tab:deb_results}
\end{table}

\begin{table}[htbp]
	\centering
	\scriptsize

	\renewcommand{\arraystretch}{0.8}
	\renewcommand{\tabcolsep}{9pt}
	\begin{tabular}{|c|c|c|c|c|}
		\hline
		\multirow{3}{*}{\textbf{PSF}} & \multicolumn{4}{c|}{\textbf{Proposed approach: Deblurring}} \\
		\cline{2-5}
		& \multirow{2}{*}{\textbf{Parameter}} & \textbf{Steepest} & \textbf{Fixed} & \multirow{2}{*}{\textbf{ADMM}} \\
		&                            & \textbf{Descent}   & \textbf{Point} &                       \\
		\hline \hline
		\multirow{6}{*}{\rotatebox[origin=c]{90}{Uniform}} & N     & 1500  & 200   & 200 \\
		& $ \sigma_f $ & $3.25 $ & $3.25 $ & $3.25 $ \\
		& $ \lambda $ & 0.02  & 0.02  & 0.02 \\
		& $ m_1 $    & --     & closed-form using FFT   & closed-form using FFT \\
		& $ m_2 $    & --     & --     & 1 \\
		& $ \beta $  & --     & --     & 0.001 \\
		\hline \hline
		\multirow{6}{*}{\rotatebox[origin=c]{90}{Gaussian}} & N     & 1500  & 200   & 200 \\
		& $ \sigma_f $ & $4.1 $ & $4.1 $ & $4.1 $ \\
		& $ \lambda $ & 0.01  & 0.01  & 0.01 \\
		& $ m_1 $    & --     & closed-form using FFT   & closed-form using FFT \\
		& $ m_2 $    & --     & --     & 1 \\
		& $ \beta $  & --     & --     & 0.001 \\
		\hline
	\end{tabular}%
	\caption{The set of parameters being used in our framework, leading to the deblurring results reported in Table~\ref{tab:deb_results} when plugging the TNRD \cite{chen2015trainable} denoiser.}
	\label{tab:params_laplacian_deblurring}%
\end{table}%

\begin{table}[htbp]
	\centering
	\scriptsize

	\renewcommand{\arraystretch}{0.8}
	\renewcommand{\tabcolsep}{9pt}
	\begin{tabular}{|c|c|c|}
		\hline
		\multirow{2}{*}{\textbf{PSF}} & \multicolumn{2}{c|}{$ P^3 $\textbf{: Deblurring}} \\
		\cline{2-3}
		& \textbf{Parameter} & \textbf{Value} \\
		\hline \hline
		\multirow{6}{*}{{\rotatebox[origin=c]{90}{Uniform}}} & N     & 200 \\
		& $ \alpha $ & 1.02 \\
		& $ \beta_0 $ & 0.0007 \\
		& $ \beta_k $ & $ \alpha^k \cdot \beta_0 $ \\
		& $ \lambda $ & $ 512 \cdot \beta_0 $ \\
		& $ \sigma_f $ & $ \sqrt{\lambda/\beta_k} $ \\
		\hline \hline
		\multirow{6}{*}{\rotatebox[origin=c]{90}{Gaussian}} & N     & 200 \\
		& $ \alpha $ & 1.02 \\
		& $ \beta_0 $ & 0.0007 \\
		& $ \beta_k $ & $ \alpha^k \cdot \beta_0 $ \\
		& $ \lambda $ & $ 320 \cdot \beta_0 $ \\
		& $ \sigma_f $ & $ \sqrt{\lambda/\beta_k} $ \\
		\hline
	\end{tabular}%
	\caption{{The set of parameters being used by the $ P^3 $ method, leading to the deblurring results reported in Table~\ref{tab:deb_results} when plugging the TNRD \cite{chen2015trainable} denoiser.}}

	\label{tab:params_ppp_deblurring}%
\end{table}%

\begin{figure}
	\begin{center}
		\subfigure[RED, Denoising Engine: Median Filter]{\includegraphics[width=0.48\linewidth]{./gauss_median_convergence_revised.pdf}}
		\subfigure[RED, Denoising Engine: TNRD]{\includegraphics[width=0.48\linewidth]{./gauss_tnrd_convergence_revised.pdf}}
	\end{center}
	\caption{{An illustration of the convergence of RED using three proposed numerical schemes -- the steepest-descent (red), the ADMM with $m_2 = 1$ (green), the ADMM with $m_2 = 3$ (blue) and the fixed-point (black). These are applied on the image \textsf{Leaves}, degraded by a Gaussian PSF, when two denoising engines are tested: (a) the median filter, and (b) TNRD \cite{chen2015trainable}.}}
	\label{fig:convergence_deb_gauss}
\end{figure}

\begin{figure*}
	\begin{center}
		\subfigure[Ground Truth]{\includegraphics[]{./Starfish_orig_c.png}}
		\subfigure[Input 20.83dB]{\includegraphics[]{./Starfish_in_c.png}}
		\subfigure[RED: SD-Median filter 25.87dB]{\includegraphics[]{./Starfish_med_c.png}}
		\subfigure[NCSR 28.39dB]{\includegraphics[]{./Starfish_ncsr_c.png}}
		\subfigure[$ P^3 $-TNRD 28.43dB]{\includegraphics[]{./Starfish_ppp_c.png}}
		\subfigure[RED: FP-TNRD 28.82dB]{\includegraphics[]{./Starfish_fp_c.png}}
	\end{center}
	\caption{Visual comparison of deblurring a cropped area from the image \textsf{Starfish}, degraded by a uniform PSF, along with the corresponding PSNR [dB] score.}
	\label{fig:comp_vis_deb_ones}
\end{figure*}

\begin{figure*}
	\begin{center}
		\subfigure[Ground Truth]{\includegraphics[]{./leaves_orig_c.png}}
		\subfigure[Input 21.40dB]{\includegraphics[]{./leaves_in_c.png}}
		\subfigure[RED:SD-Median filter 27.87dB]{\includegraphics[]{./leaves_med_c.png}}
		\subfigure[NCSR 30.03dB]{\includegraphics[]{./leaves_ncsr_c.png}}
		\subfigure[$ P^3 $-TNRD 30.36dB]{\includegraphics[]{./leaves_ppp_c.png}}
		\subfigure[RED: ADMM-TNRD 30.40dB]{\includegraphics[]{./leaves_admm_c.png}}
	\end{center}
	\caption{{Visual comparison of deblurring a cropped area from the image \textsf{Leaves}, degraded by a Gaussian PSF, along with the corresponding PSNR [dB] score.}}
	\label{fig:comp_vis_deb_gauss}
\end{figure*}

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 \cite{beck2009fast}. 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 \ref{sec:Parameters} 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~\ref{fig:convergence_deb_gauss}. 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 \ref{fig:comp_vis_deb_ones} and \ref{fig:comp_vis_deb_gauss}, 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.

\subsection{Image Super-Resolution}

Similarly to the previous subsection, we imitate the super-resolution experiments done in \cite{NCSR}. 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~\ref{tab:sr_results} 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 \cite{dong2011image} 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~\ref{tab:params_laplacian_sr} and \ref{tab:params_ppp_sr}, respectively. In the simpler case of the median filter (defined by a $ 3\times3 $ window), the number of iterations of the proposed steepest-descent algorithm is set to $ N=50 $ with a parameter $ \lambda = 0.0325 $.

\begin{table}[htbp]
	\centering
	\tiny

	\renewcommand{\arraystretch}{0.8}
	\renewcommand{\tabcolsep}{7.5pt}
	\begin{tabular}{|c||c|c|c|c|c|c|c|c|c|c|}
		\hline
		\multicolumn{11}{|c|}{\textbf{Super-Resolution, scaling = 3, $ \sigma=\sqrt{2} $}} \\
		\hline
		\textbf{Image} & \textsf{Butterfly} & \textsf{flower} & \textsf{Girl}  & \textsf{Parth.} & \textsf{Parrot} & \textsf{Raccoon} & \textsf{Bike}  & \textsf{Hat}   & \textsf{Plants} & \textbf{Average} \\
		\hline
		Bicubic & 20.74 & 24.74 & 29.61 & 24.04 & 25.43 & 26.20 & 20.75 & 27.00 & 27.62 & 25.13 \\
		ASDS-Reg \cite{dong2011image} & 26.06 & 27.83 & 31.87 & 26.22 & 29.01 & 28.01 & 23.62 & 29.61 & 31.18 & 28.16 \\
		NCSR  \cite{NCSR} & 26.86 & 28.08 & 32.03 & 26.38 & 29.51 & 28.03 & 23.80 & 29.94 & 31.73 & 28.48 \\
		$P^3$-TNRD & 27.13 & \textbf{28.23} & \textbf{32.08} & 26.50 & \textbf{29.65} & 27.95 & \textbf{24.04} & 30.30 & \textbf{31.78} & 28.63 \\
		RED: SD-Median Filter & 24.44 & 27.24 & 31.13 & 25.80 & 27.76 & 27.65 & 22.89 & 28.69 & 30.24 & 27.32 \\
		RED: SD-TNRD & \textbf{27.37} & \textbf{28.23} & \textbf{32.08} & \textbf{26.54} & 29.43 & \textbf{27.98} & \textbf{24.04} & \textbf{30.36} & \textbf{31.79} & \textbf{28.65} \\
		RED: ADMM-TNRD & 27.22	& \textbf{28.24}	& \textbf{32.08}	& 26.51	& 29.41	& \textbf{27.97}	& 23.96	& \textbf{30.35}	& 31.77	& 28.61 \\
		RED: FP-TNRD & 27.26 & \textbf{28.24} & \textbf{32.08} & 26.52 & 29.42 & \textbf{27.97} & 23.97 & \textbf{30.35} & 31.77 & 28.62 \\
		\hline
	\end{tabular}%
	\caption{Super-resolution results, measured in terms of PSNR [dB] and evaluated on the set of images provided by the authors of NCSR \cite{NCSR}. The $ P^3 $ and ours steepest-descent (SD), ADMM, and fixed-point (FP) methods build upon the TNRD \cite{chen2015trainable} as the denoising engine. We also provide the results obtained by integrating the median filter in the steepest-descent scheme. {PSNR scores being less than $ 0.01 $[dB] away from the highest result are highlighted. Similarly to the deblurring case, the $ P^3 $ does not converge here as well. Therefore, we run it for a fixed number of iterations, manually tuned to achieve the highest PSNR on average. This behavior is demonstrated in detail in Figure \ref{fig:red_vs_ppp}, and discussed in Section \ref{sec:Parameters}. }}

	\label{tab:sr_results}%
\end{table}%

\begin{table}[htbp]
	\centering
	\scriptsize

	\renewcommand{\arraystretch}{0.8}
	\renewcommand{\tabcolsep}{9pt}
	\begin{tabular}{|c|c|c|c|}
		\hline
		\multicolumn{4}{|c|}{\textbf{Proposed approach: Super-Resolution}} \\
		\hline
		\multirow{2}{*}{\textbf{Parameter}} & \textbf{Steepest} & \textbf{Fixed} & \multirow{2}{*}{\textbf{ADMM}} \\
		& \textbf{Descent}   & \textbf{Point} &                       \\
		\hline \hline
		N     & 1500  & 200   & 200 \\
		$ \sigma_f $ & $ 3 $ & $ 3 $ & $3 $ \\
		$ \lambda $ & 0.008  & 0.008  & 0.008 \\
		$ m_1 $    & --     & 200 (until convergence)   & 200 (until convergence) \\
		$ m_2 $    & --     & --     & 1 \\
		$ \beta $  & --     & --     & 0.001 \\
		\hline
	\end{tabular}%
	\caption{The set of parameters being used in our framework, leading to the super-resolution results reported in Table~\ref{tab:sr_results} when plugging the TNRD \cite{chen2015trainable} denoiser.}

	\label{tab:params_laplacian_sr}%
\end{table}%

\begin{table}[htbp]
	\centering
	\scriptsize

	\renewcommand{\arraystretch}{0.8}
	\renewcommand{\tabcolsep}{9pt}
	\begin{tabular}{|c|c|}
		\hline
		\multicolumn{2}{|c|}{$ P^3 $\textbf{: Super-Resolution}} \\
		\hline
		\textbf{Parameter} & \textbf{Value} \\
		\hline \hline
		N     & 200 \\
		$ \alpha $ & 1.02 \\
		$ \beta_0 $ & 0.001 \\
		$ \beta_k $ & $ \alpha^k \ \beta_0 $ \\
		$ \lambda $ & $ 360 \ \beta_0 $ \\
		$ \sigma_f $ & $ \sqrt{\lambda/\beta_k} $ \\
		\hline
	\end{tabular}%
	\caption{The set of parameters being used in the $ P^3 $ method, leading to the super-resolution results reported in Table~\ref{tab:sr_results}  when plugging the TNRD \cite{chen2015trainable} denoiser.}

	\label{tab:params_ppp_sr}%
\end{table}%

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 \ref{fig:comp_vis_sr_butterfly} and \ref{fig:comp_vis_sr_bike} 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 \ref{tab:sr_results}.

\begin{figure*}
		\begin{center}
		\subfigure[Bicubic 20.68dB]{\includegraphics[]{./Butterfly_in_c.png}}
		\subfigure[NCSR 26.79dB]{\includegraphics[]{./Butterfly_ncsr_c.png}}
		\subfigure[$ P^3 $-TNRD 26.61dB]{\includegraphics[]{./Butterfly_ppp_c.png}}
		\subfigure[Ours: SD-TNRD 27.39dB]{\includegraphics[]{./Butterfly_sd_c.png}}
	\end{center}
	\caption{{Visual comparison for upscaling by a factor of 3 for the image \textsf{Butterfly}, along with the corresponding PSNR [dB] score.}}
	\label{fig:comp_vis_sr_butterfly}
\end{figure*}

\begin{figure*}
	\begin{center}
		\subfigure[Bicubic 20.44dB]{\includegraphics[]{./bike_in_c.png}}
		\hfil
		\subfigure[NCSR 22.97dB]{\includegraphics[]{./bike_ncsr_c.png}}
		\\
		\subfigure[$ P^3 $-TNRD 23.25dB]{\includegraphics[]{./bike_ppp_c.png}}
		\hfil
		\subfigure[Ours: ADMM-TNRD 23.28dB]{\includegraphics[]{./bike_admm_c.png}}
	\end{center}
	\caption{{Visual comparison for upscaling by a factor of 3 for the image \textsf{Bike}, along with the corresponding PSNR [dB] score.}}
	\label{fig:comp_vis_sr_bike}
\end{figure*}

\subsection{Robustness to the Choice of Parameters} \label{sec: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 \ref{fig:red_vs_ppp} (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 \ref{fig:red_vs_ppp} (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 \ref{fig:red_vs_ppp} (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 \ref{fig:red_vs_ppp} (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 \ref{fig:algoADMM} (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 \ref{fig:red_vs_ppp} (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 \ref{fig:red_vs_ppp} (d) for the FP method\footnote{One 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.

\begin{figure*}
	\begin{center}
		\subfigure[Comparison between RED and $ P^3 $]{\includegraphics[width=0.48\linewidth]{./sr_red_vs_ppp.pdf}}
		\subfigure[$ P^3 $: Sensitivity to change in $ \beta_0 $ and $ \alpha $]{\includegraphics[width=0.48\linewidth]{./sr_ppp_parameters.pdf}}
		\subfigure[RED: Robustness to change in $ \beta $ and $ m_2 $]{\includegraphics[width=0.48\linewidth]{./sr_admm_red_parameters.pdf}}
		\subfigure[RED: The effect of a change in $ \sigma_f$]{\includegraphics[width=0.48\linewidth]{./sr_red_fp_sigmaf.pdf}}
	\end{center}
	\caption{{Image super resolution. Average PSNR as a function of the iterations of the different methods. The test images are the ones of Table \ref{tab:sr_results}.}}
	\label{fig:red_vs_ppp}
\end{figure*}

\section{Conclusions}\label{sec: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 {\em 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 \cite{P3}, 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.

redpajama_set_name: RedPajamaArXiv


For athletes with cerebral palsy (CP), understanding how well they train is extremely important. It is important that when they undergo strength training that their unaffected side isn't compensating for their affected side. The coaches and trainers play a crucial role in monitoring the athletes and ensuring that they are training correctly. However, it is difficult to identify if the athletes are training correctly. Therefore, technology can be used to measure the athletes load distributions during a strength training exercise to guarantee both limbs are being trained equally. Devices such as force plates can be extremely expensive and hard to use, stopping coaches from monitoring their athletes. A possible solution to this is the Nintendo Wii Balance Board (WBB). The WBB is capable of measuring load distributions by manipulating the weight measurements obtained from its 4 strain gauges. The portability and low cost of the Wii Balance Board make it an ideal device.