Linearized ADMM and Fast Nonlocal Denoising for Efficient Plug-and-Play Restoration

TL;DR

This paper introduces efficient linearized ADMM and a fast nonlocal means denoising algorithm to improve plug-and-play image restoration, enabling faster computations and guaranteed convergence in applications like super-resolution and photon imaging.

Contribution

It proposes a linearized ADMM approach for more efficient inversion and constraints integration, along with a significantly faster NLM denoising algorithm, ensuring convergence in plug-and-play frameworks.

Findings

Linearized ADMM reduces inversion cost.

Fast NLM is 80x faster than brute-force.

Effective in super-resolution and photon imaging.

Abstract

In plug-and-play image restoration, the regularization is performed using powerful denoisers such as nonlocal means (NLM) or BM3D. This is done within the framework of alternating direction method of multipliers (ADMM), where the regularization step is formally replaced by an off-the-shelf denoiser. Each plug-and-play iteration involves the inversion of the forward model followed by a denoising step. In this paper, we present a couple of ideas for improving the efficiency of the inversion and denoising steps. First, we propose to use linearized ADMM, which generally allows us to perform the inversion at a lower cost than standard ADMM. Moreover, we can easily incorporate hard constraints into the optimization framework as a result. Second, we develop a fast algorithm for doubly stochastic NLM, originally proposed by Sreehari et al. (IEEE TCI, 2016), which is about 80x faster than brute-force computation. This particular denoiser can be expressed as the proximal map of a convex regularizer and, as a consequence, we can guarantee convergence for linearized plug-and-play ADMM. We demonstrate the effectiveness of our proposals for super-resolution and single-photon imaging.

Tables

Table 1. Table 1 : Timings (sec) of the brute-force and fast implementations of DSG-NLM for different patch size and fixed window size N s = 21 subscript 𝑁 𝑠 21 N_{s}=21 .

Patch Size ($N_p$)	$11\times 11$	$17\times 17$	$23\times 23$	$29\times 29$
Proposed	$3.90$	$4.01$	$4.12$	$4.25$
Brute-force	$299.54$	$326.03$	$366.65$	$439.19$

Table 2. Table 2 : Primal and dual residuals and corresponding PSNRs (after 250 iterations) for image super-resolution (periodic boundary).

Error	BM3D	NLM	A-DSG-NLM	F-DSG-NLM
Primal	$3.10\times {10}^{-3}$	$7.30\times {10}^{-4}$	$1.36\times {10}^{-5}$	$2.44\times {10}^{-8}$
Dual	$1.01\times {10}^{-1}$	$4.73\times {10}^{-2}$	$1.42\times {10}^{-4}$	$2.81\times {10}^{-6}$
PSNR	$32.43$	$31.10$	$32.56$	$32.61$

Table 3. Table 3 : Same as in Table 2 but for single photon imaging.

Error	BM3D	NLM	A-DSG-NLM	F-DSG-NLM
Primal	$6.39\times {10}^{-4}$	$1.22\times {10}^{-10}$	$1.72\times {10}^{-13}$	$1.83\times {10}^{-15}$
Dual	$1.28\times {10}^{-2}$	$1.36\times {10}^{-8}$	$1.32\times {10}^{-11}$	$1.35\times {10}^{-12}$
PSNR	$32.55$	$30.98$	$32.60$	$33.24$

Equations

$$g(\mathbf{x})+\beta \Vert \mathbf{x}-\tilde{\mathbf{x}}{\Vert }^2(\beta >0)$$

$$\underset{\mathbf{x}}{\min }f(\mathbf{x})+\lambda g(\mathbf{x})(\lambda >0)$$

$$\underset{\mathbf{x}}{\min }f(\mathbf{x})+\beta \Vert \mathbf{x}-\tilde{\mathbf{x}}{\Vert }^2$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^n}{\min }f(\mathbf{x})+\lambda g(\mathbf{x})\text{s.t.}\mathbf{x}\in C,$$

$$\underset{\mathbf{x},\mathbf{v}\in {\mathbb{R}}^n}{\mathrm{argmin}}f(\mathbf{x})+\lambda g(\mathbf{v})\text{s.t.}\mathbf{x}=\mathbf{v}\text{ and }\mathbf{x}\in C.$$

$$\mathfrak{L}(\mathbf{x},\mathbf{v},\mathbf{u})=f(\mathbf{x})+\lambda g(\mathbf{v})+\frac{\rho }{2}\Vert \mathbf{x}-\mathbf{v}+\mathbf{u}{\Vert }^2-\frac{\rho }{2}\Vert \mathbf{u}{\Vert }^2,$$

$${\mathbf{x}}^{k+1}=\underset{\mathbf{x}\in C}{\mathrm{argmin}}f(\mathbf{x})+\frac{\rho }{2}\Vert \mathbf{x}-{\tilde{\mathbf{x}}}^k{\Vert }^2,$$

$${\mathbf{v}}^{k+1}=\underset{\mathbf{v}\in {\mathbb{R}}^n}{\mathrm{argmin}}\frac{\rho }{2\lambda }\Vert \mathbf{v}-{\tilde{\mathbf{v}}}^k{\Vert }^2+g(\mathbf{v}),$$

$${\mathbf{u}}^{k+1}={\mathbf{u}}^k+({\mathbf{x}}^{k+1}-{\mathbf{v}}^{k+1}).$$

$${\mathbf{v}}^{k+1}={\mathfrak{D}}_{\sigma }({\tilde{\mathbf{v}}}^k),$$

$$f({\mathbf{x}}^k)+\nabla f({\mathbf{x}}^k{)}^{\top }(\mathbf{x}-{\mathbf{x}}^k)+\frac{\alpha }{2}\Vert \mathbf{x}-{\mathbf{x}}^k{\Vert }^2,$$

$${\mathbf{x}}^{k+1}=\underset{\mathbf{x}\in C}{\mathrm{argmin}}\Vert \mathbf{x}-{\bar{\mathbf{x}}}^k{\Vert }^2={\Pi }_C({\bar{\mathbf{x}}}^k),$$

$$\bar{\bm{x}}^{k}=(\alpha+\rho)^{-1}\big{(}\alpha\bm{x}^{k}+\rho\tilde{\bm{x}}^{k}-\nabla\!f(\bm{x}^{k})\big{)}.$$

$${\hat{\mathbf{v}}}_s=\frac{{\sum }_{r\in {\Omega }_s}{k}_{s,r}{\tilde{\mathbf{v}}}_r}{{\sum }_{r\in {\Omega }_s}{k}_{s,r}},k_{s,r}=\exp \left(-\frac{\Vert {\mathbf{P}}_s-{\mathbf{P}}_r{\Vert }^2}{2N_p^2{\sigma }^2}\right),$$

$$w_{s,r}\leftarrow {\Lambda }_{s,r}{k}_{s,r},{\Lambda }_{s,r}=\Lambda \left(\frac{s-r}{N_s+1}\right),$$

$$w_{s,r}\leftarrow w_{s,r}{\left(\sum _{r\in {\Omega }_s}{w}_{s,r}\right)}^{-\frac{1}{2}}{\left(\sum _{s\in {\Omega }_r}{w}_{s,r}\right)}^{-\frac{1}{2}},$$

$$w_{s,r}\leftarrow \alpha w_{s,r}\text{and}{w}_{s,s}\leftarrow 1+w_{s,s}-\sum _{r\in {\Omega }_s}{w}_{s,r}.$$

$$\mathbf{y}=\mathbf{A}\mathbf{x}+\mathbf{n},$$

$$f(\mathbf{x})=\frac{1}{2}\Vert \mathbf{y}-\mathbf{A}\mathbf{x}{\Vert }^2.$$

$$f(\mathbf{x})=\sum _{i=1}^{n}h(x_i),$$

$$h(t)=-K_{i}^{0}\log\big{(}e^{-\frac{\eta t}{K}}\big{)}-K_{i}^{1}\log\big{(}1-e^{-\frac{\eta t}{K}}\big{)},$$

$$h^{\prime }(t)=\frac{\eta }{K}\left(K_i^0-K_i^1{e}^{-\frac{\eta t}{K}}{\left(1-e^{-\frac{\eta t}{K}}\right)}^{-1}\right).$$

Full text

Index Terms—  image restoration, ADMM, plug-and-play, nonlocal means, convergence.

1 Introduction

Following the significant progress in image denoising, researchers have experimented with the idea of using modern denoisers for image restoration problems such as deconvolution, deblurring, tomography, and compressed sensing [1, 2, 3, 4]. Generally speaking, these are based on iterative methods, where each iteration involves the inversion of the forward model followed by the application of a powerful denoiser. Despite their empirical success, it has generally been challenging to furnish guarantees on convergence and optimality. The technical hurdle in this regard stems from the fact that state-of-the-art denoisers such as NLM [5] and BM3D [6] are derived from a filtering and not an optimization perspective. To be precise, it is not known if they can be expressed as the proximal map [7] of some regularizer $g(\bm{x})$, i.e., whether the minimizer of

[TABLE]

corresponds to the output obtained by denoising $\tilde{\bm{x}}$.

It was recently shown in [8] that the existence of such a regularizer can be guaranteed for a doubly stochastic variant of NLM. In particular, it was shown that if $g(\bm{x})$ is the regularizer in question and $f(\bm{x})$ is data fidelity term associated with the forward model, then the ADMM-based solution [7] of

[TABLE]

provides a framework where the inversion and regularization steps can be decoupled, and the latter amounts to denoising using doubly stochastic NLM (see Section 2 for a detailed description). In fact, the idea of replacing the regularization step in ADMM with an arbitrary denoiser was proposed earlier in [9] under the name “plug-and-play” ADMM (PnP-ADMM), albeit without any convergence guarantees. Later, it was also shown in [10] that the iterates of PnP-ADMM are guaranteed to converge to a fixed point if the denoiser satisfies a certain boundedness property. However, the question of optimality of the fixed point was not resolved in [10]. Moreover, it is not known if any of the existing denoisers satisfy this property.

The present work has two-fold contributions. First, we propose to “linearize” the data term $f(\bm{x})$ for the proximal update

[TABLE]

used in PnP-ADMM [8]. Linearization is used in several algorithms such as ISTA [11], FISTA [12], and ADMM [13]. However, to the best of our knowledge, this has not been exploited for PnP-ADMM. Linearization allows us to perform the inversion updates at lower cost (without sacrificing convergence guarantees) for applications where (3) has to be computed in an iterative fashion. We note that techniques to cut down the inversion cost have been proposed in [10, 14]. It was shown in [14] that the inversion can be performed efficiently using primal-dual splitting, and fast inversion techniques for specific restoration problems were proposed in [10]. Our proposal requires $f(\bm{x})$ to be differentiable. Moreover, we can guarantee convergence if $\nabla\!f(\bm{x})$ is Lipschitz continuous. This is indeed the case for a wide range of problems including linear inverse problems [11]. Importantly, we can incorporate hard constraints in the optimization framework, which is difficult to do in PnP-ADMM [8]. We note that even though $\nabla\!f(\bm{x})$ fails to be Lipschitz for single-photon imaging [15], linearized ADMM is found to be stable and yields good reconstructions that are comparable with [10] (see Section 4.2).

In a different direction, we also develop a fast low-complexity algorithm for doubly stochastic NLM (DSG-NLM). In [8], the authors conceived DSG-NLM as a matrix multiplication $\mathbf{W}\bm{x}$, where $\bm{x}$ is the vectorized input image and $\mathbf{W}$ is a weight matrix computed from patches. However, $\mathbf{W}$ is derived from the weight matrix of original NLM [5] in three steps (see Section 3). As a result, it was not apparent in [8] that $\mathbf{W}\bm{x}$ can be computed using filtering (aggregation), without storing $\mathbf{W}$. We show that this can indeed be done using aggregation, however the aggregation needs to be performed not once but thrice. We also use an existing algorithm [16] for efficiently computing the patch distances. On the overall, we are able to accelerate the implementation of DSG-NLM by about $80\times$, for which PnP-ADMM (and the proposed linearization) comes with convergence and optimality guarantees.

In Section 2, we propose linearized PnP-ADMM and discuss its convergence properties. A filtering perspective of DSG-NLM is presented in Section 3, which is used to develop a fast algorithm. In Section 4, we present results for super-resolution and single-photon imaging. We conclude with a summary of the results in Section 5.

2 Linearized Plug-and-Play ADMM

The plug-and-play framework was originally proposed in [9]. Following this original work, its algorithmic properties and applications have been studied in a series of papers [8, 17, 18, 10, 19, 14]. The core idea is based on alternating direction method of multipliers (ADMM), which is tailored for solving composite optimization problems of the form in (2), e.g., see [7, 20, 13, 21]. In this work, we will consider the optimization model

[TABLE]

where the data term $f(\bm{x})$ forces consistency w.r.t. the measurements, the regularizer $g(\bm{x})$ enforces some prior, and $C\subset\mathbb{R}^{n}$ is closed and convex. Constraints of interest include $C=\mathbb{R}^{n}_{+}$ (non-negativity constraint) and $C=[0,1]^{n}$ (box constraint) [20]. We next split the variable in (4) to obtain the following equivalent problem:

[TABLE]

The important point here is that though $\bm{x}$ is constrained, $\bm{v}$ (used in the denoising step) is a free variable which is indirectly constrained by the relation $\bm{x}=\bm{v}$. The augmented Lagrangian for (5) is given by

[TABLE]

where $\bm{u}$ is the (scaled) Lagrange multiplier and $\rho>0$ is the penalty parameter [7]. The ADMM solution of (2) consists of the minimization of $\mathfrak{L}(\bm{x},\bm{v}^{k},\bm{u}^{k})$ over $\bm{x}\in C$:

[TABLE]

the minimization of $\mathfrak{L}(\bm{x}^{k+1},\bm{v},\bm{u}^{k})$ over $\bm{v}\in\mathbb{R}^{n}$:

[TABLE]

where $\tilde{\bm{x}}^{k}=\bm{v}^{k}-\bm{u}^{k}$ and $\tilde{\bm{v}}^{k}=\bm{x}^{k+1}+\bm{u}^{k}$, and the dual update

[TABLE]

Notice that (7) is the maximum-a-posteriori estimator corresponding to the observation $\tilde{\bm{v}}^{k}=\bm{v}+\bm{n}$ and prior $p(\bm{v})\propto\exp(-g(\bm{v}))$, where $\bm{n}$ is iid Gaussian with variance $\sigma^{2}=\lambda/\rho$ [22]. Different choices of $g(\bm{v})$ lead to different denoising methods, such as wavelet and total-variation denoising [20]. The key idea in [9] was to substitute step (7) with

[TABLE]

where $\mathfrak{D}_{\sigma}$ is some powerful denoiser such as NLM or BM3D (optimized to operate at noise level $\sigma$) that do not arise from any known regularizer. It was empirically demonstrated in [9, 8, 17, 18, 10, 19] that this ad-hoc modification of the ADMM algorithm produces promising results for image restoration. Using the theory of proximal maps [23], it was later demonstrated in [8] that one can associate a convex regularizer with NLM provided its weight matrix is modified to be symmetric and doubly stochastic, with eigenvalues in $[0,1]$ (discussed in detail in Section 3). Moreover, it was shown that convergence is guaranteed in this case.

In general, the optimization in (6) has to be performed iteratively, even when $C=\mathbb{R}^{n}$. That is, along with the outer iterations, we need to perform inner iterations for (6). This is particularly applicable for the applications in Section 4. Ideally, we would like to replace (6) by a simple low-complexity operation. This can be achieved using a variant of ADMM called linearized ADMM [13]. The idea is to replace $f(\bm{x})$ by its linear approximation around $\bm{x}^{k}$:

[TABLE]

where we also append a quadratic term ($\alpha>0$) for technical reasons to be discussed shortly. Substituting (9) in (6), we obtain

[TABLE]

that is, $\bm{x}^{k+1}$ is the orthogonal projection of $\bar{\bm{x}}^{k}$ onto $C$, where

[TABLE]

We note that $\Pi_{C}$ can be easily computed when $C$ is $\mathbb{R}^{n}_{+}$ or $[0,1]^{n}$. The assumption with linearization is that $f(\bm{x})$ is continuously differentiable and $\nabla\!f(\bm{x})$ is Lipschitz (and easily computable). This is true for a wide range of image restoration problems including linear inverse problems. The plug-and-play algorithm using linearized ADMM is summarized below.

If $\mathfrak{D}_{\sigma}$ is DSG-NLM, then Algorithm 1 is guaranteed to converge under some mild assumption on $f(\bm{x})$. Indeed, as shown in [8], we can associate a closed, proper, and convex regularizer $g(\bm{x})$ with $\mathfrak{D}_{\sigma}$ in this case, i.e., (8) corresponds to the ADMM update (7). Algorithm 1 is then simply an application of linearized ADMM to problem (4), and its convergence follows from existing results [13]. Specifically, if $f(\bm{x})$ is closed, proper, and convex (and satisfies some mild assumptions [13]), and $\alpha$ is larger than the Lipschitz constant of $\nabla\!f(\bm{x})$, then $f(\bm{x}^{k})+\lambda g(\bm{x}^{k})$ converges to the optimum of (4).

3 Fast DSG-NLM

We now propose a fast algorithm for DSG-NLM [8]. This denoiser is derived from NLM [5] which uses patches and nonlocal aggregation for denoising. We will follow the notation in [8] for easy comparison. In particular, we will use $\tilde{\bm{v}}$ and $\hat{\bm{v}}$ to denote the input to the denoiser and the output, both defined on a finite domain $S\subset\mathbb{Z}^{2}$. Moreover, we will use $\mathbf{P}_{s}$ (vector of length $N_{p}^{2}$) to denote a patch of size $N_{p}\times N_{p}$ centered at $s\in S$. In NLM, $\hat{\bm{v}}=\mathfrak{D}_{\sigma}(\tilde{\bm{v}})$ is given by

[TABLE]

where $\Omega_{s}$ is a search window of size $(2N_{s}+1)^{2}$ centered at $s$ [8].

Notice that we can express (10) as a linear transform $\hat{\bm{v}}=\mathbf{K}\tilde{\bm{v}}$, where $\mathbf{K}$ (derived from $\tilde{\bm{v}}$) is row stochastic, i.e., the sum of entries in each row is one. However, due to the division in (10), $\mathbf{K}$ is not symmetric in general. It was shown in [8, Section IV] that, in three simple steps, we can transform $\mathbf{K}$ into a symmetric, doubly stochastic matrix $\mathbf{W}$ whose eigenvalues are in $[0,1]$. DSG-NLM is simply the transform $\mathbf{W}\tilde{\bm{v}}$. For completeness, we recall these steps, but slightly differently from [8]. In the first step, we set $\mathbf{W}=[w_{s,r}]$, where

[TABLE]

and $\Lambda(s)=(1-\lvert s_{1}\lvert)(1-\lvert s_{2}\lvert)$ is the separable hat function. We refer the reader to [8] for the technical reason behind this and subsequent steps. The next step involves row and column normalization:

[TABLE]

We next compute the row sums of the weight matrix (at this point) and set the maximum to $1/\alpha$. The final step is given by

[TABLE]

It is proved in [8, Theorem IV.1] that the final matrix $\mathbf{W}$ is symmetric and doubly stochastic, and its eigenvalues are in $[0,1]$. As for the implementation, it is impractical to store the large matrix $\mathbf{W}$ and apply it on $\tilde{\bm{v}}$ as a matrix-vector multiplication. The present observation is that we can efficiently compute $\mathbf{W}\tilde{\bm{v}}$ without having to store $\mathbf{W}$. However, unlike NLM, we need to loop over the pixels three times instead of just once. The complete procedure is provided in Algorithm 2, where images $\bm{g}$ and $\bm{h}$ are of the size of $\tilde{\bm{v}}$. The first aggregation corresponds to (11). In the second aggregation, we compute (12) and the maximum of the row sums. In the final aggregation, we combine step (13) and the application of $\mathbf{W}$ on $\tilde{\bm{v}}$.

Similar to NLM, it is evident that the complexity of Algorithm 2 is $O(N_{s}^{2}N_{p}^{2})$ per pixel. However, similar to NLM, we can use the fast algorithm in [16] for computing patch distances, whose complexity does not depend on the patch size $N_{p}$. This significantly brings down the complexity to $O(N_{s}^{2})$. In summary, we can compute DSG-NLM via filtering and that too at a reduced complexity.

4 Experiments

We validate the performance of our restoration algorithm using super-resolution [24] and single photon imaging [15]. The reason behind this choice is that the data fidelity term is quadratic in super-resolution and non-quadratic (though convex and differentiable) in single photon imaging. Assuming (without loss of generality) that the intensity of the ground-truth is in $[0,1]$, we set the constraint set $C$ in (4) to be $C=[0,1]^{n}$ For both applications, we have used BM3D, NLM and DSG-NLM for denoising. Moreover, similar to [8], we have tried both the adaptive (A-DSG-NLM) and fixed (F-DSG-NLM) variants of DSG-NLM. In the former, the weight $\mathbf{W}$ is adapted in each iteration as discussed previously. In F-DSG-NLM, we stop adapting $\mathbf{W}$ after certain number of iterations (typically $15$), which is used for rest of the iterations. The point to note is that convergence is guaranteed for F-DSG-NLM but not for A-DSG-NLM, since the output is not a linear function of the input in the latter [8].

In Table 1, we compare the timings of DSG-NLM [8] and Algorithm 2 for a $256\times 256$ image. The simulation was done using Matlab on a quad-core $3.4$ GHz machine with $16$ GB memory. We used standard patch and window sizes [5]. The brute-force implementation of nonlocal denoising is known to be prohibitively slow [16]. In Table 1, notice that the brute-force implementation of DSG-NLM takes minutes, while Algorithm 2 takes just few seconds. The speedup is almost $100\times$ for large patch sizes. This is because the complexity of our fast algorithm does not scale with the patch size.

4.1 Single image super-resolution

The forward model for single image super-resolution is given by

[TABLE]

where $\bm{x}$ is the high-resolution image (ground truth), $\mathbf{A}$ is a linear operator, and $\bm{n}$ is iid Gaussian noise [24]. Specifically, $\mathbf{A}\bm{x}$ is the low-pass filtering of $\bm{x}$ (using blur $\bm{h}$) followed by downsampling (by factor $k$). The adjoint $\mathbf{A}^{\top}\!$ corresponds to upsampling (by factor $k$) and low-pass filtering with $\bm{h}$ (if $\bm{h}$ is symmetric) [20]. The problem is to recover $\bm{x}$ from the low-resolution image $\bm{y}$. As is well-known [22], the data fidelity term corresponding to the negative log-likelihood of $\bm{y}$ given $\bm{x}$ is

[TABLE]

The gradient and Hessian of (14) are $\mathbf{A}^{\top}\!(\bm{y}-\mathbf{A}\bm{x})$ and $\mathbf{A}^{\top}\!\mathbf{A}$. The Lipschitz constant of $\nabla\!f(\bm{x})$ is simply the largest eigenvalue of $\mathbf{A}^{\top}\!\mathbf{A}$ [20]. If we set $\alpha>\lambda_{\text{max}}(\mathbf{A}^{\top}\!\mathbf{A})$ and use A-DSG-NLM for denoising, then convergence is guaranteed for Algorithm 1.

We now present some results for image super-resolution where $\bm{h}$ (blur) is Gaussian. A typical result is reported in Figure 1, where we have also compared with the plug-and-play algorithm in [10]. For both methods, we have used F-DSG-NLM as the denoiser. In [10], the inversion step (6) is computed using fast polyphase decompositions and $\rho$ is adapted at each iteration. Notice that the reconstructions are comparable, both visually and in terms of PSNR. We note that, unlike Algorithm 1, the fast inversion in [10] works just with periodic convolutions. Following this observation, we perform a super-resolution experiment for the setup in Figure 1, but using symmetric boundary condition (for the blur). In this case, the inversion step (6) requires the solution a linear system with coefficients $\mathbf{A}^{\top}\!\mathbf{A}+\rho\mathbf{I}$. This can be done iteratively using conjugate-gradients, which is however computation intensive. The evolution of PSNR with time for standard and linearized plug-and-play ADMM are compared in Figure 2. Notice that the PSNR peaks very fast in linearized ADMM. A possible explanation is that though we use an inexact update for (6), this is partly compensated by the regularization (smoothing) step in each iteration.

4.2 Single photon imaging

We next consider single photon imaging using quanta image sensors (QIS) [15]. For an input image $\bm{x}\in\mathbb{R}^{n}$, the QIS consists of $nK$ photon detectors, which are uniformly distributed to read the incoming light ($K$ is the oversampling factor). A simplified model of the amount of photon $\bm{s}$ arriving at QIS pixels is given by $\bm{s}=\eta\mathbf{G}\bm{x}$, where $\eta$ is the sensor gain, and $\mathbf{G}=(1/K)(\mathbf{I}_{n\times n}\otimes\mathbf{1}_{K\times 1})$. More precisely, the photon count follows a Poisson distribution with parameter $\bm{s}$. The final output $\bm{y}$ is a binary image obtained by thresholding the photon count (see [10] for a detailed description). In brief, the negative log-likelihood of observing $\bm{y}$ given $\bm{x}$ is given by

[TABLE]

where

[TABLE]

and $K_{i}^{0},K_{i}^{1}$ are known design parameters. Since (15) is separable, its gradient is $\nabla\!f(\bm{x})=\left(h^{\prime}(x_{1}),\ldots,h^{\prime}(x_{n})\right)$, where

[TABLE]

It can be verified that $h^{\prime\prime}(t)$ is not bounded near the origin. Therefore, $\nabla\!f(\bm{x})$ cannot be Lipschitz. However, this is only a sufficient condition for convergence, and Algorithm 1 in fact was found to be work stably for our experiments. As shown in Figure 3, the reconstruction obtained using Algorithm 1 is comparable to that obtained using [10].

For both experiments, we computed the primal and dual residuals [7] and the PSNR for different denoisers. The results are reported in Tables 2 and 3. We observe that the residuals for BM3D, NLM, and to some extent A-DSG-NLM, do not seem to converge fully. On the other hand, the residual for F-DSG-NLM is several orders smaller than NLM and BM3D. For single photon imaging, the residual for F-DSG-NLM is $4$ orders smaller than NLM and $10$ orders smaller than BM3D. Interestingly, the PSNR for F-DSG-NLM is better than that of BM3D or NLM (a similar observation was reported in [8]).

5 Conclusion

We proposed a plug-and-play algorithm for image restoration, where both the inversion and denoising steps can be computed efficiently. The proposed algorithm comes with convergence guarantees for linear inverse problems. For super-resolution and single photon imaging, our algorithm was shown to yield reconstructions that are comparable with existing plug-and-play algorithms. Though this was not investigated due to space constraints, we can use accelerated variants of linearized ADMM to further speed up the reconstruction [13]. Another interesting possibility is to use improved variants of nonlocal means [25, 26] for the denoising step.