Averaged Deep Denoisers for Image Regularization

Classical algorithms such as FBS and DRS perform repeated inversion of the measurement model, followed by a regularization step. More precisely, the regularization is performed using the proximal operator. For a convex regularizer $h:\mathbb{R}^{n}\to\mathbb{R}$, the proximal operator is defined as

where $\|\cdot\|$ is the Euclidean norm. The proximal operator can be seen as a denoiser — it takes a noisy image $\boldsymbol{x}$ and returns an image close to $\boldsymbol{x}$ for which the penalty $h$ is small beck2009fast . This operator is at the core of FBS, DRS, and ADMM. In particular, the FBS update is given by

where $\gamma>0$ and $\boldsymbol{x}_{0}\in\mathbb{R}^{n}$ is an arbitrary initialization.

On the other hand, starting with $\boldsymbol{x}_{0}\in\mathbb{R}^{n}$ and $\gamma>0$, the DRS update is given by (e.g., see ryu2016primer )

In Plug-and-Play (PnP) regularization sreehari2016plug , instead of having to handcraft the regularizer $g$ in (1) and apply it via its proximal operator, we directly “plug” a powerful denoising operator $\mathrm{D}:\mathbb{R}^{n}\to\mathbb{R}^{n}$, such as BM3D dabov2007image or DnCNN zhang2017beyond , into (1.1) and (1.1). Thus, in the PnP variant of FBS (referred to as PnP-FBS), the update is performed by replacing the proximal operator in (1.1) by an image denoiser (Algorithm 1).

Remarkably, PnP algorithms yield state-of-the-art results across a wide range of applications, including superresolution, MRI, fusion, and tomography sreehari2016plug ; Ryu2019_PnP_trained_conv ; Sun2019_PnP_SGD ; zhang2017learning ; Tirer2019_iter_denoising ; zhang2021plug . The best performing PnP algorithms zhang2017learning ; Tirer2019_iter_denoising ; zhang2021plug use trained denoisers such as DnCNN zhang2017beyond , IRCNN zhang2017learning , and U-Net hurault2022gradient ; hurault2022proximal . This also opens up the exciting possibility of combining model-based reconstruction with deep learning zhang2017learning .

The idea of using denoisers for regularization has also been explored in Regularization-by-Denoising (RED). For example, the authors in romano2017little proposed to work with the regularizer

where $\mathrm{D}$ is some powerful denoiser as in PnP. This is solved iteratively using gradient-based methods, where the gradient of (5) is approximated using $\boldsymbol{x}-\mathrm{D}(\boldsymbol{x})$ in romano2017little . It was later shown in reehorst2018regularization that this is exact if $\mathrm{D}$ is locally homogenous and its Jacobian is symmetric.

Several variants of the original RED algorithm romano2017little have been proposed reehorst2018regularization ; Sun2019_PnP_SGD ; sun2019block ; sun2021scalable ; cohen2021regularization . One such variant called RED-PG reehorst2018regularization is described in Algorithm 2. RED algorithms have been shown to yield state-of-the-art results for phase retrieval, tomography, superresolution, deblurring, and compressed sensing metzler2018prdeep ; wu2020simba ; wu2019online ; mataev2019deepred ; sun2019block ; hu2022monotonically .

There is yet another paradigm for solving inverse problems called deep unfolding gregor2010learning ; sun2016deep ; zhang2018ista , where an iterative algorithm is unfolded into a network and trained in an end-to-end fashion. Deep unfolded networks can produce impressive results, but unlike PnP or RED (where just the denoiser needs to be trained), we have to train an unfolded network end-to-end. The present work is motivated by deep unfolding, but we propose to use them for a different purpose, namely, to construct nonexpansive denoisers by unfolding classical wavelet denoising. We then apply the trained denoiser within Algorithms 1 and 2 for image reconstruction.

We wish to point to a related work repetti2022dual , where the authors construct an unfolded denoiser by applying FBS to a sparsity-regularized image denoising problem (with box constraints on the reconstruction). However, there are important distinctions between our approach and theirs. First, splitting is applied to the Fenchel dual of the denoising problem in repetti2022dual ; the idea is to transform the nonsmooth primal (denoising) problem into a smooth, unconstrained dual problem. Second, since the sparsifying linear transform is optimized in repetti2022dual , there is no structural guarantee on the trained layers. As a result, it is difficult to come up with convergence guarantees; indeed, (repetti2022dual, , Proposition 2) assumes that the same linear transform is used across all layers, and this would generally not hold if we optimize the transforms. We present comparisons with repetti2022dual in Section 5.

The success of PnP and RED has sparked interest in the theoretical understanding of these algorithms. For example, it has been shown that for a class of linear denoisers, PnP is guaranteed to converge, even at a linear rate ACK2023-contractivity , and the limit point is the solution of a convex optimization problem sreehari2016plug ; nair2021fixed ; gavaskar2021plug . This is done by expressing the linear denoiser as the proximal operator of a convex function. Unfortunately, this does not hold for powerful nonlinear denoisers such as BM3D, DnCNN, or IRCNN. This is because, unlike the proximal operator that is nonexpansive, such denoisers have been shown to violate the nonexpansive property reehorst2018regularization .

Though developing a variational characterization for deep denoisers is challenging, it is nonetheless possible to establish iterate convergence for PnP and RED algorithms. For example, the convergence of Algorithm 1 was established for generative denoisers in jagatap2019algorithmic ; liu2021recovery ; raj2019gan . For strongly convex loss functions, PnP-FBS convergence can be guaranteed using a CNN denoiser whose residual is nonexpansive Ryu2019_PnP_trained_conv . Motivated by RED, convergence was recently established for PnP modeled on FBS, ADMM, half-quadratic splitting, and gradient descent hurault2022gradient ; cohen2021has ; hurault2022proximal . On the other hand, the convergence of RED algorithms can be guaranteed using a nonexpansive denoiser romano2017little ; reehorst2018regularization ; sun2019block ; cohen2021regularization ; liu2021recovery . For instance, the block-coordinate RED algorithm in sun2019block exhibits convergence with a nonexpansive denoiser. In liu2021recovery , the authors established convergence of RED for the compressed sensing application using a nonexpansive denoiser.

The focus of this work is on strong convergence of the PnP and RED algorithms, where the iterates converge in norm to a fixed point of a well-defined operator ryu2016primer . It is important to note that there exist other forms of iterate convergence. For example, for the denoisers in hurault2022gradient ; hurault2022proximal , the limit points of the PnP-FBS iterates are stationary points of a nonconvex objective function. In contrast, the iterates converge weakly (see (bauschke2017convex, , Section $2.5$) for a definition) to a zero of a combination of a gradient operator and a firmly nonexpansive operator in pesquet2021learning . On the other hand, the iterates of Plug-and-Play Stochastic Gradient Descent (PnP-SGD) in laumont2023maximum converge to a stationary point of a MAP problem, where the prior is derived from the denoiser. The summary of our understanding of the iterate convergence of PnP and RED is as follows:

1. 1.

   Convergence of PnP is guaranteed for averaged denoisers Sun2019_PnP_SGD ; sun2021scalable ; nair2021fixed .

2. 2.

   For nonexpansive denoisers, convergence is guaranteed for variants of RED reehorst2018regularization ; sun2019block ; cohen2021regularization , PnP-FBS applied to compressive sensing liu2021recovery , and a variant of PnP-FBS reehorst2018regularization .

For the above reasons, there has been significant research on training averaged and nonexpansive CNN denoisers. The challenge is that computing the Lipschitz constant of a neural network is NP-hard virmaux2018lipschitz . By Lipschitz constant, we mean the smallest $c$ such that $\lVert\mathrm{T}(\boldsymbol{u})-\mathrm{T}(\boldsymbol{v})\rVert\leqslant c\lVert\boldsymbol{u}-\boldsymbol{v}\rVert$ for all inputs $\boldsymbol{u}$ and $\boldsymbol{v}$, where we represent the CNN as an operator $\mathrm{T}$. In particular, while the Lipschitz constant of a convolution layer is just the largest singular value of the matrix representation of the layer (which can be computed efficiently from the FFT of its filters sedghi2018singular ) and the Lipschitz constant of common nonlinear activations (ReLU and sigmoid) is $1$, computing the Lipschitz constant for their cascaded composition is however nontrivial. Nonetheless, we can compute an upper bound for the Lipschitz constant using a direct approximation for the spectral norm of the network pesquet2021learning or the product of the spectral norms of the convolution layers Ryu2019_PnP_trained_conv . In particular, if we force the spectral norm of the convolution layers to be within $1$ at each epoch of the training process, then the final network is automatically nonexpansive. Several methods based on this idea have been proposed sedghi2018singular ; hertrich2021convolutional ; terris2020building . In particular, we will later compare with the nonexpansive CNN denoiser in Ryu2019_PnP_trained_conv , where the largest singular value of a convolutional layer (estimated using power iterations) is used to force the Lipschitz constant of the layer to be within $1$ after each gradient update.

This work was motivated by the observation that existing CNN denoisers tend to violate the nonexpansive property (see Fig. 5). In fact, existing methods for training nonexpansive denoisers either cannot guarantee the final network to be nonexpansive or are computationally expensive. For example, power iterations are used for approximate computation of the largest singular value in Ryu2019_PnP_trained_conv ; pesquet2021learning ; hence, the network is not guaranteed to be nonexpansive at each epoch. The methods in sedghi2018singular ; hertrich2021convolutional ; terris2020building for constraining the convolutional layers are computationally expensive; e.g., the time taken in sedghi2018singular is about nine-fold that required to train the network without constraints. This is because sedghi2018singular requires the SVD computation of a large matrix derived from the Fourier transform of the convolutional filters. Similarly, terris2020building requires the Fourier transform of the filters in every convolutional layer. On the other hand, the training in hertrich2021convolutional requires us to solve an optimization problem that imposes orthogonality constraints on the convolutional layers.

We also highlight a technical issue that has been sidelined in existing works Ryu2019_PnP_trained_conv ; sedghi2018singular ; hertrich2021convolutional ; terris2020building ; pesquet2021learning , namely, that the CNN filters are trained on images of fixed size but deployed on images of arbitrary size. Since the filters must be zero-padded to match the image size, their Lispchitz can change with the size. As a concrete example, consider the $3\times 3$ Sobel filter $[1,0,-1\mathchar 24635\relax\;2,0,-2\mathchar 24635\relax\;1,0,-1]$. The Lipschitz constant of this filter for an image of size $25\times 25$ is $7.9842$, whereas the Lipschitz constant for a $256\times 256$ image is $8$. Thus, a guarantee that the Lipschitz bound is independent of the image size is not available for existing CNN denoisers. To substantiate our point, we present numerical evidence (see Fig. 5 in Section 5) showing that BM3D dabov2007image , DnCNN zhang2017beyond , and N-CNN Ryu2019_PnP_trained_conv violate the nonexpansive property. Such counterexamples have also been reported in reehorst2018regularization . Our observations are summarized below:

1. 1.

   The Lipschitz constant of a CNN (trained with fixed-size filters) depends on the size of the input image and is difficult to predict.

2. 2.

   Existing CNN denoisers violate the nonexpansive property, which can result in the divergence of PnP and RED algorithms (see Tables 2 and 6 in Section 5).

3. 3.

   Developing efficient algorithms for training nonexpansive CNNs is challenging. Existing algorithms for training nonexpansive CNNs either cannot guarantee nonexpansivity or are computationally expensive.

Following the above observations, we moved away from CNNs and considered deep unfolded networks instead. In fact, from our experience, we noticed a tradeoff between nonexpansivity and the denoising performance of CNNs — a significant drop in denoising performance was observed if the denoiser was forced to be nonexpansive. Indeed, it was shown in hurault2022proximal ; hertrich2021convolutional that imposing hard Lipschitz constraints on the CNN can adversely affect its denoising performance.

In the rest of the paper, we demonstrate that it is possible to construct deep unfolded networks that come with the following desirable properties:

1. 1.

   , Unlike Ryu2019_PnP_trained_conv ; sedghi2018singular , they are averaged (or contractive) by construction and hence automatically nonexpansive.

2. 2.

   The training is not expensive compared to sedghi2018singular ; terris2020building ; pesquet2021learning . In particular, we need to only perform cheap projections onto intervals (of the learnable parameters) after every gradient step.

3. 3.

   When plugged into Algorithms 1 and 2, the reconstruction is comparable with the state-of-the-art.

Our construction is motivated by the averaged (contractivity) property of fixed-point operators arising in FBS and DRS; these are classically used to establish iterate convergence bauschke2017convex . In fact, the building blocks of our denoiser are obtained by unfolding FBS and DRS applied to wavelet denoising. These blocks are provably averaged (resp. contractive), and the proposed denoiser, being a composition of these building blocks, is automatically averaged (resp. contractive). While our denoiser is trained on small patches, we give a formula for patch aggregation that guarantees the final end-to-end image denoiser to be averaged (or contractive). The denoising capacity of our denoiser is short of BM3D and DnCNN; however, their regularization capacity for PnP and RED is comparable with BM3D and DnCNN. The source code is available in https://github.com/pravin1390/Averageddeepdenoiser.

Let $\boldsymbol{b}\in\mathbb{R}^{p}$ represent a noisy image patch, where we have stacked the pixels of a $(k\times k)$ image patch into a vector of length $p=k^{2}$. Let $\mathbf{W}$ be an orthogonal wavelet transform. Consider the wavelet denoising problem

where $\lambda$ is the regularization parameter. The above problem has the closed-form solution $\mathbf{W}^{\top}\phi_{\lambda}(\mathbf{W}\boldsymbol{b})$, where the soft-threshold function

is applied componentwise on $\mathbf{W}\boldsymbol{b}$ and $\lambda$ acts as a threshold chambolle1998nonlinear .

We consider a slightly different model with a different threshold for each wavelet coefficient. Namely, we consider the optimization problem

This also has a closed-form solution given by

where $\Lambda=(\lambda_{1},\ldots,\lambda_{p})$ are the thresholds and the thresholding operator $\mathbf{S}_{\Lambda}:\mathbb{R}^{p}\to\mathbb{R}^{p}$ is defined as

The denoising operator $\phi_{\Lambda}$ is nonexpansive, being the composition of nonexpansive and unitary operators. However, it is easy to see that $\phi_{\Lambda}$ is not contractive. To identify a contractive (resp. averaged) operator, we propose to solve the denoising problem (10) in an iterative fashion using FBS (resp. DRS). The proposed denoising network is obtained by chaining the update operator in FBS and DRS.

As a first step, we apply FBS to the optimization problem (10). Note that, a stepsize $\rho$, we can write the FBS update as

where we define the FBS block $\mathrm{U}$ to be

The FBS block is depicted in Fig. 1. Functionally, $\mathrm{U}$ is an inversion step (gradient descent on the loss) followed by wavelet regularization.

Although (10) has a closed-form solution, we are interested in the iterative solution of (10) because of the following property.

Indeed, since $\mathbf{W}$ is an orthogonal transform and the operator $\mathbf{S}_{\Lambda}$ in (12) is nonexpansive (the condition $\Lambda>0$ is required to ensure that (11) is well-defined), we have for any $\boldsymbol{x},\boldsymbol{x}^{\prime}\in\mathbb{R}^{n}$,

Since $1-\rho\in[0,1)$, we can conclude that the map $\boldsymbol{x}\mapsto\mathrm{U}(\boldsymbol{x}\mathchar 24635\relax\;\rho,\Lambda)$ is contractive.

Following Proposition 7, we chain $\mathrm{U}$ into a network. The precise definition is as follows.

Since the composition of contractive operators is again contractive, the next result follows immediately from Proposition 7.

The point is that (13) is contractive as long as $\Lambda^{(\ell)}$ and $\rho^{(\ell)}$ are in the interval stipulated in Proposition 8. We can thus tune the parameters over these intervals to optimize the denoising capacity while preserving the contractive nature.

Note that we can apply Peaceman-Rachford bauschke2017convex or DRS to problem (10). In fact, we will show that the corresponding update operator is firmly nonexpansive for DRS. Thus, the resulting denoiser obtained by chaining this operator is averaged, and we can guarantee convergence of Algorithms 1 and (2) (Corollary 6).

To arrive at the DRS update for (10), we identify (10) with (1) by setting

It follows from (1.1) that we can write the DRS update as written as

where the DRS block $\mathrm{V}$ is defined as

and $\rho>0$ is the penalty parameter of DRS.

Having defined the update operator $\mathrm{V}$, we can now establish its averaged property.

Indeed, we know that the proximal operator of a convex function is firmly nonexpansive bauschke2017convex . Thus, the operators $2\mathrm{Prox}_{\rho g}-\mathrm{I}$ and $2\mathrm{Prox}_{\rho f}-\mathrm{I}$ are nonexpansive. Thus, if we denote their composition by $\mathrm{N}$, then $\mathrm{N}$ must be nonexpansive. Since $\mathrm{V}=(1/2)(\mathrm{I}+\mathrm{N})$, $\mathrm{V}$ must be firmly nonexpansive.

We will now use Proposition 9 to construct an averaged network. As with the FBS block, the trainable parameters in (8) are $\rho$ and $\Lambda$, while $\boldsymbol{b}$ and $\mathbf{W}$ are fixed. The DRS block is depicted in Fig. 1. The DRS denoiser is obtained by chaining $\mathrm{V}$ into a network.

As a consequence of Proposition 9, we have the following structure on the denoiser.

Indeed, notice from Proposition 9 that $\mathrm{D}_{\mathrm{DRS}}$ is the composition of $L$ firmly nonexpansive operators. Now, if we compose $L$ number of $\theta$-averaged operators, then it can be shown that the resulting operator is $L\theta/(1+(L-1)\theta)$-averaged (e.g., see bauschke2017convex ). Setting $\theta=1/2$, we arrive at Proposition 10.

Recall that (13) and (16) are patch denoisers. We apply them independently on patches extracted from a noisy image and then average the denoised patches to obtain the final image zhang2017image ; dabov2007image . The overall operation defines an image-to-image denoising operator $\mathrm{D}$, which is used in Algorithms 1 and 2. In view of Corollary 3 and 6, we wish to ensure that $\mathrm{D}$ is averaged or contractive. First, we give a precise definition of the denoiser and then prove that it has the desired property.

Let $\Omega=[0,q-1]^{2}$ denote the image domain and $\mathbf{X}:\Omega\to\mathbb{R}$ be the input (noisy) image; we also view $\mathbf{X}$ as a vector in $\mathbb{R}^{n}$ where $n=q^{2}$. We assume periodic boundary condition, i.e., for $\boldsymbol{i}=(i_{1},i_{2})\in\Omega$ and $\boldsymbol{\tau}=(\tau_{1},\tau_{2})\in\mathbb{Z}^{2}$, we define

Without loss of generality, we assume that $q$ is a multiple of the patch size $k$; we can simply pad the image to satisfy this condition.

Let the patch size be $k\times k$. For $\boldsymbol{i}\in\Omega$, we define the patch operator

where $\mathrm{vec}(\cdot)$ denotes the vectorized representation of the pixels in a $(k\times k)$ window around $\boldsymbol{i}$ (in some fixed order).

We extract patches from $\mathbf{X}$ with stride $s$, where we assume that $s$ divides $k$. Thus, the starting coordinates of the patches are $\mathcal{J}=[0,(q_{s}-1)]_{s}^{2}$, where $q_{s}=q/s$. In other words, the patches are

The factor $(k/s)^{2}$ is the number of patches containing a given pixel; the assumption that $s$ divides $k$ is essential to ensure this is the same for every pixel.

Henceforth, we will assume that the patch denoiser is either $\mathrm{D}_{{\mathrm{FBS}}}$ in (13) or $\mathrm{D}_{{\mathrm{DRS}}}$ in (16). We now state the main result of this section.

The proof is given in Appendix 8. The core idea is to express the original system of overlapping patches in terms of non-overlapping patches. This allows us to use some form of orthogonality in the calculations.

We trained (13) and (16) on patches extracted from images of the BSD500 dataset. Let $\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{N}\in\mathbb{R}^{p}$ be the clean patches (normalized to $[0,1]$ intensity), where $p=k^{2}$ and $N$ is the size of traning data. The noisy patches are generated as follows:

where $\boldsymbol{\eta}\sim\mathcal{N}(\boldsymbol{0},\mathrm{I})$ and $\sigma$ is the noise level. Thus, the training examples are $({\boldsymbol{z}}_{1},\boldsymbol{x}_{1}),\ldots,({\boldsymbol{z}}_{N},\boldsymbol{x}_{N})$.

Let us use $\Theta$ to denote the parameters in (13) and (16). We will write the patch denoiser as $\mathrm{D}_{\text{patch}}(\cdot\,\mathchar 24635\relax\;\Theta)$ to account for the dependence on $\Theta$. The training problem is to optimize $\Theta$ such that $\mathrm{D}_{\text{patch}}({\boldsymbol{z}}_{i}\mathchar 24635\relax\;\Theta)\approx\boldsymbol{x}_{i}$ for all $i=1,\ldots,N$. Moreover, in view of Proposition 8, we consider the following optimization problem for training (13):

The lower bounds $\rho_{0}\in(0,1)$ and $\lambda_{0}>0$ are used to ensure that the constraint set in (19) is closed; this will be useful when we project $\Theta$ onto the intervals.

On the other hand, in view of Proposition 10, we consider the following optimization for training (16):

In Table 1, we compare the denoising performance of our denoisers with BM3D, DnCNN, and N-CNN on the Set12 dataset (see Fig. 2). Our denoisers are off by PSNR of $0.25\mbox{-}1.5$ dB and SSIM of $0.02\mbox{-}0.08$ from the best denoiser DnCNN. However, they consistently outperform N-CNN across all noise levels. Recall that the purpose of our construction was to guarantee convergence of PnP and RED. In pursuit of this, we had to make design choices that resulted in a tradeoff in denoising performance. Importantly, we will see that the reconstruction capability of our denoiser is competitive with the above denoisers. This is not surprising since there is no concrete relation between the denoising performance and the restoration capacity of a denoiser. For instance, it has been previously shown that BM3D can yield superior restoration results compared to DnCNN in specific applications Ryu2019_PnP_trained_conv , even though DnCNN generally exhibits better denoising capabilities than BM3D. This suggests that the denoising gap is compensated in PnP and RED because the denoiser is applied in every iteration. However, this is just an intuitive guess, and we do not have a precise mathematical explanation.

We conducted a comprehensive assessment of various parameter configurations for the denoiser and selected the most suitable settings. Regarding the impact of $\mathbf{W}$, our results indicate that the choice of wavelet transform does not significantly affect the denoising performance. In contrast, we observed a significant enhancement in denoising performance by increasing the number of layers $L$ from $1$ to $10$. However, increasing $L$ beyond $10$ yields marginal improvements in the generalization capability. Thus, we capped $L$ at $10$ to keep the model compact without compromising the denoising performance. Regarding stride $s$, our findings align with common trends in patch aggregation algorithms, namely, a higher stride is associated with decreased accuracy but reduced computational demands. We found that $s=8$ provides an optimal balance between accuracy and efficiency. Finally, the patch size plays a crucial role in the denoising performance. Keeping the stride constant, enlarging the patch size from $64$ to $128$ results in an accuracy improvement ranging from $0.10$ to $0.25$ dB and an enhancement in the SSIM score by $0.01$ to $0.02$. However, this adjustment also leads to an increase in the number of parameters of the denoiser.

The forward model for deblurring is $\boldsymbol{y}=\mathbf{F}\boldsymbol{\xi}+\boldsymbol{\eta}$, where $\boldsymbol{y}$ is the blurred image, $\mathbf{F}$ is a blur operator, and $\boldsymbol{\eta}$ is white Gaussian noise dong2011image . For deblurring, the loss function is $f(\boldsymbol{x})=\|\boldsymbol{y}-\mathbf{F}\boldsymbol{x}\|^{2}$, where $\mathbf{F}$ is the blur matrix. The gradient $\nabla f$ is $\beta$-Lipschitz. Specifically, we use a normalized blur for all experiments; as a result, we have $\beta=1$ ($\beta$ is the spectral norm of $\mathbf{F}$). For the deblurring experiments in Table 2 and Fig. 3, we took $\mathbf{F}$ to be a Gaussian blur of width $25\times 25$ and standard deviation $1.6$ and the noise level is $\sigma=0.04$. For deblurring using PnP-FBS, the step size is set to $\gamma=0.001$ for RED-PG; for PnP-FBS, we used $\gamma=0.5$ which satisfies the condition $\gamma\leqslant 2/\beta$ in Corollaries 3 and 6. For PnP-FBS, we used a noise level $\sigma=5$ for both our denoisers, and $\sigma=10$ for RED-PG.

In Table 2, we present numerical evidence to show that, in the absence of a theoretical guarantee, simply plugging a denoiser (in this case DnCNN) into PnP-FBS can cause the iterations to diverge. Indeed, the distance between successive iterates $\lVert\boldsymbol{x}_{k}-\boldsymbol{x}_{k-1}\rVert$ is shown to diverge for DnCNN, which means that the sequence $(\boldsymbol{x}_{k})$ does not converge. This results in spurious reconstructions, as evident from the PSNR values in Table 2. On the other hand, we notice that $\lVert\boldsymbol{x}_{k}-\boldsymbol{x}_{k-1}\rVert$ is decreasing for the proposed denoisers, thanks to their nonexpansive property.

Next, we present numerical evidence on PnP-FBS convergence using the proposed denoisers. This is done using a deblurring experiment on the Set 12 dataset (the settings are the same as in Table 2). The results are shown in Fig. 3 (a-b), where we plot $\lVert\boldsymbol{x}_{k}-\boldsymbol{x}_{k-1}\rVert$ and PSNR for different $k$. As expected, $\lVert\boldsymbol{x}_{k}-\boldsymbol{x}_{k-1}\rVert$ is decreasing and the PSNR stabilizes in less than $50$ iterations. We also demonstrate that due to the contractivity of the FBS operator $\mathrm{T}_{\mathrm{PnP}}$ using our denoiser $\mathrm{D}_{{\mathrm{FBS}}}$, the iterates converge to a unique fixed point even if we start with diverse initializations; the final reconstructions are identical, as shown in Fig 3(c), where the MSE between the different reconstructions is of the order $1\mathrm{e}\mbox{-}8$.

We next compare the regularization capacity of our denoisers with BM3D, DnCNN, and N-CNN for deblurring, where $\mathbf{F}$ is a horizontal motion blur. This is done for the RED-PG algorithm. The results are shown for a couple of images in Fig. 4. For the two images from Fig. 2, which have a mix of delicate textures, smooth regions, and broader features (top row), the reconstruction from our denoisers exhibits finer details than other denoisers.

A detailed comparison for RED-PG is provided in Table 3, where $\mathbf{F}$ is one of the blur kernels from levin2009understanding . We notice that our performance is comparable with DnCNN.

In Table 4, we compare $\mathrm{D}_{{\mathrm{DRS}}}$ with the denoiser in repetti2022dual , which we term as “DualFB”. We replicated the experimental setup used for the deblurring experiment in (repetti2022dual, , Table 1) — an asymmetric blur operation followed by the addition of noise at various levels. For a fair comparison, similar to the experiment in repetti2022dual , we use PnP-FBS as the reconstruction algorithm for our denoiser and signal-to-noise ratio (SNR) as the performance metric. We see from Table 4 that our denoiser can closely parallel the performance of DualFB. In fact, it is not surprising that DualFB yields superior results since the linear layers are optimized in repetti2022dual , whereas we use fixed handcrafted transforms.

We next compare $\mathrm{D}_{{\mathrm{DRS}}}$ with the “Prox-PnP” denoiser from hurault2022proximal ; we perform the comparison by using them as regularizers within the PnP-FBS framework. We replicated the experimental setup used for the deblurring experiment in (hurault2022proximal, , Table 3) in which the PSNR is averaged over the CBSD68 dataset for $10$ different kernels and different noise levels. As shown in Table 5, our denoiser exhibits comparable reconstruction capacity to Prox-PnP, though our gray-scale denoiser is applied channel-by-channel. It is worth noting that our denoiser is contractive by construction, whereas the denoiser in hurault2022proximal is “softly” enforced to be $1$-Lipschitz using a loss-based penalization. We assume that training our denoiser specifically for color images can boost the color deblurring performance.

In light of the above comparisons, we wish to emphasize that the key advantage of our denoiser is that it comes with convergence guarantees. On the other hand, the convergence result in (repetti2022dual, , Proposition 2) assumes that the same linear transform is used across all layers, and this is generally not true if we optimize the transforms. It remains to be seen if it is possible to optimize our construction to match the denoiser in repetti2022dual .

A widely used model for superresolution is $\boldsymbol{y}=\mathbf{S}\mathbf{B}\boldsymbol{\xi}+\boldsymbol{\eta},$ where $\mathbf{B}\in\mathbb{R}^{n\times n}$ and $\mathbf{S}\in\mathbb{R}^{m\times n}$ are the blur and subsampling operators ($m=n/k$, $k>1$), and $\boldsymbol{\eta}$ is Gaussian noise chan2017plug ; nair2021fixed . For superresolution, the loss function is $f(\boldsymbol{x})=\|\boldsymbol{y}-\mathbf{S}\mathbf{B}\boldsymbol{x}\|^{2}$, where $\mathbf{S}$ is a decimation matrix and $\mathbf{B}$ is a blur matrix. As in deblurring, we used a normalized blur. In this case, $\beta$ is bounded by the product of the spectral norms of $\mathbf{S}$ and $\mathbf{B}$, so we have $\beta\leqslant 1$. In PnP-FBS, we used the step size $\gamma=1.0$ that satisfies the condition $\gamma\leqslant 2/\beta$ in Corollaries  3 and 6. We set the noise level to $\sigma=5$ for both our denoisers for all experiments.